From b3c2a5a8f09f59cdf499e2fb82214a8bfa468c0a Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Tue, 24 Jun 2025 11:51:08 +0200 Subject: [PATCH 1/7] draft of a weighted reconstruction --- src/ExtendableFEMBase.jl | 2 +- src/reconstructionhandlers.jl | 117 ++++++++++++++++++++++++++++++--- src/reconstructionoperators.jl | 28 +++----- 3 files changed, 116 insertions(+), 31 deletions(-) diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl index df3983c5..7fbd8261 100644 --- a/src/ExtendableFEMBase.jl +++ b/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 diff --git a/src/reconstructionhandlers.jl b/src/reconstructionhandlers.jl index b27783b2..459179c8 100644 --- a/src/reconstructionhandlers.jl +++ b/src/reconstructionhandlers.jl @@ -1,3 +1,27 @@ +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 + + +""" +$(TYPEDEF) + +weighted reconstruction operator: evaluates a reconstructed version of the weighted 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 WeightedReconstruct{FETypeR, O} <: ReconstructionOperator{FETypeR, O} end + + ################## SPECIAL INTERPOLATORS #################### """ @@ -69,11 +93,21 @@ 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} 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} @@ -90,7 +124,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) 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 @@ -102,9 +136,11 @@ function ReconstructionHandler(FES::FESpace{Tv, Ti, FE1, APT}, FES_Reconst::FESp 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}(FES, FES_Reconst, xCoordinates, xFaceNodes, xFaceVolumes, xFaceNormals, xCellFaceOrientations, xCellFaces, interior_offset, interior_ndofs, coeffs) end """ @@ -113,7 +149,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] @@ -446,7 +482,7 @@ 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 @@ -474,7 +510,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 <: 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 @@ -498,7 +534,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 <: 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 +565,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 +603,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 +628,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 @@ -629,3 +665,64 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, end return nothing end + + +##### WEIGHTED + +function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: WeightedReconstruct, FE1 <: H1BR{2}, FE2 <: HDIVBDM1{2}, AT <: ON_CELLS, EG <: Union{Triangle2D, Quadrilateral2D}} + xFaceVolumes = RH.xFaceVolumes + xFaceNormals = RH.xFaceNormals + xCellFaceSigns = RH.xCellFaceOrientations + xFaceNodes = RH.xFaceNodes + xCellFaces = RH.xCellFaces + 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 + for f in 1:nfaces + face = xCellFaces[f, cell] + rmid = (xCoordinates[1, xFaceNodes[1, face]] + xCoordinates[1, xFaceNodes[2, face]]) / 2 + 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] = xFaceVolumes[face] * (1 // 6 * xCoordinates[1, xFaceNodes[n, face]] + 4 // 6 * rmid) * xFaceNormals[k, face] + # BDM1 reconstruction coefficients for P1 functions on reference element + coefficients[nfaces * (k - 1) + node, 2 * (f - 1) + 2] = xFaceVolumes[face] * xFaceNormals[k, face] * xCellFaceSigns[f, cell] * (BDM1_coeffs[n] * xCoordinates[1, xFaceNodes[n, face]] + 1 // 3 * rmid) + end + end + # RT0 reconstruction coefficients for face bubbles on reference element + coefficients[nfaces * 2 + f, 2 * (f - 1) + 1] = xFaceVolumes[face] * rmid + end + return nothing +end + +function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: WeightedReconstruct, FE1 <: H1BR{2}, FE2 <: HDIVRT0{2}, AT <: ON_CELLS, EG <: Union{Triangle2D, Quadrilateral2D}} + xFaceVolumes = RH.xFaceVolumes + xFaceNormals = RH.xFaceNormals + xFaceNodes = RH.xFaceNodes + xCellFaces = RH.xCellFaces + xCoordinates = RH.xCoordinates + face_rule = local_cellfacenodes(EG) + nnodes = size(face_rule, 1) + nfaces = size(face_rule, 2) + node = 0 + face = 0 + for f in 1:nfaces + face = xCellFaces[f, cell] + # reconstruction coefficients for P1 functions on reference element + rmid = (xCoordinates[1, xFaceNodes[1, face]] + xCoordinates[1, xFaceNodes[2, face]]) / 2 + for n in 1:nnodes + node = face_rule[n, f] + for k in 1:2 + coefficients[nfaces * (k - 1) + node, f] = xFaceVolumes[face] * (1 // 6 * xCoordinates[1, xFaceNodes[n, face]] + 1 // 3 * rmid) * xFaceNormals[k, face] + end + end + # reconstruction coefficients for face bubbles on reference element + coefficients[2 * nfaces + f, f] = xFaceVolumes[face] * rmid + end + return nothing +end diff --git a/src/reconstructionoperators.jl b/src/reconstructionoperators.jl index 42702536..24be3f95 100644 --- a/src/reconstructionoperators.jl +++ b/src/reconstructionoperators.jl @@ -1,22 +1,10 @@ -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}}) where {FETypeR, O} = "R(r" * DefaultName4Operator(O) * ")" struct FEReconstEvaluator{T, TvG, TiG, FEType, FEType2, stdop, O <: ReconstructionOperator} <: FEEvaluator{T, TvG, TiG} citem::Base.RefValue{Int} # current item @@ -24,13 +12,13 @@ struct FEReconstEvaluator{T, TvG, TiG, FEType, FEType2, stdop, O <: Reconstructi 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::ReconstructionHandler{T, TiG, O} # 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,7 +53,7 @@ function FEEvaluator( FEB = FEEvaluator(FE2, stdop, qrule, xgrid; L2G = L2G, T = T, AT = AT) ## reconstruction coefficient handler - reconst_handler = ReconstructionHandler(FE, FE2, AT, EG) + reconst_handler = ReconstructionHandler(FE, FE2, AT, EG, operator) return FEReconstEvaluator{T, TvG, TiG, FEType, FETypeReconst, stdop, operator}(FEB.citem, FE, FEB, cvals, coefficients2, reconst_handler) end From 65ef6acbe935bc2779091123e541e412c61817b3 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Wed, 16 Jul 2025 19:08:54 +0200 Subject: [PATCH 2/7] reconstruction operator with arbeitrary weight function now possible --- src/reconstructionhandlers.jl | 123 ++++++++++++++------------------- src/reconstructionoperators.jl | 13 ++-- 2 files changed, 60 insertions(+), 76 deletions(-) diff --git a/src/reconstructionhandlers.jl b/src/reconstructionhandlers.jl index 459179c8..5bb55c35 100644 --- a/src/reconstructionhandlers.jl +++ b/src/reconstructionhandlers.jl @@ -14,12 +14,16 @@ abstract type Reconstruct{FETypeR, O} <: ReconstructionOperator{FETypeR, O} end """ $(TYPEDEF) -weighted reconstruction operator: evaluates a reconstructed version of the weighted finite element function. +Weighted reconstruction operator: evaluates a reconstructed version of a finite element function, multiplied by a weight 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. +# 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). """ -abstract type WeightedReconstruct{FETypeR, O} <: ReconstructionOperator{FETypeR, O} end +abstract type WeightedReconstruct{FETypeR, O, F} <: Reconstruct{FETypeR, O} end + +weight_type(::Type{<:WeightedReconstruct{FETypeR, O, F}}) where {FETypeR, O, F} = F ################## SPECIAL INTERPOLATORS #################### @@ -103,7 +107,7 @@ $(TYPEDEF) struct for storing information needed to evaluate a reconstruction operator """ -struct ReconstructionHandler{Tv, Ti, RT, 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} @@ -115,6 +119,11 @@ struct ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG} 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 """ @@ -124,7 +133,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, RT) 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 @@ -134,13 +143,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, RT, FE1, FE2, AT, EG}(FES, FES_Reconst, xCoordinates, xFaceNodes, 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 """ @@ -487,26 +497,38 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, 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) + w[1] = weight(x1) + w[2] = weight(x2) + 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] + 1 // 3 * wmid) 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 @@ -514,22 +536,39 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, 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 @@ -666,63 +705,3 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, return nothing end - -##### WEIGHTED - -function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: WeightedReconstruct, FE1 <: H1BR{2}, FE2 <: HDIVBDM1{2}, AT <: ON_CELLS, EG <: Union{Triangle2D, Quadrilateral2D}} - xFaceVolumes = RH.xFaceVolumes - xFaceNormals = RH.xFaceNormals - xCellFaceSigns = RH.xCellFaceOrientations - xFaceNodes = RH.xFaceNodes - xCellFaces = RH.xCellFaces - 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 - for f in 1:nfaces - face = xCellFaces[f, cell] - rmid = (xCoordinates[1, xFaceNodes[1, face]] + xCoordinates[1, xFaceNodes[2, face]]) / 2 - 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] = xFaceVolumes[face] * (1 // 6 * xCoordinates[1, xFaceNodes[n, face]] + 4 // 6 * rmid) * xFaceNormals[k, face] - # BDM1 reconstruction coefficients for P1 functions on reference element - coefficients[nfaces * (k - 1) + node, 2 * (f - 1) + 2] = xFaceVolumes[face] * xFaceNormals[k, face] * xCellFaceSigns[f, cell] * (BDM1_coeffs[n] * xCoordinates[1, xFaceNodes[n, face]] + 1 // 3 * rmid) - end - end - # RT0 reconstruction coefficients for face bubbles on reference element - coefficients[nfaces * 2 + f, 2 * (f - 1) + 1] = xFaceVolumes[face] * rmid - end - return nothing -end - -function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, RT, FE1, FE2, AT, EG}, cell) where {Tv, Ti, RT <: WeightedReconstruct, FE1 <: H1BR{2}, FE2 <: HDIVRT0{2}, AT <: ON_CELLS, EG <: Union{Triangle2D, Quadrilateral2D}} - xFaceVolumes = RH.xFaceVolumes - xFaceNormals = RH.xFaceNormals - xFaceNodes = RH.xFaceNodes - xCellFaces = RH.xCellFaces - xCoordinates = RH.xCoordinates - face_rule = local_cellfacenodes(EG) - nnodes = size(face_rule, 1) - nfaces = size(face_rule, 2) - node = 0 - face = 0 - for f in 1:nfaces - face = xCellFaces[f, cell] - # reconstruction coefficients for P1 functions on reference element - rmid = (xCoordinates[1, xFaceNodes[1, face]] + xCoordinates[1, xFaceNodes[2, face]]) / 2 - for n in 1:nnodes - node = face_rule[n, f] - for k in 1:2 - coefficients[nfaces * (k - 1) + node, f] = xFaceVolumes[face] * (1 // 6 * xCoordinates[1, xFaceNodes[n, face]] + 1 // 3 * rmid) * xFaceNormals[k, face] - end - end - # reconstruction coefficients for face bubbles on reference element - coefficients[2 * nfaces + f, f] = xFaceVolumes[face] * rmid - end - return nothing -end diff --git a/src/reconstructionoperators.jl b/src/reconstructionoperators.jl index dc1a1e90..c8923933 100644 --- a/src/reconstructionoperators.jl +++ b/src/reconstructionoperators.jl @@ -6,13 +6,13 @@ Length4Operator(::Type{<:ReconstructionOperator{FETypeR, O}}, xdim, nc) where {F DefaultName4Operator(::Type{Reconstruct{FETypeR, O}}) where {FETypeR, O} = "R(" * DefaultName4Operator(O) * ")" DefaultName4Operator(::Type{WeightedReconstruct{FETypeR, O}}) where {FETypeR, O} = "R(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, O} # handler for reconstruction coefficients + reconst_handler::RH # handler for reconstruction coefficients end # constructor for reconstruction operators @@ -53,9 +53,14 @@ function FEEvaluator( FEB = FEEvaluator(FE2, stdop, qrule, xgrid; L2G = L2G, T = T, AT = AT) ## reconstruction coefficient handler - reconst_handler = ReconstructionHandler(FE, FE2, AT, EG, operator) + if operator <: WeightedReconstruct + tf = weight_type(operator) + reconst_handler = ReconstructionHandler(FE, FE2, AT, EG, operator, tf.instance) + else + reconst_handler = ReconstructionHandler(FE, FE2, AT, EG, operator) + end - return FEReconstEvaluator{T, TvG, TiG, FEType, FETypeReconst, stdop, operator}(FEB.citem, FE, FEB, cvals, coefficients2, reconst_handler) + return FEReconstEvaluator{T, TvG, TiG, FEType, FETypeReconst, stdop, typeof(reconst_handler)}(FEB.citem, FE, FEB, cvals, coefficients2, reconst_handler) end function update_basis!(FEBE::FEReconstEvaluator) From 335c67420e3a8027d837599653200c33a92b3ee6 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Wed, 16 Jul 2025 19:21:36 +0200 Subject: [PATCH 3/7] weighted BDM1 reconstruction for BR should work now --- src/reconstructionhandlers.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/reconstructionhandlers.jl b/src/reconstructionhandlers.jl index 5bb55c35..34483f06 100644 --- a/src/reconstructionhandlers.jl +++ b/src/reconstructionhandlers.jl @@ -513,8 +513,13 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, x1 = view(xCoordinates, :, xFaceNodes[1, face]) x2 = view(xCoordinates, :, xFaceNodes[2, face]) xmid .= 0.5 * (x1 .+ x2) - w[1] = weight(x1) - w[2] = weight(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] @@ -522,7 +527,7 @@ function boundary_coefficients!(coefficients, RH::ReconstructionHandler{Tv, Ti, # RT0 reconstruction coefficients for P1 functions on reference element 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] = xFaceVolumes[face] * xFaceNormals[k, face] * xCellFaceSigns[f, cell] * (BDM1_coeffs[n] * w[n] + 1 // 3 * wmid) + 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 From 630d5352b29049a9b443812087196b7f66ba4545 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Wed, 16 Jul 2025 22:34:20 +0200 Subject: [PATCH 4/7] fixed reconstruction for H1CR{2}, added warning if weighting is not available --- src/reconstructionhandlers.jl | 2 +- src/reconstructionoperators.jl | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/reconstructionhandlers.jl b/src/reconstructionhandlers.jl index 34483f06..6169680e 100644 --- a/src/reconstructionhandlers.jl +++ b/src/reconstructionhandlers.jl @@ -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, <: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 diff --git a/src/reconstructionoperators.jl b/src/reconstructionoperators.jl index c8923933..fc6ad71d 100644 --- a/src/reconstructionoperators.jl +++ b/src/reconstructionoperators.jl @@ -53,10 +53,13 @@ function FEEvaluator( FEB = FEEvaluator(FE2, stdop, qrule, xgrid; L2G = L2G, T = T, AT = AT) ## reconstruction coefficient handler - if operator <: WeightedReconstruct + 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 From 3c5022c1005f66e731b5aeaa034307331ed6e669 Mon Sep 17 00:00:00 2001 From: chmerdon Date: Wed, 16 Jul 2025 22:38:16 +0200 Subject: [PATCH 5/7] changelog and version bump --- CHANGELOG.md | 4 ++++ Project.toml | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ef22c6d..ffdd99ee 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Project.toml b/Project.toml index c55333d5..29998fc8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExtendableFEMBase" uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c" authors = ["Christian Merdon ", "Patrick Jaap "] -version = "1.3.0" +version = "1.4.0" [deps] DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" From 14729041975e42e93b36c01a7df55082b803d834 Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Thu, 17 Jul 2025 09:55:27 +0200 Subject: [PATCH 6/7] added tests, reconstruction operators can now also used with lazy_interpolate! --- src/reconstructionhandlers.jl | 9 +++---- src/reconstructionoperators.jl | 6 ++++- test/test_operators.jl | 49 ++++++++++++++++++++++++++++++++++ 3 files changed, 58 insertions(+), 6 deletions(-) diff --git a/src/reconstructionhandlers.jl b/src/reconstructionhandlers.jl index 6169680e..102645a5 100644 --- a/src/reconstructionhandlers.jl +++ b/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 - diff --git a/src/reconstructionoperators.jl b/src/reconstructionoperators.jl index fc6ad71d..91d326c8 100644 --- a/src/reconstructionoperators.jl +++ b/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) diff --git a/test/test_operators.jl b/test/test_operators.jl index d79df052..bafe31a7 100644 --- a/test/test_operators.jl +++ b/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) From 4046417784b0b849d71bdc5097a74331f2b26ffc Mon Sep 17 00:00:00 2001 From: Christian Merdon Date: Thu, 17 Jul 2025 11:53:02 +0200 Subject: [PATCH 7/7] fixed docstring, documentation should bould again --- src/reconstructionhandlers.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/reconstructionhandlers.jl b/src/reconstructionhandlers.jl index 102645a5..e5a60d67 100644 --- a/src/reconstructionhandlers.jl +++ b/src/reconstructionhandlers.jl @@ -12,9 +12,11 @@ abstract type Reconstruct{FETypeR, O} <: ReconstructionOperator{FETypeR, O} end """ -$(TYPEDEF) +WeightedReconstruct{FETypeR, O, w} <: Reconstruct{FETypeR, O} -Weighted reconstruction operator: evaluates a reconstructed version of a finite element function, multiplied by a weight function. +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).