fixes and weight function for reconstruction operators (#47)

Author: chmerdon

CHANGELOG.md

@@ -1,5 +1,9 @@
 # CHANGES
 
+## v1.4.0 July 17, 2025
+  - added new weighted reconstruction operator `WeightedReconstruct` (so far only for H1BR{2})
+  - fixed H1CR{2} reconstruction
+
 ## v1.3.0 July 09, 2025
   - some bugfixes related to new template parameter from 1.2.0
   - @show of FEVectorBlock does not crash anymore

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
 authors = ["Christian Merdon <[email protected]>", "Patrick Jaap <[email protected]>"]
-version = "1.3.0"
+version = "1.4.0"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"

src/ExtendableFEMBase.jl

@@ -127,12 +127,12 @@ export displace_mesh, displace_mesh!
 
 include("reconstructionhandlers.jl")
 export ReconstructionHandler, get_rcoefficients!
+export Reconstruct, WeightedReconstruct
 
 include("feevaluator.jl")
 export FEEvaluator, update_basis!, eval_febe!
 
 include("reconstructionoperators.jl")
-export Reconstruct
 
 include("accumvector.jl")
 export AccumulatingVector

src/reconstructionhandlers.jl

@@ -1,3 +1,33 @@
+abstract type ReconstructionOperator{FETypeR, O} <: AbstractFunctionOperator end
+
+"""
+$(TYPEDEF)
+
+reconstruction operator: evaluates a reconstructed version of the finite element function.
+
+FETypeR specifies the reconstruction space (needs to be defined for the finite element that it is applied to).
+O specifies the StandardFunctionOperator that shall be evaluated.
+"""
+abstract type Reconstruct{FETypeR, O} <: ReconstructionOperator{FETypeR, O} end
+
+
+"""
+WeightedReconstruct{FETypeR, O, w} <: Reconstruct{FETypeR, O}
+
+Weighted reconstruction operator:
+evaluates a reconstructed version of a finite element function, multiplied by a weight function.
+**Warning**: This is a prototype and currently only works for the HDIVRT0{2} and HDIVBDM1{2} reconstruction of H1BR{2}.
+
+# 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.).
+- `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, w} <: Reconstruct{FETypeR, O} end
+
+weight_type(::Type{<:WeightedReconstruct{FETypeR, O, w}}) where {FETypeR, O, w} = w
+
+
 ################## SPECIAL INTERPOLATORS ####################
 
 """
@@ -69,18 +99,33 @@ abstract type ReconstructionDofs{FE1, FE2, AT} <: AbstractGridIntegerArray2D end
 """
 $(TYPEDEF)
 
+abstract grid component type that can be used to funnel reconstruction weights into the operator
+(default is 1)
+"""
+abstract type ReconstructionWeights{AT} <: AbstractGridFloatArray2D end
+
+"""
+$(TYPEDEF)
+
 struct for storing information needed to evaluate a reconstruction operator
 """
-struct ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}
+struct ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG, F}
     FES::FESpace{Tv, Ti, FE1, ON_CELLS}
     FER::FESpace{Tv, Ti, FE2, ON_CELLS}
+    xCoordinates::Array{Tv, 2}
+    xFaceNodes::Adjacency{Ti}
     xFaceVolumes::Array{Tv, 1}
     xFaceNormals::Array{Tv, 2}
     xCellFaceOrientations::Adjacency{Ti}
     xCellFaces::Adjacency{Ti}
     interior_offset::Int
     interior_ndofs::Int
     interior_coefficients::Matrix{Tv} # coefficients for interior basis functions are precomputed
+    weight::F
+end
+
+function default_weight_function(x)
+    return 1
 end
 
 """
@@ -90,7 +135,7 @@ generates a reconstruction handler
 returns the local coefficients need to evaluate a reconstruction operator
 of one finite element space into another
 """
-function ReconstructionHandler(FES::FESpace{Tv, Ti, FE1, APT}, FES_Reconst::FESpace{Tv, Ti, FE2, APT}, AT, EG) where {Tv, Ti, FE1, FE2, APT}
+function ReconstructionHandler(FES::FESpace{Tv, Ti, FE1, APT}, FES_Reconst::FESpace{Tv, Ti, FE2, APT}, AT, EG, RT, weight = default_weight_function) where {Tv, Ti, FE1, FE2, APT}
     xgrid = FES.xgrid
     interior_offset = interior_dofs_offset(AT, FE2, EG)
     interior_ndofs = get_ndofs(AT, FE2, EG) - interior_offset
