Isoparametric FEM2D_P2 element map; positive-weight guard (v0.12.3)

FEM2D_P2 now builds the element geometry from all 7 P2+bubble node positions via
the shape functions, so displacing edge/centroid nodes genuinely curves the
element (node-varying Jacobian) instead of being silently flattened to the affine
corner triangle. The physical-derivative operators use J(node)^{-T} per node and
the quadrature weight is the signed det J(node)·R.w. Dropping abs(det) means
inverted/folded elements now produce non-positive weights.

FEM3D was already isoparametric (per-node J from D_xi/D_eta/D_zeta over all
(k+1)^3 nodes); switch its weights from abs(det J) to signed det J for the same
reason, and drop a dead duplicate z-derivative block in _fem3d_structured.

Both discretizations now refuse to build (hence to solve) when any quadrature
weight is non-positive, since the barrier method requires strictly positive
weights. Straight elements are bit-unchanged (tests green, dx/dy/dz ~1e-14).

Author: Sebastien Loisel

Project.toml

@@ -1,7 +1,7 @@
 name = "MultiGridBarrier"
 uuid = "9e2c1f1d-9131-4ad4-b32f-bd2a0b0ecd1e"
 authors = ["Sébastien Loisel"]
-version = "0.12.2"
+version = "0.12.3"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"

src/Mesh3d/Geometry.jl

@@ -72,10 +72,24 @@ function _geometric_fem3d_mg(::Type{T}=Float64; L::Int=2,
                      y_xi[i] y_eta[i] y_zeta[i];
                      z_xi[i] z_eta[i] z_zeta[i]]
 
-                detJ = abs(det(J))
+                # Signed det J: a non-positive value means the isoparametric Q_k
+                # hex map is folded / inverted at this node (guarded below).
+                detJ = det(J)
                 w_physical[start_idx + i - 1] = ref_el.weights_ref[i] * detJ
             end
         end
+
+        # The barrier method requires strictly positive quadrature weights. Refuse
+        # to build (and hence to solve) on a folded / inverted curved element.
+        if !all(>(zero(T)), w_physical)
+            bad      = findall(<=(zero(T)), w_physical)
+            badelems = sort!(unique((bad .- 1) .÷ n_nodes_per_elem .+ 1))
+            error("fem3d: non-positive quadrature weight at $(length(bad)) node(s) " *
+                  "across $(length(badelems)) element(s) (first few: $(first(badelems, 5))). " *
+                  "The isoparametric Q_k element map has det J ≤ 0 there — the curved " *
+                  "hex is folded or has inverted vertex orientation. Supply orientation-" *
+                  "preserving, non-self-intersecting elements.")
+        end
         return w_physical
     end
 
@@ -210,11 +224,6 @@ function _fem3d_structured(disc::FEM3D{T}, meshes, weights, L, k, ref_el) where
         y_xi = D_xi_dense * y_elem
         y_eta = D_eta_dense * y_elem
         y_zeta = D_zeta_dense * y_elem
-        z_xi = D_zeta_dense * z_elem  # NOTE: keep parity with original code (uses D_zeta?). See below.
-        z_eta = D_eta_dense * z_elem
-        z_zeta = D_zeta_dense * z_elem
-
-        # The original uses D_xi/eta/zeta_dense for z components. Re-implement faithfully.
         z_xi = D_xi_dense * z_elem
         z_eta = D_eta_dense * z_elem
         z_zeta = D_zeta_dense * z_elem

src/fem2d_P2.jl

@@ -134,7 +134,14 @@ legacy `geometric_mg(geom, L)` builds geometric-subdivision transfers instead.)
 # Arguments
 - `K::Array{T,3}` (`7 × N × 2`): P2+bubble per-triangle mesh; the 7 vertices per
   triangle are laid out as
