diff --git a/Project.toml b/Project.toml
index 888f80dd..6faf2c52 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FinEtools"
 uuid = "91bb5406-6c9a-523d-811d-0644c4229550"
 authors = ["Petr Krysl <pkrysl@ucsd.edu>"]
-version = "7.0.5"
+version = "7.1.0"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
diff --git a/README.md b/README.md
index 4e1f9c74..f17096f3 100644
--- a/README.md
+++ b/README.md
@@ -21,6 +21,7 @@ The package supports application packages, for instance:
 
 ## News
 
+- 10/07/2023: Update for Julia 1.10.
 - 06/19/2023: Introduce DataCache, generic linear and bilinear forms. 
 - 04/20/2023: Make all types in the library generic.
 - 04/18/2023: Implemented resizing of assembly buffers. Makematrix! resets pointers.
diff --git a/docs/issues-and-ideas.md b/docs/issues-and-ideas.md
index 316dcc9b..8f4128f2 100644
--- a/docs/issues-and-ideas.md
+++ b/docs/issues-and-ideas.md
@@ -1317,3 +1317,6 @@ end # module
    of the finite element instead of its label?
    Done. Also the quadrature point id.
    
+-- The delegation mechanism needs to be described. 09/29/23
+
+-- The finite element model machine documentation needs to be refreshed. 09/29/23
diff --git a/src/FEMMBaseModule.jl b/src/FEMMBaseModule.jl
index 65de392d..8620bca6 100644
--- a/src/FEMMBaseModule.jl
+++ b/src/FEMMBaseModule.jl
@@ -94,7 +94,7 @@ function finite_elements(self::AbstractFEMM)
 end
 
 """
-    associategeometry!(self::AbstractFEMM, geom::NodalField{FT}) where {FT}
+    associategeometry!(self::AbstractFEMM, geom::NodalField{GFT}) where {GFT}
 
 Associate geometry field with the FEMM.
 
@@ -108,44 +108,44 @@ routine is not consistent with the one for which `associategeometry!()`
 was called before, `associategeometry!()` needs to be called with
 the new geometry field.
 """
-function associategeometry!(self::AbstractFEMM, geom::NodalField{FT}) where {FT}
+function associategeometry!(self::AbstractFEMM, geom::NodalField{GFT}) where {GFT}
     return self # default is no-op
 end
 
 """
     inspectintegpoints(
         self::FEMM,
-        geom::NodalField{FT},
+        geom::NodalField{GFT},
         felist::AbstractVector{IT},
         inspector::F,
         idat,
         quantity = :Cauchy;
         context...,
-    ) where {FEMM<:AbstractFEMM, FT, IT, F<:Function}
+    ) where {FEMM<:AbstractFEMM, GFT, IT, F<:Function}
 
 Inspect integration points.
 """
 function inspectintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
+    geom::NodalField{GFT},
     felist::AbstractVector{IT},
     inspector::F,
     idat,
     quantity = :Cauchy;
     context...,
-) where {FEMM<:AbstractFEMM,FT,IT,F<:Function}
+) where {FEMM<:AbstractFEMM, GFT, IT, F<:Function}
     return idat # default is no-op
 end
 
 """
     integratefieldfunction(
         self::AbstractFEMM,
-        geom::NodalField{FT},
+        geom::NodalField{GFT},
         afield::FL,
         fh::F;
         initial::R,
         m = -1,
-    ) where {FT, FL<:NodalField{T}, T, F<:Function, R}
+    ) where {GFT, T, FL<:NodalField{T}, F<:Function,R}
 
 Integrate a nodal-field function over the discrete manifold.
 
@@ -161,12 +161,12 @@ Returns value of type `R`, which is initialized by `initial`.
 """
 function integratefieldfunction(
     self::AbstractFEMM,
-    geom::NodalField{FT},
+    geom::NodalField{GFT},
     afield::FL,
     fh::F;
     initial::R,
     m = -1,
-) where {FT,T,FL<:NodalField{T},F<:Function,R}
+) where {GFT, T, FL<:NodalField{T}, F<:Function,R}
     fes = finite_elements(self)
     if m < 0
         m = manifdim(fes)  # native  manifold dimension
@@ -192,12 +192,12 @@ end
 """
     integratefieldfunction(
         self::AbstractFEMM,
-        geom::NodalField{FT},
+        geom::NodalField{GFT},
         afield::FL,
         fh::F;
         initial::R = zero(FT),
         m = -1,
-    ) where {FT, T, FL<:ElementalField{T}, F<:Function, R}
+    ) where {GFT, T, FL<:ElementalField{T}, F<:Function,R}
 
 Integrate a elemental-field function over the discrete manifold.
 
@@ -212,12 +212,12 @@ Returns value of type `R`, which is initialized by `initial`.
 """
 function integratefieldfunction(
     self::AbstractFEMM,
-    geom::NodalField{FT},
+    geom::NodalField{GFT},
     afield::FL,
     fh::F;
     initial::R = zero(FT),
     m = -1,
-) where {FT,T,FL<:ElementalField{T},F<:Function,R}
+) where {GFT, T, FL<:ElementalField{T}, F<:Function,R}
     fes = finite_elements(self)
     if m < 0
         m = manifdim(fes)  # native  manifold dimension
