feat: Configurable Quadrature Rules for triangles/quads in ALM solver

src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp

@@ -35,6 +35,7 @@
 #include "constitutive/contact/FrictionSelector.hpp"
 #include "constitutive/solid/PorousSolid.hpp"
 #include "constitutive/solid/SolidFields.hpp"
+#include "finiteElement/FiniteElementDiscretization.hpp"
 #include "mesh/DomainPartition.hpp"
 
 namespace geos
@@ -48,10 +49,6 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta
                                                                                     Group * const parent ):
   ContactSolverBase( name, parent )
 {
-
-  m_faceTypeToFiniteElements.insert( {"Quadrilateral", std::make_unique< finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2 >()} );
-  m_faceTypeToFiniteElements.insert( {"Triangle", std::make_unique< finiteElement::H1_TriangleFace_Lagrange1_Gauss1 >()} );
-
   registerWrapper( viewKeyStruct::simultaneousString(), &m_simultaneous ).
     setInputFlag( InputFlags::OPTIONAL ).
     setApplyDefaultValue( 1 ).
@@ -167,6 +164,54 @@ void SolidMechanicsAugmentedLagrangianContact::registerDataOnMesh( dataRepositor
 
 }
 
+void SolidMechanicsAugmentedLagrangianContact::initializePostInitialConditionsPreSubGroups()
+{
+  ContactSolverBase::initializePostInitialConditionsPreSubGroups();
+
+  DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+  validateTetrahedralQuadrature( domain.getMeshBodies() );
+}
+
+void SolidMechanicsAugmentedLagrangianContact::validateTetrahedralQuadrature( Group & meshBodies )
+{
+  string const discretizationName = getDiscretizationName();
+
+  NumericalMethodsManager const & numericalMethodManager =
+    this->getGroupByPath< DomainPartition >( "/Problem/domain" ).getNumericalMethodManager();
+  FiniteElementDiscretizationManager const & feDiscretizationManager =
+    numericalMethodManager.getFiniteElementDiscretizationManager();
+  FiniteElementDiscretization const & feDiscretization =
+    feDiscretizationManager.getGroup< FiniteElementDiscretization >( discretizationName );
+
+
+  integer const useHighOrderQuadrature =
+    feDiscretization.getReference< integer >( "useHighOrderQuadratureRule" );
+
+  bool hasTetrahedra = false;
+  forDiscretizationOnMeshTargets( meshBodies, [&]( string const &,
+                                                   MeshLevel const & mesh,
+                                                   string_array const & regionNames )
+  {
+    ElementRegionManager const & elemManager = mesh.getElemManager();
+    elemManager.forElementRegions< CellElementRegion >( regionNames, [&]( localIndex const,
+                                                                          CellElementRegion const & region )
+    {
+      region.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion const & subRegion )
+      {
+        if( subRegion.getElementType() == ElementType::Tetrahedron )
+        {
+          hasTetrahedra = true;
+        }
+      } );
+    } );
+  } );
+
+  GEOS_ERROR_IF( hasTetrahedra && useHighOrderQuadrature != 1,
+                 GEOS_FMT( "{}: Tetrahedral meshes require useHighOrderQuadratureRule=\"1\" for correct integration of bubble contributions. "
+                           "Please add this attribute to your FiniteElements/{} XML block.",
+                           getName(), discretizationName ) );
+}
+
 void SolidMechanicsAugmentedLagrangianContact::setupDofs( DomainPartition const & domain,
                                                           DofManager & dofManager ) const
 {
@@ -219,6 +264,29 @@ void SolidMechanicsAugmentedLagrangianContact::setupSystem( DomainPartition & do
   PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity );
 }
 