-  `corner1, midpt(1,2), corner2, midpt(2,3), corner3, midpt(3,1), centroid`.
+  `corner1, edge(1,2), corner2, edge(2,3), corner3, edge(3,1), centroid`.
+
+The element geometry is **isoparametric**: the map from the reference triangle is
+built from all 7 node positions via the P2+bubble shape functions, so displacing
+the edge/centroid nodes off the straight midpoints/barycenter genuinely curves the
+element (with a node-varying Jacobian). Triangles must be **orientation-preserving
+and non-self-intersecting** — construction errors if any element's `det J ≤ 0` at a
+quadrature node, since the barrier method requires strictly positive weights.
 """
 function fem2d_P2(::Type{T}=Float64;
                   K::Array{T,3} = _default_K7(T),
@@ -375,18 +382,42 @@ function _fem2d_P2_structured(::Type{T}, K::Array{T,3}, K7::Array{T,3}, L::Int)
     dy_block = zeros(T, p, p, N)
     w_vec = zeros(T, n)
 
+    # Isoparametric P2+bubble map: x(ξ,η) = Σ_i N_i(ξ,η) X_i over all p geometry
+    # nodes, so the Jacobian J = [∂x/∂ξ ∂x/∂η; ∂y/∂ξ ∂y/∂η] varies node-to-node
+    # whenever the element is curved (edge nodes off the midpoints, centroid off
+    # the barycenter). The reference derivative rows give ∂N_i/∂ξ and ∂N_i/∂η at
+    # each node, so (R_dx*X)[j] = ∂x/∂ξ at node j, etc. The physical-derivative
+    # operator rows are J(node j)^{-T} applied to (R_dx, R_dy), and the quadrature
+    # weight at node j is det J(node j) · R.w[j]. For a straight triangle J is
+    # constant and this reduces to the previous affine map.
     for k in 1:N
-        u = xL[:, 1, k] - xL[:, 5, k]
-        v = xL[:, 3, k] - xL[:, 5, k]
-        A = hcat(u, v)
-        invA = inv(A)'
-        dx_block[:, :, k] = invA[1, 1] * R_dx + invA[1, 2] * R_dy
-        dy_block[:, :, k] = invA[2, 1] * R_dx + invA[2, 2] * R_dy
+        X = @view xL[1, :, k]
+        Y = @view xL[2, :, k]
+        x_xi = R_dx * X; x_eta = R_dy * X
+        y_xi = R_dx * Y; y_eta = R_dy * Y
         for j in 1:p
+            detJ = x_xi[j] * y_eta[j] - x_eta[j] * y_xi[j]
+            invdet = one(T) / detJ
+            for m in 1:p
+                dx_block[j, m, k] = ( y_eta[j] * R_dx[j, m] - y_xi[j] * R_dy[j, m]) * invdet
+                dy_block[j, m, k] = (-x_eta[j] * R_dx[j, m] + x_xi[j] * R_dy[j, m]) * invdet
+            end
             id_block[j, j, k] = one(T)
+            w_vec[(k-1)*p + j] = detJ * R.w[j]
         end
-        w_blk = abs(det(A)) * R.w
-        w_vec[(k-1)*p+1:k*p] = w_blk
+    end
+
+    # The barrier method requires strictly positive quadrature weights. A
+    # non-positive weight means det J ≤ 0 at some node: the user-supplied curved
+    # element is folded or has inverted (clockwise) orientation. Refuse to solve.
+    if !all(>(zero(T)), w_vec)
+        bad      = findall(<=(zero(T)), w_vec)
+        badelems = sort!(unique((bad .- 1) .÷ p .+ 1))
+        error("fem2d_P2: non-positive quadrature weight at $(length(bad)) node(s) " *
+              "across $(length(badelems)) element(s) (first few: $(first(badelems, 5))). " *
+              "The isoparametric element map has det J ≤ 0 there — the curved triangle " *
+              "is folded or has inverted (clockwise) vertex orientation. Supply " *
+              "orientation-preserving, non-self-intersecting elements.")
     end
 
     id = BlockDiag(id_block)

src/fem3d.jl

@@ -47,6 +47,13 @@ transfers instead.)
   `K[v, e, d]` is coordinate `d` of Lagrange node `v` of hex `e` (tensor-product
   ordering: x fastest, then y, then z).
 
+The element geometry is **isoparametric**: the map from the reference hex is built
+from all `(k+1)^3` node positions via the Q_k Lagrange basis, so displacing the
+edge/face/interior nodes genuinely curves the element (node-varying Jacobian).
+Hexes must be **orientation-preserving and non-self-intersecting** — construction
+errors if any element's `det J ≤ 0` at a quadrature node, since the barrier method
+requires strictly positive weights.
+
 The Geometry is intended for Dirichlet boundary conditions.
 """
 function fem3d(::Type{T}=Float64;