From 5eb4079b13840b5680f73cd7cd9a60c5fd65f8d8 Mon Sep 17 00:00:00 2001 From: jafranc Date: Tue, 21 Jul 2026 06:48:05 -0500 Subject: [PATCH 1/3] MACRO for tmp rather than duplication --- .../multiphysics/MultiphasePoromechanics.hpp | 6 +----- .../physicsSolvers/multiphysics/PoromechanicsSolver.hpp | 6 ++++++ .../multiphysics/SinglePhasePoromechanics.hpp | 8 ++------ .../SinglePhasePoromechanicsConformingFractures.hpp | 7 +------ .../SinglePhasePoromechanicsConformingFracturesALM.hpp | 6 +----- .../solidMechanics/contact/ContactSolverBase.hpp | 2 ++ 6 files changed, 13 insertions(+), 22 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp index de5de39883a..bf06083e779 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp @@ -101,11 +101,7 @@ class MultiphasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHANI protected: - virtual void setMGRStrategy() override - { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); - } + GEOS_MGR_STRATEGY_NOT_SUPPORTED() /** * @brief Helper function to recompute the bulk density diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp index 1b34ab94f65..1b14e979478 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp @@ -50,6 +50,12 @@ ENUM_STRINGS( StabilizationType, "Local" ); } +#define GEOS_MGR_STRATEGY_NOT_SUPPORTED() \ + virtual void setMGRStrategy() override \ + { \ + if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) \ + GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName() ) ); \ + } template< typename FLOW_SOLVER, typename MECHANICS_SOLVER = SolidMechanicsLagrangianFEM > class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER > diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp index d7e40171f79..288497421ff 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp @@ -116,12 +116,8 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN virtual void initializePostInitialConditionsPreSubGroups() override; - virtual void setMGRStrategy() override - { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); - } - + GEOS_MGR_STRATEGY_NOT_SUPPORTED() + virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override { GEOS_MARK_FUNCTION; diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp index f0ddd9bc8f1..dcbae2d46a2 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFractures.hpp @@ -81,12 +81,7 @@ class SinglePhasePoromechanicsConformingFractures : public PoromechanicsConformi * These functions provide the primary interface that is required for derived classes */ /**@{*/ - - virtual void setMGRStrategy() override - { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); - } + GEOS_MGR_STRATEGY_NOT_SUPPORTED() virtual void assembleSystem( real64 const time_n, real64 const dt, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index 4267619df0a..b4b8b8465a2 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -101,11 +101,7 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec virtual void updateState( DomainPartition & domain ) override final; - virtual void setMGRStrategy() override final - { - if( this->m_linearSolverParameters.get().preconditionerType == LinearSolverParameters::PreconditionerType::mgr ) - GEOS_ERROR( GEOS_FMT( "{}: MGR strategy is not implemented for {}", this->getName(), this->getCatalogName())); - } + GEOS_MGR_STRATEGY_NOT_SUPPORTED() /**@}*/ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp index ccfc68f6ea7..9423d8c8f2c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp @@ -59,6 +59,8 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM }; + static constexpr integer MGR_READINESS = 0; + protected: virtual void postInputInitialization() override; From bd817678c10c62bd233e88fb6f10fd39ac8ccaca Mon Sep 17 00:00:00 2001 From: jafranc Date: Tue, 21 Jul 2026 07:30:59 -0500 Subject: [PATCH 2/3] Make sequential only visible (for ALM) --- .../SinglePhasePoromechanicsConformingFracturesALM.cpp | 5 +++++ .../SinglePhasePoromechanicsConformingFracturesALM.hpp | 5 +++++ .../solidMechanics/contact/ContactSolverBase.hpp | 2 -- 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index ebc5881ee2f..863c2ce4cb2 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -40,7 +40,12 @@ template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::postInputInitialization() { Base::postInputInitialization(); + forceSequential(); +} +template< typename FLOW_SOLVER > +void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::forceSequential() +{ GEOS_WARNING_IF( this->getNonlinearSolverParameters().couplingType() == NonlinearSolverParameters::CouplingType::FullyImplicit, "FullyImplicit coupling not implemented for this solver. A sequential coupling approach will be used." ); diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index b4b8b8465a2..d62524a420a 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -179,6 +179,11 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec DofManager const & dofManager, CRSMatrix< real64, globalIndex > & localMatrix ); + /** + * @brief Force sequential mode with warning on explicit fully-implicit - temporary + */ + void forceSequential(); + virtual void mapSolutionBetweenSolvers( DomainPartition & domain, integer const solverType ) override; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp index 9423d8c8f2c..ccfc68f6ea7 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp @@ -59,8 +59,6 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM }; - static constexpr integer MGR_READINESS = 0; - protected: virtual void postInputInitialization() override; From 2a4502d84fe34950f6709489427a1d96f698dc06 Mon Sep 17 00:00:00 2001 From: jafranc Date: Thu, 23 Jul 2026 03:30:22 -0500 Subject: [PATCH 3/3] fit ALM in PM hierarchy: - normalizing ALMSinglePhase to be a PoromechanicsConformingFractures and not a SinglePhasePoromechanics anymore (as done for LM) - moving up `faceNodalArea` and `faceIntegral` to ContactSolverBase - moving up constexpr maxNumNodes - temp patching `hasStabilization/getStabilizationName` while waiting for bubble to reframe as Stabilization --- .../PoromechanicsConformingFractures.hpp | 14 +- ...asePoromechanicsConformingFracturesALM.cpp | 27 -- ...asePoromechanicsConformingFracturesALM.hpp | 30 +- .../contact/ContactSolverBase.cpp | 328 +++++++++++++++++ .../contact/ContactSolverBase.hpp | 26 ++ ...lidMechanicsAugmentedLagrangianContact.hpp | 4 + .../contact/SolidMechanicsLagrangeContact.cpp | 345 +----------------- .../contact/SolidMechanicsLagrangeContact.hpp | 30 +- 8 files changed, 406 insertions(+), 398 deletions(-) diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp index 9e0828e19d5..2024d2dbec1 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp @@ -37,11 +37,11 @@ namespace geos { -template< template< typename, typename > class POROMECHANICS_BASE, typename FLOW_SOLVER > -class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, SolidMechanicsLagrangeContact > +template< template< typename, typename > class POROMECHANICS_BASE, typename FLOW_SOLVER , typename CONTACT_SOLVER = SolidMechanicsLagrangeContact > +class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, CONTACT_SOLVER > { public: - using Base = POROMECHANICS_BASE< FLOW_SOLVER, SolidMechanicsLagrangeContact >; + using Base = POROMECHANICS_BASE< FLOW_SOLVER, CONTACT_SOLVER >; PoromechanicsConformingFractures( const string & name, dataRepository::Group * const parent ) @@ -181,6 +181,11 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager(); FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager(); + + //TODO (jafranc) - remove once ALM-bubble is frame as a stab method - tmp runtime is fine as it is tmp + if(this->solidMechanicsSolver()->hasStabilization()) + { + FluxApproximationBase const & stabilizationMethod = fvManager.getFluxApproximation( this->solidMechanicsSolver()->getStabilizationName() ); stabilizationMethod.forStencils< SurfaceElementStencil >( mesh, [&]( SurfaceElementStencil const & stencil ) @@ -224,6 +229,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, } } } ); + } } ); } @@ -431,7 +437,7 @@ class PoromechanicsConformingFractures : public POROMECHANICS_BASE< FLOW_SOLVER, } ); } ); - this->solidMechanicsSolver()->assembleContact( domain, dofManager, localMatrix, localRhs ); + this->solidMechanicsSolver()->assembleContact( time_n, dt, domain, dofManager, localMatrix, localRhs ); } virtual void assembleCouplingTerms( real64 const time_n, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp index 863c2ce4cb2..a253d104f52 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.cpp @@ -52,20 +52,6 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::forceSequent this->getNonlinearSolverParameters().m_couplingType = NonlinearSolverParameters::CouplingType::Sequential; } -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupCoupling( DomainPartition const & domain, - DofManager & dofManager ) const -{ - GEOS_MARK_FUNCTION; - - Base::setupCoupling( domain, dofManager ); - - dofManager.addCoupling( this->getFlowDofKey(), - fields::contact::traction::key(), - DofManager::Connector::Elem ); - -} - template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::setupSystem( DomainPartition & domain, @@ -142,19 +128,6 @@ void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::assembleCoup } } -template< typename FLOW_SOLVER > -void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >::updateState( DomainPartition & domain ) -{ - GEOS_MARK_FUNCTION; - - Base::updateState( domain ); - this->solidMechanicsSolver()->updateState( domain ); - - this->flowSolver()->prepareStencilWeights( domain ); - updateHydraulicApertureAndFracturePermeability( domain ); - this->flowSolver()->updateStencilWeights( domain ); -} - template< typename FLOW_SOLVER > void SinglePhasePoromechanicsConformingFracturesALM< FLOW_SOLVER >:: setUpDflux_dApertureMatrix( DomainPartition & domain, diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp index d62524a420a..19a2e643b02 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp +++ b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsConformingFracturesALM.hpp @@ -20,18 +20,21 @@ #ifndef GEOS_PHYSICSSOLVERS_MULTIPHYSICS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ #define GEOS_PHYSICSSOLVERS_MULTIPHYSICS_SINGLEPHASEPOROMECHANICSCONFORMINGFRACTURESALM_HPP_ -#include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" +#include "common/logger/Logger.hpp" +#include "physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFracturesALM.hpp" #include "physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp" +#include "physicsSolvers/multiphysics/PoromechanicsConformingFractures.hpp" +#include "physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp" namespace geos { template< typename FLOW_SOLVER = SinglePhaseBase > -class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromechanics< FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact > +class SinglePhasePoromechanicsConformingFracturesALM : public PoromechanicsConformingFractures< SinglePhasePoromechanics, FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact > { public: - using Base = SinglePhasePoromechanics< FLOW_SOLVER, SolidMechanicsAugmentedLagrangianContact >; + using Base = PoromechanicsConformingFractures< SinglePhasePoromechanics, FLOW_SOLVER , SolidMechanicsAugmentedLagrangianContact >; using Base::m_solvers; using Base::m_dofManager; using Base::m_localMatrix; @@ -82,8 +85,8 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec */ /**@{*/ - virtual void setupCoupling( DomainPartition const & domain, - DofManager & dofManager ) const override final; + // virtual void setupCoupling( DomainPartition const & domain, + // DofManager & dofManager ) const override final; virtual void setupSystem( DomainPartition & domain, DofManager & dofManager, @@ -99,12 +102,13 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) override final; - virtual void updateState( DomainPartition & domain ) override final; - - GEOS_MGR_STRATEGY_NOT_SUPPORTED() + // GEOS_MGR_STRATEGY_NOT_SUPPORTED() /**@}*/ + + + private: struct viewKeyStruct : public Base::viewKeyStruct @@ -212,6 +216,16 @@ class SinglePhasePoromechanicsConformingFracturesALM : public SinglePhasePoromec virtual void postInputInitialization() override final; + //- to be overloaded. + virtual void assembleFluidMassResidualDerivativeWrtDisplacement( MeshLevel const & GEOS_UNUSED_PARAM(mesh), + string_array const & GEOS_UNUSED_PARAM(regionNames), + DofManager const & GEOS_UNUSED_PARAM(dofManager), + CRSMatrixView< real64, globalIndex const > const & GEOS_UNUSED_PARAM(localMatrix), + arrayView1d< real64 > const & GEOS_UNUSED_PARAM(localRhs) ) override + { GEOS_ERROR("Sequential implementation only."); }; + + virtual integer numFluidComponents() const { return 1; } + }; } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp index 3faa6bdff59..c3384cb2bbd 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp @@ -27,6 +27,8 @@ #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" +#include "finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp" +#include "finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp" namespace geos { @@ -34,6 +36,8 @@ namespace geos using namespace dataRepository; using namespace constitutive; using namespace fields; +using namespace finiteElement; + ContactSolverBase::ContactSolverBase( const string & name, Group * const parent ): @@ -256,4 +260,328 @@ void ContactSolverBase::setConstitutiveNamesCallSuper( ElementSubRegionBase & su } } +void ContactSolverBase::computeFaceNodalArea( localIndex const kf0, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodeMap, + ArrayOfArraysView< localIndex const > const & faceToEdgeMap, + arrayView2d< localIndex const > const & edgeToNodeMap, + arrayView2d< real64 const > const faceCenters, + arrayView2d< real64 const > const faceNormals, + arrayView1d< real64 const > const faceAreas, + stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals ) const +{ + GEOS_MARK_FUNCTION; + localIndex const TriangularPermutation[3] = { 0, 1, 2 }; + localIndex const QuadrilateralPermutation[4] = { 0, 1, 3, 2 }; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); + + basisIntegrals.resize( numNodesPerFace ); + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + basisIntegrals[a] = 0.0; + } + localIndex const * const permutation = ( numNodesPerFace == 3 ) ? TriangularPermutation : QuadrilateralPermutation; + if( numNodesPerFace == 3 ) + { + real64 xLocal[3][3]; + for( localIndex a = 0; a < numNodesPerFace; ++a ) + { + for( localIndex j = 0; j < 3; ++j ) + { + xLocal[a][j] = nodePosition[faceToNodeMap( kf0, permutation[a] )][j]; + } + } + real64 N[3]; + for( localIndex q=0; q 4 && numNodesPerFace <= m_maxFaceNodes ) + { + // we need to L2 projector based on VEM to approximate the quadrature weights + // we need to use extra geometry information to computing L2 projector + + localIndex const MFN = m_maxFaceNodes; // Max number of face vertices. + localIndex const faceIndex = kf0; + localIndex const numFaceNodes = faceToNodeMap[ faceIndex ].size(); + + // get the face center and normal. + real64 const faceArea = faceAreas[ faceIndex ]; + localIndex faceToNodes[ MFN ]; + localIndex faceToEdges[ MFN ]; + for( localIndex i = 0; i < numFaceNodes; ++i ) + { + faceToNodes[i] = faceToNodeMap[ faceIndex ][ i ]; + faceToEdges[i] = faceToEdgeMap[ faceIndex ][ i ]; + } + // - get outward face normal and center + real64 faceNormal[3] = { faceNormals[faceIndex][0], + faceNormals[faceIndex][1], + faceNormals[faceIndex][2] }; + real64 const faceCenter[3] { faceCenters[faceIndex][0], + faceCenters[faceIndex][1], + faceCenters[faceIndex][2] }; + // - compute integrals calling auxiliary method + real64 threeDMonomialIntegrals[3] = { 0.0 }; + real64 const invCellDiameter = 0.0; + real64 const cellCenter[3] { 0.0, 0.0, 0.0 }; + computeFaceIntegrals( nodePosition, + faceToNodes, + faceToEdges, + numFaceNodes, + faceArea, + faceCenter, + faceNormal, + edgeToNodeMap, + invCellDiameter, + cellCenter, + basisIntegrals, + threeDMonomialIntegrals ); + } + else + { + GEOS_ERROR( GEOS_FMT( "Face with {} nodes. Only triangles and quadrilaterals and PEBI prisms up to 11 sides are supported.", + numNodesPerFace ), + getDataContext() ); + } +} + +void ContactSolverBase::computeFaceIntegrals( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoords, + localIndex const (&faceToNodes)[11], + localIndex const (&faceToEdges)[11], + localIndex const & numFaceVertices, + real64 const & faceArea, + real64 const (&faceCenter)[3], + real64 const (&faceNormal)[3], + arrayView2d< localIndex const > const & edgeToNodes, + real64 const & invCellDiameter, + real64 const (&cellCenter)[3], + stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals, + real64 (& threeDMonomialIntegrals)[3] ) const +{ + GEOS_MARK_FUNCTION; + localIndex const MFN = m_maxFaceNodes; // Max number of face vertices. + basisIntegrals.resize( numFaceVertices ); + // Rotate the face. + // - compute rotation matrix. + real64 faceRotationMatrix[ 3 ][ 3 ]; + computationalGeometry::RotationMatrix_3D( faceNormal, faceRotationMatrix ); + // - below we compute the diameter, the rotated vertices and the rotated center. + real64 faceRotatedVertices[ MFN ][ 2 ]; + real64 faceDiameter = 0; + + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + // apply the transpose (that is the inverse) of the rotation matrix to face vertices. + // NOTE: + // the second and third rows of the transpose of the rotation matrix rotate on the 2D face. + faceRotatedVertices[numVertex][0] = + faceRotationMatrix[ 0 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 0 ) + + faceRotationMatrix[ 1 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 1 ) + + faceRotationMatrix[ 2 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 2 ); + faceRotatedVertices[numVertex][1] = + faceRotationMatrix[ 0 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 0 ) + + faceRotationMatrix[ 1 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 1 ) + + faceRotationMatrix[ 2 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 2 ); + } + + faceDiameter = computationalGeometry::computeDiameter< 2 >( faceRotatedVertices, + numFaceVertices ); + real64 const invFaceDiameter = 1.0/faceDiameter; + // - rotate the face centroid as done for the vertices. + real64 faceRotatedCentroid[2]; + faceRotatedCentroid[0] = + faceRotationMatrix[ 0 ][ 1 ]*faceCenter[0] + + faceRotationMatrix[ 1 ][ 1 ]*faceCenter[1] + + faceRotationMatrix[ 2 ][ 1 ]*faceCenter[2]; + faceRotatedCentroid[1] = + faceRotationMatrix[ 0 ][ 2 ]*faceCenter[0] + + faceRotationMatrix[ 1 ][ 2 ]*faceCenter[1] + + faceRotationMatrix[ 2 ][ 2 ]*faceCenter[2]; + // - compute edges' lengths, outward pointing normals and local edge-to-nodes map. + real64 edgeOutwardNormals[ MFN ][ 2 ]; + real64 edgeLengths[ MFN ]; + localIndex localEdgeToNodes[ MFN ][ 2 ]; + + for( localIndex numEdge = 0; numEdge < numFaceVertices; ++numEdge ) + { + if( edgeToNodes( faceToEdges[numEdge], 0 ) == faceToNodes[ numEdge ] ) + { + localEdgeToNodes[ numEdge ][ 0 ] = numEdge; + localEdgeToNodes[ numEdge ][ 1 ] = (numEdge+1)%numFaceVertices; + } + else + { + localEdgeToNodes[ numEdge ][ 0 ] = (numEdge+1)%numFaceVertices; + localEdgeToNodes[ numEdge ][ 1 ] = numEdge; + } + real64 edgeTangent[2]; + edgeTangent[0] = faceRotatedVertices[(numEdge+1)%numFaceVertices][0] - + faceRotatedVertices[numEdge][0]; + edgeTangent[1] = faceRotatedVertices[(numEdge+1)%numFaceVertices][1] - + faceRotatedVertices[numEdge][1]; + edgeOutwardNormals[numEdge][0] = edgeTangent[1]; + edgeOutwardNormals[numEdge][1] = -edgeTangent[0]; + real64 signTestVector[2]; + signTestVector[0] = faceRotatedVertices[numEdge][0] - faceRotatedCentroid[0]; + signTestVector[1] = faceRotatedVertices[numEdge][1] - faceRotatedCentroid[1]; + if( signTestVector[0]*edgeOutwardNormals[numEdge][0] + + signTestVector[1]*edgeOutwardNormals[numEdge][1] < 0 ) + { + edgeOutwardNormals[numEdge][0] = -edgeOutwardNormals[numEdge][0]; + edgeOutwardNormals[numEdge][1] = -edgeOutwardNormals[numEdge][1]; + } + edgeLengths[numEdge] = LvArray::math::sqrt< real64 >( edgeTangent[0]*edgeTangent[0] + + edgeTangent[1]*edgeTangent[1] ); + edgeOutwardNormals[numEdge][0] /= edgeLengths[numEdge]; + edgeOutwardNormals[numEdge][1] /= edgeLengths[numEdge]; + } + + // Compute boundary quadrature weights (also equal to the integrals of basis functions on the + // boundary). + real64 boundaryQuadratureWeights[ MFN ]; + for( localIndex numWeight = 0; numWeight < numFaceVertices; ++numWeight ) + boundaryQuadratureWeights[numWeight] = 0.0; + for( localIndex numEdge = 0; numEdge < numFaceVertices; ++numEdge ) + { + boundaryQuadratureWeights[ localEdgeToNodes[ numEdge ][ 0 ] ] += 0.5*edgeLengths[numEdge]; + boundaryQuadratureWeights[ localEdgeToNodes[ numEdge ][ 1 ] ] += 0.5*edgeLengths[numEdge]; + } + + // Compute scaled monomials' integrals on edges. + real64 monomBoundaryIntegrals[3] = { 0.0 }; + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + monomBoundaryIntegrals[0] += boundaryQuadratureWeights[ numVertex ]; + monomBoundaryIntegrals[1] += (faceRotatedVertices[ numVertex ][ 0 ] - faceRotatedCentroid[0]) * + invFaceDiameter*boundaryQuadratureWeights[ numVertex ]; + monomBoundaryIntegrals[2] += (faceRotatedVertices[ numVertex ][ 1 ] - faceRotatedCentroid[1]) * + invFaceDiameter*boundaryQuadratureWeights[ numVertex ]; + } + + // Compute non constant 2D and 3D scaled monomials' integrals on the face. + real64 monomInternalIntegrals[2] = { 0.0 }; + for( localIndex numSubTriangle = 0; numSubTriangle < numFaceVertices; ++numSubTriangle ) + { + localIndex const nextVertex = (numSubTriangle+1)%numFaceVertices; + // - compute value of 2D monomials at the quadrature point on the sub-triangle (the + // barycenter). + // The result is ((v(0)+v(1)+faceCenter)/3 - faceCenter) / faceDiameter = + // = (v(0) + v(1) - 2*faceCenter)/(3*faceDiameter). + real64 monomialValues[2]; + for( localIndex i = 0; i < 2; ++i ) + { + monomialValues[i] = (faceRotatedVertices[numSubTriangle][i] + + faceRotatedVertices[nextVertex][i] - + 2.0*faceRotatedCentroid[i]) / (3.0*faceDiameter); + } + // compute value of 3D monomials at the quadrature point on the sub-triangle (the + // barycenter). The result is + // ((v(0) + v(1) + faceCenter)/3 - cellCenter)/cellDiameter. + real64 threeDMonomialValues[3]; + for( localIndex i = 0; i < 3; ++i ) + { + threeDMonomialValues[i] = ( (faceCenter[i] + + nodesCoords[faceToNodes[ numSubTriangle ]][i] + + nodesCoords[faceToNodes[ nextVertex ]][i]) / 3.0 - + cellCenter[i] ) * invCellDiameter; + } + // compute quadrature weight associated to the quadrature point (the area of the + // sub-triangle). + real64 edgesTangents[2][2]; // used to compute the area of the sub-triangle + for( localIndex i = 0; i < 2; ++i ) + { + edgesTangents[0][i] = faceRotatedVertices[numSubTriangle][i] - faceRotatedCentroid[i]; + } + for( localIndex i = 0; i < 2; ++i ) + { + edgesTangents[1][i] = faceRotatedVertices[nextVertex][i] - faceRotatedCentroid[i]; + } + real64 subTriangleArea = 0.5*LvArray::math::abs + ( edgesTangents[0][0]*edgesTangents[1][1] - + edgesTangents[0][1]*edgesTangents[1][0] ); + // compute the integrals on the sub-triangle and add it to the global integrals + for( localIndex i = 0; i < 2; ++i ) + { + monomInternalIntegrals[ i ] += monomialValues[ i ]*subTriangleArea; + } + for( localIndex i = 0; i < 3; ++i ) + { + // threeDMonomialIntegrals is assumed to be initialized to 0 by the caller + threeDMonomialIntegrals[ i ] += threeDMonomialValues[ i ]*subTriangleArea; + } + } + + // Compute integral of basis functions times normal derivative of monomials on the boundary. + real64 basisTimesMonomNormalDerBoundaryInt[ MFN ][ 2 ]; + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + for( localIndex i = 0; i < 2; ++i ) + { + basisTimesMonomNormalDerBoundaryInt[ numVertex ][ i ] = 0.0; + } + } + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + for( localIndex i = 0; i < 2; ++i ) + { + real64 thisEdgeIntTimesNormal_i = edgeOutwardNormals[numVertex][i]*edgeLengths[numVertex]; + basisTimesMonomNormalDerBoundaryInt[ localEdgeToNodes[ numVertex ][ 0 ] ][i] += thisEdgeIntTimesNormal_i; + basisTimesMonomNormalDerBoundaryInt[ localEdgeToNodes[ numVertex ][ 1 ] ][i] += thisEdgeIntTimesNormal_i; + } + } + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + for( localIndex i = 0; i < 2; ++i ) + { + basisTimesMonomNormalDerBoundaryInt[ numVertex ][ i ] *= 0.5*invFaceDiameter; + } + } + + // Compute integral mean of basis functions on this face. + real64 const invFaceArea = 1.0/faceArea; + real64 const monomialDerivativeInverse = (faceDiameter*faceDiameter)*invFaceArea; + for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) + { + real64 piNablaDofs[ 3 ]; + piNablaDofs[ 1 ] = monomialDerivativeInverse * + basisTimesMonomNormalDerBoundaryInt[ numVertex ][ 0 ]; + piNablaDofs[ 2 ] = monomialDerivativeInverse * + basisTimesMonomNormalDerBoundaryInt[ numVertex ][ 1 ]; + piNablaDofs[ 0 ] = (boundaryQuadratureWeights[ numVertex ] - + piNablaDofs[ 1 ]*monomBoundaryIntegrals[ 1 ] - + piNablaDofs[ 2 ]*monomBoundaryIntegrals[ 2 ])/monomBoundaryIntegrals[ 0 ]; + basisIntegrals[ numVertex ] = piNablaDofs[ 0 ]*faceArea + + (piNablaDofs[ 1 ]*monomInternalIntegrals[ 0 ] + + piNablaDofs[ 2 ]*monomInternalIntegrals[ 1 ]); + } +} + } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp index ccfc68f6ea7..242e52a99e6 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp @@ -49,6 +49,29 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM void synchronizeFractureState( DomainPartition & domain ) const; + void computeFaceNodalArea( localIndex const kf0, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, + ArrayOfArraysView< localIndex const > const & faceToNodeMap, + ArrayOfArraysView< localIndex const > const & faceToEdgeMap, + arrayView2d< localIndex const > const & edgeToNodeMap, + arrayView2d< real64 const > const faceCenters, + arrayView2d< real64 const > const faceNormals, + arrayView1d< real64 const > const faceAreas, + stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals ) const; + + void computeFaceIntegrals( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoords, + localIndex const (&faceToNodes)[11], + localIndex const (&faceToEdges)[11], + localIndex const & numFaceVertices, + real64 const & faceArea, + real64 const (&faceCenter)[3], + real64 const (&faceNormal)[3], + arrayView2d< localIndex const > const & edgeToNodes, + real64 const & invCellDiameter, + real64 const (&cellCenter)[3], + stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals, + real64 ( &threeDMonomialIntegrals )[3] ) const; + struct viewKeyStruct : SolidMechanicsLagrangianFEM::viewKeyStruct { constexpr static char const * fractureStateString() { return "fractureState"; } @@ -59,6 +82,9 @@ class ContactSolverBase : public SolidMechanicsLagrangianFEM }; + static constexpr localIndex m_maxFaceNodes = 11; // Maximum number of nodes on a contact face + + protected: virtual void postInputInitialization() override; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 02af0550280..4ad2fc03e7a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -114,6 +114,10 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase integer configurationLoopIter ) override final; + //TODO (jafranc) - bubble should be reframe as a stabilization + string getStabilizationName() const { return "TODO:bubble"; } + bool hasStabilization() const { return false;} + /** * @brief Loop over the finite element type on the fracture subregions of meshName and apply callback. * @tparam LAMBDA The callback function type diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp index b12a67de3ca..262fa2d528d 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp @@ -55,7 +55,6 @@ using namespace constitutive; using namespace dataRepository; using namespace fields; using namespace finiteElement; -const localIndex geos::SolidMechanicsLagrangeContact::m_maxFaceNodes = 11; SolidMechanicsLagrangeContact::SolidMechanicsLagrangeContact( const string & name, Group * const parent ): @@ -644,7 +643,7 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time, localMatrix, localRhs ); - assembleContact( domain, dofManager, localMatrix, localRhs ); + assembleContact( time, dt, domain, dofManager, localMatrix, localRhs ); // for sequential: add (fixed) pressure force contribution into residual (no derivatives) if( m_isFixedStressPoromechanicsUpdate || m_performStressInitialization ) @@ -658,7 +657,9 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time, } } -void SolidMechanicsLagrangeContact::assembleContact( DomainPartition & domain, +void SolidMechanicsLagrangeContact::assembleContact( real64 const GEOS_UNUSED_PARAM(time), + real64 const GEOS_UNUSED_PARAM(dt), + DomainPartition & domain, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ) @@ -724,9 +725,9 @@ void SolidMechanicsLagrangeContact:: Nbar[ 2 ] = faceNormal[elemsToFaces[kfe][0]][2] - faceNormal[elemsToFaces[kfe][1]][2]; LvArray::tensorOps::normalize< 3 >( Nbar ); - globalIndex rowDOF[3 * m_maxFaceNodes]; - real64 nodeRHS[3 * m_maxFaceNodes]; - stackArray1d< real64, 3 * m_maxFaceNodes > dRdP( 3*m_maxFaceNodes ); + globalIndex rowDOF[3 * ContactSolverBase::m_maxFaceNodes]; + real64 nodeRHS[3 * ContactSolverBase::m_maxFaceNodes]; + stackArray1d< real64, 3 * ContactSolverBase::m_maxFaceNodes > dRdP( 3*ContactSolverBase::m_maxFaceNodes ); for( localIndex kf=0; kf<2; ++kf ) { @@ -1000,329 +1001,7 @@ void SolidMechanicsLagrangeContact::computeRotationMatrices( DomainPartition & d } ); } -void SolidMechanicsLagrangeContact::computeFaceIntegrals( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoords, - localIndex const (&faceToNodes)[11], - localIndex const (&faceToEdges)[11], - localIndex const & numFaceVertices, - real64 const & faceArea, - real64 const (&faceCenter)[3], - real64 const (&faceNormal)[3], - arrayView2d< localIndex const > const & edgeToNodes, - real64 const & invCellDiameter, - real64 const (&cellCenter)[3], - stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals, - real64 (& threeDMonomialIntegrals)[3] ) const -{ - GEOS_MARK_FUNCTION; - localIndex const MFN = m_maxFaceNodes; // Max number of face vertices. - basisIntegrals.resize( numFaceVertices ); - // Rotate the face. - // - compute rotation matrix. - real64 faceRotationMatrix[ 3 ][ 3 ]; - computationalGeometry::RotationMatrix_3D( faceNormal, faceRotationMatrix ); - // - below we compute the diameter, the rotated vertices and the rotated center. - real64 faceRotatedVertices[ MFN ][ 2 ]; - real64 faceDiameter = 0; - - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - // apply the transpose (that is the inverse) of the rotation matrix to face vertices. - // NOTE: - // the second and third rows of the transpose of the rotation matrix rotate on the 2D face. - faceRotatedVertices[numVertex][0] = - faceRotationMatrix[ 0 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 0 ) + - faceRotationMatrix[ 1 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 1 ) + - faceRotationMatrix[ 2 ][ 1 ]*nodesCoords( faceToNodes[ numVertex ], 2 ); - faceRotatedVertices[numVertex][1] = - faceRotationMatrix[ 0 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 0 ) + - faceRotationMatrix[ 1 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 1 ) + - faceRotationMatrix[ 2 ][ 2 ]*nodesCoords( faceToNodes[ numVertex ], 2 ); - } - - faceDiameter = computationalGeometry::computeDiameter< 2 >( faceRotatedVertices, - numFaceVertices ); - real64 const invFaceDiameter = 1.0/faceDiameter; - // - rotate the face centroid as done for the vertices. - real64 faceRotatedCentroid[2]; - faceRotatedCentroid[0] = - faceRotationMatrix[ 0 ][ 1 ]*faceCenter[0] + - faceRotationMatrix[ 1 ][ 1 ]*faceCenter[1] + - faceRotationMatrix[ 2 ][ 1 ]*faceCenter[2]; - faceRotatedCentroid[1] = - faceRotationMatrix[ 0 ][ 2 ]*faceCenter[0] + - faceRotationMatrix[ 1 ][ 2 ]*faceCenter[1] + - faceRotationMatrix[ 2 ][ 2 ]*faceCenter[2]; - // - compute edges' lengths, outward pointing normals and local edge-to-nodes map. - real64 edgeOutwardNormals[ MFN ][ 2 ]; - real64 edgeLengths[ MFN ]; - localIndex localEdgeToNodes[ MFN ][ 2 ]; - - for( localIndex numEdge = 0; numEdge < numFaceVertices; ++numEdge ) - { - if( edgeToNodes( faceToEdges[numEdge], 0 ) == faceToNodes[ numEdge ] ) - { - localEdgeToNodes[ numEdge ][ 0 ] = numEdge; - localEdgeToNodes[ numEdge ][ 1 ] = (numEdge+1)%numFaceVertices; - } - else - { - localEdgeToNodes[ numEdge ][ 0 ] = (numEdge+1)%numFaceVertices; - localEdgeToNodes[ numEdge ][ 1 ] = numEdge; - } - real64 edgeTangent[2]; - edgeTangent[0] = faceRotatedVertices[(numEdge+1)%numFaceVertices][0] - - faceRotatedVertices[numEdge][0]; - edgeTangent[1] = faceRotatedVertices[(numEdge+1)%numFaceVertices][1] - - faceRotatedVertices[numEdge][1]; - edgeOutwardNormals[numEdge][0] = edgeTangent[1]; - edgeOutwardNormals[numEdge][1] = -edgeTangent[0]; - real64 signTestVector[2]; - signTestVector[0] = faceRotatedVertices[numEdge][0] - faceRotatedCentroid[0]; - signTestVector[1] = faceRotatedVertices[numEdge][1] - faceRotatedCentroid[1]; - if( signTestVector[0]*edgeOutwardNormals[numEdge][0] + - signTestVector[1]*edgeOutwardNormals[numEdge][1] < 0 ) - { - edgeOutwardNormals[numEdge][0] = -edgeOutwardNormals[numEdge][0]; - edgeOutwardNormals[numEdge][1] = -edgeOutwardNormals[numEdge][1]; - } - edgeLengths[numEdge] = LvArray::math::sqrt< real64 >( edgeTangent[0]*edgeTangent[0] + - edgeTangent[1]*edgeTangent[1] ); - edgeOutwardNormals[numEdge][0] /= edgeLengths[numEdge]; - edgeOutwardNormals[numEdge][1] /= edgeLengths[numEdge]; - } - - // Compute boundary quadrature weights (also equal to the integrals of basis functions on the - // boundary). - real64 boundaryQuadratureWeights[ MFN ]; - for( localIndex numWeight = 0; numWeight < numFaceVertices; ++numWeight ) - boundaryQuadratureWeights[numWeight] = 0.0; - for( localIndex numEdge = 0; numEdge < numFaceVertices; ++numEdge ) - { - boundaryQuadratureWeights[ localEdgeToNodes[ numEdge ][ 0 ] ] += 0.5*edgeLengths[numEdge]; - boundaryQuadratureWeights[ localEdgeToNodes[ numEdge ][ 1 ] ] += 0.5*edgeLengths[numEdge]; - } - // Compute scaled monomials' integrals on edges. - real64 monomBoundaryIntegrals[3] = { 0.0 }; - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - monomBoundaryIntegrals[0] += boundaryQuadratureWeights[ numVertex ]; - monomBoundaryIntegrals[1] += (faceRotatedVertices[ numVertex ][ 0 ] - faceRotatedCentroid[0]) * - invFaceDiameter*boundaryQuadratureWeights[ numVertex ]; - monomBoundaryIntegrals[2] += (faceRotatedVertices[ numVertex ][ 1 ] - faceRotatedCentroid[1]) * - invFaceDiameter*boundaryQuadratureWeights[ numVertex ]; - } - - // Compute non constant 2D and 3D scaled monomials' integrals on the face. - real64 monomInternalIntegrals[2] = { 0.0 }; - for( localIndex numSubTriangle = 0; numSubTriangle < numFaceVertices; ++numSubTriangle ) - { - localIndex const nextVertex = (numSubTriangle+1)%numFaceVertices; - // - compute value of 2D monomials at the quadrature point on the sub-triangle (the - // barycenter). - // The result is ((v(0)+v(1)+faceCenter)/3 - faceCenter) / faceDiameter = - // = (v(0) + v(1) - 2*faceCenter)/(3*faceDiameter). - real64 monomialValues[2]; - for( localIndex i = 0; i < 2; ++i ) - { - monomialValues[i] = (faceRotatedVertices[numSubTriangle][i] + - faceRotatedVertices[nextVertex][i] - - 2.0*faceRotatedCentroid[i]) / (3.0*faceDiameter); - } - // compute value of 3D monomials at the quadrature point on the sub-triangle (the - // barycenter). The result is - // ((v(0) + v(1) + faceCenter)/3 - cellCenter)/cellDiameter. - real64 threeDMonomialValues[3]; - for( localIndex i = 0; i < 3; ++i ) - { - threeDMonomialValues[i] = ( (faceCenter[i] + - nodesCoords[faceToNodes[ numSubTriangle ]][i] + - nodesCoords[faceToNodes[ nextVertex ]][i]) / 3.0 - - cellCenter[i] ) * invCellDiameter; - } - // compute quadrature weight associated to the quadrature point (the area of the - // sub-triangle). - real64 edgesTangents[2][2]; // used to compute the area of the sub-triangle - for( localIndex i = 0; i < 2; ++i ) - { - edgesTangents[0][i] = faceRotatedVertices[numSubTriangle][i] - faceRotatedCentroid[i]; - } - for( localIndex i = 0; i < 2; ++i ) - { - edgesTangents[1][i] = faceRotatedVertices[nextVertex][i] - faceRotatedCentroid[i]; - } - real64 subTriangleArea = 0.5*LvArray::math::abs - ( edgesTangents[0][0]*edgesTangents[1][1] - - edgesTangents[0][1]*edgesTangents[1][0] ); - // compute the integrals on the sub-triangle and add it to the global integrals - for( localIndex i = 0; i < 2; ++i ) - { - monomInternalIntegrals[ i ] += monomialValues[ i ]*subTriangleArea; - } - for( localIndex i = 0; i < 3; ++i ) - { - // threeDMonomialIntegrals is assumed to be initialized to 0 by the caller - threeDMonomialIntegrals[ i ] += threeDMonomialValues[ i ]*subTriangleArea; - } - } - - // Compute integral of basis functions times normal derivative of monomials on the boundary. - real64 basisTimesMonomNormalDerBoundaryInt[ MFN ][ 2 ]; - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - for( localIndex i = 0; i < 2; ++i ) - { - basisTimesMonomNormalDerBoundaryInt[ numVertex ][ i ] = 0.0; - } - } - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - for( localIndex i = 0; i < 2; ++i ) - { - real64 thisEdgeIntTimesNormal_i = edgeOutwardNormals[numVertex][i]*edgeLengths[numVertex]; - basisTimesMonomNormalDerBoundaryInt[ localEdgeToNodes[ numVertex ][ 0 ] ][i] += thisEdgeIntTimesNormal_i; - basisTimesMonomNormalDerBoundaryInt[ localEdgeToNodes[ numVertex ][ 1 ] ][i] += thisEdgeIntTimesNormal_i; - } - } - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - for( localIndex i = 0; i < 2; ++i ) - { - basisTimesMonomNormalDerBoundaryInt[ numVertex ][ i ] *= 0.5*invFaceDiameter; - } - } - - // Compute integral mean of basis functions on this face. - real64 const invFaceArea = 1.0/faceArea; - real64 const monomialDerivativeInverse = (faceDiameter*faceDiameter)*invFaceArea; - for( localIndex numVertex = 0; numVertex < numFaceVertices; ++numVertex ) - { - real64 piNablaDofs[ 3 ]; - piNablaDofs[ 1 ] = monomialDerivativeInverse * - basisTimesMonomNormalDerBoundaryInt[ numVertex ][ 0 ]; - piNablaDofs[ 2 ] = monomialDerivativeInverse * - basisTimesMonomNormalDerBoundaryInt[ numVertex ][ 1 ]; - piNablaDofs[ 0 ] = (boundaryQuadratureWeights[ numVertex ] - - piNablaDofs[ 1 ]*monomBoundaryIntegrals[ 1 ] - - piNablaDofs[ 2 ]*monomBoundaryIntegrals[ 2 ])/monomBoundaryIntegrals[ 0 ]; - basisIntegrals[ numVertex ] = piNablaDofs[ 0 ]*faceArea + - (piNablaDofs[ 1 ]*monomInternalIntegrals[ 0 ] + - piNablaDofs[ 2 ]*monomInternalIntegrals[ 1 ]); - } -} - -void SolidMechanicsLagrangeContact::computeFaceNodalArea( localIndex const kf0, - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, - ArrayOfArraysView< localIndex const > const & faceToNodeMap, - ArrayOfArraysView< localIndex const > const & faceToEdgeMap, - arrayView2d< localIndex const > const & edgeToNodeMap, - arrayView2d< real64 const > const faceCenters, - arrayView2d< real64 const > const faceNormals, - arrayView1d< real64 const > const faceAreas, - stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals ) const -{ - GEOS_MARK_FUNCTION; - localIndex const TriangularPermutation[3] = { 0, 1, 2 }; - localIndex const QuadrilateralPermutation[4] = { 0, 1, 3, 2 }; - localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); - - basisIntegrals.resize( numNodesPerFace ); - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - basisIntegrals[a] = 0.0; - } - localIndex const * const permutation = ( numNodesPerFace == 3 ) ? TriangularPermutation : QuadrilateralPermutation; - if( numNodesPerFace == 3 ) - { - real64 xLocal[3][3]; - for( localIndex a = 0; a < numNodesPerFace; ++a ) - { - for( localIndex j = 0; j < 3; ++j ) - { - xLocal[a][j] = nodePosition[faceToNodeMap( kf0, permutation[a] )][j]; - } - } - real64 N[3]; - for( localIndex q=0; q 4 && numNodesPerFace <= m_maxFaceNodes ) - { - // we need to L2 projector based on VEM to approximate the quadrature weights - // we need to use extra geometry information to computing L2 projector - - localIndex const MFN = m_maxFaceNodes; // Max number of face vertices. - localIndex const faceIndex = kf0; - localIndex const numFaceNodes = faceToNodeMap[ faceIndex ].size(); - - // get the face center and normal. - real64 const faceArea = faceAreas[ faceIndex ]; - localIndex faceToNodes[ MFN ]; - localIndex faceToEdges[ MFN ]; - for( localIndex i = 0; i < numFaceNodes; ++i ) - { - faceToNodes[i] = faceToNodeMap[ faceIndex ][ i ]; - faceToEdges[i] = faceToEdgeMap[ faceIndex ][ i ]; - } - // - get outward face normal and center - real64 faceNormal[3] = { faceNormals[faceIndex][0], - faceNormals[faceIndex][1], - faceNormals[faceIndex][2] }; - real64 const faceCenter[3] { faceCenters[faceIndex][0], - faceCenters[faceIndex][1], - faceCenters[faceIndex][2] }; - // - compute integrals calling auxiliary method - real64 threeDMonomialIntegrals[3] = { 0.0 }; - real64 const invCellDiameter = 0.0; - real64 const cellCenter[3] { 0.0, 0.0, 0.0 }; - computeFaceIntegrals( nodePosition, - faceToNodes, - faceToEdges, - numFaceNodes, - faceArea, - faceCenter, - faceNormal, - edgeToNodeMap, - invCellDiameter, - cellCenter, - basisIntegrals, - threeDMonomialIntegrals ); - } - else - { - GEOS_ERROR( GEOS_FMT( "Face with {} nodes. Only triangles and quadrilaterals and PEBI prisms up to 11 sides are supported.", - numNodesPerFace ), - getDataContext() ); - } -} void SolidMechanicsLagrangeContact:: assembleForceResidualDerivativeWrtTraction( MeshLevel const & mesh, @@ -1368,9 +1047,9 @@ void SolidMechanicsLagrangeContact:: localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elemsToFaces[kfe][0] ); - globalIndex rowDOF[3 * m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types - real64 nodeRHS[3 * m_maxFaceNodes]; - stackArray2d< real64, 3 * m_maxFaceNodes * 3 > dRdT( 3 * m_maxFaceNodes, 3 ); + globalIndex rowDOF[3 * ContactSolverBase::m_maxFaceNodes]; // this needs to be changed when dealing with arbitrary element types + real64 nodeRHS[3 * ContactSolverBase::m_maxFaceNodes]; + stackArray2d< real64, 3 * ContactSolverBase::m_maxFaceNodes * 3 > dRdT( 3 * ContactSolverBase::m_maxFaceNodes, 3 ); globalIndex colDOF[3]; for( localIndex i = 0; i < 3; ++i ) { @@ -1498,7 +1177,7 @@ void SolidMechanicsLagrangeContact:: if( ghostRank[kfe] < 0 ) { localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( elemsToFaces[kfe][0] ); - globalIndex nodeDOF[2 * 3 * m_maxFaceNodes]; + globalIndex nodeDOF[2 * 3 * ContactSolverBase::m_maxFaceNodes]; globalIndex elemDOF[3]; for( localIndex i = 0; i < 3; ++i ) { @@ -1508,7 +1187,7 @@ void SolidMechanicsLagrangeContact:: real64 elemRHS[3] = {0.0, 0.0, 0.0}; real64 const Ja = area[kfe]; - stackArray2d< real64, 2 * 3 * m_maxFaceNodes * 3 > dRdU( 3, 2 * 3 * m_maxFaceNodes ); + stackArray2d< real64, 2 * 3 * ContactSolverBase::m_maxFaceNodes * 3 > dRdU( 3, 2 * 3 * ContactSolverBase::m_maxFaceNodes ); stackArray2d< real64, 3 * 3 > dRdT( 3, 3 ); switch( fractureState[kfe] ) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp index 674a08876c1..8c20a6be0c5 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp @@ -102,7 +102,9 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase void updateState( DomainPartition & domain ) override final; - void assembleContact( DomainPartition & domain, + void assembleContact( real64 const time, + real64 const dt, + DomainPartition & domain, DofManager const & dofManager, CRSMatrixView< real64, globalIndex const > const & localMatrix, arrayView1d< real64 > const & localRhs ); @@ -142,32 +144,10 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase void computeTolerances( DomainPartition & domain ) const; - void computeFaceNodalArea( localIndex const kf0, - arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodePosition, - ArrayOfArraysView< localIndex const > const & faceToNodeMap, - ArrayOfArraysView< localIndex const > const & faceToEdgeMap, - arrayView2d< localIndex const > const & edgeToNodeMap, - arrayView2d< real64 const > const faceCenters, - arrayView2d< real64 const > const faceNormals, - arrayView1d< real64 const > const faceAreas, - stackArray1d< real64, FaceManager::maxFaceNodes() > & nodalArea ) const; - - void computeFaceIntegrals( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoords, - localIndex const (&faceToNodes)[11], - localIndex const (&faceToEdges)[11], - localIndex const & numFaceVertices, - real64 const & faceArea, - real64 const (&faceCenter)[3], - real64 const (&faceNormal)[3], - arrayView2d< localIndex const > const & edgeToNodes, - real64 const & invCellDiameter, - real64 const (&cellCenter)[3], - stackArray1d< real64, FaceManager::maxFaceNodes() > & basisIntegrals, - real64 ( &threeDMonomialIntegrals )[3] ) const; - real64 const machinePrecision = std::numeric_limits< real64 >::epsilon(); string getStabilizationName() const { return m_stabilizationName; } + bool hasStabilization() const { return true;} protected: @@ -186,8 +166,6 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase real64 m_stabilizationScalingCoefficient = 1.0; - static const localIndex m_maxFaceNodes; // Maximum number of nodes on a contact face - void computeFaceDisplacementJump( DomainPartition & domain ); struct viewKeyStruct : ContactSolverBase::viewKeyStruct