Conversation
Importantly this both sets the wrapper as initialized and instructs the wrapper not to destroy the PETSc matrix on exit Co-Authored-By: Claude <noreply@anthropic.com>
Co-Authored-By: Claude <noreply@anthropic.com>
Co-Authored-By: Claude <noreply@anthropic.com>
Co-Authored-By: Claude <noreply@anthropic.com>
…package Co-Authored-By: Claude <noreply@anthropic.com>
Co-Authored-By: Claude <noreply@anthropic.com>
positive_face_orientation() reports the direction in which a face's vertices run, which is a single bit. A basis whose shape functions follow the order of a face's vertices needs the whole ordering, including which vertex is the lexicographically least, so that two elements sharing a face arrive at the same shape function for a degree of freedom they share. That single bit leaves the least vertex of a face with more than three vertices open. face_orientation() carries the whole ordering, as twice the position that the least vertex takes in the face's node map plus the direction bit, so it runs over the six symmetries of a triangular face or the eight of a quadrilateral face. positive_face_orientation() now reads its low bit. edge_orientation() gives an edge's single bit the same spelling, for callers that index over edges and faces alike. The new ElemTest cases hold face_orientation() to agreeing with positive_face_orientation() and to naming the least vertex, and hold two elements sharing a face to the same ordering of that face's vertices, which is the property that lets the index select a conforming shape function. Co-Authored-By: Claude <noreply@anthropic.com>
A hierarchic basis of order three or above places degrees of freedom on the faces of a 3D element, and the two elements sharing a face have to arrive at the same shape function for each of them. Seven places did not. HEX27, in cube_indices(): face 1 case 7 and face 5 case 8 rewrote the reference coordinates with the transformation of another case, so two of the eight symmetries of those faces produced a shape function belonging to a different symmetry. PRISM*, in prism_indices(): a transposing case has to swap either the mode indices or the in-plane coordinates. The transposing cases of all three quadrilateral faces did both, which cancels, so two symmetries differing only in the direction bit produced identical shape functions. Face 2 case 1 hard-coded the offset from an edge to the second triangle's numbering, which holds only at order three, and cases 5 through 8 of that face reflected a coordinate that needs no reflection there. PRISM*, in fe_hierarchic_3D_shape(): the parity of a triangle edge mode was taken from the wrong index base, which agrees with the mode number only at order three. PRISM*, in prism_indices(): the interior basis of a triangle is not symmetric in its barycentric coordinates, so the two triangular faces need those coordinates ordered by the face's vertices, as the tetrahedral faces already are. The comment saying the TRI code handled this was mistaken; it orders the edge modes of a triangle, not the barycentric coordinates a face's interior modes read. HEX27, in cube_remap(): the flip bit of a quadrilateral face compared the least node of the face against a neighbor, which is true for any two distinct points, so the four symmetries reached through a reflection were unreachable and SIDE_HIERARCHIC was non-conforming from order three up. Rotating the least node to the origin leaves the reflection about the diagonal through it, which is settled by which of that node's two neighbors is the lesser, as orient_quad() already does. The new FEOrientationTest holds each family to two properties on a skewed mesh: that the shape functions of a degree of freedom follow the orientation index of the one entity that owns it, and that two elements sharing a side agree on the shape function of a degree of freedom they share. Every one of the fixes above is needed for a case of the second property, and the first property is what lets a cache of reference shape functions carry one entry per entity and orientation. A side-only family has no shape value at an interior point, so FEConformityTest registers the second property alone. Co-Authored-By: Claude <noreply@anthropic.com>
The one-dimensional hierarchic bubbles were (xi^p - 1)/p! for even p and (xi^p - xi)/p! for odd. The 1/p! makes a shape function shrink factorially with its order: at p = 8 it is about 2.5e-05, and a two-dimensional tensor product of two such factors carries its square. An operator assembled from that basis inherits the spread on its diagonal, because for the Laplacian the diagonal entry is the squared H1 seminorm of its shape function. On a 4x4 mesh of QUAD9 elements the diagonal of an order-8 stiffness matrix ran from 1.08e-17 to 8/3, so its smallest entry sat below the roundoff of its largest and the discretization was numerically singular in double precision however it was solved: its inverse overflowed a double from order 7 upward. The bubbles are now scaled to unit H1 seminorm on the reference interval, which replaces the factorial with a factor growing as the square root of the order. The same mesh now gives a diagonal from 0.128 to 8/3 at order 8 and a condition number of 4.4e+10, against a reciprocal machine epsilon of 4.5e+15, and the attainable residual of a solve there improves from 1.7e-12 to 4.0e-15. The two vertex functions are deliberately left unscaled. They are interpolatory, so their coefficients are values of the finite element solution at the vertices, and scaling them would change what a nodal boundary condition or a nodal output of a hierarchic variable means. This is a change of basis within one space, so the finite element solution is unaffected: the integral, L2 norm, H1 seminorm and interior point values of an order-8 solution with content in every mode all agree with the previous basis to every printed digit. Coefficient vectors change, and so does the 2-norm of a residual taken in them, which makes iteration counts and residual histories measured before and after incomparable. Two other sites carried the same factorial and move with it: the branch each of the 2D and 3D files takes when a triangle or tetrahedron edge coordinate vanishes is the limit of the general expression, so it needs the same scaling or the shape function is discontinuous there. The second derivative's general branch was also wrong for p of 8 or more, where it returned the first derivative's xi^(p-1)/(p-1)! rather than p (p-1) xi^(p-2). It is rewritten here rather than left in place, having had to be rewritten for the scaling anyway. Co-Authored-By: Claude <noreply@anthropic.com>
The coefficients a DirichletBoundary's projection determines were only ever reachable as a side effect: add the boundary, run create_dof_constraints(), read the constraint rows back, then remove the boundary and sweep again to undo the mutation. A caller that wants the coefficients rather than the constraints paid two sweeps, had to tell the rows it caused apart from the hanging-node and periodic rows already present, and mutated a DofMap that other readers share. DofMap::compute_dirichlet_values() returns them directly and constrains nothing. The projection itself is unchanged: it is the same ConstrainDirichlet over the same local element range, with a third AddConstraint implementation that records each prescribed value instead of writing a constraint row. A degree of freedom the DofMap already constrains is omitted, as AddPrimalConstraint omits it, so the values reported are those the constraint path would have added. AddConstraint no longer holds a DofMap: writing into one is a property of the two implementations that do it, not of disposing of a computed constraint, and the collecting implementation has no DofMap to write to. That leaves ConstrainDirichlet needing only a const DofMap, which is what lets the new method be const. A projected value is a coefficient in whatever basis is current, so this is how a caller obtains an essential boundary value that is correct on a modal basis such as HIERARCHIC above first order, where a degree of freedom at a node is a bubble coefficient rather than a point value of the solution. Co-Authored-By: Claude <noreply@anthropic.com>
fe_hierarchic_bubble_scaling documented the seminorm of the bubbles as sqrt(2/(2i-1))/(i-1)! for even i, which is the seminorm of the (xi^i - 1)/i! that the bubbles carried before they were normalized, not of the xi^i - 1 the same comment says they are now. The bubbles integrate to i sqrt(2/(2i-1)) and (i-1) sqrt(2/(2i-1)), and the scaling this function returns is the reciprocal of that, as the code already was. Co-Authored-By: Claude <noreply@anthropic.com>
A nodal basis whose interpolation points are the points of a quadrature rule is collocated with it: each shape function is one at its own point and zero at the others, so the mass matrix that rule assembles is diagonal, for any mapping and any element shape. That property is exact only while the basis and the rule carry the same values, so a basis that computed its own points would leave the mass matrix diagonal to roundoff rather than exactly, and any lumping that relied on it would be approximate. The table of one-dimensional points and weights is therefore reachable on its own, as QGaussLobatto::points_1D and weights_1D, keyed by the number of points rather than by the order of the rule it realizes. The two indices differ: a rule of n points is exact through degree 2n-3, so the orders 2n-4 and 2n-3 share a rule, which is why the cases of the table come in pairs and which n_points_1D inverts. A basis of degree p reads the p+1 point rule. The table's values are untouched. init_1D now copies from it, and every rule it holds is built once on first use so that reaching one needs no lock. Co-Authored-By: Claude <noreply@anthropic.com>
The basis of degree p interpolates at the p+1 Gauss-Lobatto points, which it reads from QGaussLobatto so that it is collocated with the rule of the same count. At one of its own points a shape function is one or zero and its derivatives are rows of the differentiation matrices the basis carries, so the values a collocated quadrature rule asks for cost nothing to produce. Away from them the barycentric form gives the value and its first two derivatives. The differentiation matrix is built from the ratios of the barycentric weights, with each diagonal entry the negation of the sum of its row, since the shape functions sum to one everywhere and so their derivatives sum to zero. That identity supplies the endpoint entries, -p(p+1)/4 and p(p+1)/4, without a case of their own. The second derivative matrix is the first applied twice, which is exact because the derivative of a polynomial of degree p has degree p-1 and this basis interpolates it exactly. The bases for every tabulated point count are built together on first use, a few thousand operations, so that reaching one needs no lock. Co-Authored-By: Claude <noreply@anthropic.com>
fd8d302 added meshes/quad4_tri3_hourglass.xda.gz to the data that tests/Makefile.am distributes without regenerating tests/Makefile.in, so the generated file has been a line behind its input since. Co-Authored-By: Claude <noreply@anthropic.com>
A nodal basis on the Gauss-Lobatto points, collocated with the Gauss-Lobatto rule of the same order: each shape function is one at its own point and zero at the others, so the mass matrix that rule assembles is diagonal, for any mapping and any element shape, and the values the rule asks for cost nothing to produce. The basis is a tensor product, so it lives on the tensor product element types alone. Every one of its degrees of freedom belongs to the element, which is what frees its order from the element's node count: a degree of freedom sits at a Gauss-Lobatto point, and past order two those are not nodes, so an order of eight runs on a plain QUAD4. The orders run from one to one less than the number of tabulated Gauss-Lobatto points, and asking beyond that says so. FEType::default_quadrature_rule returns the collocated rule, so a variable of this family is collocated without a caller arranging it. A rule of n points is exact through degree 2n-3, so the order yielding the p+1 points a basis of degree p interpolates at is 2p-1. That expression carries p refinement as well, since QBase raises a rule's order by twice the p_level it is initialized with and 2p-1 rises by two for each degree. Stiffness is still integrated exactly on an affine element, where the integrand has degree 2p-2 against the rule's 2p-1. The mass integrand has degree 2p, so it is the one term the rule does not reach, which is the trade that makes the matrix diagonal. Co-Authored-By: Claude <noreply@anthropic.com>
The existing FETest harness covers the family's shape functions, their first and second derivatives, and the interface reporting its degree of freedom counts, at orders one through four on the tensor product element types. The new FEGLLCollocationTest covers what is particular to this family: that the matrix of shape function values at the points of the family's own quadrature rule is the identity, and that the mass matrix that rule assembles is diagonal, carrying the weight of each point times the Jacobian there. Both hold exactly rather than to a tolerance, which is what the basis reading its points from QGaussLobatto buys, so the assertions carry none. The mesh is skewed so that a diagonal mass matrix is a property of the basis rather than of a Jacobian that happens to be constant. testPartitionOfUnity fails an unlisted family rather than assume an answer for it, so the family joins those whose shape functions sum to one. Co-Authored-By: Claude <noreply@anthropic.com>
The tensor product evaluation and the table of second derivative coordinate pairs stood in each of the three dimension files, so a fix to either wanted three applications and the C0 family would have added three more copies. They stand once now, in a header the dimension files reach through, leaving those files holding only the specializations that FE requires per dimension and family, as fe_lagrange_shape_2D.C already does for LAGRANGE and L2_LAGRANGE. What separates the families is which tensor index a shape function index names, and that alone, so it is a function of its own that each family specializes. A family without a specialization does not link. The pair table is now one table rather than one per dimension: libMesh numbers the second derivatives so that a given dimension uses the leading d(d+1)/2 of them, an edge reading the first entry, a face the first three and a cell all six. Co-Authored-By: Claude <noreply@anthropic.com>
At orders one and two on the tensor product element types the Gauss-Lobatto nodal families
and LAGRANGE are the same functions, because the Gauss-Lobatto points of those orders, {-1,1}
and {-1,0,1}, are where those elements put their nodes. They reach those functions by
different routes: LAGRANGE evaluates a closed form per order, and the Gauss-Lobatto basis
evaluates a barycentric form against its point table. Nothing held the two together, so
either could have drifted without a test noticing.
This pairs each Gauss-Lobatto shape function with the LAGRANGE one whose node sits at the
same reference point and holds them to the same values away from those points. It is what
would catch the point table, the barycentric evaluation, or the tensor index map moving away
from values that are already settled.
The serendipity types are left out. LAGRANGE spans a smaller space on QUAD8 and HEX20, eight
and twenty functions against the nine and twenty seven of a tensor product of the same order,
so there is no pairing to hold there.
Co-Authored-By: Claude <noreply@anthropic.com>
The continuous nodal basis on the Gauss-Lobatto points, collocated with the Gauss-Lobatto rule of its order as the L2 family already is, so its mass matrix is diagonal. It is diagonal globally and not only per element, since the map from a degree of freedom to the point it interpolates at is one to one. The degrees of freedom are distributed as HIERARCHIC distributes its own, one per vertex, p-1 per edge, (p-1)^2 per face and the rest on the element, so the element needs a node for each entity that owns any. That is what limits the order a type reaches: without mid-edge nodes an edge's degrees of freedom have nowhere to sit, so QUAD4 and HEX8 stop at order one, and in three dimensions a face owns degrees of freedom too, so HEX20 stops there as well while HEX27 does not. Interior degrees of freedom belong to the element, so no type needs an interior node and QUAD8 reaches any order. Two elements sharing an entity have to agree on which point each of its degrees of freedom sits at, and each numbers that entity's vertices its own way. Rather than enumerate the symmetries, the index is built from the vertices themselves: a vertex of a tensor product element sits at -1 or 1 in each coordinate, so it converts straight to a tensor index, and an edge or face degree of freedom is that index stepped along the entity from a corner both elements pick out. Elem::edge_orientation names the lesser end of an edge, and Elem::face_orientation the least vertex of a face together with the direction the rest run from it, which is an origin and two axes. So one expression covers every edge and every face of every tensor product element instead of one case apiece. extra_hanging_dofs is true, with HIERARCHIC and the other continuous high order families rather than with LAGRANGE. Past order two none of the points a mid-edge or mid-face node owns is the node itself, so that node cannot double as a vertex for the elements meeting it at a refinement interface and needs a degree of freedom of its own there. Co-Authored-By: Claude <noreply@anthropic.com>
The FETest harness covers the shape functions and their first two derivatives at orders one through four, and FEOrientationTest covers what a continuous basis needs beyond that: that the shape functions of a degree of freedom follow the orientation of the one entity owning it, and that two elements sharing a side arrive at the same shape function for a degree of freedom they share. Those are the properties the index map is built to satisfy, and they are held at orders three through five on QUAD9 and HEX27, over a skewed mesh and every permutation of an element, so each of the eight symmetries of a quadrilateral face is reached. The collocation test now runs over both Gauss-Lobatto families. They differ in the order they number their degrees of freedom rather than in where those sit, so the L2 family's matrix of shape function values is the identity while the C0 family's is a permutation of it. What both are held to is that each shape function is one at exactly one point and exactly zero at the rest and that no two claim the same point, which is what the diagonal mass matrix rests on. The mass matrix is still summed the way an assembly loop would sum it, so that the diagonal is a property of that arithmetic rather than one read back from the values just asserted. Co-Authored-By: Claude <noreply@anthropic.com>
A conforming basis on an adaptively refined mesh needs the degrees of freedom on a hanging node constrained to the coarse side they lie on. Nothing here exercised that machinery for any family, and it is invisible on a conforming mesh, so a family could carry the wrong answer for extra_hanging_dofs and every existing test would still pass. This holds the constrained space to containing the polynomials of its own order, which it must, since a polynomial that smooth is continuous across the interface and so already conforms. Projecting one onto the space therefore reproduces it, and a wrong constraint shows up as a solution that no longer does, because enforcing it moves the hanging degrees of freedom off the values that interpolate the polynomial. LAGRANGE is the control: its constraints are long settled, so a failure there would be this test's fault rather than a family's. HIERARCHIC covers the other family allocating a separate degree of freedom at a hanging node, and LAGRANGE_GLL is what the test was written for. Setting that family's extra_hanging_dofs back to false fails exactly its three cases and leaves the other four passing, so the test does discriminate. Co-Authored-By: Claude <noreply@anthropic.com>
include/include_HEADERS names the headers that are installed, and the two the Gauss-Lobatto basis adds were missing from it, so a consumer building against an installed libMesh could not reach them however well the library itself had compiled. The build tree carries its own header farm, which is why nothing inside libMesh noticed. The script that generates that list did not scan include/gpu, so the five Kokkos headers there had been added to the generated file by hand and running the script would have dropped them. It scans gpu now and those five are generated along with the rest, which leaves the list reproducible: the next person to add a header regenerates it without silently uninstalling someone else's. Co-Authored-By: Claude <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.