@@ -241,11 +241,11 @@ end
 """
     integratefunction(
         self::AbstractFEMM,
-        geom::NodalField{FT},
+        geom::NodalField{GFT},
         fh::F;
-        initial::R = zero(FT),
+        initial::R = zero(typeof(fh(zeros(ndofs(geom), 1)))),
         m = -1,
-    ) where {FT, F<:Function, R<:Number}
+    ) where {GFT<:Number, F<:Function, R<:Number}
 
 Integrate a function over the discrete manifold.
 
@@ -261,10 +261,10 @@ measure the static moment with respect to the y- axis, and so on.
 # Example:
 Compute the volume of the mesh and then its center of gravity:
 ```
-V = integratefunction(femm, geom, (x) ->  1.0)
-Sx = integratefunction(femm, geom, (x) ->  x[1])
-Sy = integratefunction(femm, geom, (x) ->  x[2])
-Sz = integratefunction(femm, geom, (x) ->  x[3])
+V = integratefunction(femm, geom, (x) ->  1.0, 0.0)
+Sx = integratefunction(femm, geom, (x) ->  x[1], 0.0)
+Sy = integratefunction(femm, geom, (x) ->  x[2], 0.0)
+Sz = integratefunction(femm, geom, (x) ->  x[3], 0.0)
 CG = vec([Sx Sy Sz]/V)
 ```
 Compute a moment of inertia of the mesh relative to the origin:
@@ -274,11 +274,11 @@ Ixx = integratefunction(femm, geom, (x) ->  x[2]^2 + x[3]^2)
 """
 function integratefunction(
     self::AbstractFEMM,
-    geom::NodalField{FT},
+    geom::NodalField{GFT},
     fh::F;
-    initial::R = zero(FT),
+    initial::R = zero(typeof(fh(zeros(ndofs(geom), 1)))),
     m = -1,
-) where {FT, F<:Function, R<:Number}
+) where {GFT<:Number, F<:Function, R<:Number}
     fes = finite_elements(self)
     if m < 0
         m = manifdim(fes)  # native  manifold dimension
@@ -620,13 +620,13 @@ end
 """
     fieldfromintegpoints(
         self::FEMM,
-        geom::NodalField{FT},
-        u::NodalField{T},
-        dT::NodalField{FT},
+        geom::NodalField{GFT},
+        u::NodalField{UFT},
+        dT::NodalField{TFT},
         quantity::Symbol,
-        component::FIntVec;
+        component::AbstractVector{IT};
         context...,
-    ) where {FEMM<:AbstractFEMM, FT<:Number, T<:Number}
+    ) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, TFT<:Number, IT<:Integer}
 
 Construct nodal field from integration points.
 
@@ -648,17 +648,18 @@ Keyword arguments
 """
 function fieldfromintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
-    u::NodalField{T},
-    dT::NodalField{FT},
+    geom::NodalField{GFT},
+    u::NodalField{UFT},
+    dT::NodalField{TFT},
     quantity::Symbol,
     component::AbstractVector{IT};
     context...,
-) where {FEMM<:AbstractFEMM,FT<:Number,T<:Number,IT<:Integer}
+) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, TFT<:Number, IT<:Integer}
     fes = finite_elements(self)
     # Constants
     nne = nodesperelem(fes) # number of nodes for element
     sdim = ndofs(geom)            # number of space dimensions
+    FT = typeof(zero(GFT) + zero(UFT) + zero(TFT))
     nodevalmethod = :invdistance
     reportat = :default
     for apair in pairs(context)
@@ -738,36 +739,36 @@ end
 
 function fieldfromintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
-    u::NodalField{T},
-    dT::NodalField{FT},
+    geom::NodalField{GFT},
+    u::NodalField{UFT},
+    dT::NodalField{TFT},
     quantity::Symbol,
     component::IT;
     context...,
-) where {FEMM<:AbstractFEMM,FT<:Number,T<:Number,IT<:Integer}
+) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, TFT<:Number, IT<:Integer}
     return fieldfromintegpoints(self, geom, u, dT, quantity, [component]; context...)
 end
 
 function fieldfromintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
-    u::NodalField{T},
+    geom::NodalField{GFT},
+    u::NodalField{UFT},
     quantity::Symbol,
     component::AbstractVector{IT};
     context...,
-) where {FEMM<:AbstractFEMM,FT<:Number,T<:Number,IT<:Integer}
+) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, IT<:Integer}
     dT = NodalField(zeros(FT, nnodes(geom), 1)) # zero difference in temperature
     return fieldfromintegpoints(self, geom, u, dT, quantity, component; context...)
 end
 
 function fieldfromintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
-    u::NodalField{T},
+    geom::NodalField{GFT},
+    u::NodalField{UFT},
     quantity::Symbol,
     component::IT;
     context...,
-) where {FEMM<:AbstractFEMM,FT<:Number,T<:Number,IT<:Integer}
+) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, IT<:Integer}
     dT = NodalField(zeros(FT, nnodes(geom), 1)) # zero difference in temperature
     return fieldfromintegpoints(self, geom, u, dT, quantity, [component]; context...)
 end
@@ -781,13 +782,13 @@ end
 """
     elemfieldfromintegpoints(
         self::FEMM,
-        geom::NodalField{FT},
-        u::NodalField{T},
-        dT::NodalField{FT},
+        geom::NodalField{GFT},
+        u::NodalField{UFT},
+        dT::NodalField{TFT},
         quantity::Symbol,
-        component::FIntVec;
+        component::AbstractVector{IT};
         context...,
-    ) where {FEMM<:AbstractFEMM, FT<:Number, T<:Number}
+    ) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, TFT<:Number, IT<:Integer}
 
 Construct elemental field from integration points.
 
@@ -805,17 +806,18 @@ Construct elemental field from integration points.
 """
 function elemfieldfromintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
-    u::NodalField{T},
-    dT::NodalField{FT},
+    geom::NodalField{GFT},
+    u::NodalField{UFT},
+    dT::NodalField{TFT},
     quantity::Symbol,
     component::AbstractVector{IT};
     context...,
-) where {FEMM<:AbstractFEMM,FT<:Number,T<:Number,IT<:Integer}
+) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, TFT<:Number, IT<:Integer}
     fes = finite_elements(self)
     # Constants
     nne = nodesperelem(fes) # number of nodes for element
     sdim = ndofs(geom)            # number of space dimensions