+void SolidMechanicsAugmentedLagrangianContact::postInputInitialization()
+{
+  ContactSolverBase::postInputInitialization();
+
+  DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
+  FiniteElementDiscretizationManager const & feDiscretizationManager =
+    numericalMethodManager.getFiniteElementDiscretizationManager();
+  FiniteElementDiscretization const & feDiscretization =
+    feDiscretizationManager.getGroup< FiniteElementDiscretization >( getDiscretizationName() );
+
+  m_faceTypeToFiniteElements.insert( {"Quadrilateral", feDiscretization.factory( ElementType::Quadrilateral )} );
+  m_faceTypeToFiniteElements.insert( {"Triangle", feDiscretization.factory( ElementType::Triangle )} );
+
+  GEOS_LOG_RANK_0( GEOS_FMT( "{} using finite element discretization {}:",
+                             getName(), getDiscretizationName() ) );
+  for( auto const & [name, fePtr] : m_faceTypeToFiniteElements )
+  {
+    GEOS_LOG_RANK_0( GEOS_FMT( "  {} face elements: {} quadrature points",
+                               name, fePtr->getNumQuadraturePoints() ) );
+  }
+}
+
 void SolidMechanicsAugmentedLagrangianContact::setSparsityPattern( DomainPartition & domain,
                                                                    DofManager & dofManager,
                                                                    CRSMatrix< real64, globalIndex > & GEOS_UNUSED_PARAM( localMatrix ),

src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp

@@ -47,8 +47,12 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase
    */
   string getCatalogName() const override { return catalogName(); }
 
+  virtual void postInputInitialization() override;
+
   virtual void registerDataOnMesh( dataRepository::Group & meshBodies ) override final;
 
+  virtual void initializePostInitialConditionsPreSubGroups() override;
+
   virtual void setupDofs( DomainPartition const & domain,
                           DofManager & dofManager ) const override;
 
@@ -213,6 +217,12 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase
 
 private:
 
+  /**
+   * @brief Validate that tetrahedral meshes use high-order quadrature rules
+   * @param meshBodies the group containing the mesh bodies
+   */
+  void validateTetrahedralQuadrature( Group & meshBodies );
+
   /**
    * @brief add the number of non-zero elements induced by the coupling between
    *   nodal and bubble displacement.

src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.cpp

@@ -27,7 +27,7 @@
 #include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsContactFaceBubbleKernels.hpp"
 #include "physicsSolvers/solidMechanics/contact/LogLevelsInfo.hpp"
 #include "physicsSolvers/LogLevelsInfo.hpp"
-
+#include "finiteElement/FiniteElementDiscretization.hpp"
 #include "constitutive/contact/FrictionSelector.hpp"
 #include "fieldSpecification/FieldSpecificationManager.hpp"
 
@@ -43,9 +43,6 @@ SolidMechanicsLagrangeContactBubbleStab::SolidMechanicsLagrangeContactBubbleStab
                                                                                   Group * const parent ):
   ContactSolverBase( name, parent )
 {
-  m_faceTypeToFiniteElements.insert( {"Quadrilateral", std::make_unique< finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2 >()} );
-  m_faceTypeToFiniteElements.insert( {"Triangle", std::make_unique< finiteElement::H1_TriangleFace_Lagrange1_Gauss1 >()} );
-
   LinearSolverParameters & linSolParams = m_linearSolverParameters.get();
   linSolParams.mgr.strategy = LinearSolverParameters::MGR::StrategyType::lagrangianContactMechanicsBubbleStab;
   linSolParams.mgr.separateComponents = true;
@@ -134,6 +131,54 @@ void SolidMechanicsLagrangeContactBubbleStab::registerDataOnMesh( Group & meshBo
   } );
 }
 
+void SolidMechanicsLagrangeContactBubbleStab::initializePostInitialConditionsPreSubGroups()
+{
+  ContactSolverBase::initializePostInitialConditionsPreSubGroups();
+
+  DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+  validateTetrahedralQuadrature( domain.getMeshBodies() );
+}
+
+void SolidMechanicsLagrangeContactBubbleStab::validateTetrahedralQuadrature( Group & meshBodies )
+{
+  string const discretizationName = getDiscretizationName();
+
+  NumericalMethodsManager const & numericalMethodManager =
+    this->getGroupByPath< DomainPartition >( "/Problem/domain" ).getNumericalMethodManager();
+  FiniteElementDiscretizationManager const & feDiscretizationManager =
+    numericalMethodManager.getFiniteElementDiscretizationManager();
+  FiniteElementDiscretization const & feDiscretization =
+    feDiscretizationManager.getGroup< FiniteElementDiscretization >( discretizationName );
+
+
+  integer const useHighOrderQuadrature =
+    feDiscretization.getReference< integer >( "useHighOrderQuadratureRule" );
+
+  bool hasTetrahedra = false;
+  forDiscretizationOnMeshTargets( meshBodies, [&]( string const &,
+                                                   MeshLevel const & mesh,
+                                                   string_array const & regionNames )
+  {
+    ElementRegionManager const & elemManager = mesh.getElemManager();
+    elemManager.forElementRegions< CellElementRegion >( regionNames, [&]( localIndex const,
+                                                                          CellElementRegion const & region )
+    {
+      region.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion const & subRegion )
+      {
+        if( subRegion.getElementType() == ElementType::Tetrahedron )
+        {
+          hasTetrahedra = true;
+        }
+      } );
+    } );
+  } );
+
+  GEOS_ERROR_IF( hasTetrahedra && useHighOrderQuadrature != 1,
+                 GEOS_FMT( "{}: Tetrahedral meshes require useHighOrderQuadratureRule=\"1\" for correct integration of bubble contributions. "
+                           "Please add this attribute to your FiniteElements/{} XML block.",
+                           getName(), discretizationName ) );
+}
+
 void SolidMechanicsLagrangeContactBubbleStab::setupDofs( DomainPartition const & domain,
                                                          DofManager & dofManager ) const
 {
@@ -192,6 +237,29 @@ void SolidMechanicsLagrangeContactBubbleStab::setupDofs( DomainPartition const &
                           meshTargets );
 }
 
+void SolidMechanicsLagrangeContactBubbleStab::postInputInitialization()
+{
+  ContactSolverBase::postInputInitialization();
+
+  DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
+  FiniteElementDiscretizationManager const & feDiscretizationManager =
+    numericalMethodManager.getFiniteElementDiscretizationManager();
+  FiniteElementDiscretization const & feDiscretization =
+    feDiscretizationManager.getGroup< FiniteElementDiscretization >( getDiscretizationName() );
+
+  m_faceTypeToFiniteElements.insert( {"Quadrilateral", feDiscretization.factory( ElementType::Quadrilateral )} );
+  m_faceTypeToFiniteElements.insert( {"Triangle", feDiscretization.factory( ElementType::Triangle )} );
+
+  GEOS_LOG_RANK_0( GEOS_FMT( "{} using finite element discretization {}:",
+                             getName(), getDiscretizationName() ) );
+  for( auto const & [name, fePtr] : m_faceTypeToFiniteElements )
+  {
+    GEOS_LOG_RANK_0( GEOS_FMT( "  {} face elements: {} quadrature points",
+                               name, fePtr->getNumQuadraturePoints() ) );
+  }
+}
+
 void SolidMechanicsLagrangeContactBubbleStab::setupSystem( DomainPartition & domain,
                                                            DofManager & dofManager,
                                                            CRSMatrix< real64, globalIndex > & localMatrix,

src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContactBubbleStab.hpp

@@ -50,6 +50,10 @@ class SolidMechanicsLagrangeContactBubbleStab : public ContactSolverBase
    */
   string getCatalogName() const override { return catalogName(); }
 
+  virtual void postInputInitialization() override;
+
+  virtual void initializePostInitialConditionsPreSubGroups() override;
+
   virtual void registerDataOnMesh( Group & MeshBodies ) override final;
 
   real64 solverStep( real64 const & time_n,
@@ -200,6 +204,12 @@ class SolidMechanicsLagrangeContactBubbleStab : public ContactSolverBase
 
 
 private:
+  /**
+   * @brief Validate that tetrahedral meshes use high-order quadrature rules
+   * @param meshBodies the group containing the mesh bodies
+   */
+  void validateTetrahedralQuadrature( Group & meshBodies );
+
   /**
    * @brief add the number of non-zero elements induced by the coupling between
    *   nodal and bubble displacement.