@@ -100,11 +145,14 @@ function ReconstructionHandler(FES::FESpace{Tv, Ti, FE1, APT}, FES_Reconst::FESp
         interior_ndofs = 0
         coeffs = zeros(Tv, 0, 0)
     end
+
     xFaceVolumes = xgrid[FaceVolumes]
     xFaceNormals = xgrid[FaceNormals]
+    xCoordinates = xgrid[Coordinates]
+    xFaceNodes = xgrid[FaceNodes]
     xCellFaceOrientations = dim_element(EG) == 2 ? xgrid[CellFaceSigns] : xgrid[CellFaceOrientations]
     xCellFaces = xgrid[CellFaces]
-    return ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}(FES, FES_Reconst, xFaceVolumes, xFaceNormals, xCellFaceOrientations, xCellFaces, interior_offset, interior_ndofs, coeffs)
+    return ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG, typeof(weight)}(FES, FES_Reconst, xCoordinates, xFaceNodes, xFaceVolumes, xFaceNormals, xCellFaceOrientations, xCellFaces, interior_offset, interior_ndofs, coeffs, weight)
 end
 
 """
@@ -113,7 +161,7 @@ $(TYPEDSIGNATURES)
 caller function to extract the coefficients of the reconstruction
 on an item
 """
-function get_rcoefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}, item) where {Tv, Ti, FE1, FE2, AT, EG}
+function get_rcoefficients!(coefficients, RH::ReconstructionHandler, item)
     boundary_coefficients!(coefficients, RH, item)
     for dof in 1:size(coefficients, 1), k in 1:RH.interior_ndofs
         coefficients[dof, RH.interior_offset + k] = RH.interior_coefficients[(dof - 1) * RH.interior_ndofs + k, item]
@@ -428,7 +476,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, <: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
@@ -446,59 +494,93 @@ end
 
 const _P1_INTO_BDM1_COEFFS = [-1 // 12, 1 // 12]
 
-function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}, cell) where {Tv, Ti, FE1 <: H1BR{2}, FE2 <: HDIVBDM1{2}, AT <: ON_CELLS, EG <: Union{Triangle2D, Quadrilateral2D}}
+function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: Reconstruct, FE1 <: H1BR{2}, FE2 <: HDIVBDM1{2}, AT <: ON_CELLS, EG <: Union{Triangle2D, Quadrilateral2D}}
     xFaceVolumes = RH.xFaceVolumes
     xFaceNormals = RH.xFaceNormals
     xCellFaceSigns = RH.xCellFaceOrientations
     xCellFaces = RH.xCellFaces
+    xFaceNodes = RH.xFaceNodes
+    xCoordinates = RH.xCoordinates
     face_rule = local_cellfacenodes(EG)
     nnodes = size(face_rule, 1)
     nfaces = size(face_rule, 2)
     node = 0
     face = 0
     BDM1_coeffs = _P1_INTO_BDM1_COEFFS
+    weight = RH.weight
+    xmid = zeros(Tv, 2)
+    w = ones(Tv, 2)
     for f in 1:nfaces
         face = xCellFaces[f, cell]
+        x1 = view(xCoordinates, :, xFaceNodes[1, face])
+        x2 = view(xCoordinates, :, xFaceNodes[2, face])
+        xmid .= 0.5 * (x1 .+ x2)
+        if xCellFaceSigns[f, cell] == -1
+            w[1] = weight(x2)
+            w[2] = weight(x1)
+        else
+            w[1] = weight(x1)
+            w[2] = weight(x2)
+        end
+        wmid = weight(xmid)
         for n in 1:nnodes
             node = face_rule[n, f]
             for k in 1:2
                 # RT0 reconstruction coefficients for P1 functions on reference element
-                coefficients[nfaces * (k - 1) + node, 2 * (f - 1) + 1] = 1 // 2 * xFaceVolumes[face] * xFaceNormals[k, face]
+                coefficients[nfaces * (k - 1) + node, 2 * (f - 1) + 1] = xFaceVolumes[face] * (1 // 6 * w[n] + 1 // 3 * wmid) * xFaceNormals[k, face]
                 # BDM1 reconstruction coefficients for P1 functions on reference element
-                coefficients[nfaces * (k - 1) + node, 2 * (f - 1) + 2] = BDM1_coeffs[n] * xFaceVolumes[face] * xFaceNormals[k, face] * xCellFaceSigns[f, cell]
+                coefficients[nfaces * (k - 1) + node, 2 * (f - 1) + 2] = xFaceVolumes[face] * xFaceNormals[k, face] * xCellFaceSigns[f, cell] * (BDM1_coeffs[n] * w[n])
             end
         end
         # RT0 reconstruction coefficients for face bubbles on reference element
-        coefficients[nfaces * 2 + f, 2 * (f - 1) + 1] = xFaceVolumes[face]
+        coefficients[nfaces * 2 + f, 2 * (f - 1) + 1] = xFaceVolumes[face] * wmid
     end
+
     return nothing
 end
 
-function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}, cell) where {Tv, Ti, FE1 <: H1BR{2}, FE2 <: HDIVRT0{2}, AT <: ON_CELLS, EG <: Union{Triangle2D, Quadrilateral2D}}
+function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: Reconstruct, FE1 <: H1BR{2}, FE2 <: HDIVRT0{2}, AT <: ON_CELLS, EG <: Union{Triangle2D, Quadrilateral2D}}
     xFaceVolumes = RH.xFaceVolumes
     xFaceNormals = RH.xFaceNormals
     xCellFaces = RH.xCellFaces
