feat: correct solid mechanics bc

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3970-15479-074f42a
+  baseline: integratedTests/baseline_integratedTests-pr3986-15734-7487221
 
 allow_fail:
   all: ''

BASELINE_NOTES.md

@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3970 (2026-02-27) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3986-15734-7487221.tar.gz>
+=====================
+Corrected traction boundary conditions
+
 PR #3970 (2026-02-11) <https://storage.googleapis.com/geosx/integratedTests/baseline_integratedTests-pr3970-15479-074f42a.tar.gz>
 =====================
 Bypass well residual calculation for closed wells

src/coreComponents/fieldSpecification/TractionBoundaryCondition.cpp

@@ -19,6 +19,8 @@
 
 #include "TractionBoundaryCondition.hpp"
 
+#include "finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp"
+#include "finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp"
 #include "functions/TableFunction.hpp"
 
 namespace geos
@@ -132,15 +134,191 @@ void TractionBoundaryCondition::initializePreSubGroups()
 }
 
 
+/**
+ * @brief Integrate a uniform traction over a face using a standard FE face formulation.
+ *
+ * @tparam FE_TYPE The finite element face type (e.g. H1_TriangleFace_Lagrange1_Gauss1_impl,
+ *                 H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl).
+ *                 Must provide: numNodes, numQuadraturePoints, calcN(), transformedQuadratureWeight().
+ * @param traction The traction vector [3] to integrate.
+ * @param xFace    The face node coordinates in FE ordering [numNodes][3].
+ * @param permutation Mapping from FE node index to mesh face-local node index.
+ *                     Accounts for ordering differences (e.g. tensor-product vs CCW).
+ * @param faceToNodeMap The face-to-node connectivity map.
+ * @param kf       The face index.
+ * @param blockLocalDofNumber Block-local DOF numbering.
+ * @param dofRankOffset The rank offset for DOFs.
+ * @param localRhs The local RHS vector to assemble into.
+ */
+template< typename FE_TYPE >
+GEOS_HOST_DEVICE
+inline
+void integrateFaceTraction( real64 const ( &traction )[3],
+                            real64 const ( &xFace )[FE_TYPE::numNodes][3],
+                            int const ( &permutation )[FE_TYPE::numNodes],
+                            ArrayOfArraysView< localIndex const > const & faceToNodeMap,
+                            localIndex const kf,
+                            arrayView1d< globalIndex const > const & blockLocalDofNumber,
+                            globalIndex const dofRankOffset,
+                            arrayView1d< real64 > const & localRhs )
+{
+  constexpr localIndex numNodes = FE_TYPE::numNodes;
+  constexpr localIndex numQuadraturePoints = FE_TYPE::numQuadraturePoints;
+
+  // Accumulate per-node integration weight: w_a = sum_q N_a(q) * |J(q)| * wq
+  real64 nodalWeight[numNodes] {};
+
+  for( localIndex q = 0; q < numQuadraturePoints; ++q )
+  {
+    real64 N[numNodes];
+    FE_TYPE::calcN( q, N );
+
+    real64 const detJxW = FE_TYPE::transformedQuadratureWeight( q, xFace );
+
+    for( localIndex a = 0; a < numNodes; ++a )
+    {
+      nodalWeight[a] += N[a] * detJxW;
+    }
+  }
+
+  // Assemble traction * nodalWeight into the RHS.
+  // Use the permutation to map FE node 'a' back to the mesh face-local node index.
+  for( localIndex a = 0; a < numNodes; ++a )
+  {
+    localIndex const meshNodeIdx = permutation[a];
+    localIndex const dof = blockLocalDofNumber[ faceToNodeMap( kf, meshNodeIdx ) ] - dofRankOffset;
+    if( dof < 0 || dof >= localRhs.size() )
+      continue;
+    RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+0], traction[0] * nodalWeight[a] );
+    RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+1], traction[1] * nodalWeight[a] );
+    RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+2], traction[2] * nodalWeight[a] );
+  }
+}
+
+/**
+ * @brief Integrate a uniform traction over a general polygon face.
+ *
+ * This is a fallback for faces that are not triangles or quadrilaterals.
+ *
+ * @param traction The traction vector [3] to integrate.
+ * @param faceArea The area of the face.
+ * @param numNodes The number of nodes of the face.
+ * @param faceToNodeMap The face-to-node connectivity map.
+ * @param kf       The face index.
+ * @param blockLocalDofNumber Block-local DOF numbering.
+ * @param dofRankOffset The rank offset for DOFs.
+ * @param localRhs The local RHS vector to assemble into.
+ */
+GEOS_HOST_DEVICE
+inline
+void integrateFaceTractionPolygon( real64 const ( &traction )[3],
+                                   real64 const faceArea,
+                                   localIndex const numNodes,
+                                   ArrayOfArraysView< localIndex const > const & faceToNodeMap,
+                                   localIndex const kf,
+                                   arrayView1d< globalIndex const > const & blockLocalDofNumber,
+                                   globalIndex const dofRankOffset,
+                                   arrayView1d< real64 > const & localRhs )
+{
+  real64 const w = faceArea / numNodes;
+  for( localIndex a = 0; a < numNodes; ++a )
+  {
+    localIndex const dof = blockLocalDofNumber[ faceToNodeMap( kf, a ) ] - dofRankOffset;
+    if( dof < 0 || dof >= localRhs.size() )
+      continue;
+    RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+0], traction[0] * w );
+    RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+1], traction[1] * w );
+    RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+2], traction[2] * w );
+  }
+}
+
+/**
+ * @brief Dispatch face traction integration to the appropriate FE formulation.
+ *
+ * Selects the FE face type based on numNodes:
+ * - 3 nodes: H1_TriangleFace_Lagrange1_Gauss1_impl
+ * - 4 nodes: H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl
+ * - Other:   Equal distribution fallback
+ *
+ * @param traction The traction vector [3].
+ * @param numNodes The number of nodes of the face.
+ * @param faceArea The area of the face (only used for polygon fallback).
+ * @param faceToNodeMap The face-to-node connectivity map.
+ * @param kf       The face index.
+ * @param nodePositions The reference positions of all nodes.
+ * @param blockLocalDofNumber Block-local DOF numbering.
+ * @param dofRankOffset The rank offset for DOFs.
+ * @param localRhs The local RHS vector to assemble into.
+ */
+GEOS_HOST_DEVICE
+inline
+void assembleFaceTraction( real64 const ( &traction )[3],
+                           localIndex const numNodes,
+                           real64 const faceArea,
+                           ArrayOfArraysView< localIndex const > const & faceToNodeMap,
+                           localIndex const kf,
+                           arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePositions,
+                           arrayView1d< globalIndex const > const & blockLocalDofNumber,
+                           globalIndex const dofRankOffset,
+                           arrayView1d< real64 > const & localRhs )
+{
+  using TriFace = finiteElement::H1_TriangleFace_Lagrange1_Gauss1_impl;
+  using QuadFace = finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2_impl;
+
+  if( numNodes == 3 )
+  {
+    // Get the permutation mapping from FE node ordering to mesh face-local node ordering.
+    int permutation[TriFace::numNodes];
+    TriFace::getPermutation( permutation );
+
+    // Gather coordinates in FE ordering: xFace[feNode] = coords of mesh node permutation[feNode]
+    real64 xFace[TriFace::numNodes][3];
+    for( localIndex a = 0; a < TriFace::numNodes; ++a )
+    {
+      localIndex const nodeIdx = faceToNodeMap( kf, permutation[a] );
+      xFace[a][0] = nodePositions( nodeIdx, 0 );
+      xFace[a][1] = nodePositions( nodeIdx, 1 );
+      xFace[a][2] = nodePositions( nodeIdx, 2 );
+    }
+    integrateFaceTraction< TriFace >( traction, xFace, permutation, faceToNodeMap, kf,
+                                      blockLocalDofNumber, dofRankOffset, localRhs );
+  }
+  else if( numNodes == 4 )
+  {
+    // Get the permutation mapping from FE node ordering to mesh face-local node ordering.
+    // For quads: FE uses Z-ordering (tensor product), mesh uses CCW ordering.
+    // getPermutation returns {0, 1, 3, 2} meaning FE node 2 → mesh node 3 and vice versa.
+    int permutation[QuadFace::numNodes];
+    QuadFace::getPermutation( permutation );
+
+    // Gather coordinates in FE ordering: xFace[feNode] = coords of mesh node permutation[feNode]
+    real64 xFace[QuadFace::numNodes][3];
+    for( localIndex a = 0; a < QuadFace::numNodes; ++a )
+    {
+      localIndex const nodeIdx = faceToNodeMap( kf, permutation[a] );
+      xFace[a][0] = nodePositions( nodeIdx, 0 );
+      xFace[a][1] = nodePositions( nodeIdx, 1 );
+      xFace[a][2] = nodePositions( nodeIdx, 2 );
+    }
+    integrateFaceTraction< QuadFace >( traction, xFace, permutation, faceToNodeMap, kf,
+                                       blockLocalDofNumber, dofRankOffset, localRhs );
+  }
+  else
+  {
+    integrateFaceTractionPolygon( traction, faceArea, numNodes, faceToNodeMap, kf,
+                                  blockLocalDofNumber, dofRankOffset, localRhs );
+  }
+}
+
+
 void TractionBoundaryCondition::launch( real64 const time,
                                         arrayView1d< globalIndex const > const blockLocalDofNumber,
                                         globalIndex const dofRankOffset,
                                         FaceManager const & faceManager,
+                                        arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePositions,
                                         SortedArrayView< localIndex const > const & targetSet,
                                         arrayView1d< real64 > const & localRhs ) const
 {
-
-
   arrayView1d< real64 const > const faceArea  = faceManager.faceArea();
   arrayView2d< real64 const > const faceNormal  = faceManager.faceNormal();
   ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst();
@@ -193,11 +371,9 @@ void TractionBoundaryCondition::launch( real64 const time,
 
       real64 const nodalScale = ( nodalScaleFlag == 1 ) ? scaleSet[i] : 1.0;
 
-      // TODO consider dispatch if appropriate
       real64 const tractionMagnitude = spatialFunction ? tractionMagnitudeArrayView[i] : (tractionMagnitude0 * nodalScale);
 
       real64 traction[3] = { 0 };
-      // TODO consider dispatch if appropriate
       if( tractionType == TractionType::vector )
       {
         traction[0] = tractionMagnitude * direction[0];
@@ -222,23 +398,11 @@ void TractionBoundaryCondition::launch( real64 const time,
           traction[1] = inputStress[5] * temp[0] + inputStress[1] * temp[1] + inputStress[3] * temp[2];
           traction[2] = inputStress[4] * temp[0] + inputStress[3] * temp[1] + inputStress[2] * temp[2];
         }
-
       }
 
-      // TODO replace with proper FEM integration
-      traction[0] *= faceArea[kf] / numNodes;
-      traction[1] *= faceArea[kf] / numNodes;
-      traction[2] *= faceArea[kf] / numNodes;
-
-      for( localIndex a=0; a<numNodes; ++a )
-      {
-        localIndex const dof = blockLocalDofNumber[ faceToNodeMap( kf, a ) ] - dofRankOffset;
-        if( dof < 0 || dof >= localRhs.size() )
-          continue;
-        RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+0], traction[0] );
-        RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+1], traction[1] );
-        RAJA::atomicAdd< parallelDeviceAtomic >( &localRhs[dof+2], traction[2] );
-      }
+      // Dispatch face traction integration to the appropriate FE formulation
+      assembleFaceTraction( traction, numNodes, faceArea[kf], faceToNodeMap, kf,
+                            nodePositions, blockLocalDofNumber, dofRankOffset, localRhs );
     } );
   }
 }

