Repository navigation
Bugfix: Dual Lagrange spaces - #218
Conversation
`assemble_local_refines!` maps every fine test cell to its parent cell and looks that cell up in the trial space. When the trial space is supported on only part of the parent mesh, the lookup returns 0 and `bels[0]` throws a BoundsError. `assemble_local_matched!` already guards against this. Added the same guard to `assemble_local_refines!`. Dual lagrange test functions hit this case because their support reaches into coarse cells that carry no trial functions.
The rule hardcoded three vertices, so it only worked for triangles. A refined curve has two-vertex charts and raised a BoundsError. Map over the chart's vertices instead. For a triangle this produces the same three arguments as before.
The three-argument restrict for degree two called numfunctions without a domain, for which there is no method, so it threw on every call. Delete it and let the generic RefSpace fallback route through interpolate, which gives the same result as the in-place five-argument form. Add a test comparing the three routes to a restriction matrix on a genuine restriction onto a child cell.
For a curve the default junction is zero-dimensional, so the call dispatched to the method that reads its second argument as a vertex list given as a mesh. The empty default junction became an empty vertex list and the constructor silently returned a space with no functions. A junction passed explicitly hit the same method, yielding dual functions on the two end vertices only instead of on every vertex. Give the one-argument form its own method and decide between junction and vertex list by comparing dimensions rather than by dispatch, which cannot separate the two readings on a curve. Surfaces are unaffected.
The three coefficients were hardcoded for a closed surface: one at the barycentre, one half at a face centre, and one over half the number of fine cells at a coarse vertex. The face-centre value assumes two coarse cells share every edge, which fails on an open surface, where a boundary edge belongs to one. The dual functions there summed to one half rather than one, so the space could not represent a constant. Replace all three with the fraction of the fine cells meeting at a vertex that belong to the coarse cell, and fold the three loops into one. On a closed surface this reproduces the previous shapes exactly.
Upstream fixed the array size, but the shape assignment itself was wrong: The construction assigned the shape at a coarse vertex only on the fine cell inside the coarse segment, and put a shape at the neighbouring segment's midpoint on the adjacent cell. The resulting functions jumped by one across every coarse vertex, so the space was not continuous. It also indexed of a surface, indexed a second adjacent cell that does not exist at the end of an open curve, reached into mesh.faces directly, and ignored the junction predicate it was passed. Assign to every fine vertex the fraction of the fine cells meeting there that belong to the coarse segment, on each incident cell, as the surface construction above now does.
Lagrange reference spaces on segments and tetrahedra had no interpolat!, so restrict! thre a MethodError for every one-dimensional cell and tetrahedral degreee one. Degree zero on a tetrahedron had no evaluator at all. Add the missing methods in the same inline style as the triangle versions. The segment nodes are the two end points followed by the interior nodes in descending order, which is the order the hand-written evaluators one to ten use. On a tetrahedron degree one interpolates at the four vertices; higher degrees have no evaluator, so nothing further is needed at the moment.
The average position was built from three literal coordinates, so every mesh whose vertices have two components raised a DimensionMismatch. Use zero(vertextype(mesh)). Makes _surface usable for curves and surfaces in the plane. Also add the testitem for dual linear Lagrange on curves. It covers a closed polygon, an open segment, a curve whose vertex numbering is inherited from a parent surface, and the junction variant, asserting continuity and the shape counts at open ends.
Pair a dual linear Lagrange space on the barycentric refinement of a circle with piecewise constants on the coarse curve, and compare the assembled single layer against the same operator with both spaces on the fine mesh. Before the refinement of a curve was tagged, this pairing fell back to conforming quadrature and produced 63 infinite entries out of 441. The test needs the CompScienceMeshes change that returns a BarycentricRefinement for curves, so it fails against released version (currently 0.12.1) until that is merged.
|
Thanks a lot for this! W.r.t. the dual lagrange functions. Have you checked that the open screen versions simply require setting the weights on the boundary to zero? For the BC functions, the situation is quite a bit more complicated. Please see the Buffa/Christiansen paper for details. I forgot to what extend this was implemented (judging by your assessment not at all for open screens). Ideally the implementation should work for mixed boundary conditions, i.e. allow for the construction of basis functions that are forced to be zero on a subset of the boundary of |
|
Thank you for the comment! I am actually not zeroing the weights on the boundary, but make sure that it is a partition of unity. If I see it correctly, this agrees with what is stated in the paper in section 4.2 (i) and is shown in figure 4 --> we have now on the edge midpoints on open boundaries a weight of 1 (automatically assigned via checking how many coarse cells meet there), where it was previously hardcoded as 1/2, which is correct for every edge midpoint on closed surfaces. Regarding construction of basis functions that are zero on a subset of the boundary, this was already implemented for surfaces before, but with the wrong 1/2 factor on the free edges. using BEAST, CompScienceMeshes
# surface case
Γ = meshrectangle(1.0, 1.0, 0.25, 3) # open square patch
bnd = boundary(Γ)
# the part of the boundary where the functions should vanish: the edges at x = 0
on_x0(m, c) = all(v -> abs(v[1]) < 1e-8, chart(m, c).vertices)
dir = submesh(on_x0, bnd) # 4 of the 16 boundary edges
X_free = duallagrangec0d1(Γ) # partition of unity up to the whole boundary
X_zero = duallagrangec0d1(Γ, bnd) # vanishes on all of the boundary
X_mixed = duallagrangec0d1(Γ, dir) # vanishes on x = 0 only
# curve case
Λ = meshsegment(1.0, 0.25, 2) # open curve, 4 segments in a 2D universe
bnd_Λ = boundary(Λ) # the two end points
dir_Λ = submesh(on_x0, bnd_Λ)
Y_free = duallagrangec0d1(Λ)
Y_zero = duallagrangec0d1(Λ, bnd_Λ)
Y_mixed = duallagrangec0d1(Λ, dir_Λ)I am less certain about the zeroing itself. Buffa and Christiansen put the homogeneous boundary conditions only on the primal spaces and construct the dual space as a partition of unity, so neither the fully nor the partly zeroed variant is the construction in the paper. The mechanism was already in the code, though, so I kept that behaviour and extended it to curves rather than questioning it. |
|
Thanks. I will merge this and see if changes are needed to support mixed boundary conditions for dual functions. |
|
Perfect, thank you! |
2 participants
Ten commits fixing bugs that were found while trying to build dual P0/P1 discretisations on curves and open surfaces
Dual lagrange spaces:
cf256c32duallagrangecxd0(curve)silently returned an empty space: the default empty zero-dimensional junction mesh dispatched to the nodes-as-mesh method. Dispatch now keys ondimension(jct).6ee61431The hardcoded1/2for face centres induallagrangec0d1(…, Val{3})assumes two coarse cells per edge, which fails on an open patch. All three coefficients are replaced by one rule: the fraction of the fine cells at a vertex that belong to the coarse cell, which reproduces the old values on a closed mesh and fixes the open case.3106d7c2duallagrangec0d1(…, Val{2})rewritten. Duallagranec0d1 error #195's array-size crash was already fixed by49213e18, but the remaining problems were not: the functions were discontinuous across cells, open end points were mishandled,SubMeshinput was unsupported, and thepredargument was ignored. It now uses the same coefficient rule asVal{3}and the generic mesh API.48f5064cunitfunctioncxd0hardcoded a three-component vertex, so it threw in a two-dimensional universe.Refinement-aware assembly and quadrature
f27e616eassemble_local_refines!indexesbels[0]when a parent cell contains no trial cells, which happens whenever the trial space does not cover the whole parent mesh.assemble_local_matched!already had this guard, which is now also included inassemble_local_refines!e193b2d7TestRefinesTrialQRulebuilt the overlap simplex fromchart.vertices[1:3]in Cartesian coordinates. Built generically from all vertices mapped to barycentric coordinates, it works in any dimension.d2a642eb,7c670438Tests for refinement-aware quadrature on a curve, usingTrialRefinesTestQRuleandTestRefinesTrialQRule.Local Lagrange interpolation
bfa8805drestrict(::LagrangeRefSpace{T,2}, dom1, dom2)threw, becausenumfunctions(f)was used without a domain and assumed a vertices-then-edges node order. Removed, so degree 2 routes through the genericinterpolatepath.ac51bfedinterpolate/interpolate!for segments (degrees 0-10) and tetrahedra (degrees 0-1, the evaluable ones), plus a{T,0,4}evaluator, following the inline node style the rest of the file now uses.Dependency on CompScienceMeshes
3106d7c2,48f5064c,d2a642eband7c670438require krcools/CompScienceMeshes.jl#76, which tags curve refinements and fixesisinclosurefor point cells.The
Val{2}rewrite above is also what keeps BEAST working after that PR: current master's version reads mesh2.faces, and the new BarycentricRefinement return type has no such field. So this and the CompScienceMeshes PR need to happen together, CompScienceMeshes first.With the CompScienceMeshes branch in the environment, the test suite passes, including the new tests.
The
[compat]lower bound here will need raising once CompScienceMeshes is tagged.