+    xCellFaceSigns = RH.xCellFaceOrientations
+    xFaceNodes = RH.xFaceNodes
+    xCoordinates = RH.xCoordinates
     face_rule = local_cellfacenodes(EG)
     nnodes = size(face_rule, 1)
     nfaces = size(face_rule, 2)
     node = 0
     face = 0
+    weight = RH.weight
+    xmid = zeros(Tv, 2)
+    w = ones(Tv, 2)
     for f in 1:nfaces
         face = xCellFaces[f, cell]
+        x1 = view(xCoordinates, :, xFaceNodes[1, face])
+        x2 = view(xCoordinates, :, xFaceNodes[2, face])
+        xmid .= 0.5 * (x1 .+ x2)
+        if xCellFaceSigns[f, cell] == -1
+            w[1] = weight(x2)
+            w[2] = weight(x1)
+        else
+            w[1] = weight(x1)
+            w[2] = weight(x2)
+        end
+        wmid = weight(xmid)
         # reconstruction coefficients for P1 functions on reference element
         for n in 1:nnodes
             node = face_rule[n, f]
             for k in 1:2
-                coefficients[nfaces * (k - 1) + node, f] = 1 // 2 * xFaceVolumes[face] * xFaceNormals[k, face]
+                coefficients[nfaces * (k - 1) + node, f] = xFaceVolumes[face] * (1 // 6 * w[n] + 1 // 3 * wmid) * xFaceNormals[k, face]
             end
         end
         # reconstruction coefficients for face bubbles on reference element
-        coefficients[2 * nfaces + f, f] = xFaceVolumes[face]
+        coefficients[2 * nfaces + f, f] = xFaceVolumes[face] * wmid
     end
     return nothing
 end
 
-function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}, cell) where {Tv, Ti, FE1 <: H1P2B{2, 2}, FE2 <: HDIVRT1{2}, AT <: ON_CELLS, EG <: Triangle2D}
+function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: Reconstruct, FE1 <: H1P2B{2, 2}, FE2 <: HDIVRT1{2}, AT <: ON_CELLS, EG <: Triangle2D}
     xFaceVolumes = RH.xFaceVolumes
     xFaceNormals = RH.xFaceNormals
     xCellFaceSigns = RH.xCellFaceOrientations
@@ -529,7 +611,7 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti,
     return nothing
 end
 
-function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}, cell) where {Tv, Ti, FE1 <: H1P2B{2, 2}, FE2 <: HDIVBDM2{2}, AT <: ON_CELLS, EG <: Triangle2D}
+function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: Reconstruct, FE1 <: H1P2B{2, 2}, FE2 <: HDIVBDM2{2}, AT <: ON_CELLS, EG <: Triangle2D}
     xFaceVolumes = RH.xFaceVolumes
     xFaceNormals = RH.xFaceNormals
     xCellFaceSigns = RH.xCellFaceOrientations
@@ -567,7 +649,7 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti,
 end
 
 
-function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}, cell) where {Tv, Ti, FE1 <: H1BR{3}, FE2 <: HDIVRT0{3}, AT <: ON_CELLS, EG <: Tetrahedron3D}
+function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: Reconstruct, FE1 <: H1BR{3}, FE2 <: HDIVRT0{3}, AT <: ON_CELLS, EG <: Tetrahedron3D}
     xFaceVolumes = RH.xFaceVolumes
     xFaceNormals = RH.xFaceNormals
     xCellFaces = RH.xCellFaces
@@ -592,7 +674,7 @@ end
 
 const _P1_INTO_BDM1_COEFFS_3D = [-1 // 36 -1 // 36 1 // 18; -1 // 36 1 // 18 -1 // 36; 1 // 18 -1 // 36 -1 // 36]
 
-function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, FE1, FE2, AT, EG}, cell) where {Tv, Ti, FE1 <: H1BR{3}, FE2 <: HDIVBDM1{3}, AT <: ON_CELLS, EG <: Tetrahedron3D}
+function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: Reconstruct, FE1 <: H1BR{3}, FE2 <: HDIVBDM1{3}, AT <: ON_CELLS, EG <: Tetrahedron3D}
     xFaceVolumes = RH.xFaceVolumes
     xFaceNormals = RH.xFaceNormals
     xCellFaces = RH.xCellFaces

