From 00157b35c45dd35ace8f1c8a416890efd7083d4f Mon Sep 17 00:00:00 2001 From: Matteo Frigo Date: Sun, 29 Sep 2024 16:00:54 -0700 Subject: [PATCH 01/18] Adding linear solving strategy --- .../linearAlgebra/CMakeLists.txt | 1 + .../interfaces/hypre/HypreMGR.cpp | 6 ++ .../AugmentedLagrangianContactMechanics.hpp | 100 ++++++++++++++++++ .../utilities/LinearSolverParameters.hpp | 2 + ...lidMechanicsAugmentedLagrangianContact.cpp | 24 ++++- ...lidMechanicsAugmentedLagrangianContact.hpp | 8 +- 6 files changed, 134 insertions(+), 7 deletions(-) create mode 100644 src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp diff --git a/src/coreComponents/linearAlgebra/CMakeLists.txt b/src/coreComponents/linearAlgebra/CMakeLists.txt index 095544ff4aa..43c75f726ef 100644 --- a/src/coreComponents/linearAlgebra/CMakeLists.txt +++ b/src/coreComponents/linearAlgebra/CMakeLists.txt @@ -102,6 +102,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/MultiphasePoromechanics.hpp interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsEmbeddedFractures.hpp interfaces/hypre/mgrStrategies/SinglePhasePoromechanicsConformingFractures.hpp diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp b/src/coreComponents/linearAlgebra/interfaces/hypre/HypreMGR.cpp index a14f0ac6b5b..73d7d64d7af 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" @@ -127,6 +128,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..28e074b5512 --- /dev/null +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp @@ -0,0 +1,100 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 precond preconditioner wrapper + * @param mgrData auxiliary MGR data + */ + void setup( LinearSolverParameters::MGR const &, + 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 ); + } +}; + +} // 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 90efd13f7ef..4a649ae9d68 100644 --- a/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp +++ b/src/coreComponents/linearAlgebra/utilities/LinearSolverParameters.hpp @@ -292,6 +292,7 @@ struct LinearSolverParameters thermalMultiphasePoromechanics, ///< thermal multiphase poromechanics with finite volume compositional multiphase flow hydrofracture, ///< hydrofracture lagrangianContactMechanics, ///< Lagrangian contact mechanics + augmentedLagrangianContactMechanics, ///< Augmented Lagrangian contact mechanics solidMechanicsEmbeddedFractures ///< Embedded fractures mechanics }; @@ -384,6 +385,7 @@ ENUM_STRINGS( LinearSolverParameters::MGR::StrategyType, "thermalMultiphasePoromechanics", "hydrofracture", "lagrangianContactMechanics", + "augmentedLagrangianContactMechanics", "solidMechanicsEmbeddedFractures" ); /// Declare strings associated with enumeration values. diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp index eb9fc81049e..32e34432e07 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -43,13 +43,27 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta m_faceTypeToFiniteElements["Quadrilateral"] = std::make_unique< finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2 >(); m_faceTypeToFiniteElements["Triangle"] = std::make_unique< finiteElement::H1_TriangleFace_Lagrange1_Gauss1 >(); - // TODO Implement the MGR strategy + 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" ); // Set the default linear solver parameters - //LinearSolverParameters & linParams = m_linearSolverParameters.get(); - //linParams.dofsPerNode = 3; - //linParams.isSymmetric = true; - //linParams.amg.separateComponents = true; + LinearSolverParameters & linSolParams = m_linearSolverParameters.get(); + // 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; + linSolParams.mgr.displacementFieldName = solidMechanics::totalDisplacement::key(); } SolidMechanicsAugmentedLagrangianContact::~SolidMechanicsAugmentedLagrangianContact() diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 54ffa40bab7..80208b9f943 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -239,6 +239,10 @@ 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"; } }; /// Tolerance for the sliding check: the tangential traction must exceed (1 + m_slidingCheckTolerance) * t_lim to activate the sliding @@ -246,11 +250,11 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase real64 const m_slidingCheckTolerance = 0.05; /// 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; }; From d9cd65e1555317a13b25bbe4fa504ba1e3e983e0 Mon Sep 17 00:00:00 2001 From: Matteo Frigo Date: Thu, 3 Oct 2024 11:41:29 -0700 Subject: [PATCH 02/18] Adding bubble functions for wedge elements --- .../H1_Wedge_Lagrange1_Gauss6.hpp | 67 ++++++++++++++++--- ...lidMechanicsAugmentedLagrangianContact.cpp | 12 ++-- 2 files changed, 67 insertions(+), 12 deletions(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp index fe349c4c393..282effab637 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp @@ -202,12 +202,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] = (1 - r - s) * s * LagrangeBasis1::valueBubble( xi ); + N[1] = (1 - r - s) * r * LagrangeBasis1::valueBubble( xi ); + N[2] = (1 - r - s) * r * s * LagrangeBasis1::value(0, xi ); + N[3] = (1 - r - s) * r * s * LagrangeBasis1::value(1, xi ); + N[4] = 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 @@ -218,8 +227,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 ); } /** @@ -609,9 +622,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] = - s * LagrangeBasis1::valueBubble( xi ); // dN0/dr + dNdXi[0][1] = (1.0 - r - 2.0 * s) * LagrangeBasis1::valueBubble( xi ); // dN0/ds + dNdXi[0][2] = (1 - r - s) * s * LagrangeBasis1::gradientBubble( xi ); // dN0/dxi + + dNdXi[1][0] = (1.0 - 2.0 * r - s) * LagrangeBasis1::valueBubble( xi ); // dN1/dr + dNdXi[1][1] = - r * LagrangeBasis1::valueBubble( xi ); // dN1/ds + dNdXi[1][2] = (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] = s * LagrangeBasis1::valueBubble( xi ); // dN4/dr + dNdXi[4][1] = r * LagrangeBasis1::valueBubble( xi ); // dN4/ds + dNdXi[4][2] = r * s * LagrangeBasis1::gradientBubble( xi ); // dN4/dxi + + for( int fi=0; fi( getUniqueFractureRegionName() ); FaceElementSubRegion const & subRegion = region.getUniqueSubRegion< FaceElementSubRegion >(); + ArrayOfArraysView< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); + array1d< localIndex > keys( subRegion.size()); array1d< localIndex > vals( subRegion.size()); array1d< localIndex > quadList; @@ -1139,8 +1141,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; @@ -1385,7 +1387,8 @@ void SolidMechanicsAugmentedLagrangianContact::addCouplingNumNonzeros( DomainPar for( localIndex kfe=0; kfe Date: Sat, 5 Oct 2024 11:32:38 -0700 Subject: [PATCH 03/18] Fixed face bubble functions for wedge elements --- .../H1_Wedge_Lagrange1_Gauss6.hpp | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp index 282effab637..55e812513ad 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp @@ -207,11 +207,11 @@ class H1_Wedge_Lagrange1_Gauss6 final : public FiniteElementBase real64 const s = pointCoord[1]; real64 const xi = pointCoord[2]; - N[0] = (1 - r - s) * s * LagrangeBasis1::valueBubble( xi ); - N[1] = (1 - r - s) * r * LagrangeBasis1::valueBubble( xi ); - N[2] = (1 - r - s) * r * s * LagrangeBasis1::value(0, xi ); - N[3] = (1 - r - s) * r * s * LagrangeBasis1::value(1, xi ); - N[4] = r * s * LagrangeBasis1::valueBubble( xi ); + 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 ); } @@ -635,13 +635,13 @@ H1_Wedge_Lagrange1_Gauss6::calcGradFaceBubbleN( localIndex const q, real64 const s = quadratureParentCoords1( q ); real64 const xi = quadratureParentCoords2( q ); - dNdXi[0][0] = - s * LagrangeBasis1::valueBubble( xi ); // dN0/dr - dNdXi[0][1] = (1.0 - r - 2.0 * s) * LagrangeBasis1::valueBubble( xi ); // dN0/ds - dNdXi[0][2] = (1 - r - s) * s * LagrangeBasis1::gradientBubble( xi ); // dN0/dxi + 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] = (1.0 - 2.0 * r - s) * LagrangeBasis1::valueBubble( xi ); // dN1/dr - dNdXi[1][1] = - r * LagrangeBasis1::valueBubble( xi ); // dN1/ds - dNdXi[1][2] = (1 - r - s) * r * LagrangeBasis1::gradientBubble( xi ); // dN1/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 @@ -651,9 +651,9 @@ H1_Wedge_Lagrange1_Gauss6::calcGradFaceBubbleN( localIndex const q, 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] = s * LagrangeBasis1::valueBubble( xi ); // dN4/dr - dNdXi[4][1] = r * LagrangeBasis1::valueBubble( xi ); // dN4/ds - dNdXi[4][2] = r * s * LagrangeBasis1::gradientBubble( xi ); // dN4/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 Date: Fri, 11 Oct 2024 16:42:41 -0700 Subject: [PATCH 04/18] Adding parameters to tune the iterative penalty coefficients --- .../SolidMechanicsAugmentedLagrangianContact.cpp | 14 ++++++++++++-- .../SolidMechanicsAugmentedLagrangianContact.hpp | 11 +++++++++++ 2 files changed, 23 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp index f949eaa5d5b..fc6d418fdeb 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -53,6 +53,16 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta 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" ); + // Set the default linear solver parameters LinearSolverParameters & linSolParams = m_linearSolverParameters.get(); // Strategy: AMG with separate displacement components @@ -1725,8 +1735,8 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio pow( rotatedInvStiffApprox[ 2 ][ 2 ], 2 )) * averageYoungModulus / 2.e+5; normalTractionTolerance[kfe] = (1.0 / 2.0) * (averageConstrainedModulus / averageBoxSize0) * (normalDisplacementTolerance[kfe]); - iterativePenalty[kfe][0] = 1e+01*averageConstrainedModulus/(averageBoxSize0*area); - iterativePenalty[kfe][1] = 1e-01*averageConstrainedModulus/(averageBoxSize0*area); + iterativePenalty[kfe][0] = m_iterPenaltyNFac*averageConstrainedModulus/(averageBoxSize0*area); + iterativePenalty[kfe][1] = m_iterPenaltyTFac*averageConstrainedModulus/(averageBoxSize0*area); } } ); } diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 80208b9f943..11ff099ea47 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -243,6 +243,11 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase constexpr static char const * simultaneousString() { return "simultaneous"; } constexpr static char const * symmetricString() { return "symmetric"; } + + constexpr static char const * iterativePenaltyNFacString() { return "epsNfac"; } + + constexpr static char const * iterativePenaltyTFacString() { return "epsTfac"; } + }; /// Tolerance for the sliding check: the tangential traction must exceed (1 + m_slidingCheckTolerance) * t_lim to activate the sliding @@ -256,6 +261,12 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase /// the normal displacement is neglected 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; + }; } /* namespace geos */ From 5a85a40e5226e7804129d35e2acc8fd1b21771ab Mon Sep 17 00:00:00 2001 From: mfrigo Date: Fri, 1 Nov 2024 13:36:49 -0700 Subject: [PATCH 05/18] Adding points for quadrature rule --- .../H1_Tetrahedron_Lagrange1_Gauss1.hpp | 236 +++++++++++++++--- .../H1_TriangleFace_Lagrange1_Gauss1.hpp | 114 ++++++++- .../physicsSolvers/contact/ContactFields.hpp | 10 +- ...lidMechanicsAugmentedLagrangianContact.cpp | 95 ++++++- ...lidMechanicsAugmentedLagrangianContact.hpp | 23 +- .../surfaceGeneration/SurfaceGenerator.cpp | 9 +- 6 files changed, 428 insertions(+), 59 deletions(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp index 92d7019857c..6f5e73843c0 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.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 * |/ \__\ | / ===== == == == @@ -63,7 +63,10 @@ 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 = 14; + //constexpr static localIndex numQuadraturePoints = 5; + //constexpr static localIndex numQuadraturePoints = 4; + //constexpr static localIndex numQuadraturePoints = 1; /// The number of sampling points per element. constexpr static int numSamplingPoints = numSamplingPointsPerDirection * numSamplingPointsPerDirection * numSamplingPointsPerDirection; @@ -221,10 +224,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 ); } @@ -321,7 +324,134 @@ 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 ) + { + + //real64 const w[numQuadraturePoints] = {-4.0/5.0, 9.0/20.0, 9.0/20.0, 9.0/20.0, 9.0/20.0 }; + //real64 const w[numQuadraturePoints] = {-4.0, 9.0/4.0, 9.0/4.0, 9.0/4.0, 9.0/4.0 }; + + + 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 ) + { + + //real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/2.0, 1.0/6.0, 1.0/6.0, 1.0/6.0 }; + + 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 ) + { + + //real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/2.0, 1.0/6.0, 1.0/6.0 }; + + 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 ) + { + + //real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/6.0, 1.0/2.0, 1.0/6.0 }; + + 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 @@ -375,10 +505,10 @@ H1_Tetrahedron_Lagrange1_Gauss1:: 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 ); } @@ -402,7 +532,6 @@ H1_Tetrahedron_Lagrange1_Gauss1:: 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] ); @@ -431,7 +560,7 @@ H1_Tetrahedron_Lagrange1_Gauss1:: } } - return detJ * weight; + return detJ * weight * quadratureWeight(q); } GEOS_HOST_DEVICE @@ -453,30 +582,31 @@ H1_Tetrahedron_Lagrange1_Gauss1::calcGradFaceBubbleN( localIndex const q, 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 ); + + //std::cout << detJ1 << " " << detJ << " " << q << " r: " << r << " s: " << s << " t: " << t << std::endl; dNdXi[0][0] = ( 1 - 2 * r - s - t ) * t; // dN0/dr dNdXi[0][1] = -r * t; // dN0/ds @@ -484,7 +614,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 @@ -501,7 +631,50 @@ 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; +/////////////////////////////////////////////////////////////////////////////////////////////////// +/* + real64 dNdXi1[numNodes][3] = {{0}}; + + dNdXi1[0][0] = -1; // dN0/dr + dNdXi1[0][1] = -1; // dN0/ds + dNdXi1[0][2] = -1; // dN0/dt + + dNdXi1[1][0] = 1; // dN1/dr + dNdXi1[1][1] = 0; // dN1/ds + dNdXi1[1][2] = 0; // dN1/dt + + dNdXi1[2][0] = 0; // dN2/dr + dNdXi1[2][1] = 1; // dN2/ds + dNdXi1[2][2] = 0; // dN2/dt + + dNdXi1[3][0] = 0; // dN3/dr + dNdXi1[3][1] = 0; // dN3/ds + dNdXi1[3][2] = 1; // dN3/dt + + real64 gradN1[numNodes][3] = {{0}}; + + for( int fi=0; fi, + 0, + LEVEL_0, + NO_WRITE, + "Tangential traction" ); DECLARE_FIELD( deltaDispJump, "deltaDisplacementJump", diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp index fc6d418fdeb..8f5f0198938 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -63,6 +63,26 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta 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( 1.e-05 ). + setDescription( "Tolerance for the sliding check" ); + // Set the default linear solver parameters LinearSolverParameters & linSolParams = m_linearSolverParameters.get(); // Strategy: AMG with separate displacement components @@ -1311,6 +1331,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 ) @@ -1318,6 +1351,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()); @@ -1730,11 +1814,16 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio LvArray::tensorOps::scale< 3, 3 >( 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]); + + GEOS_LOG_LEVEL( 2, + GEOS_FMT( "kfe: {}, normalDisplacementTolerance: {}, slidingTolerance: {}, normalTractionTolerance: {}", + kfe, normalDisplacementTolerance[kfe], slidingTolerance[kfe], normalTractionTolerance[kfe])); + iterativePenalty[kfe][0] = m_iterPenaltyNFac*averageConstrainedModulus/(averageBoxSize0*area); iterativePenalty[kfe][1] = m_iterPenaltyTFac*averageConstrainedModulus/(averageBoxSize0*area); } diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 11ff099ea47..2fbd601b6c5 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -244,15 +244,23 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase constexpr static char const * symmetricString() { return "symmetric"; } - constexpr static char const * iterativePenaltyNFacString() { return "epsNfac"; } + constexpr static char const * iterativePenaltyNFacString() { return "iterPenaltyN"; } - constexpr static char const * iterativePenaltyTFacString() { return "epsTfac"; } + 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 = 1.e-05; /// Flag to update the Lagrange multiplier at each Newton iteration (true), or only after the Newton loop has converged (false) int m_simultaneous = 1; @@ -267,6 +275,15 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase /// 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; + }; } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp index 9bcce6125f0..8645032eb2b 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp @@ -721,23 +721,16 @@ int SurfaceGenerator::separationDriver( DomainPartition & domain, for( localIndex kfe=0; kfe Date: Mon, 4 Nov 2024 15:39:18 -0800 Subject: [PATCH 06/18] Scaling dispJump by area --- .../physicsSolvers/contact/ContactFields.hpp | 4 +-- .../contact/ContactSolverBase.cpp | 2 ++ .../SolidMechanicsALMJumpUpdateKernels.hpp | 5 ++-- .../contact/SolidMechanicsALMKernels.hpp | 5 ++++ .../contact/SolidMechanicsALMKernelsBase.hpp | 11 ++++---- .../SolidMechanicsALMSimultaneousKernels.hpp | 16 ++++++------ ...lidMechanicsAugmentedLagrangianContact.cpp | 25 ++++++++++++++++--- 7 files changed, 48 insertions(+), 20 deletions(-) diff --git a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp index 33fc27fff22..191e8bba68d 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp +++ b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp @@ -79,7 +79,7 @@ DECLARE_FIELD( slip, array1d< real64 >, 0, LEVEL_0, - NO_WRITE, + WRITE_AND_READ, "Slip" ); DECLARE_FIELD( tangentialTraction, @@ -87,7 +87,7 @@ DECLARE_FIELD( tangentialTraction, array1d< real64 >, 0, LEVEL_0, - NO_WRITE, + WRITE_AND_READ, "Tangential traction" ); DECLARE_FIELD( deltaDispJump, diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp index 185a4c798c9..13ed1393343 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp @@ -91,6 +91,8 @@ void ContactSolverBase::registerDataOnMesh( dataRepository::Group & meshBodies ) subRegion.registerField< fields::contact::oldFractureState >( getName() ); subRegion.registerField< fields::contact::slip >( getName() ); + + subRegion.registerField< fields::contact::tangentialTraction >( getName() ); } ); } ); diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMJumpUpdateKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMJumpUpdateKernels.hpp index 1f434b3a13e..7bf92d2a644 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMJumpUpdateKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMJumpUpdateKernels.hpp @@ -54,6 +54,7 @@ class ALMJumpUpdate : using Base::m_faceToNodes; using Base::m_rotationMatrix; using Base::m_dispJump; + using Base::m_area; /** * @brief Constructor @@ -234,8 +235,8 @@ class ALMJumpUpdate : // Store the results for( int i=0; i<3; ++i ) { - m_dispJump[ k ][ i ] = stack.dispJumpLocal[ i ]; - m_deltaDispJump[ k ][ i ] = stack.deltaDispJumpLocal[ i ]; + m_dispJump[ k ][ i ] = stack.dispJumpLocal[ i ]/m_area[k]; + m_deltaDispJump[ k ][ i ] = stack.deltaDispJumpLocal[ i ]/m_area[k]; } return 0.0; diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp index 9f5c53f5132..539db85eb18 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernels.hpp @@ -63,6 +63,7 @@ class ALM : using Base::m_oldDispJump; using Base::m_matrix; using Base::m_rhs; + using Base::m_area; /** * @brief Constructor @@ -276,6 +277,10 @@ class ALM : tractionNew, fractureState ); + // Divide localPenalty by area + real64 const fac = 1.0/m_area[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 ); diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernelsBase.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernelsBase.hpp index 4f3504e3d5f..229f8b854d0 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernelsBase.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMKernelsBase.hpp @@ -98,7 +98,8 @@ class ALMKernelsBase : m_rotationMatrix( elementSubRegion.getField< fields::contact::rotationMatrix >().toViewConst()), m_dispJump( elementSubRegion.getField< fields::contact::dispJump >().toView() ), m_oldDispJump( elementSubRegion.getField< fields::contact::oldDispJump >().toViewConst() ), - m_penalty( elementSubRegion.getField< fields::contact::iterativePenalty >().toViewConst() ) + m_penalty( elementSubRegion.getField< fields::contact::iterativePenalty >().toViewConst() ), + m_area( elementSubRegion.getElementArea().toViewConst() ) {} //*************************************************************************** @@ -264,6 +265,8 @@ class ALMKernelsBase : /// The array containing the penalty coefficients for each element. arrayView2d< real64 const > const m_penalty; + /// The array containing the element area + arrayView1d< real64 const > const m_area; }; /** @@ -325,7 +328,6 @@ struct ConstraintCheckKernel * @param[in] normalDisplacementTolerance Check tolerance (compenetration) * @param[in] slidingTolerance Check tolerance (sliding) * @param[in] slidingCheckTolerance Check tolerance (if shear strass exceeds tauLim) - * @param[in] area interface element area * @param[in] fractureState the array containing the fracture state * @param[out] condConv the array containing the convergence flag: * 0: Constraint conditions satisfied @@ -346,7 +348,6 @@ struct ConstraintCheckKernel arrayView1d< real64 const > const & normalDisplacementTolerance, arrayView1d< real64 const > const & slidingTolerance, real64 const slidingCheckTolerance, - arrayView1d< real64 const > const & area, arrayView1d< integer const > const & fractureState, arrayView1d< integer > const & condConv ) { @@ -360,8 +361,8 @@ struct ConstraintCheckKernel traction[k], fractureState[k], normalTractionTolerance[k], - normalDisplacementTolerance[k]*area[k], - slidingTolerance[k]*area[k], + normalDisplacementTolerance[k], + slidingTolerance[k], slidingCheckTolerance, condConv[k] ); } diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp index f9cd307b6fe..1b9129f6561 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsALMSimultaneousKernels.hpp @@ -61,6 +61,7 @@ class ALMSimultaneous : using Base::m_oldDispJump; using Base::m_matrix; using Base::m_rhs; + using Base::m_area; /** * @brief Constructor @@ -223,18 +224,19 @@ class ALMSimultaneous : } // The minus sign is consistent with the sign of the Jacobian - stack.localPenalty[0][0] = -m_penalty( k, 0 ); + real64 const fac = 1.0 / m_area[k]; + stack.localPenalty[0][0] = -m_penalty( k, 0 )*fac; - stack.localPenalty[1][1] = -m_penalty( k, 2 ); - stack.localPenalty[2][2] = -m_penalty( k, 3 ); - stack.localPenalty[1][2] = -m_penalty( k, 4 ); - stack.localPenalty[2][1] = -m_penalty( k, 4 ); + stack.localPenalty[1][1] = -m_penalty( k, 2 )*fac; + stack.localPenalty[2][2] = -m_penalty( k, 3 )*fac; + stack.localPenalty[1][2] = -m_penalty( k, 4 )*fac; + stack.localPenalty[2][1] = -m_penalty( k, 4 )*fac; for( int i=0; i 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]; + + + } ); } ); @@ -829,7 +847,7 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit arrayView1d< real64 const > const & normalTractionTolerance = subRegion.getReference< array1d< real64 > >( viewKeyStruct::normalTractionToleranceString() ); - arrayView1d< real64 const > const area = subRegion.getElementArea().toViewConst(); + //arrayView1d< real64 const > const area = subRegion.getElementArea().toViewConst(); std::ptrdiff_t const sizes[ 2 ] = {subRegion.size(), 3}; traction_new.resize( 2, sizes ); @@ -885,7 +903,6 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit normalDisplacementTolerance, slidingTolerance, slidingCheckTolerance, - area, fractureState, condConv_v ); } ); @@ -1824,8 +1841,8 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio GEOS_FMT( "kfe: {}, normalDisplacementTolerance: {}, slidingTolerance: {}, normalTractionTolerance: {}", kfe, normalDisplacementTolerance[kfe], slidingTolerance[kfe], normalTractionTolerance[kfe])); - iterativePenalty[kfe][0] = m_iterPenaltyNFac*averageConstrainedModulus/(averageBoxSize0*area); - iterativePenalty[kfe][1] = m_iterPenaltyTFac*averageConstrainedModulus/(averageBoxSize0*area); + iterativePenalty[kfe][0] = m_iterPenaltyNFac*averageConstrainedModulus/(averageBoxSize0); + iterativePenalty[kfe][1] = m_iterPenaltyTFac*averageConstrainedModulus/(averageBoxSize0); } } ); } From 69ab77acf944b6eb873a85e4782bbc8c6c28da98 Mon Sep 17 00:00:00 2001 From: mfrigo Date: Wed, 13 Nov 2024 15:06:48 -0800 Subject: [PATCH 07/18] Added var to set the number of quadrature points for tetrahedron and triangle --- .../H1_Tetrahedron_Lagrange1_Gauss1.hpp | 103 ++++++++---------- .../H1_TriangleFace_Lagrange1_Gauss1.hpp | 75 ++++++++----- 2 files changed, 94 insertions(+), 84 deletions(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp index 6f5e73843c0..a9d1b7010de 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp @@ -29,6 +29,8 @@ namespace geos namespace finiteElement { +#define TETRAHEDRON_QUADRATURE_POINTS 14 + /** * This class contains the kernel accessible functions specific to the * H1-conforming nodal linear tetrahedron finite element with a 1-point @@ -63,10 +65,13 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static localIndex maxSupportPoints = numNodes; /// The number of quadrature points per element. +#if TETRAHEDRON_QUADRATURE_POINTS == 14 constexpr static localIndex numQuadraturePoints = 14; - //constexpr static localIndex numQuadraturePoints = 5; - //constexpr static localIndex numQuadraturePoints = 4; - //constexpr static localIndex numQuadraturePoints = 1; +#elif TETRAHEDRON_QUADRATURE_POINTS == 5 + constexpr static localIndex numQuadraturePoints = 5; +#elif TETRAHEDRON_QUADRATURE_POINTS == 1 + constexpr static localIndex numQuadraturePoints = 1; +#endif /// The number of sampling points per element. constexpr static int numSamplingPoints = numSamplingPointsPerDirection * numSamplingPointsPerDirection * numSamplingPointsPerDirection; @@ -337,10 +342,7 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureWeight( localIndex const q ) { - //real64 const w[numQuadraturePoints] = {-4.0/5.0, 9.0/20.0, 9.0/20.0, 9.0/20.0, 9.0/20.0 }; - //real64 const w[numQuadraturePoints] = {-4.0, 9.0/4.0, 9.0/4.0, 9.0/4.0, 9.0/4.0 }; - - +#if TETRAHEDRON_QUADRATURE_POINTS == 14 real64 const w[numQuadraturePoints] = { 0.073493043116361949544, 0.073493043116361949544, 0.073493043116361949544, @@ -356,6 +358,14 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase 0.042546020777081466438, 0.042546020777081466438 }; +#elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; + +#elif TETRAHEDRON_QUADRATURE_POINTS == 1 + real64 const w[numQuadraturePoints] = { 1.0 }; + +#endif + return w[q]; } @@ -371,8 +381,8 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureParentCoords0( localIndex const q ) { - //real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/2.0, 1.0/6.0, 1.0/6.0, 1.0/6.0 }; +#if TETRAHEDRON_QUADRATURE_POINTS == 14 real64 const qCoords[numQuadraturePoints] = { 0.72179424906732632079, 0.092735250310891226402, 0.092735250310891226402, @@ -388,6 +398,13 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase 0.45449629587435035051, 0.45449629587435035051 }; +#elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; + +#elif TETRAHEDRON_QUADRATURE_POINTS == 1 + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + +#endif return qCoords[q]; } @@ -403,8 +420,7 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureParentCoords1( localIndex const q ) { - //real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/2.0, 1.0/6.0, 1.0/6.0 }; - +#if TETRAHEDRON_QUADRATURE_POINTS == 14 real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402, 0.72179424906732632079, 0.092735250310891226402, @@ -420,6 +436,14 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase 0.045503704125649649492, 0.45449629587435035051 }; +#elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; + +#elif TETRAHEDRON_QUADRATURE_POINTS == 1 + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + +#endif + return qCoords[q]; } @@ -434,8 +458,7 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureParentCoords2( localIndex const q ) { - //real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/6.0, 1.0/2.0, 1.0/6.0 }; - +#if TETRAHEDRON_QUADRATURE_POINTS == 14 real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402, 0.092735250310891226402, 0.72179424906732632079, @@ -450,6 +473,15 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase 0.045503704125649649492, 0.45449629587435035051, 0.045503704125649649492 }; + +#elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; + +#elif TETRAHEDRON_QUADRATURE_POINTS == 1 + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + +#endif + return qCoords[q]; } @@ -606,8 +638,6 @@ H1_Tetrahedron_Lagrange1_Gauss1::calcGradFaceBubbleN( localIndex const q, real64 const s = quadratureParentCoords1( q ); real64 const t = quadratureParentCoords2( q ); - //std::cout << detJ1 << " " << detJ << " " << q << " r: " << r << " s: " << s << " t: " << t << std::endl; - dNdXi[0][0] = ( 1 - 2 * r - s - t ) * t; // dN0/dr dNdXi[0][1] = -r * t; // dN0/ds dNdXi[0][2] = ( 1 - r - s - 2 * t ) * r; // dN0/dt @@ -631,49 +661,6 @@ 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]; } -/////////////////////////////////////////////////////////////////////////////////////////////////// -/* - real64 dNdXi1[numNodes][3] = {{0}}; - - dNdXi1[0][0] = -1; // dN0/dr - dNdXi1[0][1] = -1; // dN0/ds - dNdXi1[0][2] = -1; // dN0/dt - - dNdXi1[1][0] = 1; // dN1/dr - dNdXi1[1][1] = 0; // dN1/ds - dNdXi1[1][2] = 0; // dN1/dt - - dNdXi1[2][0] = 0; // dN2/dr - dNdXi1[2][1] = 1; // dN2/ds - dNdXi1[2][2] = 0; // dN2/dt - - dNdXi1[3][0] = 0; // dN3/dr - dNdXi1[3][1] = 0; // dN3/ds - dNdXi1[3][2] = 1; // dN3/dt - - real64 gradN1[numNodes][3] = {{0}}; - - for( int fi=0; fi Date: Fri, 10 Jan 2025 02:51:39 -0800 Subject: [PATCH 08/18] Cleaning up SolidMechanicsALMSimultaneousKernels.hpp --- .../SolidMechanicsALMSimultaneousKernels.hpp | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp index 5a22f2155b4..fb8ea11b687 100644 --- a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp @@ -233,19 +233,18 @@ class ALMSimultaneous : } // The minus sign is consistent with the sign of the Jacobian - real64 const fac = 1.0 / m_faceArea[k]; - stack.localPenalty[0][0] = -m_penalty( k, 0 )*fac; + stack.localPenalty[0][0] = -m_penalty( k, 0 ); - stack.localPenalty[1][1] = -m_penalty( k, 2 )*fac; - stack.localPenalty[2][2] = -m_penalty( k, 3 )*fac; - stack.localPenalty[1][2] = -m_penalty( k, 4 )*fac; - stack.localPenalty[2][1] = -m_penalty( k, 4 )*fac; + stack.localPenalty[1][1] = -m_penalty( k, 2 ); + stack.localPenalty[2][2] = -m_penalty( k, 3 ); + stack.localPenalty[1][2] = -m_penalty( k, 4 ); + stack.localPenalty[2][1] = -m_penalty( k, 4 ); for( int i=0; i( 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 ); @@ -331,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 ); From e5ff26189ad41577eb2a9b580ecf25b75f07ca70 Mon Sep 17 00:00:00 2001 From: mfrigo Date: Fri, 10 Jan 2025 08:25:37 -0800 Subject: [PATCH 09/18] uncrustify_style --- .../H1_Tetrahedron_Lagrange1_Gauss1.hpp | 36 +++++++++---------- .../H1_TriangleFace_Lagrange1_Gauss1.hpp | 14 ++++---- .../H1_Wedge_Lagrange1_Gauss6.hpp | 16 ++++----- ...lidMechanicsAugmentedLagrangianContact.cpp | 32 ++++++++--------- .../kernels/SolidMechanicsALMKernels.hpp | 2 +- .../SolidMechanicsALMSimultaneousKernels.hpp | 4 +-- 6 files changed, 52 insertions(+), 52 deletions(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp index ee4a6a36df3..f1ae8908671 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp @@ -348,10 +348,10 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase 0.073493043116361949544, 0.073493043116361949544, 0.073493043116361949544, - 0.11268792571801585080, - 0.11268792571801585080, - 0.11268792571801585080, - 0.11268792571801585080, + 0.11268792571801585080, + 0.11268792571801585080, + 0.11268792571801585080, + 0.11268792571801585080, 0.042546020777081466438, 0.042546020777081466438, 0.042546020777081466438, @@ -422,23 +422,23 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase { #if TETRAHEDRON_QUADRATURE_POINTS == 14 - real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402, + real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402, 0.72179424906732632079, - 0.092735250310891226402, - 0.092735250310891226402, + 0.092735250310891226402, + 0.092735250310891226402, 0.31088591926330060980, 0.067342242210098170608, 0.31088591926330060980, 0.31088591926330060980, - 0.045503704125649649492, - 0.45449629587435035051, - 0.45449629587435035051, - 0.045503704125649649492, - 0.045503704125649649492, + 0.045503704125649649492, + 0.45449629587435035051, + 0.45449629587435035051, + 0.045503704125649649492, + 0.045503704125649649492, 0.45449629587435035051 }; #elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0, 1.0/6.0, 1.0/2.0, 1.0/6.0, 1.0/6.0 }; #elif TETRAHEDRON_QUADRATURE_POINTS == 1 real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; @@ -471,9 +471,9 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase 0.45449629587435035051, 0.045503704125649649492, 0.45449629587435035051, - 0.045503704125649649492, + 0.045503704125649649492, 0.45449629587435035051, - 0.045503704125649649492 }; + 0.045503704125649649492 }; #elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; @@ -593,7 +593,7 @@ H1_Tetrahedron_Lagrange1_Gauss1:: } } - return detJ * weight * quadratureWeight(q); + return detJ * weight * quadratureWeight( q ); } GEOS_HOST_DEVICE @@ -662,7 +662,7 @@ 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 * quadratureWeight(q); + return detJ * weight * quadratureWeight( q ); } //************************************************************************************************* @@ -677,7 +677,7 @@ H1_Tetrahedron_Lagrange1_Gauss1:: real64 detJ = determinantJacobianTransformation( X ); - return detJ * weight * quadratureWeight(q); + return detJ * weight * quadratureWeight( q ); } /// @endcond diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp index 6b3809ca26f..7ab18e27e0c 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp @@ -182,7 +182,7 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase real64 (& N)[1] ) { - real64 const qCoords[2] = {quadratureParentCoords0(q), quadratureParentCoords1(q) }; + real64 const qCoords[2] = {quadratureParentCoords0( q ), quadratureParentCoords1( q ) }; calcBubbleN( qCoords, N ); } @@ -259,9 +259,9 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase #elif TRIANGLE_QUADRATURE_POINTS == 4 real64 const w[numQuadraturePoints] = {-0.562500000000000, - 0.520833333333333, - 0.520833333333333, - 0.520833333333333 }; + 0.520833333333333, + 0.520833333333333, + 0.520833333333333 }; #elif TRIANGLE_QUADRATURE_POINTS == 1 @@ -381,9 +381,9 @@ H1_TriangleFace_Lagrange1_Gauss1:: real64 (& N)[numNodes] ) { - real64 const qCoords[2] = {quadratureParentCoords0(q), quadratureParentCoords1(q) }; + real64 const qCoords[2] = {quadratureParentCoords0( q ), quadratureParentCoords1( q ) }; - calcN(qCoords, N); + calcN( qCoords, N ); } @@ -411,7 +411,7 @@ H1_TriangleFace_Lagrange1_Gauss1:: ( 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 * quadratureWeight(q); + return sqrt( n[0] * n[0] + n[1] * n[1] + n[2] * n[2] ) * weight * quadratureWeight( q ); } /// @endcond diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp index 88dedd4623b..d901a13572b 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Wedge_Lagrange1_Gauss6.hpp @@ -210,8 +210,8 @@ class H1_Wedge_Lagrange1_Gauss6 final : public FiniteElementBase 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[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 ); } @@ -644,13 +644,13 @@ H1_Wedge_Lagrange1_Gauss6::calcGradFaceBubbleN( localIndex const q, 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[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[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 diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 600a03d9fcb..7a162920649 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -85,7 +85,7 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta setInputFlag( InputFlags::OPTIONAL ). setApplyDefaultValue( 1.e-05 ). setDescription( "Tolerance for the sliding check" ); - + // Set the default linear solver parameters LinearSolverParameters & linSolParams = m_linearSolverParameters.get(); addLogLevel< logInfo::Configuration >(); @@ -1355,7 +1355,7 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti arrayView2d< real64 const > const elemsCenter = cellElementSubRegion.getElementCenter().toViewConst(); arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const nodePosition = nodeManager.referencePosition(); -*/ + */ ////////////////////////////////////////////////////////////////////////////////////////// forAll< parallelDevicePolicy<> >( nBubElems, @@ -1375,9 +1375,9 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti 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]; + 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; @@ -1385,36 +1385,36 @@ void SolidMechanicsAugmentedLagrangianContact::createBubbleCellList( DomainParti 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) + if (faceElemsList_v[k][1] == 0) { - std::cout << elemsToNodeMap(kfe, 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) + else if (faceElemsList_v[k][1] == 1) { - std::cout << elemsToNodeMap(kfe, 0) << " "; + 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) + else if (faceElemsList_v[k][1] == 2) { - std::cout << elemsToNodeMap(kfe, 0) << " "; + 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) + else if (faceElemsList_v[k][1] == 3) { - std::cout << elemsToNodeMap(kfe, 1) << " "; + 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()); @@ -1835,8 +1835,8 @@ void SolidMechanicsAugmentedLagrangianContact::computeTolerances( DomainPartitio (normalDisplacementTolerance[kfe]); GEOS_LOG_LEVEL( 2, - GEOS_FMT( "kfe: {}, normalDisplacementTolerance: {}, slidingTolerance: {}, normalTractionTolerance: {}", - kfe, normalDisplacementTolerance[kfe], slidingTolerance[kfe], normalTractionTolerance[kfe])); + 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/kernels/SolidMechanicsALMKernels.hpp b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernels.hpp index f14bb5f0d12..0154a36173d 100644 --- a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMKernels.hpp @@ -275,7 +275,7 @@ class ALM : // Divide localPenalty by area real64 const fac = 1.0/m_faceArea[k]; - LvArray::tensorOps::scale< 3, 3 >( stack.localPenalty, fac ); + LvArray::tensorOps::scale< 3, 3 >( stack.localPenalty, fac ); // transp(R) * Atu LvArray::tensorOps::Rij_eq_AkiBkj< 3, numUdofs, 3 >( matRRtAtu, stack.localRotationMatrix, diff --git a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp index fb8ea11b687..07f526cdece 100644 --- a/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp +++ b/src/coreComponents/physicsSolvers/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp @@ -309,8 +309,8 @@ class ALMSimultaneous : matRRtAtb ); real64 const fac = 1.0 / m_faceArea[k]; - LvArray::tensorOps::scale< 3, numUdofs >( matDRtAtu, fac); - LvArray::tensorOps::scale< 3, numBdofs >( matDRtAtb, fac); + 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, From 459aef6f0a63a692820adfc17350de736f7288e5 Mon Sep 17 00:00:00 2001 From: mfrigo Date: Mon, 13 Jan 2025 21:47:38 -0800 Subject: [PATCH 10/18] Cleaning up H1_TriangleFace_Lagrange1_Gauss.hpp and H1_Tetrahedron_Lagrange1_Gauss.hpp --- .../finiteElement/CMakeLists.txt | 4 +- .../finiteElement/FiniteElementDispatch.hpp | 4 +- ...hpp => H1_Tetrahedron_Lagrange1_Gauss.hpp} | 329 ++++++++++-------- ...pp => H1_TriangleFace_Lagrange1_Gauss.hpp} | 221 ++++++------ .../testH1_Tetrahedron_Lagrange1_Gauss1.cpp | 2 +- .../testH1_TriangleFace_Lagrange1_Gauss1.cpp | 2 +- .../mesh/utilities/ComputationalGeometry.hpp | 2 +- .../contact/SolidMechanicsLagrangeContact.cpp | 2 +- 8 files changed, 316 insertions(+), 250 deletions(-) rename src/coreComponents/finiteElement/elementFormulations/{H1_Tetrahedron_Lagrange1_Gauss1.hpp => H1_Tetrahedron_Lagrange1_Gauss.hpp} (67%) rename src/coreComponents/finiteElement/elementFormulations/{H1_TriangleFace_Lagrange1_Gauss1.hpp => H1_TriangleFace_Lagrange1_Gauss.hpp} (65%) 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/FiniteElementDispatch.hpp b/src/coreComponents/finiteElement/FiniteElementDispatch.hpp index 5db1789d6e8..2cc1590b61d 100644 --- a/src/coreComponents/finiteElement/FiniteElementDispatch.hpp +++ b/src/coreComponents/finiteElement/FiniteElementDispatch.hpp @@ -24,13 +24,13 @@ #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 \ diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp similarity index 67% rename from src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp rename to src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp index f1ae8908671..bc4246cf1e0 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" @@ -29,8 +29,6 @@ namespace geos namespace finiteElement { -#define TETRAHEDRON_QUADRATURE_POINTS 14 - /** * This class contains the kernel accessible functions specific to the * H1-conforming nodal linear tetrahedron finite element with a 1-point @@ -48,7 +46,8 @@ namespace finiteElement * 0 1 * */ -class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase +template< int NUM_Q_POINTS > +class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase { public: @@ -65,19 +64,13 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static localIndex maxSupportPoints = numNodes; /// The number of quadrature points per element. -#if TETRAHEDRON_QUADRATURE_POINTS == 14 - constexpr static localIndex numQuadraturePoints = 14; -#elif TETRAHEDRON_QUADRATURE_POINTS == 5 - constexpr static localIndex numQuadraturePoints = 5; -#elif TETRAHEDRON_QUADRATURE_POINTS == 1 - constexpr static localIndex numQuadraturePoints = 1; -#endif + constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS; /// 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 @@ -342,32 +335,40 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase inline constexpr static real64 quadratureWeight( localIndex const q ) { - -#if TETRAHEDRON_QUADRATURE_POINTS == 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 }; - -#elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; - -#elif TETRAHEDRON_QUADRATURE_POINTS == 1 - real64 const w[numQuadraturePoints] = { 1.0 }; - -#endif - - return w[q]; + if constexpr (NUM_Q_POINTS == 1) + { + real64 const w[numQuadraturePoints] = { 1.0 }; + return w[q]; + } + else if constexpr (NUM_Q_POINTS == 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 (NUM_Q_POINTS == 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]; + } + else + { + GEOS_UNUSED_VAR( q ); + GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); + return 0; + } } @@ -382,32 +383,41 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureParentCoords0( localIndex const q ) { + if constexpr (NUM_Q_POINTS == 1) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + return qCoords[q]; + } + else if constexpr (NUM_Q_POINTS == 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 (NUM_Q_POINTS == 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]; + } + else + { + GEOS_UNUSED_VAR( q ); + GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); + return 0; + } -#if TETRAHEDRON_QUADRATURE_POINTS == 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 }; - -#elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; - -#elif TETRAHEDRON_QUADRATURE_POINTS == 1 - real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; - -#endif - - return qCoords[q]; } /** @@ -421,31 +431,41 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureParentCoords1( localIndex const q ) { -#if TETRAHEDRON_QUADRATURE_POINTS == 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 }; - -#elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; - -#elif TETRAHEDRON_QUADRATURE_POINTS == 1 - real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; - -#endif - - return qCoords[q]; + if constexpr (NUM_Q_POINTS == 1) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + return qCoords[q]; + } + else if constexpr (NUM_Q_POINTS == 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 (NUM_Q_POINTS == 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]; + } + else + { + GEOS_UNUSED_VAR( q ); + GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); + return 0; + } + } /** @@ -459,31 +479,42 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureParentCoords2( localIndex const q ) { -#if TETRAHEDRON_QUADRATURE_POINTS == 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 }; - -#elif TETRAHEDRON_QUADRATURE_POINTS == 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 }; - -#elif TETRAHEDRON_QUADRATURE_POINTS == 1 - real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; - -#endif - - return qCoords[q]; + if constexpr (NUM_Q_POINTS == 1) + { + real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; + return qCoords[q]; + } + else if constexpr (NUM_Q_POINTS == 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 (NUM_Q_POINTS == 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]; + } + else + { + GEOS_UNUSED_VAR( q ); + GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); + return 0; + } + } /** @@ -499,11 +530,12 @@ class H1_Tetrahedron_Lagrange1_Gauss1 final : public FiniteElementBase /// @cond Doxygen_Suppress +template< int 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] ) ) @@ -512,12 +544,13 @@ H1_Tetrahedron_Lagrange1_Gauss1:: //************************************************************************************************* +template< int 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]; @@ -531,12 +564,13 @@ H1_Tetrahedron_Lagrange1_Gauss1:: } +template< int 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] ) { real64 const pointCoord[3] = {quadratureParentCoords0( q ), quadratureParentCoords1( q ), @@ -545,25 +579,27 @@ H1_Tetrahedron_Lagrange1_Gauss1:: calcN( pointCoord, N ); } +template< int 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< int 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] ) { 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] ); @@ -596,23 +632,25 @@ H1_Tetrahedron_Lagrange1_Gauss1:: return detJ * weight * quadratureWeight( q ); } +template< int 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< int 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] ) { //real64 detJ = determinantJacobianTransformation( X ); @@ -667,12 +705,13 @@ H1_Tetrahedron_Lagrange1_Gauss1::calcGradFaceBubbleN( localIndex const q, //************************************************************************************************* +template< int 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] ) { real64 detJ = determinantJacobianTransformation( X ); @@ -682,9 +721,11 @@ H1_Tetrahedron_Lagrange1_Gauss1:: /// @endcond -#undef TETRAHEDRON_QUADRATURE_POINTS +using H1_Tetrahedron_Lagrange1_Gauss1 = H1_Tetrahedron_Lagrange1_Gauss< 1 >; +using H1_Tetrahedron_Lagrange1_Gauss5 = H1_Tetrahedron_Lagrange1_Gauss< 5 >; +using H1_Tetrahedron_Lagrange1_Gauss14 = H1_Tetrahedron_Lagrange1_Gauss< 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 65% rename from src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss1.hpp rename to src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp index 7ab18e27e0c..58fa12ce640 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" @@ -29,8 +29,6 @@ namespace geos namespace finiteElement { -#define TRIANGLE_QUADRATURE_POINTS 4 - /** * This class contains the kernel accessible functions specific to the * H1-conforming nodal linear triangular face finite element with a @@ -49,7 +47,8 @@ namespace finiteElement * 0 1 * */ -class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase +template< int NUM_Q_POINTS > +class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase { public: @@ -62,16 +61,10 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static localIndex maxSupportPoints = numNodes; /// The number of quadrature points per element. -#if TRIANGLE_QUADRATURE_POINTS == 6 - constexpr static localIndex numQuadraturePoints = 6; -#elif TRIANGLE_QUADRATURE_POINTS == 4 - constexpr static localIndex numQuadraturePoints = 4; -#elif TRIANGLE_QUADRATURE_POINTS == 1 - constexpr static localIndex numQuadraturePoints = 1; -#endif + constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS; GEOS_HOST_DEVICE - virtual ~H1_TriangleFace_Lagrange1_Gauss1() override + virtual ~H1_TriangleFace_Lagrange1_Gauss() override {} GEOS_HOST_DEVICE @@ -248,28 +241,35 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureWeight( localIndex const q ) { -#if TRIANGLE_QUADRATURE_POINTS == 6 - real64 const w[numQuadraturePoints] = { 1.0/6.0, - 1.0/6.0, - 1.0/6.0, - 1.0/6.0, - 1.0/6.0, - 1.0/6.0 }; - -#elif TRIANGLE_QUADRATURE_POINTS == 4 - - real64 const w[numQuadraturePoints] = {-0.562500000000000, - 0.520833333333333, - 0.520833333333333, - 0.520833333333333 }; - -#elif TRIANGLE_QUADRATURE_POINTS == 1 - - real64 const w[numQuadraturePoints] = { 1.0 }; -#endif - - return w[q]; - + if constexpr (NUM_Q_POINTS == 1) + { + real64 const w[numQuadraturePoints] = { 1.0 }; + return w[q]; + } + else if constexpr (NUM_Q_POINTS == 4) + { + real64 const w[numQuadraturePoints] = {-0.562500000000000, + 0.520833333333333, + 0.520833333333333, + 0.520833333333333 }; + return w[q]; + } + else if constexpr (NUM_Q_POINTS == 6) + { + real64 const w[numQuadraturePoints] = { 1.0/6.0, + 1.0/6.0, + 1.0/6.0, + 1.0/6.0, + 1.0/6.0, + 1.0/6.0 }; + return w[q]; + } + else + { + GEOS_UNUSED_VAR( q ); + GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); + return 0; + } } /** @@ -283,27 +283,35 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureParentCoords0( localIndex const q ) { -#if TRIANGLE_QUADRATURE_POINTS == 6 - real64 const qCoords[numQuadraturePoints] = { 0.659027622374092, - 0.109039009072877, - 0.231933368553031, - 0.659027622374092, - 0.109039009072877, - 0.231933368553031 }; - -#elif TRIANGLE_QUADRATURE_POINTS == 4 - real64 const qCoords[numQuadraturePoints] = { 0.333333333333333, - 0.600000000000000, - 0.200000000000000, - 0.200000000000000 }; - -#elif TRIANGLE_QUADRATURE_POINTS == 1 - real64 const qCoords[numQuadraturePoints] = { 1/3 }; - -#endif - - - return qCoords[q]; + if constexpr (NUM_Q_POINTS == 1) + { + real64 const qCoords[numQuadraturePoints] = { 1/3 }; + return qCoords[q]; + } + else if constexpr (NUM_Q_POINTS == 4) + { + real64 const qCoords[numQuadraturePoints] = { 0.333333333333333, + 0.600000000000000, + 0.200000000000000, + 0.200000000000000 }; + return qCoords[q]; + } + else if constexpr (NUM_Q_POINTS == 6) + { + real64 const qCoords[numQuadraturePoints] = { 0.659027622374092, + 0.109039009072877, + 0.231933368553031, + 0.659027622374092, + 0.109039009072877, + 0.231933368553031 }; + return qCoords[q]; + } + else + { + GEOS_UNUSED_VAR( q ); + GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); + return 0; + } } /** @@ -317,53 +325,65 @@ class H1_TriangleFace_Lagrange1_Gauss1 final : public FiniteElementBase constexpr static real64 quadratureParentCoords1( localIndex const q ) { -#if TRIANGLE_QUADRATURE_POINTS == 6 - real64 const qCoords[numQuadraturePoints] = { 0.231933368553031, - 0.659027622374092, - 0.109039009072877, - 0.109039009072877, - 0.231933368553031, - 0.659027622374092 }; - -#elif TRIANGLE_QUADRATURE_POINTS == 4 - real64 const qCoords[numQuadraturePoints] = { 0.333333333333333, - 0.200000000000000, - 0.600000000000000, - 0.200000000000000 }; - -#elif TRIANGLE_QUADRATURE_POINTS == 1 - real64 const qCoords[numQuadraturePoints] = { 1/3 }; - -#endif - - return qCoords[q]; + if constexpr (NUM_Q_POINTS == 1) + { + real64 const qCoords[numQuadraturePoints] = { 1/3 }; + return qCoords[q]; + } + else if constexpr (NUM_Q_POINTS == 4) + { + real64 const qCoords[numQuadraturePoints] = { 0.333333333333333, + 0.200000000000000, + 0.600000000000000, + 0.200000000000000 }; + return qCoords[q]; + } + else if constexpr (NUM_Q_POINTS == 6) + { + real64 const qCoords[numQuadraturePoints] = { 0.231933368553031, + 0.659027622374092, + 0.109039009072877, + 0.109039009072877, + 0.231933368553031, + 0.659027622374092 }; + return qCoords[q]; + } + else + { + GEOS_UNUSED_VAR( q ); + GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); + return 0; + } } }; /// @cond Doxygen_Suppress + +template< int 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< int 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]; @@ -373,12 +393,13 @@ H1_TriangleFace_Lagrange1_Gauss1:: N[2] = s; } +template< int 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] ) { real64 const qCoords[2] = {quadratureParentCoords0( q ), quadratureParentCoords1( q ) }; @@ -387,24 +408,26 @@ H1_TriangleFace_Lagrange1_Gauss1:: } +template< int 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< int 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 ); 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] ), @@ -416,8 +439,10 @@ H1_TriangleFace_Lagrange1_Gauss1:: /// @endcond -#undef TRIANGLE_QUADRATURE_POINTS +using H1_TriangleFace_Lagrange1_Gauss1 = H1_TriangleFace_Lagrange1_Gauss< 1 >; +using H1_TriangleFace_Lagrange1_Gauss4 = H1_TriangleFace_Lagrange1_Gauss< 4 >; +using H1_TriangleFace_Lagrange1_Gauss6 = H1_TriangleFace_Lagrange1_Gauss< 6 >; } } -#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS1_HPP_ +#endif //GEOS_FINITEELEMENT_ELEMENTFORMULATIONS_H1TRIANGLEFACELAGRANGE1GAUSS_HPP_ diff --git a/src/coreComponents/finiteElement/unitTests/testH1_Tetrahedron_Lagrange1_Gauss1.cpp b/src/coreComponents/finiteElement/unitTests/testH1_Tetrahedron_Lagrange1_Gauss1.cpp index 6ff8c6ee6ec..ee6f25dfc5f 100644 --- a/src/coreComponents/finiteElement/unitTests/testH1_Tetrahedron_Lagrange1_Gauss1.cpp +++ b/src/coreComponents/finiteElement/unitTests/testH1_Tetrahedron_Lagrange1_Gauss1.cpp @@ -17,7 +17,7 @@ * @file testH1_Tetrahedron_Lagrange1_Gauss1.cpp */ -#include "finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss1.hpp" +#include "finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "gtest/gtest.h" diff --git a/src/coreComponents/finiteElement/unitTests/testH1_TriangleFace_Lagrange1_Gauss1.cpp b/src/coreComponents/finiteElement/unitTests/testH1_TriangleFace_Lagrange1_Gauss1.cpp index 058e20d2db2..d3b36e40df8 100644 --- a/src/coreComponents/finiteElement/unitTests/testH1_TriangleFace_Lagrange1_Gauss1.cpp +++ b/src/coreComponents/finiteElement/unitTests/testH1_TriangleFace_Lagrange1_Gauss1.cpp @@ -22,7 +22,7 @@ #include "gtest/gtest.h" #include -#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/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/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 ) From 5cd89d9b0cd4f3e74034160b850a6ff14124a8a4 Mon Sep 17 00:00:00 2001 From: mfrigo Date: Tue, 14 Jan 2025 13:47:48 -0800 Subject: [PATCH 11/18] Fixing some review comments --- .../FiniteElementDiscretization.cpp | 24 +++++- .../FiniteElementDiscretization.hpp | 4 + .../finiteElement/FiniteElementDispatch.hpp | 4 +- .../H1_Tetrahedron_Lagrange1_Gauss.hpp | 83 ++++++++----------- .../H1_TriangleFace_Lagrange1_Gauss.hpp | 69 +++++++-------- .../AugmentedLagrangianContactMechanics.hpp | 1 + .../kernels/SolidMechanicsKernels.cmake | 1 + 7 files changed, 97 insertions(+), 89 deletions(-) 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 2cc1590b61d..f7abaf86252 100644 --- a/src/coreComponents/finiteElement/FiniteElementDispatch.hpp +++ b/src/coreComponents/finiteElement/FiniteElementDispatch.hpp @@ -37,6 +37,7 @@ 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_Gauss.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp index bc4246cf1e0..2736955bff6 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_Tetrahedron_Lagrange1_Gauss.hpp @@ -46,11 +46,17 @@ namespace finiteElement * 0 1 * */ -template< int NUM_Q_POINTS > +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; @@ -64,7 +70,7 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase constexpr static localIndex maxSupportPoints = numNodes; /// The number of quadrature points per element. - constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS; + constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS::value; /// The number of sampling points per element. constexpr static int numSamplingPoints = numSamplingPointsPerDirection * numSamplingPointsPerDirection * numSamplingPointsPerDirection; @@ -335,17 +341,17 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase inline constexpr static real64 quadratureWeight( localIndex const q ) { - if constexpr (NUM_Q_POINTS == 1) + if constexpr (numQuadraturePoints == 1) { real64 const w[numQuadraturePoints] = { 1.0 }; return w[q]; } - else if constexpr (NUM_Q_POINTS == 5) + 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 (NUM_Q_POINTS == 14) + else if constexpr (numQuadraturePoints == 14) { real64 const w[numQuadraturePoints] = { 0.073493043116361949544, 0.073493043116361949544, @@ -363,12 +369,6 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase 0.042546020777081466438 }; return w[q]; } - else - { - GEOS_UNUSED_VAR( q ); - GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); - return 0; - } } @@ -383,17 +383,17 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase constexpr static real64 quadratureParentCoords0( localIndex const q ) { - if constexpr (NUM_Q_POINTS == 1) + if constexpr (numQuadraturePoints == 1) { real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; return qCoords[q]; } - else if constexpr (NUM_Q_POINTS == 5) + 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 (NUM_Q_POINTS == 14) + else if constexpr (numQuadraturePoints == 14) { real64 const qCoords[numQuadraturePoints] = { 0.72179424906732632079, 0.092735250310891226402, @@ -411,12 +411,6 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase 0.45449629587435035051 }; return qCoords[q]; } - else - { - GEOS_UNUSED_VAR( q ); - GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); - return 0; - } } @@ -431,17 +425,17 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase constexpr static real64 quadratureParentCoords1( localIndex const q ) { - if constexpr (NUM_Q_POINTS == 1) + if constexpr (numQuadraturePoints == 1) { real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; return qCoords[q]; } - else if constexpr (NUM_Q_POINTS == 5) + 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 (NUM_Q_POINTS == 14) + else if constexpr (numQuadraturePoints == 14) { real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402, 0.72179424906732632079, @@ -459,12 +453,6 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase 0.45449629587435035051 }; return qCoords[q]; } - else - { - GEOS_UNUSED_VAR( q ); - GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); - return 0; - } } @@ -479,17 +467,17 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase constexpr static real64 quadratureParentCoords2( localIndex const q ) { - if constexpr (NUM_Q_POINTS == 1) + if constexpr (numQuadraturePoints == 1) { real64 const qCoords[numQuadraturePoints] = { 1.0/4.0 }; return qCoords[q]; } - else if constexpr (NUM_Q_POINTS == 5) + 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 (NUM_Q_POINTS == 14) + else if constexpr (numQuadraturePoints == 14) { real64 const qCoords[numQuadraturePoints] = { 0.092735250310891226402, 0.092735250310891226402, @@ -508,12 +496,6 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase return qCoords[q]; } - else - { - GEOS_UNUSED_VAR( q ); - GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); - return 0; - } } @@ -530,7 +512,7 @@ class H1_Tetrahedron_Lagrange1_Gauss final : public FiniteElementBase /// @cond Doxygen_Suppress -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 @@ -544,7 +526,7 @@ determinantJacobianTransformation( real64 const (&X)[numNodes][3] ) //************************************************************************************************* -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void @@ -564,7 +546,7 @@ calcN( real64 const (&coords)[3], } -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void @@ -579,7 +561,7 @@ calcN( localIndex const q, calcN( pointCoord, N ); } -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: @@ -592,7 +574,7 @@ calcN( localIndex const q, //************************************************************************************************* -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 @@ -632,7 +614,7 @@ calcGradN( localIndex const q, return detJ * weight * quadratureWeight( q ); } -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >:: @@ -644,7 +626,7 @@ calcGradN( localIndex const q, return calcGradN( q, X, gradN ); } -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 @@ -705,7 +687,7 @@ H1_Tetrahedron_Lagrange1_Gauss< NUM_Q_POINTS >::calcGradFaceBubbleN( localIndex //************************************************************************************************* -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 @@ -721,9 +703,12 @@ transformedQuadratureWeight( localIndex const q, /// @endcond -using H1_Tetrahedron_Lagrange1_Gauss1 = H1_Tetrahedron_Lagrange1_Gauss< 1 >; -using H1_Tetrahedron_Lagrange1_Gauss5 = H1_Tetrahedron_Lagrange1_Gauss< 5 >; -using H1_Tetrahedron_Lagrange1_Gauss14 = H1_Tetrahedron_Lagrange1_Gauss< 14 >; +/// @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 > >; } } diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp index 58fa12ce640..7e564375518 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp @@ -47,11 +47,17 @@ namespace finiteElement * 0 1 * */ -template< int NUM_Q_POINTS > +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; @@ -61,7 +67,7 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase constexpr static localIndex maxSupportPoints = numNodes; /// The number of quadrature points per element. - constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS; + constexpr static localIndex numQuadraturePoints = NUM_Q_POINTS::value; GEOS_HOST_DEVICE virtual ~H1_TriangleFace_Lagrange1_Gauss() override @@ -241,12 +247,12 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase constexpr static real64 quadratureWeight( localIndex const q ) { - if constexpr (NUM_Q_POINTS == 1) + if constexpr (numQuadraturePoints == 1) { real64 const w[numQuadraturePoints] = { 1.0 }; return w[q]; } - else if constexpr (NUM_Q_POINTS == 4) + else if constexpr (numQuadraturePoints == 4) { real64 const w[numQuadraturePoints] = {-0.562500000000000, 0.520833333333333, @@ -254,7 +260,7 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase 0.520833333333333 }; return w[q]; } - else if constexpr (NUM_Q_POINTS == 6) + else if constexpr (numQuadraturePoints == 6) { real64 const w[numQuadraturePoints] = { 1.0/6.0, 1.0/6.0, @@ -264,12 +270,7 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase 1.0/6.0 }; return w[q]; } - else - { - GEOS_UNUSED_VAR( q ); - GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); - return 0; - } + } /** @@ -283,12 +284,12 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase constexpr static real64 quadratureParentCoords0( localIndex const q ) { - if constexpr (NUM_Q_POINTS == 1) + if constexpr (numQuadraturePoints == 1) { real64 const qCoords[numQuadraturePoints] = { 1/3 }; return qCoords[q]; } - else if constexpr (NUM_Q_POINTS == 4) + else if constexpr (numQuadraturePoints == 4) { real64 const qCoords[numQuadraturePoints] = { 0.333333333333333, 0.600000000000000, @@ -296,7 +297,7 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase 0.200000000000000 }; return qCoords[q]; } - else if constexpr (NUM_Q_POINTS == 6) + else if constexpr (numQuadraturePoints == 6) { real64 const qCoords[numQuadraturePoints] = { 0.659027622374092, 0.109039009072877, @@ -306,12 +307,7 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase 0.231933368553031 }; return qCoords[q]; } - else - { - GEOS_UNUSED_VAR( q ); - GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); - return 0; - } + } /** @@ -325,12 +321,12 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase constexpr static real64 quadratureParentCoords1( localIndex const q ) { - if constexpr (NUM_Q_POINTS == 1) + if constexpr (numQuadraturePoints == 1) { real64 const qCoords[numQuadraturePoints] = { 1/3 }; return qCoords[q]; } - else if constexpr (NUM_Q_POINTS == 4) + else if constexpr (numQuadraturePoints == 4) { real64 const qCoords[numQuadraturePoints] = { 0.333333333333333, 0.200000000000000, @@ -338,7 +334,7 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase 0.200000000000000 }; return qCoords[q]; } - else if constexpr (NUM_Q_POINTS == 6) + else if constexpr (numQuadraturePoints == 6) { real64 const qCoords[numQuadraturePoints] = { 0.231933368553031, 0.659027622374092, @@ -348,12 +344,7 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase 0.659027622374092 }; return qCoords[q]; } - else - { - GEOS_UNUSED_VAR( q ); - GEOS_ERROR( "NUM_Q_POINTS not available for this element type" ); - return 0; - } + } }; @@ -361,7 +352,7 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase /// @cond Doxygen_Suppress -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > template< localIndex NUMDOFSPERTRIALSUPPORTPOINT, bool UPPER > GEOS_HOST_DEVICE inline @@ -377,7 +368,7 @@ addGradGradStabilization( StackVariables const & stack, GEOS_UNUSED_VAR( scaleFactor ); } -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void @@ -393,7 +384,7 @@ calcN( real64 const (&coords)[2], N[2] = s; } -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void @@ -408,7 +399,7 @@ calcN( localIndex const q, } -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline void H1_TriangleFace_Lagrange1_Gauss< NUM_Q_POINTS >:: @@ -421,7 +412,7 @@ calcN( localIndex const q, //************************************************************************************************* -template< int NUM_Q_POINTS > +template< typename NUM_Q_POINTS > GEOS_HOST_DEVICE inline real64 @@ -439,9 +430,13 @@ transformedQuadratureWeight( localIndex const q, /// @endcond -using H1_TriangleFace_Lagrange1_Gauss1 = H1_TriangleFace_Lagrange1_Gauss< 1 >; -using H1_TriangleFace_Lagrange1_Gauss4 = H1_TriangleFace_Lagrange1_Gauss< 4 >; -using H1_TriangleFace_Lagrange1_Gauss6 = H1_TriangleFace_Lagrange1_Gauss< 6 >; +/// @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 > >; + } } diff --git a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp index 98a605cdb4b..4f6186ca84c 100644 --- a/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp +++ b/src/coreComponents/linearAlgebra/interfaces/hypre/mgrStrategies/AugmentedLagrangianContactMechanics.hpp @@ -76,6 +76,7 @@ class AugmentedLagrangianContactMechanics : public MGRStrategyBase< 1 > /** * @brief Setup the MGR strategy. + * @param mgrParams MGR configuration parameters * @param precond preconditioner wrapper * @param mgrData auxiliary MGR data */ 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 From 907dcae713bd486bc5bbe3313b143021821530c6 Mon Sep 17 00:00:00 2001 From: mfrigo Date: Tue, 14 Jan 2025 16:08:44 -0800 Subject: [PATCH 12/18] Fixing compilation issues --- .../poromechanicsKernels/PoromechanicsKernels.cmake | 2 ++ 1 file changed, 2 insertions(+) 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 From ab95729e3d09504fb0db1b7de0d9437679a29e1c Mon Sep 17 00:00:00 2001 From: mfrigo Date: Tue, 14 Jan 2025 17:43:32 -0800 Subject: [PATCH 13/18] Fixed bug (int-real type casting) --- .../H1_TriangleFace_Lagrange1_Gauss.hpp | 66 +++++++++---------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp index 7e564375518..a68c7bc6f2a 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp @@ -249,25 +249,25 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase if constexpr (numQuadraturePoints == 1) { - real64 const w[numQuadraturePoints] = { 1.0 }; + constexpr real64 w[numQuadraturePoints] = { 1.0 }; return w[q]; } else if constexpr (numQuadraturePoints == 4) { - real64 const w[numQuadraturePoints] = {-0.562500000000000, - 0.520833333333333, - 0.520833333333333, - 0.520833333333333 }; + constexpr real64 w[numQuadraturePoints] = {-0.562500000000000, + 0.520833333333333, + 0.520833333333333, + 0.520833333333333 }; return w[q]; } else if constexpr (numQuadraturePoints == 6) { - real64 const w[numQuadraturePoints] = { 1.0/6.0, - 1.0/6.0, - 1.0/6.0, - 1.0/6.0, - 1.0/6.0, - 1.0/6.0 }; + real64 const w[numQuadraturePoints] = { 0.166666666666666, + 0.166666666666666, + 0.166666666666666, + 0.166666666666666, + 0.166666666666666, + 0.166666666666666 }; return w[q]; } @@ -286,25 +286,25 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase if constexpr (numQuadraturePoints == 1) { - real64 const qCoords[numQuadraturePoints] = { 1/3 }; + constexpr real64 qCoords[numQuadraturePoints] = { 1.0/3.0 }; return qCoords[q]; } else if constexpr (numQuadraturePoints == 4) { - real64 const qCoords[numQuadraturePoints] = { 0.333333333333333, - 0.600000000000000, - 0.200000000000000, - 0.200000000000000 }; + constexpr real64 qCoords[numQuadraturePoints] = { 0.333333333333333, + 0.600000000000000, + 0.200000000000000, + 0.200000000000000 }; return qCoords[q]; } else if constexpr (numQuadraturePoints == 6) { - real64 const qCoords[numQuadraturePoints] = { 0.659027622374092, - 0.109039009072877, - 0.231933368553031, - 0.659027622374092, - 0.109039009072877, - 0.231933368553031 }; + constexpr real64 qCoords[numQuadraturePoints] = { 0.659027622374092, + 0.109039009072877, + 0.231933368553031, + 0.659027622374092, + 0.109039009072877, + 0.231933368553031 }; return qCoords[q]; } @@ -323,25 +323,25 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase if constexpr (numQuadraturePoints == 1) { - real64 const qCoords[numQuadraturePoints] = { 1/3 }; + constexpr real64 qCoords[numQuadraturePoints] = { 1.0/3.0 }; return qCoords[q]; } else if constexpr (numQuadraturePoints == 4) { - real64 const qCoords[numQuadraturePoints] = { 0.333333333333333, - 0.200000000000000, - 0.600000000000000, - 0.200000000000000 }; + constexpr real64 qCoords[numQuadraturePoints] = { 0.333333333333333, + 0.200000000000000, + 0.600000000000000, + 0.200000000000000 }; return qCoords[q]; } else if constexpr (numQuadraturePoints == 6) { - real64 const qCoords[numQuadraturePoints] = { 0.231933368553031, - 0.659027622374092, - 0.109039009072877, - 0.109039009072877, - 0.231933368553031, - 0.659027622374092 }; + constexpr real64 qCoords[numQuadraturePoints] = { 0.231933368553031, + 0.659027622374092, + 0.109039009072877, + 0.109039009072877, + 0.231933368553031, + 0.659027622374092 }; return qCoords[q]; } From 53eed41a895a192ffe44272321af688bd4a0984b Mon Sep 17 00:00:00 2001 From: mfrigo Date: Tue, 14 Jan 2025 19:45:09 -0800 Subject: [PATCH 14/18] uncrustify --- .../elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp index a68c7bc6f2a..61a3086598b 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_TriangleFace_Lagrange1_Gauss.hpp @@ -255,9 +255,9 @@ class H1_TriangleFace_Lagrange1_Gauss final : public FiniteElementBase else if constexpr (numQuadraturePoints == 4) { constexpr real64 w[numQuadraturePoints] = {-0.562500000000000, - 0.520833333333333, - 0.520833333333333, - 0.520833333333333 }; + 0.520833333333333, + 0.520833333333333, + 0.520833333333333 }; return w[q]; } else if constexpr (numQuadraturePoints == 6) From aba0a41778fb8cce14fda44dd7034fe9774ff1f5 Mon Sep 17 00:00:00 2001 From: mfrigo Date: Fri, 17 Jan 2025 13:33:54 -0800 Subject: [PATCH 15/18] Modified the default input parameters --- .../contact/SolidMechanicsAugmentedLagrangianContact.cpp | 2 +- .../contact/SolidMechanicsAugmentedLagrangianContact.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 7a162920649..090f551ca81 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -83,7 +83,7 @@ SolidMechanicsAugmentedLagrangianContact::SolidMechanicsAugmentedLagrangianConta registerWrapper( viewKeyStruct::tolTauLimitString(), &m_slidingCheckTolerance ). setInputFlag( InputFlags::OPTIONAL ). - setApplyDefaultValue( 1.e-05 ). + setApplyDefaultValue( 5.e-02 ). setDescription( "Tolerance for the sliding check" ); // Set the default linear solver parameters diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 543346449ed..ce7f43edad0 100644 --- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -260,7 +260,7 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase /// Tolerance for the sliding check: the tangential traction must exceed (1 + m_slidingCheckTolerance) * t_lim to activate the sliding /// condition - real64 m_slidingCheckTolerance = 1.e-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) int m_simultaneous = 1; From 31f145ecf744c8be2dcb01954b232103e7833d6f Mon Sep 17 00:00:00 2001 From: mfrigo Date: Wed, 22 Jan 2025 15:29:41 -0800 Subject: [PATCH 16/18] Update baseline --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 622c1da3d1b..932707643e2 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-9813-e80036d 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. From f642064786701fdc02b0e3578889009b03af8f4c Mon Sep 17 00:00:00 2001 From: mfrigo Date: Wed, 22 Jan 2025 17:04:27 -0800 Subject: [PATCH 17/18] Update .integrated_tests.yaml --- .integrated_tests.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 932707643e2..ea7e4b5b516 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3395-9813-e80036d + baseline: integratedTests/baseline_integratedTests-pr3395-9832-31f145e allow_fail: all: '' From a87c9f2a2e7f4766f35f451ce86c36e5a936284f Mon Sep 17 00:00:00 2001 From: mfrigo Date: Thu, 23 Jan 2025 13:46:30 -0800 Subject: [PATCH 18/18] Fix xml files --- inputFiles/inducedSeismicity/SCEC_BP6_QD_S_explicit.xml | 4 ++-- inputFiles/inducedSeismicity/SCEC_BP6_QD_S_implicit.xml | 4 ++-- inputFiles/inducedSeismicity/SCEC_BP7_QD_base.xml | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) 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 +