added tests, reconstruction operators can now also used with lazy_interpolate!

Author: chmerdon

src/reconstructionhandlers.jl

@@ -19,11 +19,11 @@ Weighted reconstruction operator: evaluates a reconstructed version of a finite
 # Parameters
 - `FETypeR`: The reconstruction finite element space type (target space for reconstruction).
 - `O`: The standard function operator to be evaluated (e.g., identity, gradient, etc.).
-- `F`: The type of the weight function (should be callable, e.g., a function or functor).
+- `w`: The type of the weight function (should be callable, e.g., a function or functor of type w(x)).
 """
-abstract type WeightedReconstruct{FETypeR, O, F} <: Reconstruct{FETypeR, O} end
+abstract type WeightedReconstruct{FETypeR, O, w} <: Reconstruct{FETypeR, O} end
 
-weight_type(::Type{<:WeightedReconstruct{FETypeR, O, F}}) where {FETypeR, O, F} = F
+weight_type(::Type{<:WeightedReconstruct{FETypeR, O, w}}) where {FETypeR, O, w} = w
 
 
 ################## SPECIAL INTERPOLATORS ####################
@@ -474,7 +474,7 @@ boundary_coefficients!(coefficients, RH::ReconstructionHandler, cell)
 
 generates the coefficients for the facial dofs of the reconstruction operator on the cell
 """
-function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, <: Reconstruct, <:H1CR{ncomponents}, <:HDIVRT0{ncomponents}, <:ON_CELLS, EG}, cell) where {Tv, Ti, ncomponents, EG}
+function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, <:Reconstruct, <:H1CR{ncomponents}, <:HDIVRT0{ncomponents}, <:ON_CELLS, EG}, cell) where {Tv, Ti, ncomponents, EG}
     xFaceVolumes = RH.xFaceVolumes
     xFaceNormals = RH.xFaceNormals
     xCellFaces = RH.xCellFaces
@@ -709,4 +709,3 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti,
     end
     return nothing
 end
-

src/reconstructionoperators.jl

@@ -4,7 +4,7 @@ ReconstructionSpace(::Type{<:ReconstructionOperator{FETypeR, O}}) where {FETypeR
 NeededDerivative4Operator(::Type{<:ReconstructionOperator{FETypeR, O}}) where {FETypeR, O} = NeededDerivative4Operator(O)
 Length4Operator(::Type{<:ReconstructionOperator{FETypeR, O}}, xdim, nc) where {FETypeR, O} = Length4Operator(O, xdim, nc)
 DefaultName4Operator(::Type{Reconstruct{FETypeR, O}}) where {FETypeR, O} = "R(" * DefaultName4Operator(O) * ")"
-DefaultName4Operator(::Type{WeightedReconstruct{FETypeR, O}}) where {FETypeR, O} = "R(r" * DefaultName4Operator(O) * ")"
+DefaultName4Operator(::Type{WeightedReconstruct{FETypeR, O, w}}) where {FETypeR, O, w} = "wR(r" * DefaultName4Operator(O) * ")"
 
 struct FEReconstEvaluator{T, TvG, TiG, FEType, FEType2, stdop, RH} <: FEEvaluator{T, TvG, TiG}
     citem::Base.RefValue{Int}                       # current item
@@ -66,6 +66,10 @@ function FEEvaluator(
     return FEReconstEvaluator{T, TvG, TiG, FEType, FETypeReconst, stdop, typeof(reconst_handler)}(FEB.citem, FE, FEB, cvals, coefficients2, reconst_handler)
 end
 
+function relocate_xref!(FEB::FEReconstEvaluator, new_xref)
+    return relocate_xref!(FEB.FEB, new_xref)
+end
+
 function update_basis!(FEBE::FEReconstEvaluator)
     ## evaluate standard operator in reconstruction basis (->cvals_reconst)
     update_basis!(FEBE.FEB)

test/test_operators.jl

@@ -8,9 +8,58 @@ function run_operator_tests()
         @test error < 1.0e-14
         error = test_derivatives3D()
         @test error < 1.0e-14
+        test_reconstructions()
     end
 end
 
+function test_reconstructions()
+    ## divergence-free axisymmetric velocity field u(r,z) = (r,-2z) in cylindrical coordinates
+    function u!(result, qpinfo)
+        x = qpinfo.x
+        result[1] = x[1]
+        return result[2] = -2 * x[2]
+    end
+
+    ## vector field times r, should have zero divergence div(ru) = (d/dr, d/dz) \cdot (ru) = 0
+    function ru!(result, qpinfo)
+        x = qpinfo.x
+        result[1] = x[1]^2
+        return result[2] = -2 * x[1] * x[2]
+    end
+
+    xgrid = testgrid(Triangle2D)
+
+    ## interpolate u into H1BR{2} (inf-sup stable Stokes element)
+    FES = FESpace{H1BR{2}}(xgrid)
+    uh = FEVector(FES)
+    interpolate!(uh[1], u!; bonus_quadorder = 2)
+
+    for FETypeR in [HDIVRT0{2}, HDIVBDM1{2}]
+        ## interpolate ru into HDIVRTO{2} or HDIVBDM1{2}
+        FES2 = FESpace{FETypeR}(xgrid)
+        Πur = FEVector(FES2)
+        interpolate!(Πur[1], ru!; bonus_quadorder = 2)
+
+        ## test if interpolate of ru is divergence-free by interpolating into P0 function and checking its coefficients
+        FES3 = FESpace{L2P0{1}}(xgrid)
+        divΠur = FEVector(FES3)
+        lazy_interpolate!(divΠur[1], Πur, [(1, Divergence)])
+        @test sqrt(sum((divΠur.entries .^ 2))) < 1.0e-14
+
+        ## test if r-weighted reconstruction of uh is divergence-free by interpolating into P0 function and checking its coefficients
+        weight = (x) -> (x[1])
+        lazy_interpolate!(divΠur[1], uh, [(1, WeightedReconstruct{FETypeR, Divergence, typeof(weight)})])
+        @test sqrt(sum((divΠur.entries .^ 2))) < 1.0e-14
+
+        ## test if weighted reconstruction of uh and interpolation of ru are identical
+        FES4 = FESpace{L2P1{1}}(xgrid)
+        diff = FEVector(FES4)
+        lazy_interpolate!(diff[1], [uh[1], Πur[1]], [(1, WeightedReconstruct{FETypeR, Identity, typeof(weight)}), (2, Identity)]; postprocess = (result, args, qpinfo) -> (result[1] = (args[1] - args[3])^2 + (args[2] - args[4])^2))
+        @test sqrt(sum((diff.entries))) < 1.0e-14
+    end
+    return
+end
+
 function test_derivatives2D()
     ## define test function and expected operator evals
     function testf(result, qpinfo)