src/reconstructionoperators.jl

@@ -1,36 +1,24 @@
-abstract type ReconstructionOperator <: AbstractFunctionOperator end
+StandardFunctionOperator(::Type{<:ReconstructionOperator{FETypeR, O}}) where {FETypeR, O} = O
+ReconstructionSpace(::Type{<:ReconstructionOperator{FETypeR, O}}) where {FETypeR, O} = FETypeR
 
-"""
-$(TYPEDEF)
-
-reconstruction operator: evaluates a reconstructed version of the finite element function.
-
-FETypeR specifies the reconstruction space (needs to be defined for the finite element that it is applied to).
-O specifies the StandardFunctionOperator that shall be evaluated.
-"""
-abstract type Reconstruct{FETypeR, O} <: ReconstructionOperator where {FETypeR <: AbstractFiniteElement, O <: StandardFunctionOperator} end # calculates the jump between both sides of the face
-
-
-StandardFunctionOperator(::Type{Reconstruct{FETypeR, O}}) where {FETypeR, O} = O
-ReconstructionSpace(::Type{Reconstruct{FETypeR, O}}) where {FETypeR, O} = FETypeR
-
-NeededDerivative4Operator(::Type{Reconstruct{FETypeR, O}}) where {FETypeR, O} = NeededDerivative4Operator(O)
-Length4Operator(::Type{Reconstruct{FETypeR, O}}, xdim, nc) where {FETypeR, O} = Length4Operator(O, xdim, nc)
+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, w}}) where {FETypeR, O, w} = "wR(r" * DefaultName4Operator(O) * ")"
 
-struct FEReconstEvaluator{T, TvG, TiG, FEType, FEType2, stdop, O <: ReconstructionOperator} <: FEEvaluator{T, TvG, TiG}
+struct FEReconstEvaluator{T, TvG, TiG, FEType, FEType2, stdop, RH} <: FEEvaluator{T, TvG, TiG}
     citem::Base.RefValue{Int}                       # current item
     FE::FESpace{TvG, TiG, FEType}                     # link to full FE (e.g. for coefficients)
     FEB::SingleFEEvaluator{T, TvG, TiG, stdop, FEType2}                # FEBasisEvaluator for stdop in reconstruction space
     cvals::Array{T, 3}                               # current operator vals on item (reconstruction)
     coefficients::Array{TvG, 2}                      # additional coefficients for reconstruction
-    reconst_handler::ReconstructionHandler{T, TiG}  # handler for reconstruction coefficients
+    reconst_handler::RH                             # handler for reconstruction coefficients
 end
 
 # constructor for reconstruction operators
 function FEEvaluator(
         FE::FESpace{TvG, TiG, FEType, FEAPT},
-        operator::Type{<:Reconstruct{FETypeReconst, stdop}},
+        operator::Type{<:ReconstructionOperator{FETypeReconst, stdop}},
         qrule::QuadratureRule{TvR, EG},
         xgrid = FE.xgrid;
         L2G = nothing,
@@ -65,9 +53,21 @@ function FEEvaluator(
     FEB = FEEvaluator(FE2, stdop, qrule, xgrid; L2G = L2G, T = T, AT = AT)
 
     ## reconstruction coefficient handler
-    reconst_handler = ReconstructionHandler(FE, FE2, AT, EG)
+    if operator <: WeightedReconstruct && FEType in [H1BR{2}]
+        tf = weight_type(operator)
+        reconst_handler = ReconstructionHandler(FE, FE2, AT, EG, operator, tf.instance)
+    else
+        if operator <: WeightedReconstruct
+            @warn "weighted reconstruction not available for $FEType, ignoring the weight function"
+        end
+        reconst_handler = ReconstructionHandler(FE, FE2, AT, EG, operator)
+    end
+
+    return FEReconstEvaluator{T, TvG, TiG, FEType, FETypeReconst, stdop, typeof(reconst_handler)}(FEB.citem, FE, FEB, cvals, coefficients2, reconst_handler)
+end
 
-    return FEReconstEvaluator{T, TvG, TiG, FEType, FETypeReconst, stdop, operator}(FEB.citem, FE, FEB, cvals, coefficients2, reconst_handler)
+function relocate_xref!(FEB::FEReconstEvaluator, new_xref)
+    return relocate_xref!(FEB.FEB, new_xref)
 end
 
 function update_basis!(FEBE::FEReconstEvaluator)