src/coreComponents/fieldSpecification/TractionBoundaryCondition.hpp

@@ -23,6 +23,7 @@
 
 #include "FieldSpecificationBase.hpp"
 #include "mesh/FaceManager.hpp"
+#include "mesh/NodeManager.hpp"
 
 namespace geos
 {
@@ -71,13 +72,15 @@ class TractionBoundaryCondition : public FieldSpecificationBase
    * @param blockLocalDofNumber Array of block local DOF numbers for the displacement.
    * @param dofRankOffset The rank offset for the DOF.
    * @param faceManager Reference to the face manager (Tractions are applied on faces)
+   * @param nodePositions The reference position of all nodes (used for proper FEM integration on faces)
    * @param targetSet The set of faces to apply the BC to.
    * @param localRhs The RHS of the system to add contributions to.
    */
   void launch( real64 const time,
                arrayView1d< globalIndex const > const blockLocalDofNumber,
                globalIndex const dofRankOffset,
                FaceManager const & faceManager,
+               arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePositions,
                SortedArrayView< localIndex const > const & targetSet,
                arrayView1d< real64 > const & localRhs ) const;
 

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp

@@ -499,6 +499,18 @@ real64 SolidMechanicsLagrangianFEM::solverStep( real64 const & time_n,
   {
     int const maxNumResolves = m_maxNumResolves;
     int globallyFractured = 0;
+
+    // Check if mesh was modified by an external event (e.g., SurfaceGen) before this solver step
+    // This ensures setupSystem is called before implicitStepSetup when topology changes occur
+    Timestamp const initialMeshModificationTimestamp = getMeshModificationTimestamp( domain );
+    MeshLevel & meshLevel = domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName );
+    if( initialMeshModificationTimestamp > getSystemSetupTimestamp() || meshLevel.getModificationTimestamp() > getSystemSetupTimestamp() )
+    {
+      m_dofManager.clear();
+      setupSystem( domain, m_dofManager, m_localMatrix, m_rhs, m_solution, true );
+      setSystemSetupTimestamp( LvArray::math::max( initialMeshModificationTimestamp, meshLevel.getModificationTimestamp() ) );
+    }
+
     implicitStepSetup( time_n, dt, domain );
     for( int solveIter=0; solveIter<maxNumResolves+1; ++solveIter )
     {
@@ -507,10 +519,10 @@ real64 SolidMechanicsLagrangianFEM::solverStep( real64 const & time_n,
       Timestamp const meshModificationTimestamp = getMeshModificationTimestamp( domain );
 
       // Only build the sparsity pattern if the mesh has changed
-      if( meshModificationTimestamp > getSystemSetupTimestamp() || globallyFractured )
+      if( meshModificationTimestamp > getSystemSetupTimestamp() || meshLevel.getModificationTimestamp() > getSystemSetupTimestamp() || globallyFractured )
       {
         setupSystem( domain, m_dofManager, m_localMatrix, m_rhs, m_solution );
-        setSystemSetupTimestamp( meshModificationTimestamp );
+        setSystemSetupTimestamp( LvArray::math::max( meshModificationTimestamp, meshLevel.getModificationTimestamp() ) );
       }
 
       dtReturn = nonlinearImplicitStep( time_n,
@@ -798,6 +810,7 @@ void SolidMechanicsLagrangianFEM::applyTractionBC( real64 const time,
                  blockLocalDofNumber,
                  dofRankOffset,
                  faceManager,
+                 nodeManager.referencePosition(),
                  targetSet,
                  localRhs );
     } );