+    FT = typeof(zero(GFT) + zero(UFT) + zero(TFT))
     # Container of intermediate results
     idat = MeanValueInspectorData(
         zeros(IT, count(fes)),
@@ -859,36 +861,36 @@ end
 
 function elemfieldfromintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
-    u::NodalField{T},
-    dT::NodalField{FT},
+    geom::NodalField{GFT},
+    u::NodalField{UFT},
+    dT::NodalField{TFT},
     quantity::Symbol,
     component::IT;
     context...,
-) where {FEMM<:AbstractFEMM,FT<:Number,T<:Number,IT<:Integer}
+) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, TFT<:Number, IT<:Integer}
     return elemfieldfromintegpoints(self, geom, u, dT, quantity, [component]; context...)
 end
 
 function elemfieldfromintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
-    u::NodalField{T},
+    geom::NodalField{GFT},
+    u::NodalField{UFT},
     quantity::Symbol,
     component::IT;
     context...,
-) where {FEMM<:AbstractFEMM,FT<:Number,T<:Number,IT<:Integer}
+) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, IT<:Integer}
     dT = NodalField(zeros(FT, nnodes(geom), 1)) # zero difference in temperature
     return elemfieldfromintegpoints(self, geom, u, dT, quantity, [component]; context...)
 end
 
 function elemfieldfromintegpoints(
     self::FEMM,
-    geom::NodalField{FT},
-    u::NodalField{T},
+    geom::NodalField{GFT},
+    u::NodalField{UFT},
     quantity::Symbol,
     component::AbstractVector{IT};
     context...,
-) where {FEMM<:AbstractFEMM,FT<:Number,T<:Number,IT<:Integer}
+) where {FEMM<:AbstractFEMM, GFT<:Number, UFT<:Number, IT<:Integer}
     dT = NodalField(zeros(FT, nnodes(geom), 1)) # zero difference in temperature
     return elemfieldfromintegpoints(self, geom, u, dT, quantity, component; context...)
 end
diff --git a/src/FinEtools.jl b/src/FinEtools.jl
index cdcf691d..2e470a99 100644
--- a/src/FinEtools.jl
+++ b/src/FinEtools.jl
@@ -17,10 +17,6 @@ include("allmodules.jl")
 ###########################################################################
 # General facilities
 ###########################################################################
-using .FTypesModule:
-    FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
-# Exported: basic numerical types, type of data dictionary
-export FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
 
 using .DataCacheModule: DataCache
 # Exported: vector-cache type and methods to invoke the update callback
diff --git a/test/test_basics.jl b/test/test_basics.jl
index 6c0e4683..afc1c1de 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -315,7 +315,7 @@ function test()
         0.8320502943378436 0.5547001962252293 0.0
         0.0 0.0 1.0
     ]
-    # function compute!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, feid::FInt)
+    # function compute!(csmatout, XYZ, tangents, feid)
     #     # Cylindrical coordinate system
     #     xyz = XYZ[:] .- center
     #     xyz[3] = 0.0
@@ -346,7 +346,7 @@ function test()
     XYZ = reshape([0.2, 0.3, 0.4], 1, 3)
     tangents = rand(3, 3)
     rfcsmat = Matrix(1.0 * I, 3, 3)
-    # function compute!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, feid::FInt)
+    # function compute!(csmatout, XYZ, tangents, feid)
     #     # Cylindrical coordinate system
     #     xyz = XYZ[:] .- center
     #     xyz[3] = 0.0
@@ -442,59 +442,59 @@ using Test
 function test()
     let
         N = 8
