From 79c04227e7a96d0bf0ee7aba616055e63abb3337 Mon Sep 17 00:00:00 2001 From: Matteo Frigo <43966286+matteofrigo5@users.noreply.github.com> Date: Thu, 23 Jan 2025 17:21:34 -0800 Subject: [PATCH] feat: MGR Strategy for ALM - Face Bubble Functions for Wedge Elements - Multipoint Integration Rules for Tetrahedra and Triangles (#3395) * Adding linear solving strategy * Adding bubble functions for wedge elements * Fixed face bubble functions for wedge elements * Adding parameters to tune the iterative penalty coefficients * Adding points for quadrature rule * Scaling dispJump by area * Added var to set the number of quadrature points for tetrahedron and triangle * Cleaning up SolidMechanicsALMSimultaneousKernels.hpp * Cleaning up H1_TriangleFace_Lagrange1_Gauss.hpp and H1_Tetrahedron_Lagrange1_Gauss.hpp * Fixing some review comments * Fixing compilation issues * Fixed bug (int-real type casting) * Modified the default input parameters --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 + .../SCEC_BP6_QD_S_explicit.xml | 4 +- .../SCEC_BP6_QD_S_implicit.xml | 4 +- .../inducedSeismicity/SCEC_BP7_QD_base.xml | 4 +- .../constitutive/contact/CoulombFriction.hpp | 16 +- .../constitutive/contact/FrictionBase.hpp | 6 +- .../finiteElement/CMakeLists.txt | 4 +- .../FiniteElementDiscretization.cpp | 24 +- .../FiniteElementDiscretization.hpp | 4 + .../finiteElement/FiniteElementDispatch.hpp | 8 +- ...hpp => H1_Tetrahedron_Lagrange1_Gauss.hpp} | 317 ++++++++++++++---- ...pp => H1_TriangleFace_Lagrange1_Gauss.hpp} | 205 +++++++++-- .../H1_Wedge_Lagrange1_Gauss6.hpp | 67 +++- .../testH1_Tetrahedron_Lagrange1_Gauss1.cpp | 2 +- .../testH1_TriangleFace_Lagrange1_Gauss1.cpp | 2 +- .../linearAlgebra/CMakeLists.txt | 1 + .../interfaces/hypre/HypreMGR.cpp | 6 + .../AugmentedLagrangianContactMechanics.hpp | 101 ++++++ .../utilities/LinearSolverParameters.hpp | 2 + .../mesh/utilities/ComputationalGeometry.hpp | 2 +- .../physicsSolvers/contact/ContactFields.hpp | 12 +- .../contact/ContactSolverBase.cpp | 2 + ...lidMechanicsAugmentedLagrangianContact.cpp | 173 ++++++++-- ...lidMechanicsAugmentedLagrangianContact.hpp | 38 ++- .../contact/SolidMechanicsLagrangeContact.cpp | 2 +- .../kernels/SolidMechanicsALMKernels.hpp | 8 +- .../kernels/SolidMechanicsALMKernelsBase.hpp | 2 - .../SolidMechanicsALMSimultaneousKernels.hpp | 18 +- .../PoromechanicsKernels.cmake | 2 + .../kernels/SolidMechanicsKernels.cmake | 1 + .../surfaceGeneration/SurfaceGenerator.cpp | 9 +- 32 files changed, 863 insertions(+), 189 deletions(-) rename src/coreComponents/finiteElement/elementFormulations/{H1_Tetrahedron_Lagrange1_Gauss1.hpp => H1_Tetrahedron_Lagrange1_Gauss.hpp} (57%) rename src/coreComponents/finiteElement/elementFormulations/{H1_TriangleFace_Lagrange1_Gauss1.hpp => H1_TriangleFace_Lagrange1_Gauss.hpp} (55%) create mode 100644 src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 622c1da3d1b..ea7e4b5b516 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3416-9790-5bdb1fa + baseline: integratedTests/baseline_integratedTests-pr3395-9832-31f145e allow_fail: all: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 4304e1b9671..fe597c2c6c8 100644 --- a/BASELINE_NOTES.md +++ b/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 #3395 (2024-01-22) +===================== +Add new fields and change the default input for some tests. + PR #3416 (2024-01-21) ===================== Refactoring of induced seismicity EQ solvers to add coupling. diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml index fbbfccd5b38..ffab615c0f9 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml @@ -64,7 +64,7 @@ @@ -136,4 +136,4 @@ targetExactTimestep="0" target="/Outputs/timeHistoryOutput"/> - \ No newline at end of file + diff --git a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml index 1b74521b1ad..87d8a3e8455 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml @@ -65,7 +65,7 @@ @@ -137,4 +137,4 @@ targetExactTimestep="0" target="/Outputs/timeHistoryOutput"/> - \ No newline at end of file + diff --git a/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml index 08fc7cb9930..8df6ff94139 100644 --- a/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml +++ b/inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml @@ -45,7 +45,7 @@ @@ -248,4 +248,4 @@ targetExactTimestep="0" target="/Outputs/restart"/> --> - \ No newline at end of file + diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp index d64453fe31a..9041ac93ec0 100644 --- a/src/coreComponents/constitutive/contact/CoulombFriction.hpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.hpp @@ -103,7 +103,6 @@ class CoulombFrictionUpdates : public FrictionBaseUpdates arraySlice1d< real64 const > const & dispJump, arraySlice1d< real64 const > const & penalty, arraySlice1d< real64 const > const & traction, - real64 const faceArea, bool const symmetric, bool const fixedLimitTau, real64 const normalTractionTolerance, @@ -119,7 +118,6 @@ class CoulombFrictionUpdates : public FrictionBaseUpdates arraySlice1d< real64 const > const & deltaDispJump, arraySlice1d< real64 const > const & penalty, arraySlice1d< real64 const > const & traction, - real64 const faceArea, arraySlice1d< real64 > const & tractionNew ) const override final; GEOS_HOST_DEVICE @@ -363,7 +361,6 @@ inline void CoulombFrictionUpdates::updateTraction( arraySlice1d< real64 const > arraySlice1d< real64 const > const & dispJump, arraySlice1d< real64 const > const & penalty, arraySlice1d< real64 const > const & traction, - real64 const faceArea, bool const symmetric, bool const fixedLimitTau, real64 const normalTractionTolerance, @@ -380,9 +377,9 @@ inline void CoulombFrictionUpdates::updateTraction( arraySlice1d< real64 const > // Compute the trial traction real64 tractionTrial[ 3 ]; - tractionTrial[ 0 ] = traction[0] + penalty[0] * dispJump[0] * faceArea; - tractionTrial[ 1 ] = traction[1] + penalty[1] * (dispJump[1] - oldDispJump[1]) * faceArea; - tractionTrial[ 2 ] = traction[2] + penalty[1] * (dispJump[2] - oldDispJump[2]) * faceArea; + tractionTrial[ 0 ] = traction[0] + penalty[0] * dispJump[0]; + tractionTrial[ 1 ] = traction[1] + penalty[1] * (dispJump[1] - oldDispJump[1]); + tractionTrial[ 2 ] = traction[2] + penalty[1] * (dispJump[2] - oldDispJump[2]); // Compute tangential trial traction norm real64 const tau[2] = { tractionTrial[1], @@ -502,16 +499,15 @@ inline void CoulombFrictionUpdates::updateTractionOnly( arraySlice1d< real64 con arraySlice1d< real64 const > const & deltaDispJump, arraySlice1d< real64 const > const & penalty, arraySlice1d< real64 const > const & traction, - real64 const faceArea, arraySlice1d< real64 > const & tractionNew ) const { // TODO: Pass this tol as an argument or define a new class member real64 const zero = LvArray::NumericLimits< real64 >::epsilon; - tractionNew[0] = traction[0] + penalty[0] * dispJump[0] * faceArea; - tractionNew[1] = traction[1] + penalty[1] * deltaDispJump[1] * faceArea; - tractionNew[2] = traction[2] + penalty[1] * deltaDispJump[2] * faceArea; + tractionNew[0] = traction[0] + penalty[0] * dispJump[0]; + tractionNew[1] = traction[1] + penalty[1] * deltaDispJump[1]; + tractionNew[2] = traction[2] + penalty[1] * deltaDispJump[2]; real64 const tau[2] = { tractionNew[1], tractionNew[2] }; diff --git a/src/coreComponents/constitutive/contact/FrictionBase.hpp b/src/coreComponents/constitutive/contact/FrictionBase.hpp index d44a9fae7d4..b7583429e9b 100644 --- a/src/coreComponents/constitutive/contact/FrictionBase.hpp +++ b/src/coreComponents/constitutive/contact/FrictionBase.hpp @@ -125,7 +125,6 @@ class FrictionBaseUpdates arraySlice1d< real64 const > const & dispJump, arraySlice1d< real64 const > const & penalty, arraySlice1d< real64 const > const & traction, - real64 const faceArea, bool const symmetric, bool const fixedLimitTau, real64 const normalTractionTolerance, @@ -134,7 +133,7 @@ class FrictionBaseUpdates real64 ( & tractionNew )[3], integer & fractureState ) const { - GEOS_UNUSED_VAR( oldDispJump, dispJump, penalty, traction, faceArea, symmetric, fixedLimitTau, + GEOS_UNUSED_VAR( oldDispJump, dispJump, penalty, traction, symmetric, fixedLimitTau, normalTractionTolerance, tangentialTractionTolerance, dTraction_dDispJump, tractionNew, fractureState ); } @@ -153,9 +152,8 @@ class FrictionBaseUpdates arraySlice1d< real64 const > const & deltaDispJump, arraySlice1d< real64 const > const & penalty, arraySlice1d< real64 const > const & traction, - real64 const faceArea, arraySlice1d< real64 > const & tractionNew ) const - { GEOS_UNUSED_VAR( dispJump, deltaDispJump, penalty, traction, faceArea, tractionNew ); } + { GEOS_UNUSED_VAR( dispJump, deltaDispJump, penalty, traction, tractionNew ); } /** * @brief Check for the constraint satisfaction diff --git a/src/coreComponents/finiteElement/CMakeLists.txt b/src/coreComponents/finiteElement/CMakeLists.txt index 1fd77c86b1a..53ef7404158 100644 --- a/src/coreComponents/finiteElement/CMakeLists.txt +++ b/src/coreComponents/finiteElement/CMakeLists.txt @@ -32,8 +32,8 @@ set( finiteElement_headers elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp elementFormulations/H1_Pyramid_Lagrange1_Gauss5.hpp elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp - elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp - elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp + elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp + elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp elementFormulations/ConformingVirtualElementOrder1.hpp elementFormulations/ConformingVirtualElementOrder1_impl.hpp diff --git a/src/coreComponents/finiteElement/FiniteElementDiscretization.cpp b/src/coreComponents/finiteElement/FiniteElementDiscretization.cpp index 85f18da5efe..e809e456021 100644 --- a/src/coreComponents/finiteElement/FiniteElementDiscretization.cpp +++ b/src/coreComponents/finiteElement/FiniteElementDiscretization.cpp @@ -49,6 +49,11 @@ FiniteElementDiscretization::FiniteElementDiscretization( string const & name, G setInputFlag( InputFlags::OPTIONAL ). setApplyDefaultValue( 0 ). setDescription( "Specifier to indicate whether to force the use of VEM" ); + + registerWrapper( viewKeyStruct::useHighOrderQuadratureRuleString(), &m_useHighOrderQuadratureRule ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "Specifier to indicate whether to use a high order quadrature rule" ); } FiniteElementDiscretization::~FiniteElementDiscretization() @@ -68,7 +73,15 @@ FiniteElementDiscretization::factory( ElementType const parentElementShape ) con { switch( parentElementShape ) { - case ElementType::Triangle: return std::make_unique< H1_TriangleFace_Lagrange1_Gauss1 >(); + case ElementType::Triangle: + if( m_useHighOrderQuadratureRule == 1 ) + { + return std::make_unique< H1_TriangleFace_Lagrange1_Gauss4 >(); + } + else + { + return std::make_unique< H1_TriangleFace_Lagrange1_Gauss1 >(); + } case ElementType::Quadrilateral: return std::make_unique< H1_QuadrilateralFace_Lagrange1_GaussLegendre2 >(); // On polyhedra where FEM are available, we use VEM only if useVirtualElements is set to 1 in // the input file. @@ -80,7 +93,14 @@ FiniteElementDiscretization::factory( ElementType const parentElementShape ) con } else { - return std::make_unique< H1_Tetrahedron_Lagrange1_Gauss1 >(); + if( m_useHighOrderQuadratureRule == 1 ) + { + return std::make_unique< H1_Tetrahedron_Lagrange1_Gauss14 >(); + } + else + { + return std::make_unique< H1_Tetrahedron_Lagrange1_Gauss1 >(); + } } } case ElementType::Pyramid: diff --git a/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp b/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp index 0ee428550a9..21d7e9f71d5 100644 --- a/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp +++ b/src/coreComponents/finiteElement/FiniteElementDiscretization.hpp @@ -97,6 +97,7 @@ class FiniteElementDiscretization : public dataRepository::Group static constexpr char const * orderString() { return "order"; } static constexpr char const * formulationString() { return "formulation"; } static constexpr char const * useVemString() { return "useVirtualElements"; } + static constexpr char const * useHighOrderQuadratureRuleString() { return "useHighOrderQuadratureRule"; } }; /// The order of the finite element basis @@ -108,6 +109,9 @@ class FiniteElementDiscretization : public dataRepository::Group /// Optional parameter indicating if the class should use Virtual Elements. int m_useVem; + /// Optional parameter indicating if the class should use a high order quadrature rule. + int m_useHighOrderQuadratureRule; + void postInputInitialization() override final; }; diff --git a/src/coreComponents/finiteElement/FiniteElementDispatch.hpp b/src/coreComponents/finiteElement/FiniteElementDispatch.hpp index 5db1789d6e8..f7abaf86252 100644 --- a/src/coreComponents/finiteElement/FiniteElementDispatch.hpp +++ b/src/coreComponents/finiteElement/FiniteElementDispatch.hpp @@ -24,19 +24,20 @@ #include "elementFormulations/ConformingVirtualElementOrder1.hpp" #include "elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp" #include "elementFormulations/H1_Pyramid_Lagrange1_Gauss5.hpp" -#include "elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp" +#include "elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp" #include "elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp" #if !defined( GEOS_USE_HIP ) #include "elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp" #endif #include "elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp" -#include "elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp" +#include "elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp" #include "LvArray/src/system.hpp" #define FE_1_TYPES \ finiteElement::H1_Hexahedron_Lagrange1_GaussLegendre2, \ finiteElement::H1_Wedge_Lagrange1_Gauss6, \ finiteElement::H1_Tetrahedron_Lagrange1_Gauss1, \ + finiteElement::H1_Tetrahedron_Lagrange1_Gauss14, \ finiteElement::H1_Pyramid_Lagrange1_Gauss5 #define GL_FE_TYPES \ @@ -88,7 +89,8 @@ #define FE_TYPES_2D \ finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2, \ - finiteElement::H1_TriangleFace_Lagrange1_Gauss1 + finiteElement::H1_TriangleFace_Lagrange1_Gauss1, \ + finiteElement::H1_TriangleFace_Lagrange1_Gauss4 #define BASE_FE_TYPES_2D FE_TYPES_2D diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp similarity index 57% rename from src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp rename to src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp index 0750d7a196b..2736955bff6 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp @@ -14,11 +14,11 @@ */ /** - * @file H1_Tetrahedron_Lagrange1_Gauss1.hpp + * @file H1_Tetrahedron_Lagrange1_Gauss.hpp */ -#ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS1_HPP_ -#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS1_HPP_ +#ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS_HPP_ +#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS_HPP_ #include "FiniteElementBase.hpp" #include "LagrangeBasis1.hpp" @@ -38,7 +38,7 @@ namespace finiteElement * + Node r s t * |\\_ ===== == == == * || \_ 0 0 0 0 - * | \ \_ t 1 1 0 1 + * | \ \_ t 1 1 0 0 * | +__ \_ | s 2 0 1 0 * | /2 \__ \_ | / 3 0 0 1 * |/ \__\ | / ===== == == == @@ -46,10 +46,17 @@ namespace finiteElement * 0 1 * */ -class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase +template< typename NUM_Q_POINTS > +class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase { public: + /// Check that the number of quadrature points is valid. + static_assert( ( NUM_Q_POINTS::value == 1 || + NUM_Q_POINTS::value == 5 || + NUM_Q_POINTS::value == 14 ), + "NUM_Q_POINTS::value must be 1, 5, or 14!" ); + /// The type of basis used for this element using BASIS = LagrangeBasis1; @@ -63,13 +70,13 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static localIndex maxSupportPoints = numNodes; /// The number of quadrature points per element. - constexpr static localIndex numQuadraturePoints = 1; + constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS::value; /// The number of sampling points per element. constexpr static int numSamplingPoints = numSamplingPointsPerDirection * numSamplingPointsPerDirection * numSamplingPointsPerDirection; GEOS_HOST_DEVICE - virtual ~H1_Tetrahedron_Lagrange1_Gauss1() override + virtual ~H1_Tetrahedron_Lagrange1_Gauss() override {} GEOS_HOST_DEVICE @@ -222,10 +229,10 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase static void calcFaceBubbleN( localIndex const q, real64 (& N)[numFaces] ) { - GEOS_UNUSED_VAR( q ); - // single quadrature point (centroid), i.e. r = s = t = 1/4 - real64 const pointCoord[3] = {0.25, 0.25, 0.25}; + real64 const pointCoord[3] = {quadratureParentCoords0( q ), + quadratureParentCoords1( q ), + quadratureParentCoords2( q )}; calcFaceBubbleN( pointCoord, N ); } @@ -322,7 +329,175 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 parentVolume = 1.0 / 6.0; /// The weight of each quadrature point. - constexpr static real64 weight = parentVolume / numQuadraturePoints; + constexpr static real64 weight = parentVolume; + //constexpr static real64 weight = parentVolume / numQuadraturePoints; + + /** + * @brief Calculate the weights + * @param q + * @return weight. + */ + GEOS_HOST_DEVICE + inline + constexpr static real64 quadratureWeight( localIndex const q ) + { + if constexpr (numQuadraturePoints == 1) + { + real64 const w[numQuadraturePoints] = { 1.0 }; + return w[q]; + } + else if constexpr (numQuadraturePoints == 5) + { + real64 const w[numQuadraturePoints] = {-4.0/5.0, 9.0/20.0, 9.0/20.0, 9.0/20.0, 9.0/20.0 }; + return w[q]; + } + else if constexpr (numQuadraturePoints == 14) + { + real64 const w[numQuadraturePoints] = { 0.073493043116361949544, + 0.073493043116361949544, + 0.073493043116361949544, + 0.073493043116361949544, + 0.11268792571801585080, + 0.11268792571801585080, + 0.11268792571801585080, + 0.11268792571801585080, + 0.042546020777081466438, + 0.042546020777081466438, + 0.042546020777081466438, + 0.042546020777081466438, + 0.042546020777081466438, + 0.042546020777081466438 }; + return w[q]; + } + + } + + /** + * @brief Calculate the parent coordinates for the r direction, given the + * linear index of a quadrature point. + * @param a The linear index of quadrature point + * @return parent coordinate in the r direction. + */ + GEOS_HOST_DEVICE + inline + constexpr static real64 quadratureParentCoords0( localIndex const q ) + { + + if constexpr (numQuadraturePoints == 1) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 5) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/2.0, 1.0/6.0, 1.0/6.0, 1.0/6.0 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 14) + { + real64 const qCoords[numQuadraturePoints] = { 0.72179424906732632079, + 0.092735250310891226402, + 0.092735250310891226402, + 0.092735250310891226402, + 0.067342242210098170608, + 0.31088591926330060980, + 0.31088591926330060980, + 0.31088591926330060980, + 0.045503704125649649492, + 0.045503704125649649492, + 0.045503704125649649492, + 0.45449629587435035051, + 0.45449629587435035051, + 0.45449629587435035051 }; + return qCoords[q]; + } + + } + + /** + * @brief Calculate the parent coordinates for the s direction, given the + * linear index of a quadrature point. + * @param q The linear index of quadrature point + * @return parent coordinate in the s direction. + */ + GEOS_HOST_DEVICE + inline + constexpr static real64 quadratureParentCoords1( localIndex const q ) + { + + if constexpr (numQuadraturePoints == 1) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 5) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/2.0, 1.0/6.0, 1.0/6.0 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 14) + { + real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402, + 0.72179424906732632079, + 0.092735250310891226402, + 0.092735250310891226402, + 0.31088591926330060980, + 0.067342242210098170608, + 0.31088591926330060980, + 0.31088591926330060980, + 0.045503704125649649492, + 0.45449629587435035051, + 0.45449629587435035051, + 0.045503704125649649492, + 0.045503704125649649492, + 0.45449629587435035051 }; + return qCoords[q]; + } + + } + + /** + * @brief Calculate the parent coordinates for the t direction, given the + * linear index of a quadrature point. + * @param q The linear index of quadrature point + * @return parent coordinate in the t direction. + */ + GEOS_HOST_DEVICE + inline + constexpr static real64 quadratureParentCoords2( localIndex const q ) + { + + if constexpr (numQuadraturePoints == 1) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 5) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/6.0, 1.0/2.0, 1.0/6.0 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 14) + { + real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402, + 0.092735250310891226402, + 0.72179424906732632079, + 0.092735250310891226402, + 0.31088591926330060980, + 0.31088591926330060980, + 0.067342242210098170608, + 0.31088591926330060980, + 0.45449629587435035051, + 0.045503704125649649492, + 0.45449629587435035051, + 0.045503704125649649492, + 0.45449629587435035051, + 0.045503704125649649492 }; + + return qCoords[q]; + } + + } /** * @brief Calculates the determinant of the Jacobian of the isoparametric @@ -337,11 +512,12 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase /// @cond Doxygen_Suppress +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 -H1_Tetrahedron_Lagrange1_Gauss1:: - determinantJacobianTransformation( real64 const (&X)[numNodes][3] ) +H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: +determinantJacobianTransformation( real64 const (&X)[numNodes][3] ) { return ( X[1][0] - X[0][0] )*( ( X[2][1] - X[0][1] )*( X[3][2] - X[0][2] ) - ( X[3][1] - X[0][1] )*( X[2][2] - X[0][2] ) ) + ( X[1][1] - X[0][1] )*( ( X[3][0] - X[0][0] )*( X[2][2] - X[0][2] ) - ( X[2][0] - X[0][0] )*( X[3][2] - X[0][2] ) ) @@ -350,12 +526,13 @@ H1_Tetrahedron_Lagrange1_Gauss1:: //************************************************************************************************* +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void -H1_Tetrahedron_Lagrange1_Gauss1:: - calcN( real64 const (&coords)[3], - real64 (& N)[numNodes] ) +H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcN( real64 const (&coords)[3], + real64 (& N)[numNodes] ) { real64 const r = coords[0]; real64 const s = coords[1]; @@ -369,41 +546,43 @@ H1_Tetrahedron_Lagrange1_Gauss1:: } +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void -H1_Tetrahedron_Lagrange1_Gauss1:: - calcN( localIndex const q, - real64 (& N)[numNodes] ) +H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcN( localIndex const q, + real64 (& N)[numNodes] ) { - GEOS_UNUSED_VAR( q ); + real64 const pointCoord[3] = {quadratureParentCoords0( q ), + quadratureParentCoords1( q ), + quadratureParentCoords2( q )}; - // single quadrature point (centroid), i.e. r = s = t = 1/4 - real64 const pointCoord[3] = {0.25, 0.25, 0.25}; calcN( pointCoord, N ); } +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline -void H1_Tetrahedron_Lagrange1_Gauss1:: - calcN( localIndex const q, - StackVariables const & GEOS_UNUSED_PARAM( stack ), - real64 ( & N )[numNodes] ) +void H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcN( localIndex const q, + StackVariables const & GEOS_UNUSED_PARAM( stack ), + real64 ( & N )[numNodes] ) { return calcN( q, N ); } //************************************************************************************************* +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 -H1_Tetrahedron_Lagrange1_Gauss1:: - calcGradN( localIndex const q, - real64 const (&X)[numNodes][3], - real64 (& gradN)[numNodes][3] ) +H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcGradN( localIndex const q, + real64 const (&X)[numNodes][3], + real64 (& gradN)[numNodes][3] ) { - GEOS_UNUSED_VAR( q ); gradN[0][0] = X[1][1]*( X[3][2] - X[2][2] ) - X[2][1]*( X[3][2] - X[1][2] ) + X[3][1]*( X[2][2] - X[1][2] ); gradN[0][1] = -X[1][0]*( X[3][2] - X[2][2] ) + X[2][0]*( X[3][2] - X[1][2] ) - X[3][0]*( X[2][2] - X[1][2] ); @@ -432,52 +611,53 @@ H1_Tetrahedron_Lagrange1_Gauss1:: } } - return detJ * weight; + return detJ * weight * quadratureWeight( q ); } +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline -real64 H1_Tetrahedron_Lagrange1_Gauss1:: - calcGradN( localIndex const q, - real64 const (&X)[numNodes][3], - StackVariables const & GEOS_UNUSED_PARAM( stack ), - real64 ( & gradN )[numNodes][3] ) +real64 H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcGradN( localIndex const q, + real64 const (&X)[numNodes][3], + StackVariables const & GEOS_UNUSED_PARAM( stack ), + real64 ( & gradN )[numNodes][3] ) { return calcGradN( q, X, gradN ); } +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 -H1_Tetrahedron_Lagrange1_Gauss1::calcGradFaceBubbleN( localIndex const q, - real64 const (&X)[numNodes][3], - real64 (& gradN)[numFaces][3] ) +H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::calcGradFaceBubbleN( localIndex const q, + real64 const (&X)[numNodes][3], + real64 (& gradN)[numFaces][3] ) { - GEOS_UNUSED_VAR( q ); - - real64 detJ = determinantJacobianTransformation( X ); - real64 factor = 1.0 / ( detJ ); + //real64 detJ = determinantJacobianTransformation( X ); real64 J[3][3] = {{0}}; - J[0][0] = (-X[0][1]*( X[3][2] - X[2][2] ) + X[2][1]*( X[3][2] - X[0][2] ) - X[3][1]*( X[2][2] - X[0][2] ))*factor; - J[0][1] = ( X[0][0]*( X[3][2] - X[2][2] ) - X[2][0]*( X[3][2] - X[0][2] ) + X[3][0]*( X[2][2] - X[0][2] ))*factor; - J[0][2] = (-X[0][0]*( X[3][1] - X[2][1] ) + X[2][0]*( X[3][1] - X[0][1] ) - X[3][0]*( X[2][1] - X[0][1] ))*factor; + J[0][0] = X[1][0] - X[0][0]; + J[0][1] = X[2][0] - X[0][0]; + J[0][2] = X[3][0] - X[0][0]; - J[1][0] = ( X[0][1]*( X[3][2] - X[1][2] ) - X[1][1]*( X[3][2] - X[0][2] ) + X[3][1]*( X[1][2] - X[0][2] ))*factor; - J[1][1] = (-X[0][0]*( X[3][2] - X[1][2] ) + X[1][0]*( X[3][2] - X[0][2] ) - X[3][0]*( X[1][2] - X[0][2] ))*factor; - J[1][2] = ( X[0][0]*( X[3][1] - X[1][1] ) - X[1][0]*( X[3][1] - X[0][1] ) + X[3][0]*( X[1][1] - X[0][1] ))*factor; + J[1][0] = X[1][1] - X[0][1]; + J[1][1] = X[2][1] - X[0][1]; + J[1][2] = X[3][1] - X[0][1]; - J[2][0] = (-X[0][1]*( X[2][2] - X[1][2] ) + X[1][1]*( X[2][2] - X[0][2] ) - X[2][1]*( X[1][2] - X[0][2] ))*factor; - J[2][1] = ( X[0][0]*( X[2][2] - X[1][2] ) - X[1][0]*( X[2][2] - X[0][2] ) + X[2][0]*( X[1][2] - X[0][2] ))*factor; - J[2][2] = (-X[0][0]*( X[2][1] - X[1][1] ) + X[1][0]*( X[2][1] - X[0][1] ) - X[2][0]*( X[1][1] - X[0][1] ))*factor; + J[2][0] = X[1][2] - X[0][2]; + J[2][1] = X[2][2] - X[0][2]; + J[2][2] = X[3][2] - X[0][2]; + + real64 const detJ = LvArray::tensorOps::invert< 3 >( J ); real64 dNdXi[numFaces][3] = {{0}}; - // single quadrature point (centroid), i.e. r = s = t = 1/4 - real64 const r = 1.0 / 4.0; - real64 const s = 1.0 / 4.0; - real64 const t = 1.0 / 4.0; + + real64 const r = quadratureParentCoords0( q ); + real64 const s = quadratureParentCoords1( q ); + real64 const t = quadratureParentCoords2( q ); dNdXi[0][0] = ( 1 - 2 * r - s - t ) * t; // dN0/dr dNdXi[0][1] = -r * t; // dN0/ds @@ -485,7 +665,7 @@ H1_Tetrahedron_Lagrange1_Gauss1::calcGradFaceBubbleN( localIndex const q, dNdXi[1][0] = ( 1 - 2 * r - s - t ) * s; // dN1/dr dNdXi[1][1] = ( 1 - r - 2 * s - t ) * r; // dN1/ds - dNdXi[1][2] = -r * s; // dN1/dt + dNdXi[1][2] = -r * s; // dN1/dt dNdXi[2][0] = -t * s; // dN2/dr dNdXi[2][1] = ( 1 - r - 2 * s - t ) * t; // dN2/ds @@ -502,28 +682,35 @@ H1_Tetrahedron_Lagrange1_Gauss1::calcGradFaceBubbleN( localIndex const q, gradN[fi][2] = dNdXi[fi][0] * J[0][2] + dNdXi[fi][1] * J[1][2] + dNdXi[fi][2] * J[2][2]; } - return detJ * weight; + return detJ * weight * quadratureWeight( q ); } //************************************************************************************************* +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 -H1_Tetrahedron_Lagrange1_Gauss1:: - transformedQuadratureWeight( localIndex const q, - real64 const (&X)[numNodes][3] ) +H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: +transformedQuadratureWeight( localIndex const q, + real64 const (&X)[numNodes][3] ) { - GEOS_UNUSED_VAR( q ); real64 detJ = determinantJacobianTransformation( X ); - return detJ * weight; + return detJ * weight * quadratureWeight( q ); } /// @endcond +/// @brief Instantiate the H1_Tetrahedron_Lagrange1_Gauss class for the 1-point Gaussian quadrature rule. +using H1_Tetrahedron_Lagrange1_Gauss1 = H1_Tetrahedron_Lagrange1_Gauss< std::integral_constant< int, 1 > >; +/// @brief Instantiate the H1_Tetrahedron_Lagrange1_Gauss class for the 5-point Gaussian quadrature rule. +using H1_Tetrahedron_Lagrange1_Gauss5 = H1_Tetrahedron_Lagrange1_Gauss< std::integral_constant< int, 5 > >; +/// @brief Instantiate the H1_Tetrahedron_Lagrange1_Gauss class for the 14-point Gaussian quadrature rule. +using H1_Tetrahedron_Lagrange1_Gauss14 = H1_Tetrahedron_Lagrange1_Gauss< std::integral_constant< int, 14 > >; + } } -#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS1_HPP_ +#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TETRAHEDRONLAGRANGE1GAUSS_HPP_ diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp similarity index 55% rename from src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp rename to src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp index 9d87f889397..61a3086598b 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp @@ -14,11 +14,11 @@ */ /** - * @file H1_TriangleFace_Lagrange1_Gauss1.hpp + * @file H1_TriangleFace_Lagrange1_Gauss.hpp */ -#ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS1_HPP_ -#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS1_HPP_ +#ifndef GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS_HPP_ +#define GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS_HPP_ #include "FiniteElementBase.hpp" #include "LagrangeBasis1.hpp" @@ -47,10 +47,17 @@ namespace finiteElement * 0 1 * */ -class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase +template< typename NUM_Q_POINTS > +class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase { public: + /// Check that the number of quadrature points is valid. + static_assert( ( NUM_Q_POINTS::value == 1 || + NUM_Q_POINTS::value == 4 || + NUM_Q_POINTS::value == 6 ), + "NUM_Q_POINTS::value must be 1, 4, or 6!" ); + /// The type of basis used for this element using BASIS = LagrangeBasis1; @@ -60,10 +67,10 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static localIndex maxSupportPoints = numNodes; /// The number of quadrature points per element. - constexpr static localIndex numQuadraturePoints = 1; + constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS::value; GEOS_HOST_DEVICE - virtual ~H1_TriangleFace_Lagrange1_Gauss1() override + virtual ~H1_TriangleFace_Lagrange1_Gauss() override {} GEOS_HOST_DEVICE @@ -173,10 +180,9 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase static void calcBubbleN( localIndex const q, real64 (& N)[1] ) { - GEOS_UNUSED_VAR( q ); - // single quadrature point (centroid), i.e. r = s = 1/3 - real64 const qCoords[2] = { 1.0 / 3.0, 1.0 / 3.0}; + real64 const qCoords[2] = {quadratureParentCoords0( q ), quadratureParentCoords1( q ) }; + calcBubbleN( qCoords, N ); } @@ -228,33 +234,147 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 parentArea = 0.5; /// The weight of each quadrature point. - constexpr static real64 weight = parentArea / numQuadraturePoints; + //constexpr static real64 weight = parentArea / numQuadraturePoints; + constexpr static real64 weight = parentArea; + + /** + * @brief Calculate the weights + * @param q + * @return weight. + */ + GEOS_HOST_DEVICE + inline + constexpr static real64 quadratureWeight( localIndex const q ) + { + + if constexpr (numQuadraturePoints == 1) + { + constexpr real64 w[numQuadraturePoints] = { 1.0 }; + return w[q]; + } + else if constexpr (numQuadraturePoints == 4) + { + constexpr real64 w[numQuadraturePoints] = {-0.562500000000000, + 0.520833333333333, + 0.520833333333333, + 0.520833333333333 }; + return w[q]; + } + else if constexpr (numQuadraturePoints == 6) + { + real64 const w[numQuadraturePoints] = { 0.166666666666666, + 0.166666666666666, + 0.166666666666666, + 0.166666666666666, + 0.166666666666666, + 0.166666666666666 }; + return w[q]; + } + + } + + /** + * @brief Calculate the parent coordinates for the r direction, given the + * linear index of a quadrature point. + * @param a The linear index of quadrature point + * @return parent coordinate in the r direction. + */ + GEOS_HOST_DEVICE + inline + constexpr static real64 quadratureParentCoords0( localIndex const q ) + { + + if constexpr (numQuadraturePoints == 1) + { + constexpr real64 qCoords[numQuadraturePoints] = { 1.0/3.0 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 4) + { + constexpr real64 qCoords[numQuadraturePoints] = { 0.333333333333333, + 0.600000000000000, + 0.200000000000000, + 0.200000000000000 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 6) + { + constexpr real64 qCoords[numQuadraturePoints] = { 0.659027622374092, + 0.109039009072877, + 0.231933368553031, + 0.659027622374092, + 0.109039009072877, + 0.231933368553031 }; + return qCoords[q]; + } + + } + + /** + * @brief Calculate the parent coordinates for the s direction, given the + * linear index of a quadrature point. + * @param q The linear index of quadrature point + * @return parent coordinate in the s direction. + */ + GEOS_HOST_DEVICE + inline + constexpr static real64 quadratureParentCoords1( localIndex const q ) + { + + if constexpr (numQuadraturePoints == 1) + { + constexpr real64 qCoords[numQuadraturePoints] = { 1.0/3.0 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 4) + { + constexpr real64 qCoords[numQuadraturePoints] = { 0.333333333333333, + 0.200000000000000, + 0.600000000000000, + 0.200000000000000 }; + return qCoords[q]; + } + else if constexpr (numQuadraturePoints == 6) + { + constexpr real64 qCoords[numQuadraturePoints] = { 0.231933368553031, + 0.659027622374092, + 0.109039009072877, + 0.109039009072877, + 0.231933368553031, + 0.659027622374092 }; + return qCoords[q]; + } + + } }; /// @cond Doxygen_Suppress + +template< typename NUM_Q_POINTS > template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, bool UPPER > GEOS_HOST_DEVICE inline -void H1_TriangleFace_Lagrange1_Gauss1:: - addGradGradStabilization( StackVariables const & stack, - real64 ( & matrix ) - [maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT] - [maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT], - real64 const & scaleFactor ) +void H1_TriangleFace_Lagrange1_Gauss< NUM_Q_POINTS >:: +addGradGradStabilization( StackVariables const & stack, + real64 ( & matrix ) + [maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT] + [maxSupportPoints * NUMDOFSPERTRIALSUPPORTPOINT], + real64 const & scaleFactor ) { GEOS_UNUSED_VAR( stack ); GEOS_UNUSED_VAR( matrix ); GEOS_UNUSED_VAR( scaleFactor ); } +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void -H1_TriangleFace_Lagrange1_Gauss1:: - calcN( real64 const (&coords)[2], - real64 ( & N )[numNodes] ) +H1_TriangleFace_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcN( real64 const (&coords)[2], + real64 ( & N )[numNodes] ) { real64 const r = coords[0]; real64 const s = coords[1]; @@ -264,49 +384,60 @@ H1_TriangleFace_Lagrange1_Gauss1:: N[2] = s; } +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void -H1_TriangleFace_Lagrange1_Gauss1:: - calcN( localIndex const q, - real64 (& N)[numNodes] ) +H1_TriangleFace_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcN( localIndex const q, + real64 (& N)[numNodes] ) { - GEOS_UNUSED_VAR( q ); - // single quadrature point (centroid), i.e. r = s = 1/3 - N[0] = 1.0 / 3.0; // N0 = 1 - r - s - N[1] = N[0]; // N1 = r - N[2] = N[0]; // N2 = s + real64 const qCoords[2] = {quadratureParentCoords0( q ), quadratureParentCoords1( q ) }; + + calcN( qCoords, N ); + } +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline -void H1_TriangleFace_Lagrange1_Gauss1:: - calcN( localIndex const q, - StackVariables const & GEOS_UNUSED_PARAM( stack ), - real64 ( & N )[numNodes] ) +void H1_TriangleFace_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcN( localIndex const q, + StackVariables const & GEOS_UNUSED_PARAM( stack ), + real64 ( & N )[numNodes] ) { return calcN( q, N ); } //************************************************************************************************* +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 -H1_TriangleFace_Lagrange1_Gauss1:: - transformedQuadratureWeight( localIndex const q, - real64 const (&X)[numNodes][3] ) +H1_TriangleFace_Lagrange1_Gauss< NUM_Q_POINTS >:: +transformedQuadratureWeight( localIndex const q, + real64 const (&X)[numNodes][3] ) { - GEOS_UNUSED_VAR( q ); + //GEOS_UNUSED_VAR( q ); real64 n[3] = { ( X[1][1] - X[0][1] ) * ( X[2][2] - X[0][2] ) - ( X[2][1] - X[0][1] ) * ( X[1][2] - X[0][2] ), ( X[2][0] - X[0][0] ) * ( X[1][2] - X[0][2] ) - ( X[1][0] - X[0][0] ) * ( X[2][2] - X[0][2] ), ( X[1][0] - X[0][0] ) * ( X[2][1] - X[0][1] ) - ( X[2][0] - X[0][0] ) * ( X[1][1] - X[0][1] )}; - return sqrt( n[0] * n[0] + n[1] * n[1] + n[2] * n[2] ) * weight; + + return sqrt( n[0] * n[0] + n[1] * n[1] + n[2] * n[2] ) * weight * quadratureWeight( q ); } /// @endcond +/// @brief Istanciation of the class with 1 quadrature points. +using H1_TriangleFace_Lagrange1_Gauss1 = H1_TriangleFace_Lagrange1_Gauss< std::integral_constant< int, 1 > >; +/// @brief Istanciation of the class with 4 quadrature points. +using H1_TriangleFace_Lagrange1_Gauss4 = H1_TriangleFace_Lagrange1_Gauss< std::integral_constant< int, 4 > >; +/// @brief Istanciation of the class with 6 quadrature points. +using H1_TriangleFace_Lagrange1_Gauss6 = H1_TriangleFace_Lagrange1_Gauss< std::integral_constant< int, 6 > >; + + } } -#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS1_HPP_ +#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS_HPP_ diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp index 3afcee034fc..d901a13572b 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp @@ -203,12 +203,21 @@ class H1_Wedge_Lagrange1_Gauss6 final : public FiniteElementBase static void calcFaceBubbleN( real64 const (&pointCoord)[3], real64 (& N)[numFaces] ) { - GEOS_UNUSED_VAR( pointCoord, N ); - GEOS_ERROR( "Unsupported bubble functions for wedge elements" ); + + real64 const r = pointCoord[0]; + real64 const s = pointCoord[1]; + real64 const xi = pointCoord[2]; + + N[0] = 4.0 * (1.0 - r - s) * s * LagrangeBasis1::valueBubble( xi ); + N[1] = 4.0 * (1.0 - r - s) * r * LagrangeBasis1::valueBubble( xi ); + N[2] = (1.0 - r - s) * r * s * LagrangeBasis1::value( 0, xi ); + N[3] = (1.0 - r - s) * r * s * LagrangeBasis1::value( 1, xi ); + N[4] = 4.0 * r * s * LagrangeBasis1::valueBubble( xi ); + } /** - * @brief Calculate shape functions values for each support face at a + * @brief Calculate face bubble functions values for each face at a * quadrature point. * @param q Index of the quadrature point. * @param N An array to pass back the shape function values for each support @@ -219,8 +228,12 @@ class H1_Wedge_Lagrange1_Gauss6 final : public FiniteElementBase static void calcFaceBubbleN( localIndex const q, real64 (& N)[numFaces] ) { - GEOS_UNUSED_VAR( q, N ); - GEOS_ERROR( "Unsupported bubble functions for wedge elements" ); + + real64 const pointCoord[3] = {quadratureParentCoords0( q ), + quadratureParentCoords1( q ), + quadratureParentCoords2( q )}; + + calcFaceBubbleN( pointCoord, N ); } /** @@ -610,9 +623,47 @@ H1_Wedge_Lagrange1_Gauss6::calcGradFaceBubbleN( localIndex const q, real64 const (&X)[numNodes][3], real64 (& gradN)[numFaces][3] ) { - GEOS_UNUSED_VAR( q, X, gradN ); - GEOS_ERROR( "Unsupported bubble functions for wedge elements" ); - return 0.0; + + real64 J[3][3] = {{0}}; + + jacobianTransformation( q, X, J ); + + real64 const detJ = LvArray::tensorOps::invert< 3 >( J ); + + real64 dNdXi[numFaces][3] = {{0}}; + + real64 const r = quadratureParentCoords0( q ); + real64 const s = quadratureParentCoords1( q ); + real64 const xi = quadratureParentCoords2( q ); + + dNdXi[0][0] = -4.0 * s * LagrangeBasis1::valueBubble( xi ); // dN0/dr + dNdXi[0][1] = 4.0 * (1.0 - r - 2.0 * s) * LagrangeBasis1::valueBubble( xi ); // dN0/ds + dNdXi[0][2] = 4.0 *(1 - r - s) * s * LagrangeBasis1::gradientBubble( xi ); // dN0/dxi + + dNdXi[1][0] = 4.0 * (1.0 - 2.0 * r - s) * LagrangeBasis1::valueBubble( xi ); // dN1/dr + dNdXi[1][1] = -4.0 * r * LagrangeBasis1::valueBubble( xi ); // dN1/ds + dNdXi[1][2] = 4.0 * (1 - r - s) * r * LagrangeBasis1::gradientBubble( xi ); // dN1/dxi + + dNdXi[2][0] = (1.0 - 2.0 * r - s) * s * LagrangeBasis1::value( 0, xi ); // dN2/dr + dNdXi[2][1] = (1.0 - r - 2.0 * s) * r * LagrangeBasis1::value( 0, xi ); // dN2/ds + dNdXi[2][2] = (1 - r - s) * r * s * LagrangeBasis1::gradient( 0, xi ); // dN2/dxi + + dNdXi[3][0] = (1.0 - 2.0 * r - s) * s * LagrangeBasis1::value( 1, xi ); // dN3/dr + dNdXi[3][1] = (1.0 - r - 2.0 * s) * r * LagrangeBasis1::value( 1, xi ); // dN3/ds + dNdXi[3][2] = (1 - r - s) * r * s * LagrangeBasis1::gradient( 1, xi ); // dN3/dxi + + dNdXi[4][0] = 4.0 * s * LagrangeBasis1::valueBubble( xi ); // dN4/dr + dNdXi[4][1] = 4.0 * r * LagrangeBasis1::valueBubble( xi ); // dN4/ds + dNdXi[4][2] = 4.0 * r * s * LagrangeBasis1::gradientBubble( xi ); // dN4/dxi + + for( int fi=0; fi -#include "finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp" +#include "finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp" using namespace geos; using namespace finiteElement; diff --git a/src/coreComponents/linearAlgebra/CMakeLists.txt b/src/coreComponents/linearAlgebra/CMakeLists.txt index b99b0697b03..9208438c9ba 100644 --- a/src/coreComponents/linearAlgebra/CMakeLists.txt +++ b/src/coreComponents/linearAlgebra/CMakeLists.txt @@ -128,6 +128,7 @@ if( ENABLE_HYPRE ) interfaces/hypre/mgrStrategies/HybridSinglePhasePoromechanics.hpp interfaces/hypre/mgrStrategies/Hydrofracture.hpp interfaces/hypre/mgrStrategies/LagrangianContactMechanics.hpp + interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp interfaces/hypre/mgrStrategies/LagrangianContactMechanicsBubbleStabilization.hpp interfaces/hypre/mgrStrategies/MultiphasePoromechanics.hpp interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp index 466fc6ff04d..5bfa3f08bf4 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp @@ -20,6 +20,7 @@ #include "HypreMGR.hpp" +#include "linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/CompositionalMultiphaseFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/CompositionalMultiphaseHybridFVM.hpp" #include "linearAlgebra/interfaces/hypre/mgrStrategies/CompositionalMultiphaseReservoirFVM.hpp" @@ -134,6 +135,11 @@ void hypre::mgr::createMGR( LinearSolverParameters const & params, setStrategy< Hydrofracture >( params.mgr, numComponentsPerField, precond, mgrData ); break; } + case LinearSolverParameters::MGR::StrategyType::augmentedLagrangianContactMechanics: + { + setStrategy< AugmentedLagrangianContactMechanics >( params.mgr, numComponentsPerField, precond, mgrData ); + break; + } case LinearSolverParameters::MGR::StrategyType::lagrangianContactMechanics: { setStrategy< LagrangianContactMechanics >( params.mgr, numComponentsPerField, precond, mgrData ); diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp new file mode 100644 index 00000000000..4f6186ca84c --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp @@ -0,0 +1,101 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file AugmentedLagrangianContactMechanics.hpp + */ + +#ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRALMCONTACTMECHANICS_HPP_ +#define GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRALMCONTACTMECHANICS_HPP_ + +#include "linearAlgebra/interfaces/hypre/HypreMGR.hpp" + +namespace geos +{ + +namespace hypre +{ + +namespace mgr +{ + +/** + * @brief AugmentedLagrangianContactMechanics strategy + * + * Augmented Lagrangina contact mechanics with bubble functions + * + * dofLabel: 0 = displacement, x-component + * dofLabel: 1 = displacement, y-component + * dofLabel: 2 = displacement, z-component + * dofLabel: 3 = displacement bubble function, x-component + * dofLabel: 4 = displacement bubble function, y-component + * dofLabel: 5 = displacement bubble function, z-component + * + * Ingredients: + * 1. F-points w (3 - 4 - 5), C-points displacement (0 - 1 - 2) + * 2. F-points smoother: block Jacobi + * 3. C-points coarse-grid/Schur complement solver: boomer AMG + * 4. Global smoother: none + * + */ +class AugmentedLagrangianContactMechanics : public MGRStrategyBase< 1 > +{ +public: + + /** + * @brief Constructor. + */ + explicit AugmentedLagrangianContactMechanics( arrayView1d< int const > const & ) + : MGRStrategyBase( 6 ) + { + // Level 0: we keep all three displacement + m_labels[0] = { 0, 1, 2 }; + + setupLabels(); + + // Level 0 + m_levelFRelaxType[0] = MGRFRelaxationType::l1jacobi; + m_levelFRelaxIters[0] = 1; + m_levelInterpType[0] = MGRInterpolationType::blockJacobi; + m_levelRestrictType[0] = MGRRestrictionType::injection; + m_levelCoarseGridMethod[0] = MGRCoarseGridMethod::galerkin; + m_levelGlobalSmootherType[0] = MGRGlobalSmootherType::none; + } + + /** + * @brief Setup the MGR strategy. + * @param mgrParams MGR configuration parameters + * @param precond preconditioner wrapper + * @param mgrData auxiliary MGR data + */ + void setup( LinearSolverParameters::MGR const & mgrParams, + HyprePrecWrapper & precond, + HypreMGRData & mgrData ) + { + setReduction( precond, mgrData ); + + // Configure the BoomerAMG solver used as mgr coarse solver for the displacement reduced system + // (note that no separate displacement component approach is used here) + setDisplacementAMG( mgrData.coarseSolver, mgrParams.separateComponents ); + } +}; + +} // namespace mgr + +} // namespace hypre + +} // namespace geos + +#endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREMGRALMCONTACTMECHANICS_HPP_*/ diff --git a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp index 4d43955a0e1..f46d85a57cd 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -295,6 +295,7 @@ struct LinearSolverParameters thermalMultiphasePoromechanics, ///< thermal multiphase poromechanics with finite volume compositional multiphase flow hydrofracture, ///< hydrofracture lagrangianContactMechanics, ///< Lagrangian contact mechanics + augmentedLagrangianContactMechanics, ///< Augmented Lagrangian contact mechanics lagrangianContactMechanicsBubbleStab, ///< Lagrangian contact mechanics with bubble stabilization solidMechanicsEmbeddedFractures ///< Embedded fractures mechanics }; @@ -388,6 +389,7 @@ ENUM_STRINGS( LinearSolverParameters::MGR::StrategyType, "thermalMultiphasePoromechanics", "hydrofracture", "lagrangianContactMechanics", + "augmentedLagrangianContactMechanics", "lagrangianContactMechanicsBubbleStab", "solidMechanicsEmbeddedFractures" ); diff --git a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp index d9c054d3401..ef6f7750ff1 100644 --- a/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp +++ b/src/coreComponents/mesh/utilities/ComputationalGeometry.hpp @@ -24,7 +24,7 @@ #include "common/DataLayouts.hpp" #include "finiteElement/elementFormulations/H1_Hexahedron_Lagrange1_GaussLegendre2.hpp" #include "finiteElement/elementFormulations/H1_Pyramid_Lagrange1_Gauss5.hpp" -#include "finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp" +#include "finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp" #include "finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp" #include "LvArray/src/output.hpp" #include "LvArray/src/tensorOps.hpp" diff --git a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp index 2c61762010c..7528f02da14 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp +++ b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp @@ -87,8 +87,16 @@ DECLARE_FIELD( slip, array1d< real64 >, 0, LEVEL_0, - NO_WRITE, - "Slip." ); + WRITE_AND_READ, + "Slip" ); + +DECLARE_FIELD( tangentialTraction, + "tangentialTraction", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Tangential traction" ); DECLARE_FIELD( deltaSlip, "deltaSlip", diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp index c1c886aa3dd..edd3bd4ec24 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp @@ -99,6 +99,8 @@ void ContactSolverBase::registerDataOnMesh( dataRepository::Group & meshBodies ) subRegion.registerField< fields::contact::slip >( getName() ); + subRegion.registerField< fields::contact::tangentialTraction >( getName() ); + subRegion.registerField< fields::contact::deltaSlip >( getName() ). setDimLabels( 1, labelsTangent ).reference().resizeDimension< 1 >( 2 ); diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp index c940cb45289..090f551ca81 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -46,14 +46,59 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta m_faceTypeToFiniteElements["Quadrilateral"] = std::make_unique< finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2 >(); m_faceTypeToFiniteElements["Triangle"] = std::make_unique< finiteElement::H1_TriangleFace_Lagrange1_Gauss1 >(); - LinearSolverParameters & linParams = m_linearSolverParameters.get(); + registerWrapper( viewKeyStruct::simultaneousString(), &m_simultaneous ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1 ). + setDescription( "Flag to update the Lagrange multiplier at each Newton iteration (true), or only after the Newton loop has converged (false)" ); + + registerWrapper( viewKeyStruct::symmetricString(), &m_symmetric ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1 ). + setDescription( "Flag to neglect the non-symmetric contribution in the tangential matrix" ); + + registerWrapper( viewKeyStruct::iterativePenaltyNFacString(), &m_iterPenaltyNFac ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 10.0 ). + setDescription( "Factor for tuning the iterative penalty coefficient for normal traction" ); + + registerWrapper( viewKeyStruct::iterativePenaltyTFacString(), &m_iterPenaltyTFac ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0.1 ). + setDescription( "Factor for tuning the iterative penalty coefficient for tangential traction" ); + + registerWrapper( viewKeyStruct::tolJumpDispNFacString(), &m_tolJumpDispNFac ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1.e-07 ). + setDescription( "Factor to adjust the tolerance for normal jump" ); + + registerWrapper( viewKeyStruct::tolJumpDispTFacString(), &m_tolJumpDispTFac ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 1.e-05 ). + setDescription( "Factor to adjust the tolerance for tangential jump" ); + + registerWrapper( viewKeyStruct::tolNormalTracFacString(), &m_tolNormalTracFac ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0.5 ). + setDescription( "Factor to adjust the tolerance for normal traction" ); + + registerWrapper( viewKeyStruct::tolTauLimitString(), &m_slidingCheckTolerance ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 5.e-02 ). + setDescription( "Tolerance for the sliding check" ); + + // Set the default linear solver parameters + LinearSolverParameters & linSolParams = m_linearSolverParameters.get(); addLogLevel< logInfo::Configuration >(); - linParams.isSymmetric = true; - linParams.dofsPerNode = 3; - linParams.mgr.separateComponents = true; - // TODO Implement the MGR strategy - //linParams.mgr.strategy = LinearSolverParameters::MGR::StrategyType::solidMechanicsAugumentedLagrangianContact; + // Strategy: AMG with separate displacement components + linSolParams.dofsPerNode = 3; + linSolParams.isSymmetric = true; + linSolParams.amg.separateComponents = true; + + // Strategy: static condensation of bubble dofs using MGR + linSolParams.mgr.strategy = LinearSolverParameters::MGR::StrategyType::augmentedLagrangianContactMechanics; + linSolParams.mgr.separateComponents = true; + } SolidMechanicsAugmentedLagrangianContact::~SolidMechanicsAugmentedLagrangianContact() @@ -531,16 +576,34 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepComplete( real64 cons arrayView2d< real64 > const oldDispJump = subRegion.getField< contact::oldDispJump >(); arrayView2d< real64 > const deltaDispJump = subRegion.getField< contact::deltaDispJump >(); + arrayView2d< real64 > const traction = subRegion.getField< contact::traction >(); + arrayView1d< integer const > const fractureState = subRegion.getField< contact::fractureState >(); arrayView1d< integer > const oldFractureState = subRegion.getField< contact::oldFractureState >(); + arrayView1d< real64 > const slip = subRegion.getField< contact::slip >(); + arrayView1d< real64 > const tangentialTraction = subRegion.getField< contact::tangentialTraction >(); + forAll< parallelDevicePolicy<> >( subRegion.size(), [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) { + // Compute the slip + real64 const deltaDisp[2] = { deltaDispJump[kfe][1], + deltaDispJump[kfe][2] }; + slip[kfe] = LvArray::tensorOps::l2Norm< 2 >( deltaDisp ); + + // Compute current Tau and limit Tau + real64 const tau[2] = { traction[kfe][1], + traction[kfe][2] }; + tangentialTraction[kfe] = LvArray::tensorOps::l2Norm< 2 >( tau ); + LvArray::tensorOps::fill< 3 >( deltaDispJump[kfe], 0.0 ); LvArray::tensorOps::copy< 3 >( oldDispJump[kfe], dispJump[kfe] ); oldFractureState[kfe] = fractureState[kfe]; + + + } ); } ); @@ -783,8 +846,6 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit arrayView1d< real64 const > const & normalTractionTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::normalTractionToleranceString() ); - arrayView1d< real64 const > const faceElementArea = subRegion.getElementArea().toViewConst(); - std::ptrdiff_t const sizes[ 2 ] = {subRegion.size(), 3}; traction_new.resize( 2, sizes ); arrayView2d< real64 > const traction_new_v = traction_new.toView(); @@ -806,7 +867,6 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit traction, dispJump, deltaDispJump, - faceElementArea, traction_new_v ); } else @@ -818,7 +878,6 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit traction, dispJump, deltaDispJump, - faceElementArea, traction_new_v ); } } ); @@ -964,8 +1023,6 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit arrayView2d< real64 const > const deltaDispJump = subRegion.getField< contact::deltaDispJump >(); - arrayView1d< real64 const > const faceElementArea = subRegion.getField< fields::elementArea >(); - constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) { using FrictionType = TYPEOFREF( castedFrictionLaw ); @@ -977,7 +1034,6 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit oldDispJump, dispJump, iterativePenalty, - faceElementArea, m_symmetric, normalTractionTolerance, traction, @@ -1113,6 +1169,8 @@ void SolidMechanicsAugmentedLagrangianContact::createFaceTypeList( DomainPartiti SurfaceElementRegion const & region = elemManager.getRegion< SurfaceElementRegion >( getUniqueFractureRegionName() ); FaceElementSubRegion const & subRegion = region.getUniqueSubRegion< FaceElementSubRegion >(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); + array1d< localIndex > keys( subRegion.size()); array1d< localIndex > vals( subRegion.size()); array1d< localIndex > quadList; @@ -1127,8 +1185,8 @@ void SolidMechanicsAugmentedLagrangianContact::createFaceTypeList( DomainPartiti forAll< parallelDevicePolicy<> >( subRegion.size(), [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) { - - localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kfe ); + localIndex const kf0 = elemsToFaces[kfe][0]; + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kf0 ); if( numNodesPerFace == 3 ) { keys_v[kfe]=0; @@ -1287,6 +1345,19 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti arrayView2d< localIndex > const faceElemsList_v = faceElemsList.toView(); +////////////////////////////////////////////////////////////////////////////////////////// +/* + FaceManager const & faceManager = mesh.getFaceManager(); + ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); + NodeManager const & nodeManager = mesh.getNodeManager(); + + arrayView2d< localIndex const > const elemsToNodeMap = cellElementSubRegion.nodeList().toViewConst(); + arrayView2d< real64 const > const elemsCenter = cellElementSubRegion.getElementCenter().toViewConst(); + + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = nodeManager.referencePosition(); + */ +////////////////////////////////////////////////////////////////////////////////////////// + forAll< parallelDevicePolicy<> >( nBubElems, [ = ] GEOS_HOST_DEVICE ( localIndex const k ) @@ -1294,6 +1365,57 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti localIndex const kfe = bubbleElemsList_v[k]; faceElemsList_v[k][0] = elemsToFaces[kfe][keys_v[k]]; faceElemsList_v[k][1] = keys_v[k]; +////////////////////////////////////////////////////////////////////////////////////////// +/* + std::cout << "GlobID: " << faceElemsList_v[k][0] << " LocalID: " << faceElemsList_v[k][1] << std::endl; + std::cout << "FACE NODES: " << std::endl; + real64 x = 0; + real64 y = 0; + real64 z = 0; + for (int nod=0; nod<3; ++nod) + { + std::cout << faceToNodeMap( faceElemsList_v[k][0], nod) << " "; + x += nodePosition[faceToNodeMap( faceElemsList_v[k][0], nod)][0]; + y += nodePosition[faceToNodeMap( faceElemsList_v[k][0], nod)][1]; + z += nodePosition[faceToNodeMap( faceElemsList_v[k][0], nod)][2]; + } + std::cout << std::endl; + std::cout << "FACE CENTERS: " << std::endl; + std::cout << x/3 << " " << y/3 << " " << z/3 << std::endl; + std::cout << "ELEMENT CENTERS: " << kfe << std::endl; + std::cout << elemsCenter[kfe][0] << " " << elemsCenter[kfe][1] << " " << elemsCenter[kfe][2] << std::endl; + std::cout << "ELMENT NODES: " << std::endl; + if (faceElemsList_v[k][1] == 0) + { + std::cout << elemsToNodeMap(kfe, 0) << " "; + std::cout << elemsToNodeMap(kfe, 1) << " "; + std::cout << elemsToNodeMap(kfe, 3) << " "; + std::cout << std::endl; + } + else if (faceElemsList_v[k][1] == 1) + { + std::cout << elemsToNodeMap(kfe, 0) << " "; + std::cout << elemsToNodeMap(kfe, 2) << " "; + std::cout << elemsToNodeMap(kfe, 1) << " "; + std::cout << std::endl; + } + else if (faceElemsList_v[k][1] == 2) + { + std::cout << elemsToNodeMap(kfe, 0) << " "; + std::cout << elemsToNodeMap(kfe, 3) << " "; + std::cout << elemsToNodeMap(kfe, 2) << " "; + std::cout << std::endl; + } + else if (faceElemsList_v[k][1] == 3) + { + std::cout << elemsToNodeMap(kfe, 1) << " "; + std::cout << elemsToNodeMap(kfe, 2) << " "; + std::cout << elemsToNodeMap(kfe, 3) << " "; + std::cout << std::endl; + } + std::cout << std::endl; + */ +////////////////////////////////////////////////////////////////////////////////////////// } ); cellElementSubRegion.setFaceElementsList( faceElemsList.toViewConst()); @@ -1373,7 +1495,8 @@ void SolidMechanicsAugmentedLagrangianContact::addCouplingNumNonzeros( DomainPar for( localIndex kfe=0; kfe( rotatedInvStiffApprox, area ); // Finally, compute tolerances for the given fracture element - normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus / 2.e+7; + normalDisplacementTolerance[kfe] = rotatedInvStiffApprox[ 0 ][ 0 ] * averageYoungModulus * m_tolJumpDispNFac; slidingTolerance[kfe] = sqrt( pow( rotatedInvStiffApprox[ 1 ][ 1 ], 2 ) + - pow( rotatedInvStiffApprox[ 2 ][ 2 ], 2 )) * averageYoungModulus / 2.e+5; - normalTractionTolerance[kfe] = (1.0 / 2.0) * (averageConstrainedModulus / averageBoxSize0) * + pow( rotatedInvStiffApprox[ 2 ][ 2 ], 2 )) * averageYoungModulus * m_tolJumpDispTFac; + normalTractionTolerance[kfe] = m_tolNormalTracFac * (averageConstrainedModulus / averageBoxSize0) * (normalDisplacementTolerance[kfe]); - iterativePenalty[kfe][0] = 1e+01*averageConstrainedModulus/(averageBoxSize0*area); - iterativePenalty[kfe][1] = 1e-01*averageConstrainedModulus/(averageBoxSize0*area); + + GEOS_LOG_LEVEL( 2, + GEOS_FMT( "kfe: {}, normalDisplacementTolerance: {}, slidingTolerance: {}, normalTractionTolerance: {}", + kfe, normalDisplacementTolerance[kfe], slidingTolerance[kfe], normalTractionTolerance[kfe] )); + + iterativePenalty[kfe][0] = m_iterPenaltyNFac*averageConstrainedModulus/(averageBoxSize0); + iterativePenalty[kfe][1] = m_iterPenaltyTFac*averageConstrainedModulus/(averageBoxSize0); } } ); } diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 0ed5d059b8a..ce7f43edad0 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -239,18 +239,50 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase constexpr static char const * slidingToleranceString() { return "slidingTolerance"; } constexpr static char const * dispJumpUpdPenaltyString() { return "dispJumpUpdPenalty"; } + + constexpr static char const * simultaneousString() { return "simultaneous"; } + + constexpr static char const * symmetricString() { return "symmetric"; } + + constexpr static char const * iterativePenaltyNFacString() { return "iterPenaltyN"; } + + constexpr static char const * iterativePenaltyTFacString() { return "iterPenaltyT"; } + + constexpr static char const * tolJumpDispNFacString() { return "tolJumpN"; } + + constexpr static char const * tolJumpDispTFacString() { return "tolJumpT"; } + + constexpr static char const * tolNormalTracFacString() { return "tolNormalTrac"; } + + constexpr static char const * tolTauLimitString() { return "tolTauLimit"; } + }; /// Tolerance for the sliding check: the tangential traction must exceed (1 + m_slidingCheckTolerance) * t_lim to activate the sliding /// condition - real64 const m_slidingCheckTolerance = 0.05; + real64 m_slidingCheckTolerance = 5.e-02; /// Flag to update the Lagrange multiplier at each Newton iteration (true), or only after the Newton loop has converged (false) - bool m_simultaneous = true; + int m_simultaneous = 1; /// Flag to neglect the non-symmetric contribution in the tangential matrix, i.e., the derivative of tangential traction with respect to /// the normal displacement is neglected - bool m_symmetric = true; + int m_symmetric = 1; + + /// Factor for tuning the iterative penalty coefficient for normal traction + real64 m_iterPenaltyNFac = 10.0; + + /// Factor for tuning the iterative penalty coefficient for tangential traction + real64 m_iterPenaltyTFac = 0.1; + + /// Factor to adjust the tolerance for normal jump + real64 m_tolJumpDispNFac = 1.e-07; + + /// Factor to adjust the tolerance for tangential jump + real64 m_tolJumpDispTFac = 1.e-05; + + /// Factor to adjust the tolerance for normal traction + real64 m_tolNormalTracFac = 0.5; }; diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp index 9a76b95632b..41f39281a2d 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp @@ -42,7 +42,7 @@ #include "linearAlgebra/solvers/PreconditionerBlockJacobi.hpp" #include "linearAlgebra/solvers/BlockPreconditioner.hpp" #include "linearAlgebra/solvers/SeparateComponentPreconditioner.hpp" -#include "finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp" +#include "finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp" #include "finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp" #if defined( __INTEL_COMPILER ) diff --git a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernels.hpp b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernels.hpp index b93406cf877..0154a36173d 100644 --- a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernels.hpp @@ -265,7 +265,6 @@ class ALM : m_dispJump[k], m_penalty[k], m_traction[k], - m_faceArea[k], m_symmetric, m_symmetric, zero, @@ -274,6 +273,10 @@ class ALM : tractionNew, fractureState ); + // Divide localPenalty by area + real64 const fac = 1.0/m_faceArea[k]; + LvArray::tensorOps::scale< 3, 3 >( stack.localPenalty, fac ); + // transp(R) * Atu LvArray::tensorOps::Rij_eq_AkiBkj< 3, numUdofs, 3 >( matRRtAtu, stack.localRotationMatrix, stack.localAtu ); @@ -412,7 +415,6 @@ struct ComputeTractionKernel arrayView2d< real64 const > const & traction, arrayView2d< real64 const > const & dispJump, arrayView2d< real64 const > const & deltaDispJump, - arrayView1d< real64 const > const & faceArea, arrayView2d< real64 > const & tractionNew ) { @@ -420,7 +422,7 @@ struct ComputeTractionKernel { contactWrapper.updateTractionOnly( dispJump[k], deltaDispJump[k], - penalty[k], traction[k], faceArea[k], tractionNew[k] ); + penalty[k], traction[k], tractionNew[k] ); } ); } diff --git a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernelsBase.hpp b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernelsBase.hpp index b3f96f9f182..8ba57f7a05d 100644 --- a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernelsBase.hpp +++ b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernelsBase.hpp @@ -117,7 +117,6 @@ struct UpdateStateKernel arrayView2d< real64 const > const & oldDispJump, arrayView2d< real64 const > const & dispJump, arrayView2d< real64 > const & penalty, - arrayView1d< real64 const > const & faceArea, bool const symmetric, arrayView1d< real64 const > const & normalTractionTolerance, arrayView2d< real64 > const & traction, @@ -135,7 +134,6 @@ struct UpdateStateKernel dispJump[k], penalty[k], traction[k], - faceArea[k], symmetric, false, normalTractionTolerance[k], diff --git a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp index 2769a101963..07f526cdece 100644 --- a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp @@ -276,9 +276,9 @@ class ALMSimultaneous : // Compute the trial traction real64 dispJump[ 3 ]; - dispJump[0] = stack.dispJumpLocal[0] * m_faceArea[k]; - dispJump[1] = ( stack.dispJumpLocal[1] - stack.oldDispJumpLocal[1] ) * m_faceArea[k]; - dispJump[2] = ( stack.dispJumpLocal[2] - stack.oldDispJumpLocal[2] ) * m_faceArea[k]; + dispJump[0] = stack.dispJumpLocal[0]; + dispJump[1] = ( stack.dispJumpLocal[1] - stack.oldDispJumpLocal[1] ); + dispJump[2] = ( stack.dispJumpLocal[2] - stack.oldDispJumpLocal[2] ); LvArray::tensorOps::scaledCopy< 3 >( tractionNew, stack.tLocal, -1.0 ); LvArray::tensorOps::Ri_add_AijBj< 3, 3 >( tractionNew, stack.localPenalty, dispJump ); @@ -308,6 +308,10 @@ class ALMSimultaneous : LvArray::tensorOps::Rij_eq_AikBkj< 3, numBdofs, 3 >( matDRtAtb, stack.localPenalty, matRRtAtb ); + real64 const fac = 1.0 / m_faceArea[k]; + LvArray::tensorOps::scale< 3, numUdofs >( matDRtAtu, fac ); + LvArray::tensorOps::scale< 3, numBdofs >( matDRtAtb, fac ); + // R*RtAtu LvArray::tensorOps::Rij_eq_AikBkj< 3, numUdofs, 3 >( matRRtAtu, stack.localRotationMatrix, matDRtAtu ); @@ -330,6 +334,7 @@ class ALMSimultaneous : LvArray::tensorOps::Rij_eq_AkiBkj< numUdofs, numBdofs, 3 >( stack.localAutAtb, stack.localAtu, matRRtAtb ); + // Compute the local residuals LvArray::tensorOps::scaledAdd< numUdofs >( stack.localRu, tractionR, 1 ); @@ -423,17 +428,16 @@ struct ComputeTractionSimultaneousKernel arrayView2d< real64 const > const & traction, arrayView2d< real64 const > const & dispJump, arrayView2d< real64 const > const & deltaDispJump, - arrayView1d< real64 const > const & faceElementArea, arrayView2d< real64 > const & tractionNew ) { forAll< POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const kfe ) { - tractionNew[kfe][0] = traction[kfe][0] + penalty[kfe][0] * dispJump[kfe][0] * faceElementArea[kfe]; + tractionNew[kfe][0] = traction[kfe][0] + penalty[kfe][0] * dispJump[kfe][0]; tractionNew[kfe][1] = traction[kfe][1] + ( penalty[kfe][2] * deltaDispJump[kfe][1]+ - penalty[kfe][4] * deltaDispJump[kfe][2] ) * faceElementArea[kfe]; + penalty[kfe][4] * deltaDispJump[kfe][2] ); tractionNew[kfe][2] = traction[kfe][2] + ( penalty[kfe][3] * deltaDispJump[kfe][2] + - penalty[kfe][4] * deltaDispJump[kfe][1] ) * faceElementArea[kfe]; + penalty[kfe][4] * deltaDispJump[kfe][1] ); } ); } diff --git a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsKernels.cmake b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsKernels.cmake index c3aa48ddb39..4c0e705c4bb 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsKernels.cmake +++ b/src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/PoromechanicsKernels.cmake @@ -31,6 +31,7 @@ set( porousSolidDispatch PorousSolid set( finiteElementDispatch H1_Hexahedron_Lagrange1_GaussLegendre2 H1_Wedge_Lagrange1_Gauss6 H1_Tetrahedron_Lagrange1_Gauss1 + H1_Tetrahedron_Lagrange1_Gauss14 H1_Pyramid_Lagrange1_Gauss5 H1_Tetrahedron_VEM_Gauss1 H1_Prism5_VEM_Gauss1 @@ -75,6 +76,7 @@ set( porousSolidDispatch PorousSolid ) set( finiteElementDispatch H1_Hexahedron_Lagrange1_GaussLegendre2 H1_Wedge_Lagrange1_Gauss6 H1_Tetrahedron_Lagrange1_Gauss1 + H1_Tetrahedron_Lagrange1_Gauss14 H1_Pyramid_Lagrange1_Gauss5 H1_Tetrahedron_VEM_Gauss1 H1_Prism5_VEM_Gauss1 diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsKernels.cmake b/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsKernels.cmake index c71c31bacd2..1110dfe768b 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsKernels.cmake +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsKernels.cmake @@ -35,6 +35,7 @@ set( solidBaseDispatch DamageSpectral set( finiteElementDispatch H1_Hexahedron_Lagrange1_GaussLegendre2 H1_Wedge_Lagrange1_Gauss6 H1_Tetrahedron_Lagrange1_Gauss1 + H1_Tetrahedron_Lagrange1_Gauss14 H1_Pyramid_Lagrange1_Gauss5 H1_Tetrahedron_VEM_Gauss1 H1_Prism5_VEM_Gauss1 diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp index aaf9430df2b..788c078ada3 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp @@ -761,23 +761,16 @@ int SurfaceGenerator::separationDriver( DomainPartition & domain, for( localIndex kfe=0; kfe