diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ef22c6..ffdd99e 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 c55333d..29998fc 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" diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl index f2f4e14..6d3cef6 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 b27783b..e5a60d6 100644 --- a/src/reconstructionhandlers.jl +++ b/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,11 +99,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, 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} @@ -81,6 +121,11 @@ struct ReconstructionHandler{Tv, Ti, 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 """ @@ -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 diff --git a/src/reconstructionoperators.jl b/src/reconstructionoperators.jl index 853564c..91d326c 100644 --- a/src/reconstructionoperators.jl +++ b/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) diff --git a/test/test_operators.jl b/test/test_operators.jl index d79df05..bafe31a 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)