-        Ke::FFltMat, gradN::FFltMat, mult::FFlt = fill(0.0, N, N), rand(N, 2), 3.0
-        add_mggt_ut_only!(Ke::FFltMat, gradN::FFltMat, mult::FFlt)
-        complete_lt!(Ke::FFltMat)
+        Ke, gradN, mult = fill(0.0, N, N), rand(N, 2), 3.0
+        add_mggt_ut_only!(Ke, gradN, mult)
+        complete_lt!(Ke)
         @test norm(Ke - (mult .* (gradN * gradN'))) / norm(Ke) <= 1.0e-9
     end
 
 
     let
         N = 8
-        Ke::FFltMat, gradN::FFltMat, Jac_w::FFlt = fill(0.0, N, N), rand(N, 3), 0.33
-        kappa_bar::FFltMat = rand(3, 3)
+        Ke, gradN, Jac_w = fill(0.0, N, N), rand(N, 3), 0.33
+        kappa_bar = rand(3, 3)
         kappa_bar = kappa_bar + kappa_bar'
-        kappa_bargradNT::FFltMat = fill(0.0, 3, N)
+        kappa_bargradNT = fill(0.0, 3, N)
         add_gkgt_ut_only!(
-            Ke::FFltMat,
-            gradN::FFltMat,
-            Jac_w::FFlt,
-            kappa_bar::FFltMat,
-            kappa_bargradNT::FFltMat,
+            Ke,
+            gradN,
+            Jac_w,
+            kappa_bar,
+            kappa_bargradNT,
         )
-        complete_lt!(Ke::FFltMat)
+        complete_lt!(Ke)
         @test norm(Ke - (Jac_w .* (gradN * kappa_bar * gradN'))) / norm(Ke) <= 1.0e-9
     end
 
 
     let
         N = 12
-        Ke::FFltMat, B::FFltMat, Jac_w::FFlt = fill(0.0, N, N), rand(3, N), 0.33
-        D::FFltMat = rand(3, 3)
+        Ke, B, Jac_w = fill(0.0, N, N), rand(3, N), 0.33
+        D = rand(3, 3)
         D = D + D'
-        DB::FFltMat = fill(0.0, 3, N)
-        add_btdb_ut_only!(Ke::FFltMat, B::FFltMat, Jac_w::FFlt, D::FFltMat, DB::FFltMat)
-        complete_lt!(Ke::FFltMat)
+        DB = fill(0.0, 3, N)
+        add_btdb_ut_only!(Ke, B, Jac_w, D, DB)
+        complete_lt!(Ke)
         @test norm(Ke - (Jac_w .* (B' * D * B))) / norm(Ke) <= 1.0e-9
     end
 
 
     let
         N = 16
-        Fe::FFltVec, B::FFltMat, coefficient::FFlt, sigma::FFltVec =
+        Fe, B, coefficient, sigma =
             fill(0.0, N), rand(3, N), 0.33, rand(3)
-        D::FFltMat = rand(3, 3)
+        D = rand(3, 3)
         D = D + D'
-        add_btsigma!(Fe::FFltVec, B::FFltMat, coefficient::FFlt, sigma::FFltVec)
+        add_btsigma!(Fe, B, coefficient, sigma)
         @test norm(Fe - (coefficient * (B' * sigma))) / norm(Fe) <= 1.0e-9
     end
 
 
     let
         M = 21
-        Ke::FFltMat, Nn::FFltMat, Jac_w_coeff::FFlt = fill(0.0, M, M), rand(M, 1), 0.533
-        add_nnt_ut_only!(Ke::FMat{FFlt}, Nn::FFltMat, Jac_w_coeff::FFlt)
-        complete_lt!(Ke::FFltMat)
+        Ke, Nn, Jac_w_coeff = fill(0.0, M, M), rand(M, 1), 0.533
+        add_nnt_ut_only!(Ke, Nn, Jac_w_coeff)
+        complete_lt!(Ke)
         @test norm(Ke - ((Nn * Nn') * (Jac_w_coeff))) / norm(Ke) <= 1.0e-9
     end
 
@@ -583,10 +583,10 @@ using FinEtools.MeshExportModule: VTK
 using LinearAlgebra
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         100.0, 200.0, 3, 6, pi / 3, :a
     fens, fes =
-        T3annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol)
+        T3annulus(rin, rex, nr, nc, Angl, orientation::Symbol)
     # File = "mesh.vtk"
     # VTK.vtkexportmesh(File, fens, fes)
     geom = NodalField(fens.xyz)
@@ -1082,7 +1082,7 @@ function test()
     gradN, gradNparams, redJ =
         fill(0.0, nodesperelem(fes), 2), bfundpar(fes, [0.57, 0.57]), Matrix(1.0 * I, 2, 2)
 
-    gradN!(fes, gradN::FFltMat, gradNparams::FFltMat, redJ::FFltMat)
+    gradN!(fes, gradN, gradNparams, redJ)
     @test norm(
         gradN - [
             -0.10750000000000001 -0.10750000000000001
@@ -1115,7 +1115,7 @@ function test()
     bfundpar(fes, [0.57, 0.57, -0.57]),
     Matrix(1.0 * I, 3, 3)
 
-    gradN!(fes, gradN::FFltMat, gradNparams::FFltMat, redJ::FFltMat)
+    gradN!(fes, gradN, gradNparams, redJ)
     # @show gradN
     @test norm(
         gradN - [
@@ -1184,14 +1184,14 @@ using Test
 function test()
     ndimensions = 2
     tangents, feid, qpid = reshape([1.0; 1.0], 2, 1), 0, 0
-    n = SurfaceNormal(ndimensions::FInt)
-    normal = updatenormal!(n, [0.0 0.0 0.0], tangents::FFltMat, feid::FInt, qpid)
+    n = SurfaceNormal(ndimensions)
+    normal = updatenormal!(n, [0.0 0.0 0.0], tangents, feid, qpid)
     @test norm(normal - [0.7071067811865475, -0.7071067811865475]) <= 1.0e-5
 
     ndimensions = 3
     tangents, feid = reshape([1.0 0.0; 1.0 0.0; 0.0 1.0], 3, 2), 0
-    n = SurfaceNormal(ndimensions::FInt)
-    normal = updatenormal!(n, [0.0 0.0 0.0], tangents::FFltMat, feid::FInt, qpid)
+    n = SurfaceNormal(ndimensions)
+    normal = updatenormal!(n, [0.0 0.0 0.0], tangents, feid, qpid)
     @test norm(normal - [0.7071067811865475, -0.7071067811865475, 0.0]) <= 1.0e-5
     true
 end
@@ -2500,7 +2500,7 @@ using FinEtools.SurfaceNormalModule: SurfaceNormal, updatenormal!
 using LinearAlgebra
 using Test
 
-function __computenormal!(normalout::FFltVec, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+function __computenormal!(normalout, XYZ, tangents, feid, qpid)
     fill!(normalout, 0.0)
     # We are assuming a surface element here!
     if (size(tangents,1) == 3) && (size(tangents,2) == 2)# surface in three dimensions
diff --git a/test/test_meshing.jl b/test/test_meshing.jl
index 27afe3c1..5e6a9482 100644
--- a/test/test_meshing.jl
+++ b/test/test_meshing.jl
@@ -2418,9 +2418,9 @@ module mpointmm1
 using FinEtools
 using Test
 function test()
-    x::FFltMat = [i == j ? one(FFlt) : zero(FFlt) for i in 1:4, j in 1:3]
+    x = [i == j ? one(Float64) : zero(Float64) for i in 1:4, j in 1:3]
     fes = FESetP1(reshape([1 2 4 3], 4, 1))
-    pt::FFltVec = x[4, :]
+    pt = x[4, :]
     pc, success = map2parametric(fes, reshape(x[4, :], 1, 3), pt)
     @test success
     @test pc[1] == 00.0
@@ -3542,7 +3542,7 @@ using FinEtools
 using FinEtools.MeshExportModule
 using Test
 function test()
-    Radius::FFlt, Length::FFlt, nperradius, nL = 1.0, 2.0, 7, 9
+    Radius, Length, nperradius, nL = 1.0, 2.0, 7, 9
 
     fens, fes = H8cylindern(Radius, Length, nperradius, nL)
 
@@ -3562,7 +3562,7 @@ using FinEtools
 using FinEtools.MeshExportModule
 using Test
 function test()
-    Radius::FFlt, Length::FFlt, nperradius, nL = 1.0, 2.0, 17, 3
+    Radius, Length, nperradius, nL = 1.0, 2.0, 17, 3
 
     fens, fes = H8cylindern(Radius, Length, nperradius, nL)
 
@@ -4271,7 +4271,7 @@ using FinEtools
 using Test
 function test()
     img = fill(1, 30, 20, 10)
-    fens, fes = T4voximg(img, FFlt[1.0, 1.1, 1.3], [1])
+    fens, fes = T4voximg(img, Float64[1.0, 1.1, 1.3], [1])
     @test count(fes) == 30 * 20 * 10 * 5
     true
 end
@@ -4309,10 +4309,10 @@ using FinEtools
 using FinEtools.MeshExportModule: VTK
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         100.0, 200.0, 3, 6, pi / 3, :a
     fens, fes =
-        T6annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol)
+        T6annulus(rin, rex, nr, nc, Angl, orientation::Symbol)
     # @show count(fens), count(fes)
     @test (count(fens), count(fes)) == (91, 36)
     fens.xyz = xyz3(fens)
@@ -4331,10 +4331,10 @@ using FinEtools
 using FinEtools.MeshExportModule: VTK
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         100.0, 200.0, 3, 6, pi / 3, :a
     fens, fes =
-        T3annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol)
+        T3annulus(rin, rex, nr, nc, Angl, orientation::Symbol)
     # @show count(fens), count(fes)
     @test (count(fens), count(fes)) == (28, 36)
     fens.xyz = xyz3(fens)
@@ -4353,10 +4353,10 @@ using FinEtools
 using FinEtools.MeshExportModule: VTK
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         100.0, 200.0, 3, 6, pi / 3, :a
     fens, fes =
-        T3annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol)
+        T3annulus(rin, rex, nr, nc, Angl, orientation::Symbol)
     fens, fes = T3refine(fens, fes)
     # @show count(fens), count(fes)
     @test (count(fens), count(fes)) == (91, 144)
@@ -4426,7 +4426,7 @@ using FinEtools.MeshExportModule: VTK
 using Test
 function test()
     img = fill(1, 30, 20, 10)
-    fens, fes = H8voximg(img, FFlt[1.0, 1.1, 1.3], [1])
+    fens, fes = H8voximg(img, Float64[1.0, 1.1, 1.3], [1])
     @test (count(fens), count(fes)) == (7161, 6000)
     # File = "mesh.vtk"
     # VTK.vtkexportmesh(File, fens, fes)
@@ -4470,10 +4470,10 @@ using FinEtools
 using FinEtools.MeshExportModule: VTK
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         100.0, 200.0, 13, 16, pi / 3, :a
     fens, fes =
-        T3annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol)
+        T3annulus(rin, rex, nr, nc, Angl, orientation::Symbol)
     d = (fens.xyz[:, 1] .^ 2 + fens.xyz[:, 2] .^ 2)
     File = "mesh.vtk"
     # Export of multiple scalar fields
@@ -4491,10 +4491,10 @@ using FinEtools
 using FinEtools.MeshExportModule: STL
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         1.0, 2.0, 7, 16, pi / 3, :a
     fens, fes =
-        T3annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol)
+        T3annulus(rin, rex, nr, nc, Angl, orientation::Symbol)
     fens.xyz = xyz3(fens)
     filename = "mesh.stl"
     e = STLExporter(filename::AbstractString)
@@ -4519,10 +4519,10 @@ using FinEtools
 using FinEtools.MeshExportModule: VTK
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         100.0, 200.0, 13, 16, pi / 3, :a
     fens, fes =
-        T3annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol)
+        T3annulus(rin, rex, nr, nc, Angl, orientation::Symbol)
     d = (fens.xyz[:, 1] .^ 2 + fens.xyz[:, 2] .^ 2)
     File = "mesh.vtk"
     # Export of multiple scalar fields
@@ -4545,10 +4545,10 @@ using FinEtools
 using FinEtools.MeshExportModule: VTK
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         100.0, 200.0, 3, 6, pi / 3, :a
     fens, fes =
-        T3annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol)
+        T3annulus(rin, rex, nr, nc, Angl, orientation::Symbol)
     fens.xyz = xyz3(fens)
     d = (fens.xyz[:, 1] .^ 2 + fens.xyz[:, 2] .^ 2)
     File = "mesh.vtk"
@@ -4775,7 +4775,7 @@ module mtetvt1
 using FinEtools
 using FinEtools.MeshTetrahedronModule: tetv
 using Test
-function reftetv(X::FFltMat)
+function reftetv(X)
     local one6th = 1.0 / 6
     # @assert size(X, 1) == 4
     # @assert size(X, 2) == 3
@@ -4830,25 +4830,25 @@ using FinEtools
 using FinEtools.MeshExportModule: VTK
 using Test
 function test()
-    xradius::FFlt,
-    yradius::FFlt,
-    L::FFlt,
-    H::FFlt,
-    T::FFlt,
-    nL::FInt,
-    nH::FInt,
-    nW::FInt,
-    nT::FInt = 0.3, 0.45, 0.9, 0.8, 0.1, 5, 4, 3, 2
+    xradius,
+    yradius,
+    L,
+    H,
+    T,
+    nL,
+    nH,
+    nW,
+    nT = 0.3, 0.45, 0.9, 0.8, 0.1, 5, 4, 3, 2
     fens, fes = H8elliphole(
-        xradius::FFlt,
-        yradius::FFlt,
-        L::FFlt,
-        H::FFlt,
-        T::FFlt,
-        nL::FInt,
-        nH::FInt,
-        nW::FInt,
-        nT::FInt,
+        xradius,
+        yradius,
+        L,
+        H,
+        T,
+        nL,
+        nH,
+        nW,
+        nT,
     )
     # @show count(fens), count(fes)
     @test (count(fens), count(fes)) == (120, 54)
@@ -4932,9 +4932,9 @@ using FinEtools
 using FinEtools.MeshExportModule: VTK
 using Test
 function test()
-    rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt, orientation::Symbol =
+    rin, rex, nr, nc, Angl, orientation::Symbol =
         100.0, 200.0, 13, 16, pi / 3, :a
-    fens, fes = Q8annulus(rin::FFlt, rex::FFlt, nr::FInt, nc::FInt, Angl::FFlt)
+    fens, fes = Q8annulus(rin, rex, nr, nc, Angl)
     d = (fens.xyz[:, 1] .^ 2 + fens.xyz[:, 2] .^ 2)
     File = "mesh.vtk"
     # Export of multiple scalar fields
diff --git a/test/test_miscellaneous.jl b/test/test_miscellaneous.jl
index 839190d5..4da347be 100644
--- a/test/test_miscellaneous.jl
+++ b/test/test_miscellaneous.jl
@@ -850,7 +850,7 @@ function test()
     fens, fes = L2block(L, nl)
     fens.xyz = xyz3(fens)
     fens.xyz[2, 3] += L / 2
-    csmatout = zeros(FFlt, 3, 1)
+    csmatout = zeros(Float64, 3, 1)
     gradNparams = FESetModule.bfundpar(fes, vec([0.0]))
     J = transpose(fens.xyz) * gradNparams
     CSysModule.gen_iso_csmat!(csmatout, mean(fens.xyz, dims = 1), J, 0, 0)
@@ -879,7 +879,7 @@ function test()
     fens.xyz[2, 3] += L / 2
     File = "mesh.vtk"
     MeshExportModule.VTK.vtkexportmesh(File, fens, fes)
-    csmatout = zeros(FFlt, 3, 2)
+    csmatout = zeros(Float64, 3, 2)
     gradNparams = FESetModule.bfundpar(fes, vec([0.0 0.0]))
     J = zeros(3, 2)
     for i in eachindex(fes.conn[1])
@@ -1297,7 +1297,7 @@ function test()
         0.6324
         0.0975
     ]
-    x0 = fill(zero(FFlt), 6)
+    x0 = fill(zero(Float64), 6)
     maxiter = 20
     x = conjugategradient(A, b, x0, maxiter)
 end
@@ -2037,7 +2037,7 @@ function test()
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
     c = DataCache([10.0])
-    c(XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    c(XYZ, tangents, feid, qpid)
     @test c._cache == [10.0]
 end
 end
@@ -2051,11 +2051,11 @@ function test()
     XYZ = reshape([0.0, 0.0], 2, 1)
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
-    setvector!(v, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid) = begin
+    setvector!(v, XYZ, tangents, feid, qpid) = begin
         v .= [10.0]
     end
-    c = DataCache(FFlt[1], setvector!)
-    c(XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    c = DataCache(Float64[1], setvector!)
+    c(XYZ, tangents, feid, qpid)
     @test c._cache == [10.0]
 end
 end
@@ -2070,7 +2070,7 @@ function test()
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
     t = Ref(0.0)
-    setvector!(v, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid) = begin
+    setvector!(v, XYZ, tangents, feid, qpid) = begin
         if t[] < 5.0
             v .= [10.0]
         else
@@ -2081,11 +2081,11 @@ function test()
 
 
 
-    c = DataCache(FFlt[1], (cacheout, XYZ, tangents, feid, qpid) -> setvector!(cacheout, XYZ, tangents, feid, qpid))
-    c(XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    c = DataCache(Float64[1], (cacheout, XYZ, tangents, feid, qpid) -> setvector!(cacheout, XYZ, tangents, feid, qpid))
+    c(XYZ, tangents, feid, qpid)
     @test c._cache == [10.0]
     t[] = 6.0
-    c(XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    c(XYZ, tangents, feid, qpid)
     @test c._cache == [0.0]
 end
 end
@@ -2099,12 +2099,12 @@ function test()
     XYZ = reshape([0.0, 0.0], 2, 1)
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
-    setvector!(v, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid) = begin
+    setvector!(v, XYZ, tangents, feid, qpid) = begin
         v .= [10.0]
         return v
     end
-    c = DataCache(FFlt[1], setvector!)
-    c(XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    c = DataCache(Float64[1], setvector!)
+    c(XYZ, tangents, feid, qpid)
     @test c._cache == [10.0]
 end
 end
@@ -2120,7 +2120,7 @@ function test()
     feid = 0; qpid = 1
     vector = [10.0]
     fi = ForceIntensity(vector)
-    v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    v = updateforce!(fi, XYZ, tangents, feid, qpid)
     @test v == [10.0]
 end
 end
@@ -2134,12 +2134,12 @@ function test()
     XYZ = reshape([0.0, 0.0], 2, 1)
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
-    setvector!(v, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid) = begin
+    setvector!(v, XYZ, tangents, feid, qpid) = begin
         v .= [10.0]
     end
     vector = [10.0]
-    fi = ForceIntensity(FFlt, length(vector), setvector!)
-    v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    fi = ForceIntensity(Float64, length(vector), setvector!)
+    v = updateforce!(fi, XYZ, tangents, feid, qpid)
     @test v == [10.0]
 end
 end
@@ -2154,7 +2154,7 @@ function test()
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
     t = Ref(0.0)
-    setvector!(v, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid) = begin
+    setvector!(v, XYZ, tangents, feid, qpid) = begin
         if t[] < 5.0
             v .= [10.0]
         else
@@ -2164,12 +2164,12 @@ function test()
     end
     vector = [10.0]
     t[] = 0.0
-    fi = ForceIntensity(FFlt, length(vector), (out, XYZ, tangents, feid, qpid) -> setvector!(out, XYZ, tangents, feid, qpid))
-    v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    fi = ForceIntensity(Float64, length(vector), (out, XYZ, tangents, feid, qpid) -> setvector!(out, XYZ, tangents, feid, qpid))
+    v = updateforce!(fi, XYZ, tangents, feid, qpid)
     @test v == [10.0]
     t[] = 6.0
-    fi = ForceIntensity(FFlt, length(vector), (out, XYZ, tangents, feid, q) -> setvector!(out, XYZ, tangents, feid, qpid))
-    v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    fi = ForceIntensity(Float64, length(vector), (out, XYZ, tangents, feid, q) -> setvector!(out, XYZ, tangents, feid, qpid))
+    v = updateforce!(fi, XYZ, tangents, feid, qpid)
     @test v == [0.0]
 end
 end
@@ -2183,15 +2183,15 @@ function test()
     XYZ = reshape([0.0, 0.0], 2, 1)
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
-    setvector!(v, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid) = begin
+    setvector!(v, XYZ, tangents, feid, qpid) = begin
         v .= [10.0]
         return v
     end
     vector = [10.0]
-    fi = ForceIntensity(FFlt, length(vector), setvector!)
-    v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    fi = ForceIntensity(Float64, length(vector), setvector!)
+    v = updateforce!(fi, XYZ, tangents, feid, qpid)
     @test v == [10.0]
-    v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    v = updateforce!(fi, XYZ, tangents, feid, qpid)
     @test v == [10.0]
 end
 end
@@ -2207,7 +2207,7 @@ function test()
     feid = 0; qpid = 1
     vector = [10.0, -3.0]
     fi = SurfaceNormal(vector)
-    v = updatenormal!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    v = updatenormal!(fi, XYZ, tangents, feid, qpid)
     @test v == [10.0, -3.0]
 end
 end
@@ -2221,11 +2221,11 @@ function test()
     XYZ = reshape([0.0, 0.0], 2, 1)
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
-    setvector!(v, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid) = begin
+    setvector!(v, XYZ, tangents, feid, qpid) = begin
         v .= [10.0, -3.13]
     end
     fi = SurfaceNormal(2, setvector!)
-    v = updatenormal!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    v = updatenormal!(fi, XYZ, tangents, feid, qpid)
     @test v == [10.0, -3.13]
 end
 end
@@ -2240,7 +2240,7 @@ function test()
     tangents = reshape([-1.0, 1.0], 2, 1)
     feid = 0; qpid = 1
     fi = SurfaceNormal(2)
-    v = updatenormal!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    v = updatenormal!(fi, XYZ, tangents, feid, qpid)
     @test v == [0.7071067811865475, 0.7071067811865475]
 end
 end
@@ -2254,11 +2254,11 @@ function test()
     XYZ = reshape([0.0, 0.0], 2, 1)
     tangents = reshape([0.0, 1.0], 2, 1)
     feid = 0; qpid = 1
-    setvector!(v, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid) = begin
+    setvector!(v, XYZ, tangents, feid, qpid) = begin
         v .= [10.0, -3.13]
     end
-    fi = SurfaceNormal(2, zero(FFlt), setvector!)
-    v = updatenormal!(fi, XYZ::FFltMat, tangents::FFltMat, feid::FInt, qpid)
+    fi = SurfaceNormal(2, zero(Float64), setvector!)
+    v = updatenormal!(fi, XYZ, tangents, feid, qpid)
     @test v == [10.0, -3.13]
 end
 end
@@ -2625,25 +2625,29 @@ th = 7.0
 f(x) = cos(0.93 * pi * x[1] / Lx) + sin(1.7 * pi * x[2] / Ly)
 g(x) = -cos(0.93 * pi * x[1] / Lx) + 2 * sin(1.7 * pi * x[2] / Ly)
 
+mutable struct FEMMAdhocBase{ID<:IntegDomain} <: AbstractFEMM
+    integdomain::ID # domain data
+end
+
 function inspectintegpoints(
-    self::FEMM,
-    geom::NodalField{FFlt},
-    u::NodalField{FFlt},
-    dT::NodalField{FFlt},
-    felist::FIntVec,
+    self::FEMMAdhocBase,
+    geom::NodalField{Float64},
+    u::NodalField{Float64},
+    dT::NodalField{Float64},
+    felist,
     inspector::F,
     idat,
     quantity = :Cauchy;
     context...,
-) where {FEMM<:FEMMBase,F<:Function}
+) where {F<:Function}
     fes = self.integdomain.fes
     npts, Ns, gradNparams, w, pc = integrationdata(self.integdomain)
     nne = nodesperelem(fes)
     sdim = ndofs(geom)   # number of space dimensions
     mdim = manifdim(fes) # manifold dimension of the element
-    ecoords = fill(zero(FFlt), nne, ndofs(geom)) # array of Element coordinates
-    loc = fill(zero(FFlt), 1, sdim) # buffer
-    J = fill(zero(FFlt), sdim, mdim) # buffer
+    ecoords = fill(zero(Float64), nne, ndofs(geom)) # array of Element coordinates
+    loc = fill(zero(Float64), 1, sdim) # buffer
+    J = fill(zero(Float64), sdim, mdim) # buffer
     for ilist in eachindex(felist) # Loop over elements
         i = felist[ilist]
         gathervalues_asmat!(geom, ecoords, fes.conn[i])
@@ -2659,7 +2663,7 @@ end
 function test()
 
     fens, fes = Q4block(Lx, Ly, nx, ny) # Mesh
-    femm = FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
+    femm = FEMMAdhocBase(IntegDomain(fes, GaussRule(2, 2), th))
     geom = NodalField(fens.xyz)
     u = NodalField(0.0 .* fens.xyz)
     dT = NodalField(fill(0.0, count(fens), 1))
@@ -2706,25 +2710,29 @@ th = 7.0
 f(x) = cos(0.93 * pi * x[1] / Lx) + sin(1.7 * pi * x[2] / Ly)
 g(x) = -cos(0.93 * pi * x[1] / Lx) + 2 * sin(1.7 * pi * x[2] / Ly)
 
+mutable struct FEMMAdhocBase{ID<:IntegDomain} <: AbstractFEMM
+    integdomain::ID # domain data
+end
+
 function inspectintegpoints(
-    self::FEMM,
-    geom::NodalField{FFlt},
-    u::NodalField{FFlt},
-    dT::NodalField{FFlt},
-    felist::FIntVec,
+    self::FEMMAdhocBase,
+    geom::NodalField{Float64},
+    u::NodalField{Float64},
+    dT::NodalField{Float64},
+    felist,
     inspector::F,
     idat,
     quantity = :Cauchy;
     context...,
-) where {FEMM<:FEMMBase,F<:Function}
+) where {F<:Function}
     fes = self.integdomain.fes
     npts, Ns, gradNparams, w, pc = integrationdata(self.integdomain)
     nne = nodesperelem(fes)
     sdim = ndofs(geom)   # number of space dimensions
     mdim = manifdim(fes) # manifold dimension of the element
-    ecoords = fill(zero(FFlt), nne, ndofs(geom)) # array of Element coordinates
-    loc = fill(zero(FFlt), 1, sdim) # buffer
-    J = fill(zero(FFlt), sdim, mdim) # buffer
+    ecoords = fill(zero(Float64), nne, ndofs(geom)) # array of Element coordinates
+    loc = fill(zero(Float64), 1, sdim) # buffer
+    J = fill(zero(Float64), sdim, mdim) # buffer
     for ilist in eachindex(felist) # Loop over elements
         i = felist[ilist]
         gathervalues_asmat!(geom, ecoords, fes.conn[i])
@@ -2740,7 +2748,7 @@ end
 function test()
 
     fens, fes = Q4block(Lx, Ly, nx, ny) # Mesh
-    femm = FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
+    femm = FEMMAdhocBase(IntegDomain(fes, GaussRule(2, 2), th))
     geom = NodalField(fens.xyz)
     u = NodalField(0.0 .* fens.xyz)
     dT = NodalField(fill(0.0, count(fens), 1))
@@ -2779,14 +2787,14 @@ using Test
 function test()
     ndimensions = 2
     tangents, feid, qpid = reshape([1.0; 1.0], 2, 1), 0, 0
-    n = SurfaceNormal(ndimensions::FInt)
-    normal = updatenormal!(n, [0.0 0.0 0.0], tangents::FFltMat, feid::FInt, qpid)
+    n = SurfaceNormal(ndimensions)
+    normal = updatenormal!(n, [0.0 0.0 0.0], tangents, feid, qpid)
     @test norm(normal - [0.7071067811865475, -0.7071067811865475]) <= 1.0e-5
 
     ndimensions = 3
     tangents, feid = reshape([1.0 0.0; 1.0 0.0; 0.0 1.0], 3, 2), 0
-    n = SurfaceNormal(ndimensions::FInt)
-    normal = updatenormal!(n, [0.0 0.0 0.0], tangents::FFltMat, feid::FInt, qpid)
+    n = SurfaceNormal(ndimensions)
+    normal = updatenormal!(n, [0.0 0.0 0.0], tangents, feid, qpid)
     @test norm(normal - [0.7071067811865475, -0.7071067811865475, 0.0]) <= 1.0e-5
     true
 end
@@ -2811,25 +2819,29 @@ th = 7.0
 f(x) = cos(0.93 * pi * x[1] / Lx) + sin(1.7 * pi * x[2] / Ly)
 g(x) = -cos(0.93 * pi * x[1] / Lx) + 2 * sin(1.7 * pi * x[2] / Ly)
 
+mutable struct FEMMAdhocBase{ID<:IntegDomain} <: AbstractFEMM
+    integdomain::ID # domain data
+end
+
 function inspectintegpoints(
-    self::FEMM,
-    geom::NodalField{FFlt},
-    u::NodalField{FFlt},
-    dT::NodalField{FFlt},
-    felist::FIntVec,
+    self::FEMMAdhocBase,
+    geom::NodalField{Float64},
+    u::NodalField{Float64},
+    dT::NodalField{Float64},
+    felist,
     inspector::F,
     idat,
     quantity = :Cauchy;
     context...,
-) where {FEMM<:FEMMBase,F<:Function}
+) where {F<:Function}
     fes = self.integdomain.fes
     npts, Ns, gradNparams, w, pc = integrationdata(self.integdomain)
     nne = nodesperelem(fes)
     sdim = ndofs(geom)   # number of space dimensions
     mdim = manifdim(fes) # manifold dimension of the element
-    ecoords = fill(zero(FFlt), nne, ndofs(geom)) # array of Element coordinates
-    loc = fill(zero(FFlt), 1, sdim) # buffer
-    J = fill(zero(FFlt), sdim, mdim) # buffer
+    ecoords = fill(zero(Float64), nne, ndofs(geom)) # array of Element coordinates
+    loc = fill(zero(Float64), 1, sdim) # buffer
+    J = fill(zero(Float64), sdim, mdim) # buffer
     for ilist in eachindex(felist) # Loop over elements
         i = felist[ilist]
         gathervalues_asmat!(geom, ecoords, fes.conn[i])
@@ -2845,7 +2857,7 @@ end
 function test()
 
     fens, fes = Q4block(Lx, Ly, nx, ny) # Mesh
-    femm = FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
+    femm = FEMMAdhocBase(IntegDomain(fes, GaussRule(2, 2), th))
     geom = NodalField(fens.xyz)
     u = NodalField(0.0 .* fens.xyz)
     dT = NodalField(fill(0.0, count(fens), 1))
@@ -2895,13 +2907,14 @@ using DelimitedFiles
 using SparseArrays
 using FinEtools
 using FinEtools.MatrixUtilityModule: export_sparse, import_sparse
+using LinearAlgebra
 using Test
 function test()
     N = 10
     A = sprand(N, N, 0.1)
     export_sparse("A.txt", A)
     B = import_sparse("A.txt")
-    @test (A - B) == spzeros(size(B)...)
+    @test maximum(abs.(A[:] - B[:]))  < 1.0e-5
     try
         rm("A.txt")
     catch
@@ -3135,3 +3148,24 @@ function test()
 end
 test()
 end
+
+# struct T
+# end
+
+# f(a::T) = "Main"
+
+# module M1
+# import Main: T, f
+# f(a::T) = "M1"
+# end
+
+# module M2
+# import Main: T, f
+# f(a::T) = "M2"
+# end
+
+# using .M1
+# using .M2
+
+# M1.f(T())
+# M2.f(T())