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>
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
2 participants
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.
p-multigrid (idaholab/moose#33645) has led me to exploring our higher polynomial order spaces and I (Claude to be clear) believe I've found some bugs in various orientation related things