changes suggested by pjaap's review

Author: chmerdon

CHANGELOG.md

@@ -1,6 +1,11 @@
 # CHANGES
 
-## v1.0.0
+## v1.1.0 May 15, 2025
+  - overhaul of interpolators: each finite element
+    interpolation is composed of the three new
+    structures NodalInterpolator, MomentInterpolator, FunctionalInterpolator
+  
+## v1.0.0 April 7, 2025
   - doc improvements
   - several bugfixes regarding finite elements defined on subgrids
   - integrate now reacts on regions parameter and has proper docstrings

src/interpolators.jl

@@ -32,6 +32,7 @@ function NodalInterpolator(
     broken = FES.broken,
     component_offset = FES.coffset,
     kwargs...) where {T}
+
     FEType = eltype(FES)
     ncomponents = get_ncomponents(FEType)
     xCoordinates = grid[Coordinates]
@@ -115,7 +116,20 @@ are solved with a mass matrix of interior basis functions (evaluated with operat
 basis functions of the moments of type FEType_ref and with moments_operator. In the bestapprox mode
 the mass matrix instead is formed from the scalar product of the interior basis functions.
 """
-function MomentInterpolator(FE::FESpace{Tv, Ti, FEType, APT}, AT::Type{<:AssemblyType}, xgrid = FE.dofgrid; operator = Identity, FEType_ref = :auto, FEType_moments = :auto, moments_operator = operator, moments_dofs = Int[], bestapprox = false, order = 0, coffset::Int = -1, componentwise = true, kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
+function MomentInterpolator(
+    FE::FESpace{Tv, Ti, FEType, APT},
+    AT::Type{<:AssemblyType},
+    xgrid = FE.dofgrid;
+    operator = Identity,
+    FEType_ref = :auto,
+    FEType_moments = :auto,
+    moments_operator = operator,
+    moments_dofs = Int[],
+    bestapprox = false,
+    order = 0,
+    coffset::Int = -1,
+    componentwise = true,
+    kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
 
     itemvolumes = xgrid[GridComponentVolumes4AssemblyType(AT)]
     itemnodes = xgrid[GridComponentNodes4AssemblyType(AT)]
@@ -174,7 +188,10 @@ function MomentInterpolator(FE::FESpace{Tv, Ti, FEType, APT}, AT::Type{<:Assembl
     end
 
     ## check if mass matrix can be computed once or needs to be recomputed on every mesh
-    if FEType_ref <: AbstractH1FiniteElement && !(FEType_ref <: AbstractH1FiniteElementWithCoefficients) && FEType_moments <: AbstractH1FiniteElement && moments_operator == Identity
+    if FEType_ref <: AbstractH1FiniteElement &&
+        !(FEType_ref <: AbstractH1FiniteElementWithCoefficients) &&
+        FEType_moments <: AbstractH1FiniteElement &&
+         moments_operator == Identity
         fixed_mass_matrix = true
     else
         fixed_mass_matrix = false
@@ -444,7 +461,11 @@ function FunctionalInterpolator(
     AT::Type{<:AssemblyType} = ON_FACES,
     xgrid = FE.dofgrid;
     bonus_quadorder = 0,
-    operator = NormalFlux, nfluxes = 0, dofs = [], mean = false, kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
+    operator = NormalFlux,
+    nfluxes = 0,
+    dofs = [],
+    mean = false,
+    kwargs...) where {Tv, Ti, FEType <: AbstractFiniteElement, APT}
 
     itemvolumes = xgrid[GridComponentVolumes4AssemblyType(AT)]
     itemnodes = xgrid[GridComponentNodes4AssemblyType(AT)]

test/test_interpolators.jl

@@ -2,62 +2,59 @@ function run_interpolator_tests()
 
     # list of FETypes that should be tested
     TestCatalog1D = [
-        L2P0{1},
-        H1P1{1},
-        H1P2{1, 1},
-        H1P3{1, 1},
-        H1Pk{1, 1, 3},
-        H1Pk{1, 1, 4},
-        H1Pk{1, 1, 5},
+        L2P0{1} => 0,
+        H1P1{1} => 1,
+        H1P2{1, 1} => 2,
+        H1P3{1, 1} => 3,
+        H1Pk{1, 1, 3} => 3,
+        H1Pk{1, 1, 4} => 4,
+        H1Pk{1, 1, 5} => 5
     ]
-    ExpectedOrders1D = [0, 1, 2, 3, 3, 4, 5]
 
     TestCatalog2D = [
-        HCURLN0{2},
-        HCURLN1{2},
-        HDIVRT0{2},
-        HDIVRTk{2, 0},
-        HDIVBDM1{2},
-        HDIVRT1{2},
-        HDIVRTk{2, 1},
-        HDIVBDM2{2},
-        HDIVRTk{2, 2},
-        HDIVRTk{2, 3},
-        HDIVRTk{2, 4},
-        L2P0{2},
-        L2P1{2},
-        H1P1{2},
-        H1Q1{2},
-        H1CR{2},
-        H1MINI{2, 2},
-        H1P1TEB{2},
-        H1BR{2},
-        H1P2{2, 2},
-        H1P2B{2, 2},
-        H1Q2{2, 2},
-        H1P3{2, 2},
-        H1Pk{2, 2, 3},
-        H1Pk{2, 2, 4},
-        H1Pk{2, 2, 5},
+        HCURLN0{2} => 0,
+        HCURLN1{2} => 1,
+        HDIVRT0{2} => 0,
+        HDIVRTk{2, 0} => 0,
+        HDIVBDM1{2} => 1,
+        HDIVRT1{2} => 1,
+        HDIVRTk{2, 1} => 1,
+        HDIVBDM2{2} => 2,
+        HDIVRTk{2, 2} => 2,
+        HDIVRTk{2, 3} => 3,
+        HDIVRTk{2, 4} => 4,
+        L2P0{2} => 0,
+        L2P1{2} => 1,
+        H1P1{2} => 1,
+        H1Q1{2} => 1,
+        H1CR{2} => 1,
+        H1MINI{2, 2} => 1,
+        H1P1TEB{2} => 1,
+        H1BR{2} => 1,
+        H1P2{2, 2} => 2,
+        H1P2B{2, 2} => 2,
+        H1Q2{2, 2} => 2,
+        H1P3{2, 2} => 3,
+        H1Pk{2, 2, 3} => 3,
+        H1Pk{2, 2, 4} => 4,
+        H1Pk{2, 2, 5} => 5
     ]
-    ExpectedOrders2D = [0, 1, 0, 0, 1, 1, 1, 2, 2, 3, 4, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 4, 5]
 
     TestCatalog3D = [
-        HCURLN0{3},
-        HDIVRT0{3},
-        HDIVBDM1{3},
-        HDIVRT1{3},
-        L2P0{3},
-        H1P1{3},
-        H1Q1{3},
-        H1CR{3},
-        H1MINI{3, 3},
-        H1P1TEB{3},
-        H1BR{3},
-        H1P2{3, 3},
-        H1P3{3, 3},
+        HCURLN0{3} => 0,
+        HDIVRT0{3} => 0,
+        HDIVBDM1{3} => 1,
+        HDIVRT1{3} => 1,
+        L2P0{3} => 0,
+        H1P1{3} => 1,
+        H1Q1{3} => 1,
+        H1CR{3} => 1,
+        H1MINI{3, 3} => 1,
+        H1P1TEB{3} => 1,
+        H1BR{3} => 1,
+        H1P2{3, 3} => 2,
+        H1P3{3, 3} => 3
     ]
-    ExpectedOrders3D = [0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 2, 3]
 
     ## function that computes errors at enough quadrature points for polynomial of degree order
     function compute_error(uh::FEVectorBlock, u::Function, order = get_polynomialorder(get_FEType(uh), uh.FES.xgrid[CellGeometries][1]))
@@ -123,8 +120,8 @@ function run_interpolator_tests()
         println("============================")
         xgrid = testgrid(Edge1D)
         for n in 1:length(TestCatalog1D)
-            test_interpolation(xgrid, TestCatalog1D[n], ExpectedOrders1D[n])
-            test_interpolation(xgrid, TestCatalog1D[n], ExpectedOrders1D[n], true)
+            test_interpolation(xgrid, TestCatalog1D[n].first, TestCatalog1D[n].second)
+            test_interpolation(xgrid, TestCatalog1D[n].first, TestCatalog1D[n].second, true)
         end
         println("\n")
         println("============================")
@@ -134,8 +131,8 @@ function run_interpolator_tests()
             xgrid = uniform_refine(reference_domain(EG), 1)
             println("EG = $EG")
             for n in 1:length(TestCatalog2D), broken in (false, true)
-                if ExtendableFEMBase.isdefined(TestCatalog2D[n], EG, broken)
-                    test_interpolation(xgrid, TestCatalog2D[n], ExpectedOrders2D[n], broken)
+                if ExtendableFEMBase.isdefined(TestCatalog2D[n].first, EG, broken)
+                    test_interpolation(xgrid, TestCatalog2D[n].first, TestCatalog2D[n].second, broken)
                 else
                     @warn "$(TestCatalog2D[n]) (broken = $broken) not defined on $EG (skipping test case)"
                 end
@@ -148,8 +145,8 @@ function run_interpolator_tests()
         for EG in [Tetrahedron3D, Parallelepiped3D]
             xgrid = uniform_refine(reference_domain(EG), 1)
             for n in 1:length(TestCatalog3D), broken in (false, true)
-                if ExtendableFEMBase.isdefined(TestCatalog3D[n], EG, broken)
-                    test_interpolation(xgrid, TestCatalog3D[n], ExpectedOrders3D[n], broken)
+                if ExtendableFEMBase.isdefined(TestCatalog3D[n].first, EG, broken)
+                    test_interpolation(xgrid, TestCatalog3D[n].first, TestCatalog3D[n].second, broken)
                 else
                     @warn "$(TestCatalog3D[n]) (broken = $broken) not defined on $EG (skipping test case)"
                 end