diff --git a/inputFiles/mortar/TwoDomainBeam.xml b/inputFiles/mortar/TwoDomainBeam.xml new file mode 100644 index 00000000000..05a9448c8e4 --- /dev/null +++ b/inputFiles/mortar/TwoDomainBeam.xml @@ -0,0 +1,266 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/constitutive/PVTPackage b/src/coreComponents/constitutive/PVTPackage index a2cdca31546..374b72962f7 160000 --- a/src/coreComponents/constitutive/PVTPackage +++ b/src/coreComponents/constitutive/PVTPackage @@ -1 +1 @@ -Subproject commit a2cdca315464da11f2bdcc902c894ed1340b8572 +Subproject commit 374b72962f7e605e74a1943c541e94d7b3f492a3 diff --git a/src/coreComponents/finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp b/src/coreComponents/finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp index 6aa3001b291..02174b9bf74 100644 --- a/src/coreComponents/finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp +++ b/src/coreComponents/finiteElement/elementFormulations/H1_QuadrilateralFace_Lagrange1_GaussLegendre2.hpp @@ -142,6 +142,18 @@ class H1_QuadrilateralFace_Lagrange1_GaussLegendre2 final : public FiniteElement StackVariables const & stack, real64 ( &N )[numNodes] ); + /** + * @brief Calculate shape functions gradients for each support point at a + * quadrature point. + * @param coords The parent coordinates at which to evaluate the shape function gradient (unused) + * @param gradN An array to pass back the shape function gradients for each support + * point. + */ + GEOS_HOST_DEVICE + inline + static void calcGradN( real64 const (&coords)[2], + real64 ( &gradN )[numNodes][2] ); + /** * @brief Calculate shape bubble functions values at a given point in the parent space. * @param pointCoord coordinates of the given point. @@ -298,8 +310,8 @@ H1_QuadrilateralFace_Lagrange1_GaussLegendre2:: for( localIndex a=0; a +GEOS_HOST_DEVICE +inline +void H1_TriangleFace_Lagrange1_Gauss< NUM_Q_POINTS >:: +calcGradN( real64 const (&coords)[2], + real64 (& gradN)[numNodes][2] ) +{ + GEOS_UNUSED_VAR( coords ); + + gradN[0][0] = -1.0; + gradN[0][1] = -1.0; + gradN[1][0] = 1.0; + gradN[1][1] = 0.0; + gradN[2][0] = 0.0; + gradN[2][1] = 1.0; +} + //************************************************************************************************* template< typename NUM_Q_POINTS > diff --git a/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt b/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt index ced15c47ba5..bf8e61c58d0 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt @@ -45,6 +45,7 @@ set( solidMechanicsSolvers_headers contact/SolidMechanicsLagrangeContact.hpp contact/SolidMechanicsLagrangeContactBubbleStab.hpp contact/SolidMechanicsAugmentedLagrangianContact.hpp + contact/SolidMechanicsMortarContact.hpp contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp contact/kernels/SolidMechanicsDisplacementJumpUpdateKernels.hpp contact/kernels/SolidMechanicsEFEMKernelsBase.hpp @@ -55,7 +56,8 @@ set( solidMechanicsSolvers_headers contact/kernels/SolidMechanicsALMKernels.hpp contact/kernels/SolidMechanicsConformingContactKernelsHelper.hpp contact/kernels/SolidMechanicsContactFaceBubbleKernels.hpp - contact/kernels/SolidMechanicsLagrangeContactKernels.hpp + contact/kernels/SolidMechanicsLagrangeContactKernels.hpp + contact/kernels/SolidMechanicsMortarContactKernels.hpp contact/LogLevelsInfo.hpp ) # Specify solver sources @@ -68,7 +70,8 @@ set( solidMechanicsSolvers_sources contact/SolidMechanicsEmbeddedFractures.cpp contact/SolidMechanicsLagrangeContact.cpp contact/SolidMechanicsLagrangeContactBubbleStab.cpp - contact/SolidMechanicsAugmentedLagrangianContact.cpp ) + contact/SolidMechanicsAugmentedLagrangianContact.cpp + contact/SolidMechanicsMortarContact.cpp ) #include( kernels/SolidMechanicsKernels.cmake) @@ -112,4 +115,4 @@ install( TARGETS solidMechanicsSolvers LIBRARY DESTINATION ${CMAKE_INSTALL_PREFI if( externalComponentDeps ) target_include_directories( solidMechanicsSolvers PUBLIC ${CMAKE_SOURCE_DIR}/externalComponents ) -endif() \ No newline at end of file +endif() diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/CMakeLists.txt b/src/coreComponents/physicsSolvers/solidMechanics/contact/CMakeLists.txt index 6d6457f7597..d7a31331209 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/CMakeLists.txt @@ -8,6 +8,7 @@ set( solidMechanicsSolvers_headers contact/SolidMechanicsLagrangeContact.hpp contact/SolidMechanicsLagrangeContactBubbleStab.hpp contact/SolidMechanicsAugmentedLagrangianContact.hpp + contact/SolidMechanicsMortarContact.hpp contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp contact/kernels/SolidMechanicsDisplacementJumpUpdateKernels.hpp contact/kernels/SolidMechanicsEFEMKernelsBase.hpp @@ -18,7 +19,8 @@ set( solidMechanicsSolvers_headers contact/kernels/SolidMechanicsALMKernels.hpp contact/kernels/SolidMechanicsConformingContactKernelsHelper.hpp contact/kernels/SolidMechanicsContactFaceBubbleKernels.hpp - contact/kernels/SolidMechanicsLagrangeContactKernels.hpp + contact/kernels/SolidMechanicsLagrangeContactKernels.hpp + conact/kernels/SolidMechanicsMortarContactKernels.hpp PARENT_SCOPE ) @@ -30,4 +32,5 @@ set( solidMechanicsSolvers_sources contact/SolidMechanicsLagrangeContact.cpp contact/SolidMechanicsLagrangeContactBubbleStab.cpp contact/SolidMechanicsAugmentedLagrangianContact.cpp - PARENT_SCOPE ) \ No newline at end of file + contact/SolidMechanicsMortarContact.cpp + PARENT_SCOPE ) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp index 469aa255c3d..94fa4649363 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/ContactSolverBase.cpp @@ -27,6 +27,7 @@ #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "common/GEOS_RAJA_Interface.hpp" #include "fieldSpecification/FieldSpecificationManager.hpp" +#include "physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.hpp" namespace geos { @@ -44,9 +45,6 @@ ContactSolverBase::ContactSolverBase( const string & name, this->getWrapper< string >( viewKeyStruct::surfaceGeneratorNameString() ). setInputFlag( dataRepository::InputFlags::FALSE ); - - addLogLevel< logInfo::ConfigurationStatistics >(); - addLogLevel< logInfo::ContactTolerance >(); } void ContactSolverBase::postInputInitialization() @@ -70,10 +68,15 @@ void ContactSolverBase::registerDataOnMesh( dataRepository::Group & meshBodies ) string const labels[3] = { "normal", "tangent1", "tangent2" }; string const labelsTangent[2] = { "tangent1", "tangent2" }; - forFractureRegionOnMeshTargets( meshBodies, [&] ( SurfaceElementRegion & fractureRegion ) + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + string_array const & ) { - fractureRegion.forElementSubRegions< SurfaceElementSubRegion >( [&]( SurfaceElementSubRegion & subRegion ) + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) { + subRegion.registerField< contact::dispJump >( getName() ). setDimLabels( 1, labels ). reference().resizeDimension< 1 >( 3 ); @@ -125,10 +128,14 @@ void ContactSolverBase::setFractureRegions( dataRepository::Group const & meshBo } ); // TODO remove once multiple regions are fully supported - GEOS_THROW_IF( m_fractureRegionNames.size() > 1, - GEOS_FMT( "{} {}: The number of fracture regions can not be more than one", - this->getCatalogName(), this->getName() ), - InputError ); + // Disable this check for mortar contact solver + if( m_fractureRegionNames.size() > 1 && !dynamic_cast< SolidMechanicsMortarContact * >( this ) ) + { + GEOS_THROW_IF( m_fractureRegionNames.size() > 1, + GEOS_FMT( "{} {}: The number of fracture regions can not be more than one", + this->getCatalogName(), this->getName() ), + InputError ); + } } void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh, @@ -200,21 +207,24 @@ void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh, void ContactSolverBase::outputConfigurationStatistics( DomainPartition const & domain ) const { - globalIndex numStick = 0; - globalIndex numNewSlip = 0; - globalIndex numSlip = 0; - globalIndex numOpen = 0; - - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, - MeshLevel const & mesh, - string_array const & ) + if( getLogLevel() >=1 ) { - computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen ); + globalIndex numStick = 0; + globalIndex numNewSlip = 0; + globalIndex numSlip = 0; + globalIndex numOpen = 0; + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel const & mesh, + string_array const & ) + { + computeFractureStateStatistics( mesh, numStick, numNewSlip, numSlip, numOpen ); - GEOS_LOG_RANK_0( GEOS_FMT( " Number of element for each fracture state:" - " stick = {:12} | new slip = {:12} | slip = {:12} | open = {:12}", - numStick, numNewSlip, numSlip, numOpen ) ); - } ); + GEOS_LOG_RANK_0( GEOS_FMT( " Number of element for each fracture state:" + " stick: {:12} | new slip: {:12} | slip: {:12} | open: {:12}", + numStick, numNewSlip, numSlip, numOpen ) ); + } ); + } } real64 ContactSolverBase::explicitStep( real64 const & GEOS_UNUSED_PARAM( time_n ), diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.cpp new file mode 100644 index 00000000000..abd3e4ae905 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.cpp @@ -0,0 +1,2294 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +/* + * SolidMechanicsMortarContact.cpp + */ + +#include "mesh/DomainPartition.hpp" +#include "SolidMechanicsMortarContact.hpp" + + +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsConformingContactKernelsBase.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernels.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMKernelsBase.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsALMSimultaneousKernels.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsDisplacementJumpUpdateKernels.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsContactFaceBubbleKernels.hpp" +#include "physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsMortarContactKernels.hpp" +#include "physicsSolvers/solidMechanics/contact/LogLevelsInfo.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" +#include "physicsSolvers/LogLevelsInfo.hpp" +#include "constitutive/contact/FrictionSelector.hpp" +#include "mesh/utilities/ComputationalGeometry.hpp" +#include "denseLinearAlgebra/denseLASolvers.hpp" + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; +using namespace fields; + +SolidMechanicsMortarContact::SolidMechanicsMortarContact( const string & name, + Group * const parent ): + ContactSolverBase( name, parent ) +{ + + m_faceTypeToMortarFiniteElements.get_inserted( ElementShape::Quadrilateral ) = std::make_unique< finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2 >(); + m_faceTypeToMortarFiniteElements.get_inserted( ElementShape::Triangle ) = std::make_unique< finiteElement::H1_TriangleFace_Lagrange1_Gauss4 >(); + + registerWrapper( viewKeyStruct::masterString(), &m_masterName ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the master surface" ); + + registerWrapper( viewKeyStruct::slaveString(), &m_slaveName ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the slave surface" ); + +} + +SolidMechanicsMortarContact::~SolidMechanicsMortarContact() +{} + +void SolidMechanicsMortarContact::registerDataOnMesh( dataRepository::Group & meshBodies ) +{ + GEOS_MARK_FUNCTION; + ContactSolverBase::registerDataOnMesh( meshBodies ); + +} + +void SolidMechanicsMortarContact::registerMortarDataOnMesh( ) +{ + + MeshLevel & mesh = *m_mortarSide.at( MortarSide::Slave ).mesh; + FaceElementSubRegion & subRegion = + m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + + FaceManager & faceManager = mesh.getFaceManager(); + + // Register the total bubble displacement + faceManager.registerField< contact::totalBubbleDisplacement >( this->getName() ). + reference().resizeDimension< 1 >( 3 ); + + // Register the incremental bubble displacement + faceManager.registerField< contact::incrementalBubbleDisplacement >( this->getName() ). + reference().resizeDimension< 1 >( 3 ); + + // Register the rotation matrix + subRegion.registerField< contact::rotationMatrix >( this->getName() ). + reference().resizeDimension< 1, 2 >( 3, 3 ); + + subRegion.registerField< contact::deltaTraction >( getName() ). + reference().resizeDimension< 1 >( 3 ); + + subRegion.registerField< contact::targetIncrementalJump >( getName() ). + reference().resizeDimension< 1 >( 3 ); + +} + + +void SolidMechanicsMortarContact::setupDofs( DomainPartition const & domain, + DofManager & dofManager ) const +{ + GEOS_MARK_FUNCTION; + + SolidMechanicsLagrangianFEM::setupDofs( domain, dofManager ); + + map< std::pair< string, string >, string_array > meshTargets; + + // bubble and tractions are defined on the slave side only + string_array regions; + + MeshLevel & meshSlave = *m_mortarSide.at( MortarSide::Slave ).mesh; + ElementRegionManager const & elementRegionManager = meshSlave.getElemManager(); + elementRegionManager.forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion const & region ) + { + regions.emplace_back( region.getName() ); + } ); + + meshTargets[std::make_pair( m_meshSlaveName, meshSlave.getName())] = std::move( regions ); + + + dofManager.addField( contact::traction::key(), + FieldLocation::Elem, + 3, + meshTargets ); + + dofManager.addField( contact::totalBubbleDisplacement::key(), + FieldLocation::Face, + 3, + meshTargets ); + + + // Add coupling between bubble + dofManager.addCoupling( contact::totalBubbleDisplacement::key(), + contact::totalBubbleDisplacement::key(), + DofManager::Connector::Elem ); + + dofManager.addCoupling( contact::traction::key(), + contact::traction::key(), + DofManager::Connector::Elem ); + + dofManager.addCoupling( solidMechanics::totalDisplacement::key(), + contact::traction::key(), + DofManager::Connector::Elem, + meshTargets ); + + dofManager.addCoupling( contact::totalBubbleDisplacement::key(), + contact::traction::key(), + DofManager::Connector::Elem, + meshTargets ); + +} + + +void SolidMechanicsMortarContact::setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity ) +{ + + ElementShape shapes[2] = { ElementShape::Triangle, ElementShape::Quadrilateral }; + + setMortarSurfaces( domain ); + + registerMortarDataOnMesh(); + + createFaceTypeListMortar( MortarSide::Master ); + createFaceTypeListMortar( MortarSide::Slave ); + + // create list of bubbles for the slave side + this->createBubbleCellList(); + + // perform rough screen to find potential connections between master and slave faces + connectivityMapType connectivityMap; + getConnectivityMap( connectivityMap ); + + // preprocess mortar integration data + computeMortarInterpolation( connectivityMap ); + + // setup dofs for the problem + dofManager.setDomain( domain ); + setupDofs( domain, dofManager ); + dofManager.reorderByRank(); + + // Set the sparsity pattern without the Abu and Aub blocks. + SparsityPattern< globalIndex > patternDiag; + dofManager.setSparsityPattern( patternDiag ); + + // Manually define the sparsity pattern induced by the slave traction <-> master displacements coupling + // Get the original row lengths (diagonal blocks only) + array1d< localIndex > rowLengths( patternDiag.numRows() ); + for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) + { + rowLengths[localRow] = patternDiag.numNonZeros( localRow ); + } + + // Add the number of nonzeros induced by coupling + this->addBubbleCouplingNumNonzeros( dofManager, rowLengths.toView() ); + + for( auto slaveShape : shapes ) + { + for( auto masterShape : shapes ) + { + this->addMortarCouplingNumNonzeros( dofManager, + slaveShape, + masterShape, + connectivityMap, + rowLengths.toView()); + } + } + + + // Create a new pattern with enough capacity for coupled matrix + SparsityPattern< globalIndex > pattern; + pattern.resizeFromRowCapacities< parallelHostPolicy >( patternDiag.numRows(), patternDiag.numColumns(), rowLengths.data() ); + + // Copy the original nonzeros + for( localIndex localRow = 0; localRow < patternDiag.numRows(); ++localRow ) + { + globalIndex const * cols = patternDiag.getColumns( localRow ).dataIfContiguous(); + pattern.insertNonZeros( localRow, cols, cols + patternDiag.numNonZeros( localRow ) ); + } + + // Add the nonzeros from coupling + this->addBubbleCouplingSparsityPattern( dofManager, pattern.toView() ); + + for( auto slaveShape : shapes ) + { + for( auto masterShape : shapes ) + { + this->addMortarCouplingSparsityPattern( dofManager, + slaveShape, + masterShape, + pattern.toView() ); + } + } + + // // Finally, steal the pattern into a CRS matrix + localMatrix.assimilate< parallelDevicePolicy<> >( std::move( pattern ) ); + localMatrix.setName( this->getName() + "/localMatrix" ); + + rhs.setName( this->getName() + "/rhs" ); + rhs.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + + solution.setName( this->getName() + "/solution" ); + solution.create( dofManager.numLocalDofs(), MPI_COMM_GEOS ); + + computeRotationMatrices( ); + + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( setSparsity ); + +} + +void SolidMechanicsMortarContact::implicitStepSetup( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) +{ + + GEOS_MARK_FUNCTION; + SolidMechanicsLagrangianFEM::implicitStepSetup( time_n, dt, domain ); + + + +} + +void SolidMechanicsMortarContact::assembleSystem( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + + GEOS_MARK_FUNCTION; + synchronizeFractureState( domain ); + + SolidMechanicsLagrangianFEM::assembleSystem( time, + dt, + domain, + dofManager, + localMatrix, + localRhs ); + + assembleBubbles( dt, domain, dofManager, localMatrix, localRhs ); + + assembleMortar( dt, dofManager, localMatrix, localRhs ); + +} + + +void SolidMechanicsMortarContact::assembleBubbles( real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + + // Loop for assembling contributes of bubble elements (Abb, Abu, Aub) + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, + MeshLevel & mesh, + string_array const & regionNames ) + { + + if( meshName != m_meshSlaveName ) + { + return; + } + else + { + NodeManager const & nodeManager = mesh.getNodeManager(); + FaceManager const & faceManager = mesh.getFaceManager(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & bubbleDofKey = dofManager.getKey( contact::totalBubbleDisplacement::key() ); + + arrayView1d< globalIndex const > const dispDofNumber = nodeManager.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + + real64 const gravityVectorData[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( gravityVector() ); + + + solidMechanicsConformingContactKernels::FaceBubbleFactory kernelFactory( dispDofNumber, + bubbleDofNumber, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + gravityVectorData ); + + real64 maxTraction = finiteElement:: + regionBasedKernelApplication + < parallelDevicePolicy< >, + ElasticIsotropic, + CellElementSubRegion >( mesh, + regionNames, + getDiscretizationName(), + SolidMechanicsLagrangianFEM::viewKeyStruct::solidMaterialNamesString(), + kernelFactory ); + + GEOS_UNUSED_VAR( maxTraction ); + + + } + + } ); + + // assembling interface bubble contribution: Abt, Atb, inactive Abb + assembleMortarBubbles< ElementShape::Triangle >( dofManager, localMatrix, localRhs ); + assembleMortarBubbles< ElementShape::Quadrilateral >( dofManager, localMatrix, localRhs ); + +} + + +template< ElementShape shape > +void SolidMechanicsMortarContact::assembleMortarBubbles( DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + + arrayView1d< localIndex const > elementList = m_faceTypeToElementList.at( MortarSide::Slave ).at( shape ); + + if( elementList.size() == 0 ) + return; + + using FEType = typename selectFE< shape >::type; + localIndex constexpr numQuadraturePoints = FEType::numQuadraturePoints; + localIndex constexpr numNodesPerElem = FEType::numNodes; + + MeshLevel & mesh = *m_mortarSide.at( MortarSide::Slave ).mesh; + NodeManager const & nodeManager = mesh.getNodeManager(); + FaceManager const & faceManager = mesh.getFaceManager(); + FaceElementSubRegion const & surfRegion = m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + arrayView2d< localIndex const > const elemsToFaces = surfRegion.faceList().toViewConst(); + + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const coords = nodeManager.referencePosition(); + ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); + + + globalIndex const rankOffset = dofManager.rankOffset(); + string const bubbleDofKey = dofManager.getKey( contact::totalBubbleDisplacement::key() ); + string const tractionDofKey = dofManager.getKey( contact::traction::key() ); + + arrayView2d< real64 const > const traction = surfRegion.getField< fields::contact::traction >().toViewConst(); + arrayView2d< real64 const > const bubbleDisplacement = faceManager.getField< fields::contact::totalBubbleDisplacement >().toViewConst(); + arrayView3d< real64 const > const rotationMatrix = surfRegion.getField< fields::contact::rotationMatrix >().toViewConst(); + + arrayView1d< globalIndex const > const bubbleDofNumber = faceManager.getReference< globalIndex_array >( bubbleDofKey ); + arrayView1d< globalIndex const > const tractionDofNumber = surfRegion.getReference< globalIndex_array >( tractionDofKey ); + + int permutation[numNodesPerElem]; + FEType::getPermutation( permutation ); + + forAll< parallelDevicePolicy<> >( elementList.size(), [=] GEOS_HOST_DEVICE ( localIndex const kk ) + { + localIndex const k = elementList[kk]; + localIndex const kf0 = elemsToFaces[k][0]; + localIndex const kf1 = elemsToFaces[k][1]; + + globalIndex bDof1[3], bDof2[3], tDof[3]; + real64 bLoc[3], tLoc[3]; + for( localIndex i=0; i<3; ++i ) + { + bDof1[i] = bubbleDofNumber[kf0]+i; + bDof2[i] = bubbleDofNumber[kf1]+i; + tDof[i] = tractionDofNumber[k]+i; + bLoc[i] = bubbleDisplacement[kf0][i]; + tLoc[i] = traction[k][i]; + } + + real64 X[numNodesPerElem][3]; + + for( localIndex a=0; a( RtAtb, R, localAtb ); + LvArray::tensorOps::copy< 3, 3 >( localAtb, RtAtb ); + LvArray::tensorOps::transpose< 3, 3 >( localAbt, localAtb ); + LvArray::tensorOps::Ri_eq_AijBj< 3, 3 >( rhsB, localAbt, tLoc ); + LvArray::tensorOps::Ri_eq_AijBj< 3, 3 >( rhsT, localAtb, bLoc ); + + // assemble matrix and rhs contribution of Abt + for( localIndex i = 0; i < 3; ++i ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( bDof1[i] - rankOffset ); + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRow< parallelHostAtomic >( localRow, + tDof, + localAbt[i], + 3 ); + RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], rhsB[i] ); + } + } + + // assemble matrix and rhs contribution of Atb + for( localIndex i = 0; i < 3; ++i ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( tDof[i] - rankOffset ); + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + localMatrix.addToRow< parallelHostAtomic >( localRow, + bDof1, + localAtb[i], + 3 ); + RAJA::atomicAdd( parallelHostAtomic{}, &localRhs[localRow], rhsT[i] ); + } + } + + real64 diagVal[1]; + globalIndex colDof[1]; + diagVal[0] = 1.0; + + // assemble diagonal contribution of inactive bubbles (dirichlet bcs) + for( localIndex i = 0; i < 3; ++i ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( bDof2[i] - rankOffset ); + if( localRow >= 0 && localRow < localMatrix.numRows() ) + { + colDof[0] = bDof2[i]; + localMatrix.addToRow< parallelHostAtomic >( localRow, + colDof, + diagVal, + 1 ); + } + } + + + } ); + +} + + +void SolidMechanicsMortarContact::assembleMortar( real64 const dt, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + // Call the templated version of assembleMortar with the appropriate element shapes + assembleMortar< ElementShape::Triangle, ElementShape::Triangle >( dt, dofManager, localMatrix, localRhs ); + assembleMortar< ElementShape::Triangle, ElementShape::Quadrilateral >( dt, dofManager, localMatrix, localRhs ); + assembleMortar< ElementShape::Quadrilateral, ElementShape::Triangle >( dt, dofManager, localMatrix, localRhs ); + assembleMortar< ElementShape::Quadrilateral, ElementShape::Quadrilateral >( dt, dofManager, localMatrix, localRhs ); +} + +template< ElementShape shpS, ElementShape shpM > +void SolidMechanicsMortarContact::assembleMortar( real64 const dt, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + + if( m_triCellsDet.find( std::make_pair( shpS, shpM )) == m_triCellsDet.end()) + { + // no connections for the given element shapes + return; + } + + MeshLevel & meshSlave = *m_mortarSide.at( MortarSide::Slave ).mesh; + FaceElementSubRegion & subRegionSlave = + m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + NodeManager const & nodeManagerSlave = meshSlave.getNodeManager(); + //FaceManager const & faceManagerSlave = meshSlave.getFaceManager(); + + MeshLevel & meshMaster = *m_mortarSide.at( MortarSide::Master ).mesh; + FaceElementSubRegion const & subRegionMaster = + m_mortarSide.at( MortarSide::Master ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + NodeManager const & nodeManagerMaster = meshMaster.getNodeManager(); + //FaceManager const & faceManagerMaster = meshMaster.getFaceManager(); + + string const & dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const & tractionDofKey = dofManager.getKey( contact::traction::key() ); + + arrayView1d< globalIndex const > const dispDofNumberSlave = nodeManagerSlave.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const dispDofNumberMaster = nodeManagerMaster.getReference< globalIndex_array >( dispDofKey ); + + finiteElement::FiniteElementBase const & feSlave = *(m_faceTypeToMortarFiniteElements.at( shpS )); + finiteElement::FiniteElementBase const & feMaster = *(m_faceTypeToMortarFiniteElements.at( shpM )); + + arrayView1d< real64 const > detPlus = m_triCellsDet.at( std::make_pair( shpS, shpM )).toView(); + localIndex const nTri = detPlus.size(); + array1d< real64 > detMinus( nTri ); + + for( localIndex i=0; i listSlave = m_triCells.at( MortarSide::Slave ).at( std::make_pair( shpS, shpM )).toView(); + arrayView1d< localIndex const > listMaster = m_triCells.at( MortarSide::Master ).at( std::make_pair( shpS, shpM )).toView(); + arrayView3d< real64 const > gpLocalCoordsSlave = m_gpLocalCoords.at( MortarSide::Slave ).at( std::make_pair( shpS, shpM )).toView(); + arrayView3d< real64 const > gpLocalCoordsMaster = m_gpLocalCoords.at( MortarSide::Master ).at( std::make_pair( shpS, shpM )).toView(); + + solidMechanicsMortarContactKernels::MortarFactory + kernelFactorySlave( dispDofNumberSlave, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + subRegionSlave, + subRegionSlave, + listSlave, + listSlave, + detPlus, + gpLocalCoordsSlave, + tractionDofKey ); + + real64 maxTraction1 = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy< >, + FrictionBase >( meshSlave, + m_slaveName, + listSlave, + feSlave, + viewKeyStruct::frictionLawNameString(), + kernelFactorySlave ); + + + solidMechanicsMortarContactKernels::MortarFactory + kernelFactoryMaster( dispDofNumberMaster, + dofManager.rankOffset(), + localMatrix, + localRhs, + dt, + subRegionSlave, + subRegionMaster, + listSlave, + listMaster, + detMinus.toView(), + gpLocalCoordsMaster, + tractionDofKey ); + + real64 maxTraction2 = finiteElement:: + interfaceBasedKernelApplication + < parallelDevicePolicy< >, + FrictionBase >( meshMaster, + m_masterName, + listSlave, + feMaster, + viewKeyStruct::frictionLawNameString(), + kernelFactoryMaster ); + + GEOS_UNUSED_VAR( maxTraction1 ); + GEOS_UNUSED_VAR( maxTraction2 ); + +} + + +void SolidMechanicsMortarContact::implicitStepComplete( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) +{ + + SolidMechanicsLagrangianFEM::implicitStepComplete( time_n, dt, domain ); + + FaceElementSubRegion & subRegion = + m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView2d< real64 > const deltaTraction = subRegion.getField< contact::deltaTraction >(); + // arrayView2d< real64 > const deltaDispJump = subRegion.getField< contact::deltaDispJump >(); + // arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + // arrayView2d< real64 > const oldDispJump = subRegion.getField< contact::oldDispJump >() + + forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const kfe ) + { + //LvArray::tensorOps::fill< 3 >( deltaDispJump[kfe], 0.0 ); + LvArray::tensorOps::fill< 3 >( deltaTraction[kfe], 0.0 ); + //LvArray::tensorOps::copy< 3 >( oldDispJump[kfe], dispJump[kfe] ); + } ); + +} + +real64 SolidMechanicsMortarContact::calculateResidualNorm( real64 const & time, + real64 const & dt, + DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localRhs ) +{ + GEOS_MARK_FUNCTION; + + real64 const solidResidual = SolidMechanicsLagrangianFEM::calculateResidualNorm( time, dt, domain, dofManager, localRhs ); + + real64 const contactResidual = calculateContactResidualNorm( domain, dofManager, localRhs ); + + return sqrt( solidResidual * solidResidual + contactResidual * contactResidual ); +} + +real64 SolidMechanicsMortarContact::calculateContactResidualNorm( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localRhs ) +{ + + GEOS_UNUSED_VAR( domain ); + + string const & dofKey = dofManager.getKey( contact::traction::key() ); + globalIndex const rankOffset = dofManager.rankOffset(); + + real64 stickResidual = 0.0; + + FaceElementSubRegion const & subRegion = + m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + + arrayView1d< globalIndex const > const & dofNumber = subRegion.getReference< array1d< globalIndex > >( dofKey ); + arrayView1d< integer const > const & ghostRank = subRegion.ghostRank(); + arrayView1d< real64 const > const & area = subRegion.getElementArea(); + RAJA::ReduceSum< parallelHostReduce, real64 > stickSum( 0.0 ); + forAll< parallelHostPolicy >( subRegion.size(), [=] ( localIndex const k ) + { + if( ghostRank[k] < 0 ) + { + localIndex const localRow = LvArray::integerConversion< localIndex >( dofNumber[k] - rankOffset ); + for( localIndex dim = 0; dim < 3; ++dim ) + { + real64 const norm = localRhs[localRow + dim] / area[k]; + stickSum += norm * norm; + } + } + } ); + stickResidual += stickSum.get(); + + stickResidual = MpiWrapper::sum( stickResidual ); + stickResidual = sqrt( stickResidual ); + + if( getLogLevel() >= 1 && logger::internal::rank==0 ) + { + std::cout << GEOS_FMT( " ( Rt ) = ( {:15.6e} )", stickResidual ); + } + + return sqrt( stickResidual * stickResidual ); +} + +void SolidMechanicsMortarContact::applySystemSolution( DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution, + real64 const scalingFactor, + real64 const dt, + DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + + SolidMechanicsLagrangianFEM::applySystemSolution( dofManager, localSolution, scalingFactor, dt, domain ); + + dofManager.addVectorToField( localSolution, + contact::traction::key(), + contact::deltaTraction::key(), + scalingFactor ); + + dofManager.addVectorToField( localSolution, + contact::traction::key(), + contact::traction::key(), + scalingFactor ); + + dofManager.addVectorToField( localSolution, + contact::totalBubbleDisplacement::key(), + contact::totalBubbleDisplacement::key(), + scalingFactor ); + + dofManager.addVectorToField( localSolution, + contact::totalBubbleDisplacement::key(), + contact::incrementalBubbleDisplacement::key(), + scalingFactor ); + + MeshLevel & mesh = *m_mortarSide.at( MortarSide::Slave ).mesh; + + FieldIdentifiers fieldsToBeSync; + + fieldsToBeSync.addFields( FieldLocation::Face, + { contact::incrementalBubbleDisplacement::key(), + contact::totalBubbleDisplacement::key() } ); + fieldsToBeSync.addElementFields( { contact::traction::key(), + contact::deltaTraction::key(), + contact::dispJump::key() }, + { getUniqueFractureRegionName() } ); + + CommunicationTools::getInstance().synchronizeFields( fieldsToBeSync, + mesh, + domain.getNeighbors(), + true ); + +} + + + +void SolidMechanicsMortarContact::addBubbleCouplingNumNonzeros( DofManager & dofManager, + arrayView1d< localIndex > const & rowLengths ) const +{ + + MeshLevel & mesh = *m_mortarSide.at( MortarSide::Slave ).mesh; + + ElementRegionManager const & elemManagerSlave = mesh.getElemManager(); + NodeManager const & nodeManagerSlave = mesh.getNodeManager(); + FaceManager const & faceManagerSlave = mesh.getFaceManager(); + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const bubbleDofKey = dofManager.getKey( contact::totalBubbleDisplacement::key() ); + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + + arrayView1d< globalIndex const > const bubbleDofNumber = faceManagerSlave.getReference< globalIndex_array >( bubbleDofKey ); + arrayView1d< globalIndex const > const dispSlaveDofNumber = nodeManagerSlave.getReference< globalIndex_array >( dispDofKey ); + + // add coupling between bubble and displacement in the slave side + elemManagerSlave.forElementSubRegions< CellElementSubRegion >( [&]( const CellElementSubRegion & cellElementSubRegion ) + { + + arrayView1d< localIndex const > const bubbleElemsList = cellElementSubRegion.bubbleElementsList(); + arrayView2d< localIndex const > const faceElemsList = cellElementSubRegion.faceElementsList(); + + localIndex const numDispDof = 3*cellElementSubRegion.numNodesPerElement(); + + for( localIndex bi=0; bi( bubbleDofNumber[k] - rankOffset ); + + if( localRow >= 0 && localRow < rowLengths.size() ) + { + for( localIndex i=0; i<3; ++i ) + { + rowLengths[localRow + i] += numDispDof; + } + } + + for( localIndex a=0; a( dispSlaveDofNumber[node] - rankOffset ); + + if( localDispRow >= 0 && localDispRow < rowLengths.size() ) + { + for( int d=0; d<3; ++d ) + { + rowLengths[localDispRow + d] += 3; + } + } + } + } + + } ); + +} + + +void SolidMechanicsMortarContact::addMortarCouplingNumNonzeros( DofManager & dofManager, + ElementShape const & slaveShape, + ElementShape const & masterShape, + connectivityMapType const & connectivityMap, + arrayView1d< localIndex > const & rowLengths ) const +{ + + FaceElementSubRegion const & surfRegionSlave = + m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + + MeshLevel & meshMaster = *m_mortarSide.at( MortarSide::Master ).mesh; + FaceElementSubRegion const & surfRegionMaster = + m_mortarSide.at( MortarSide::Master ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + NodeManager const & nodeManagerMaster = meshMaster.getNodeManager(); + FaceManager const & faceManagerMaster = meshMaster.getFaceManager(); + ArrayOfArraysView< localIndex const > const faceToNodeMapMaster = faceManagerMaster.nodeList().toViewConst(); + arrayView2d< localIndex const > const elemsToFacesMaster = surfRegionMaster.faceList().toViewConst(); + + ArrayOfArraysView< localIndex const > const connections = connectivityMap.at( {slaveShape, masterShape} ).toViewConst(); + arrayView1d< localIndex const > const masterElementList = m_faceTypeToElementList.at( MortarSide::Master ).at( masterShape ).toViewConst(); + arrayView1d< localIndex const > const slaveElementList = m_faceTypeToElementList.at( MortarSide::Slave ).at( slaveShape ).toViewConst(); + + if( slaveElementList.size() == 0 || masterElementList.size() == 0 ) + { + return; + } + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const tractionDofKey = dofManager.getKey( contact::traction::key() ); + arrayView1d< globalIndex const > const dispMasterDofNumber = nodeManagerMaster.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const tractionDofNumber = surfRegionSlave.getReference< globalIndex_array >( tractionDofKey ); + + + // add coupling between tractions and master displacements + for( localIndex im=0; im( dispMasterDofNumber[node] - rankOffset ); + + if( localDispRow >= 0 && localDispRow < rowLengths.size() ) + { + for( int d=0; d<3; ++d ) + { + rowLengths[localDispRow + d] += 3 * nSlave; + } + } + } + + // loop over connected slave elements + for( localIndex is=0; is( tractionDofNumber[kSlave] - rankOffset ); + if( localRow >= 0 && localRow < rowLengths.size() ) + { + for( localIndex i=0; i<3; ++i ) + { + rowLengths[localRow + i] += numDispDof; + } + } + } + + } +} + + +void SolidMechanicsMortarContact::addBubbleCouplingSparsityPattern( DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const +{ + + MeshLevel & mesh = *m_mortarSide.at( MortarSide::Slave ).mesh; + ElementRegionManager const & elemManagerSlave = mesh.getElemManager(); + NodeManager const & nodeManagerSlave = mesh.getNodeManager(); + FaceManager const & faceManagerSlave = mesh.getFaceManager(); + + //FaceElementSubRegion const & surfRegionSlave = m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion + // >(); + //arrayView2d< localIndex const > const elemsToFace = surfRegionSlave.faceList().toViewConst(); + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const bubbleDofKey = dofManager.getKey( contact::totalBubbleDisplacement::key() ); + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + arrayView1d< globalIndex const > const bubbleDofNumber = faceManagerSlave.getReference< globalIndex_array >( bubbleDofKey ); + arrayView1d< globalIndex const > const dispDofNumber = nodeManagerSlave.getReference< globalIndex_array >( dispDofKey ); + + static constexpr int maxNumDispDof = 3 * 8; + + elemManagerSlave.forElementSubRegions< CellElementSubRegion >( [&]( const CellElementSubRegion & cellElementSubRegion ) + { + arrayView1d< localIndex const > const bubbleElemsList = cellElementSubRegion.bubbleElementsList(); + arrayView2d< localIndex const > const faceElemsList = cellElementSubRegion.faceElementsList(); + localIndex const numDispDof = 3*cellElementSubRegion.numNodesPerElement(); + for( localIndex bi=0; bi eqnRowIndicesDisp ( numDispDof ); + stackArray1d< globalIndex, 3 > eqnRowIndicesBubble( 3 ); + stackArray1d< globalIndex, maxNumDispDof > dofColIndicesDisp ( numDispDof ); + stackArray1d< globalIndex, 3 > dofColIndicesBubble( 3 ); + for( localIndex idof = 0; idof < 3; ++idof ) + { + eqnRowIndicesBubble[idof] = bubbleDofNumber[k] + idof - rankOffset; + dofColIndicesBubble[idof] = bubbleDofNumber[k] + idof; + } + for( localIndex a=0; a= 0 && eqnRowIndicesDisp[i] < pattern.numRows() ) + { + for( localIndex j = 0; j < dofColIndicesBubble.size(); ++j ) + { + pattern.insertNonZero( eqnRowIndicesDisp[i], dofColIndicesBubble[j] ); + } + } + } + for( localIndex i = 0; i < eqnRowIndicesBubble.size(); ++i ) + { + if( eqnRowIndicesBubble[i] >= 0 && eqnRowIndicesBubble[i] < pattern.numRows() ) + { + for( localIndex j=0; j < dofColIndicesDisp.size(); ++j ) + { + pattern.insertNonZero( eqnRowIndicesBubble[i], dofColIndicesDisp[j] ); + } + } + } + } + } ); + +} + + +void SolidMechanicsMortarContact::addMortarCouplingSparsityPattern( DofManager & dofManager, + ElementShape const & slaveShape, + ElementShape const & masterShape, + SparsityPatternView< globalIndex > const & pattern ) const +{ + + FaceElementSubRegion const & surfRegionSlave = m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + + MeshLevel & meshMaster = *m_mortarSide.at( MortarSide::Master ).mesh; + FaceElementSubRegion const & surfRegionMaster = m_mortarSide.at( MortarSide::Master ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + NodeManager const & nodeManagerMaster = meshMaster.getNodeManager(); + FaceManager const & faceManagerMaster = meshMaster.getFaceManager(); + ArrayOfArraysView< localIndex const > const faceToNodeMapMaster = faceManagerMaster.nodeList().toViewConst(); + arrayView2d< localIndex const > const elemsToFacesMaster = surfRegionMaster.faceList().toViewConst(); + + //ArrayOfArraysView< localIndex const > const connections = connectivityMap.at({slaveShape, masterShape}).toViewConst(); + //arrayView1d< localIndex const > const masterElementList = m_faceTypeToElementList.at(MortarSide::Master).at(masterShape).toViewConst(); + //arrayView1d< localIndex const > const slaveElementList = m_faceTypeToElementList.at(MortarSide::Slave).at(slaveShape).toViewConst(); + + if( m_triCellsDet.find( std::make_pair( slaveShape, masterShape )) == m_triCellsDet.end()) + { + // no connections for the given element shapes + return; + } + + arrayView1d< localIndex const > const masterElementList = + m_triCells.at( MortarSide::Master ).at( std::make_pair( slaveShape, masterShape )).toViewConst(); + + arrayView1d< localIndex const > const slaveElementList = + m_triCells.at( MortarSide::Slave ).at( std::make_pair( slaveShape, masterShape )).toViewConst(); + + + globalIndex const rankOffset = dofManager.rankOffset(); + + string const dispDofKey = dofManager.getKey( solidMechanics::totalDisplacement::key() ); + string const tractionDofKey = dofManager.getKey( contact::traction::key() ); + arrayView1d< globalIndex const > const dispMasterDofNumber = nodeManagerMaster.getReference< globalIndex_array >( dispDofKey ); + arrayView1d< globalIndex const > const tractionDofNumber = surfRegionSlave.getReference< globalIndex_array >( tractionDofKey ); + + // populate sparsity pattern for coupling between tractions and master elements + static constexpr int maxNumDispFaceDof = 3 * 4; + + for( localIndex k=0; k eqnRowIndicesDisp ( numDispDof ); + stackArray1d< globalIndex, 3 > eqnRowIndicesTraction( 3 ); + stackArray1d< globalIndex, maxNumDispFaceDof > dofColIndicesDisp( numDispDof ); + stackArray1d< globalIndex, 3 > dofColIndicesTraction( 3 ); + for( localIndex idof = 0; idof < 3; ++idof ) + { + eqnRowIndicesTraction[idof] = tractionDofNumber[kElSlave] + idof - rankOffset; + dofColIndicesTraction[idof] = tractionDofNumber[kElSlave] + idof; + } + for( localIndex a=0; a= 0 && eqnRowIndicesDisp[i] < pattern.numRows() ) + { + for( localIndex j = 0; j < dofColIndicesTraction.size(); ++j ) + { + pattern.insertNonZero( eqnRowIndicesDisp[i], dofColIndicesTraction[j] ); + } + } + } + for( localIndex i = 0; i < eqnRowIndicesTraction.size(); ++i ) + { + if( eqnRowIndicesTraction[i] >= 0 && eqnRowIndicesTraction[i] < pattern.numRows() ) + { + for( localIndex j=0; j < dofColIndicesDisp.size(); ++j ) + { + pattern.insertNonZero( eqnRowIndicesTraction[i], dofColIndicesDisp[j] ); + } + } + } + } + + +} + + +void SolidMechanicsMortarContact::computeRotationMatrices( ) +{ + + MeshLevel & mesh = *m_mortarSide.at( MortarSide::Slave ).mesh; + FaceElementSubRegion & subRegion = + m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + + FaceManager & faceManager = mesh.getFaceManager(); + //ElementRegionManager & elemManager = mesh.getElemManager(); + + arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); + + // remember, view modify the original object + arrayView3d< real64 > const rotationMatrix = + subRegion.getField< contact::rotationMatrix >().toView(); + + arrayView2d< real64 > const unitNormal = subRegion.getNormalVector(); + arrayView2d< real64 > const unitTangent1 = subRegion.getTangentVector1(); + arrayView2d< real64 > const unitTangent2 = subRegion.getTangentVector2(); + + // Compute rotation matrices + solidMechanicsConformingContactKernels::ComputeRotationMatricesKernel::launch< parallelDevicePolicy<> >( subRegion.size(), + faceNormal, + elemsToFaces, + rotationMatrix, + unitNormal, + unitTangent1, + unitTangent2 ); +} + + + +void SolidMechanicsMortarContact::updateState( DomainPartition & domain ) +{ + GEOS_UNUSED_VAR( domain ); +} + + + +void SolidMechanicsMortarContact::createFaceTypeListMortar( MortarSide side ) +{ + // Generate lists containing elements of various face types + MeshLevel & mesh = *m_mortarSide.at( side ).mesh; + FaceElementSubRegion const & surfRegion = m_mortarSide.at( side ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + + FaceManager const & faceManager = mesh.getFaceManager(); + //ElementRegionManager const & elemManager = mesh.getElemManager(); + ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst(); + + array1d< localIndex > keys( surfRegion.size()); + array1d< localIndex > vals( surfRegion.size()); + array1d< localIndex > quadList; + array1d< localIndex > triList; + RAJA::ReduceSum< ReducePolicy< parallelDevicePolicy<> >, localIndex > nTri_r( 0 ); + RAJA::ReduceSum< ReducePolicy< parallelDevicePolicy<> >, localIndex > nQuad_r( 0 ); + + arrayView1d< localIndex > const keys_v = keys.toView(); + arrayView1d< localIndex > const vals_v = vals.toView(); + // Determine the size of the lists and generate the vector keys and vals for parallel indexing into lists. + // (With RAJA, parallelizing this operation seems the most viable approach.) + forAll< parallelDevicePolicy<> >( surfRegion.size(), + [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) + { + + localIndex const numNodesPerFace = faceToNodeMap.sizeOfArray( kfe ); + if( numNodesPerFace == 3 ) + { + keys_v[kfe]=0; + vals_v[kfe]=kfe; + nTri_r += 1; + GEOS_ERROR( "SolidMechanicsMortarContact:: triangular face type not yet available" ); + } + else if( numNodesPerFace == 4 ) + { + keys_v[kfe]=1; + vals_v[kfe]=kfe; + nQuad_r += 1; + } + else + { + GEOS_ERROR( "SolidMechanicsMortarContact:: invalid face type" ); + } + } ); + + localIndex nQuad = static_cast< localIndex >(nQuad_r.get()); + localIndex nTri = static_cast< localIndex >(nTri_r.get()); + + // Sort vals according to keys to ensure that + // elements of the same type are adjacent in the vals list. + // This arrangement allows for efficient copying into the container + // by leveraging parallelism. + RAJA::sort_pairs< parallelDevicePolicy<> >( keys_v, vals_v ); + + quadList.resize( nQuad ); + triList.resize( nTri ); + arrayView1d< localIndex > const quadList_v = quadList.toView(); + arrayView1d< localIndex > const triList_v = triList.toView(); + + forAll< parallelDevicePolicy<> >( nTri, [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) + { + triList_v[kfe] = vals_v[kfe]; + } ); + + forAll< parallelDevicePolicy<> >( nQuad, [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) + { + quadList_v[kfe] = vals_v[nTri+kfe]; + } ); + + stdMap< string, array1d< localIndex > > faceTypeList; + + m_faceTypeToElementList.get_inserted( side ).get_inserted( ElementShape::Quadrilateral ) = quadList; + m_faceTypeToElementList.get_inserted( side ).get_inserted( ElementShape::Triangle ) = triList; + +} + + +void SolidMechanicsMortarContact::createBubbleCellList( ) const +{ + MeshLevel & mesh = *m_mortarSide.at( MortarSide::Slave ).mesh; + FaceElementSubRegion const & surfRegion = m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + + ElementRegionManager & elemManager = mesh.getElemManager(); + + // Array to store face indexes + array1d< localIndex > tmpSpace( 2*surfRegion.size()); + SortedArray< localIndex > faceIdList; + + arrayView1d< localIndex > const tmpSpace_v = tmpSpace.toView(); + // Store indexes of faces in the temporany array. + { + arrayView2d< localIndex const > const elemsToFaces = surfRegion.faceList().toViewConst(); + + forAll< parallelDevicePolicy<> >( surfRegion.size(), [ = ] GEOS_HOST_DEVICE ( localIndex const kfe ) + { + + localIndex const kf0 = elemsToFaces[kfe][0], kf1 = elemsToFaces[kfe][1]; + tmpSpace_v[2*kfe] = kf0, tmpSpace_v[2*kfe+1] = kf1; + + } ); + } + + // Sort indexes to enable efficient searching using binary search. + RAJA::stable_sort< parallelDevicePolicy<> >( tmpSpace_v ); + faceIdList.insert( tmpSpace_v.begin(), tmpSpace_v.end()); + + // Search for bubble element on each CellElementSubRegion and + // store element indexes, global and local face indexes. + elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & cellElementSubRegion ) + { + + arrayView2d< localIndex const > const elemsToFaces = cellElementSubRegion.faceList().toViewConst(); + + RAJA::ReduceSum< ReducePolicy< parallelDevicePolicy<> >, localIndex > nBubElems_r( 0 ); + + localIndex const n_max = cellElementSubRegion.size() * elemsToFaces.size( 1 ); + array1d< localIndex > keys( n_max ); + array1d< localIndex > perms( n_max ); + array1d< localIndex > vals( n_max ); + array1d< localIndex > localFaceIds( n_max ); + + arrayView1d< localIndex > const keys_v = keys.toView(); + arrayView1d< localIndex > const perms_v = perms.toView(); + arrayView1d< localIndex > const vals_v = vals.toView(); + arrayView1d< localIndex > const localFaceIds_v = localFaceIds.toView(); + SortedArrayView< localIndex const > const faceIdList_v = faceIdList.toViewConst(); + + forAll< parallelDevicePolicy<> >( cellElementSubRegion.size(), + [ = ] + GEOS_HOST_DEVICE ( localIndex const kfe ) + { + for( int i=0; i < elemsToFaces.size( 1 ); ++i ) + { + perms_v[kfe*elemsToFaces.size( 1 )+i] = kfe*elemsToFaces.size( 1 )+i; + if( faceIdList_v.contains( elemsToFaces[kfe][i] )) + { + keys_v[kfe*elemsToFaces.size( 1 )+i] = 0; + vals_v[kfe*elemsToFaces.size( 1 )+i] = kfe; + localFaceIds_v[kfe*elemsToFaces.size( 1 )+i] = i; + nBubElems_r += 1; + } + else + { + keys_v[kfe*elemsToFaces.size( 1 )+i] = 1; + vals_v[kfe*elemsToFaces.size( 1 )+i] = -1; + localFaceIds_v[kfe*elemsToFaces.size( 1 )+i] = -1; + } + } + } ); + + // Sort perms according to keys to ensure that bubble elements are adjacent + // and occupy the first positions of the list. + // This arrangement allows for efficient copying into the container + // by leveraging parallelism. + localIndex nBubElems = static_cast< localIndex >(nBubElems_r.get()); + RAJA::sort_pairs< parallelDevicePolicy<> >( keys_v, perms_v ); + + array1d< localIndex > bubbleElemsList; + bubbleElemsList.resize( nBubElems ); + + arrayView1d< localIndex > const bubbleElemsList_v = bubbleElemsList.toView(); + + forAll< parallelDevicePolicy<> >( n_max, [ = ] GEOS_HOST_DEVICE ( localIndex const k ) + { + keys_v[k] = vals_v[perms_v[k]]; + } ); + + forAll< parallelDevicePolicy<> >( nBubElems, [ = ] GEOS_HOST_DEVICE ( localIndex const k ) + { + bubbleElemsList_v[k] = keys_v[k]; + } ); + cellElementSubRegion.setBubbleElementsList( bubbleElemsList.toViewConst()); + + forAll< parallelDevicePolicy<> >( n_max, [ = ] GEOS_HOST_DEVICE ( localIndex const k ) + { + keys_v[k] = localFaceIds_v[perms_v[k]]; + } ); + + array2d< localIndex > faceElemsList; + faceElemsList.resize( nBubElems, 2 ); + + arrayView2d< localIndex > const faceElemsList_v = faceElemsList.toView(); + + forAll< parallelDevicePolicy<> >( nBubElems, + [ = ] + GEOS_HOST_DEVICE ( localIndex const k ) + { + localIndex const kfe = bubbleElemsList_v[k]; + faceElemsList_v[k][0] = elemsToFaces[kfe][keys_v[k]]; + faceElemsList_v[k][1] = keys_v[k]; + } ); + cellElementSubRegion.setFaceElementsList( faceElemsList.toViewConst()); + + } ); + +} + +void SolidMechanicsMortarContact::setMortarSurfaces( DomainPartition & domain ) +{ + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, + MeshLevel & mesh, + string_array const & ) + { + + GEOS_UNUSED_VAR( meshName ); + ElementRegionManager & elemManager = mesh.getElemManager(); + elemManager.forElementRegions< SurfaceElementRegion >( [&, this]( SurfaceElementRegion & region ) + { + + string surfacePath = region.getPath(); + // check if region is master or slave and populate member maps + if( surfacePath.find( m_slaveName ) != std::string::npos ) + { + MortarSurface surfaceSlave; + surfaceSlave.mesh = &mesh; + surfaceSlave.surface = ®ion; + m_mortarSide.get_inserted( MortarSide::Slave ) = surfaceSlave; + m_meshSlaveName = meshName; + } + else if( surfacePath.find( m_masterName ) != std::string::npos ) + { + MortarSurface surfaceMaster; + surfaceMaster.mesh = &mesh; + surfaceMaster.surface = ®ion; + m_mortarSide.get_inserted( MortarSide::Master ) = surfaceMaster; + } + } ); + } ); + + // debug log surfaces path + // string pathMaster; + // string pathMeshMaster; + // pathMaster = m_mortarSide[MortarSide::Master].surface->getPath(); + // pathMeshMaster = m_mortarSide[MortarSide::Master].mesh->getPath(); + // std::cout << "Path of master surface: " << pathMaster << std::endl; + // std::cout << "Path of master mesh level: " << pathMeshMaster << std::endl; + // string pathSlave; + // string pathMeshSlave; + // pathSlave = m_mortarSide[MortarSide::Slave].surface->getPath(); + // pathMeshSlave = m_mortarSide[MortarSide::Slave].mesh->getPath(); + // std::cout << "Path of slave surface: " << pathSlave << std::endl; + // std::cout << "Path of slave mesh level: " << pathMeshSlave << std::endl; + +} + +void SolidMechanicsMortarContact::computeMortarInterpolation ( + connectivityMapType & connectivityMap ) +{ + computeMortarInterpolation< ElementShape::Triangle, ElementShape::Triangle >( + connectivityMap[{ElementShape::Triangle, ElementShape::Triangle}] ); + + computeMortarInterpolation< ElementShape::Triangle, ElementShape::Quadrilateral >( + connectivityMap[{ElementShape::Triangle, ElementShape::Quadrilateral}] ); + + computeMortarInterpolation< ElementShape::Quadrilateral, ElementShape::Triangle >( + connectivityMap[{ElementShape::Quadrilateral, ElementShape::Triangle}] ); + + computeMortarInterpolation< ElementShape::Quadrilateral, ElementShape::Quadrilateral >( + connectivityMap[{ElementShape::Quadrilateral, ElementShape::Quadrilateral}] ); +} + +template< ElementShape slaveShape, ElementShape masterShape > +void SolidMechanicsMortarContact::computeMortarInterpolation( ArrayOfArrays< localIndex > const & connections ) +{ + + arrayView1d< localIndex const > faceListMaster = m_faceTypeToElementList.at( MortarSide::Master ).at( masterShape ).toViewConst(); + arrayView1d< localIndex const > faceListSlave = m_faceTypeToElementList.at( MortarSide::Slave ).at( slaveShape ).toViewConst(); + + localIndex numbConnections; + numbConnections = 0; + for( localIndex i = 0; i < connections.size(); ++i ) + { + numbConnections += connections.sizeOfArray( i ); + } + + // allocate mortar quantities + // auto const & slaveFE = getFE< slaveShape >(); + // auto const & masterFE = getFE< masterShape >(); + using FETypeMaster = typename selectFE< masterShape >::type; + using FETypeSlave = typename selectFE< slaveShape >::type; + localIndex constexpr nNodeMaster = FETypeMaster::numNodes; + localIndex constexpr nNodeSlave = FETypeSlave::numNodes; + // maximum number of expected subtriangles for each mortar pair + localIndex constexpr maxSubTri = nNodeMaster + nNodeSlave - 2; + localIndex nSubTriangles = maxSubTri * numbConnections; + + array1d< real64 > numbTrianglesInPairs( numbConnections ); + array2d< localIndex > cellPairs( nSubTriangles, 2 ); + array1d< real64 > subTriDeterminants( nSubTriangles ); + array3d< real64 > localCoordsSlave( nSubTriangles, nGPtri, 2 ); + array3d< real64 > localCoordsMaster( nSubTriangles, nGPtri, 2 ); + // get views + //arrayView1d< real64 > numbTrianglesInPairs_v = numbTrianglesInPairs.toView(); + arrayView2d< localIndex > cellPairs_v = cellPairs.toView(); + arrayView1d< real64 > subTriDeterminants_v = subTriDeterminants.toView(); + arrayView3d< real64 > localCoordsSlave_v = localCoordsSlave.toView(); + arrayView3d< real64 > localCoordsMaster_v = localCoordsMaster.toView(); + + if( numbConnections==0 ) + { + return; + } + + // get mesh manager objects for both sides + MeshLevel & meshMaster = *m_mortarSide.at( MortarSide::Master ).mesh; + FaceElementSubRegion const & surfMaster = m_mortarSide.at( MortarSide::Master ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + FaceManager const & faceManagerMaster = meshMaster.getFaceManager(); + NodeManager const & nodeManagerMaster = meshMaster.getNodeManager(); + arrayView2d< localIndex const > const elemsToFacesMaster = surfMaster.faceList().toViewConst(); + ArrayOfArraysView< localIndex const > const & faceToNodeMapMaster = faceManagerMaster.nodeList().toViewConst(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const coordsMaster = nodeManagerMaster.referencePosition(); + + MeshLevel & meshSlave = *m_mortarSide.at( MortarSide::Slave ).mesh; + FaceElementSubRegion const & surfSlave = m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(); + FaceManager const & faceManagerSlave = meshSlave.getFaceManager(); + NodeManager const & nodeManagerSlave = meshSlave.getNodeManager(); + arrayView2d< localIndex const > const elemsToFacesSlave = surfSlave.faceList().toViewConst(); + ArrayOfArraysView< localIndex const > const & faceToNodeMapSlave = faceManagerSlave.nodeList().toViewConst(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const coordsSlave = nodeManagerSlave.referencePosition(); + + array1d< localIndex > numConnectionsPerMaster( faceListMaster.size()); + for( localIndex im = 0; im < faceListMaster.size(); ++im ) + { + numConnectionsPerMaster[im] = connections.sizeOfArray( im ); + } + + array1d< localIndex > pairOffsets( faceListMaster.size()); + pairOffsets[0] = 0; + for( localIndex i = 1; i < faceListMaster.size(); ++i ) + { + pairOffsets[i] = pairOffsets[i-1] + numConnectionsPerMaster[i-1]; + } + + RAJA::ReduceSum< ReducePolicy< parallelHostPolicy >, localIndex > totTriSum( 0 ); + localIndex totTri = 0; + + + forAll< parallelHostPolicy >( faceListMaster.size(), [&]( localIndex im ) + { + localIndex masterFaceId = faceListMaster[im]; + localIndex k = pairOffsets[im]; + + for( localIndex is = 0; is < connections.sizeOfArray( im ); ++is ) + { + localIndex slaveFaceId = faceListSlave[connections( im, is )]; + + localIndex kfeM = elemsToFacesMaster[masterFaceId][0]; + localIndex kfeS = elemsToFacesSlave[slaveFaceId][0]; + + arraySlice1d< localIndex const > nodesMaster = faceToNodeMapMaster[kfeM]; + arraySlice1d< localIndex const > nodesSlave = faceToNodeMapSlave[kfeS]; + + localIndex nTriPair = processMortarPair< slaveShape, masterShape >( slaveFaceId, + masterFaceId, + nodesSlave, + nodesMaster, + coordsSlave, + coordsMaster, + cellPairs_v, + subTriDeterminants_v, + localCoordsSlave_v, + localCoordsMaster_v, + k + ); + + numbTrianglesInPairs[k] = nTriPair; + totTriSum += nTriPair; + ++k; + } + } ); + + totTri += totTriSum.get(); + + // allocate and populate member maps with segment based preprocessed info + array1d< localIndex > triCellsListSlave( totTri ); + array1d< localIndex > triCellsListMaster( totTri ); + array1d< real64 > triCellsDetList( totTri ); + array3d< real64 > gpLocalCoordsSlave( totTri, nGPtri, 2 ); + array3d< real64 > gpLocalCoordsMaster( totTri, nGPtri, 2 ); + + localIndex kp = 0; + for( localIndex i=0; i +localIndex SolidMechanicsMortarContact::processMortarPair( localIndex const slaveFaceId, + localIndex const masterFaceId, + arraySlice1d< localIndex const > const & nodesSlave, + arraySlice1d< localIndex const > const & nodesMaster, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & coordsSlave, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & coordsMaster, + arrayView2d< localIndex > & cellPairs, + arrayView1d< real64 > & subTriDeterminants, + arrayView3d< real64 > & localCoordsSlave, + arrayView3d< real64 > & localCoordsMaster, + localIndex const & kPair ) +{ + + std::cout << "Processing mortar pair slave/master: " << slaveFaceId << " - " << masterFaceId << std::endl; + + using FETypeMaster = typename selectFE< masterShape >::type; + using FETypeSlave = typename selectFE< slaveShape >::type; + + constexpr static localIndex numNodeMaster = FETypeMaster::numNodes; + constexpr static localIndex numNodeSlave = FETypeSlave::numNodes; + + // STEPS FOR THE SEGMENT BASED SCHEME + // 1. find auxiliary plane on the slave element (use computational geometry) + real64 planeCentroid[3]; + real64 planeNormal[3]; + real64 area = computationalGeometry::centroid_3DPolygon( nodesSlave, + coordsSlave, + planeCentroid, + planeNormal ); + + GEOS_UNUSED_VAR( area ); + + // 2. project slave and master nodes on the auxiliary plane (use computational geometry) + array2d< real64 > projSlave( numNodeSlave, 2 ); + array2d< real64 > projMaster( numNodeMaster, 2 ); + + for( localIndex i=0; i( coord3d, planeNormal, planeCentroid, projected2d ); + projSlave[i][0] = projected2d[0]; + projSlave[i][1] = projected2d[1]; + } + + for( localIndex i=0; i( coord3d, planeNormal, planeCentroid, projected2d ); + projMaster[i][0] = projected2d[0]; + projMaster[i][1] = projected2d[1]; + } + + + // 3. intersect slave and master element and return the resulting intersecting polygon + array2d< real64 > clipPoly( numNodeSlave, 2 ); + LvArray::tensorOps::copy< numNodeSlave, 2 >( clipPoly, projSlave ); + polygonClipping< numNodeSlave, numNodeMaster >( clipPoly, projMaster ); + + + if( !validateClip( clipPoly )) + { + std::cout << "Invalid clip polygon discarded" << std::endl; + return 0; + } + + localIndex clipSize = clipPoly.size( 0 ); + + std::cout << "Clipped polygon has " << clipSize << " vertices:" << std::endl; + for( localIndex i = 0; i < clipSize; ++i ) + { + std::cout << " (" << clipPoly( i, 0 ) << ", " << clipPoly( i, 1 ) << ")" << std::endl; + } + + // 4. project gauss point for each local triangle into master and slave side + for( localIndex i = 1; i < clipSize - 1; ++i ) + { + real64 coordsTri[3][2]; + // get coordinates of sub-triangle + coordsTri[0][0] = clipPoly( 0, 0 ); + coordsTri[0][1] = clipPoly( 0, 1 ); + coordsTri[1][0] = clipPoly( i, 0 ); + coordsTri[1][1] = clipPoly( i, 1 ); + coordsTri[2][0] = clipPoly( i+1, 0 ); + coordsTri[2][1] = clipPoly( i+1, 1 ); + + real64 const subTriDeterminant = (coordsTri[1][0] - coordsTri[0][0]) * (coordsTri[2][1] - coordsTri[0][1]) - + (coordsTri[2][0] - coordsTri[0][0]) * (coordsTri[1][1] - coordsTri[0][1]); + + if( subTriDeterminant < 1e-10 ) + { + GEOS_ERROR( "Found degenerate triangle" ); + } + + real64 xiSlave[nGPtri][2] = {{}}; + real64 xiMaster[nGPtri][2] = {{}}; + + projectGP< slaveShape >( coordsTri, projSlave.toViewConst(), xiSlave ); + + projectGP< masterShape >( coordsTri, projMaster.toViewConst(), xiMaster ); + + // 5. populate output lists + localIndex maxTri = numNodeSlave+numNodeMaster-2; + cellPairs( kPair*maxTri+i-1, 0 ) = slaveFaceId; + cellPairs( kPair*maxTri+i-1, 1 ) = masterFaceId; + + subTriDeterminants[kPair*maxTri+i-1] = subTriDeterminant; + + for( localIndex j = 0; j < nGPtri; ++j ) + { + localCoordsSlave( kPair*maxTri+i-1, j, 0 ) = xiSlave[j][0]; + localCoordsSlave( kPair*maxTri+i-1, j, 1 ) = xiSlave[j][1]; + + localCoordsMaster( kPair*maxTri+i-1, j, 0 ) = xiMaster[j][0]; + localCoordsMaster( kPair*maxTri+i-1, j, 1 ) = xiMaster[j][1]; + } + + std::cout << std::endl; + + } + + return clipSize - 2; // return number of triangles for the mortar pair + +} + + +template< MortarSide side > +void SolidMechanicsMortarContact::projectPointInPlane( real64 const (&coord3d)[3], + real64 const (&normal)[3], + real64 const (&origin)[3], + real64 (& proj2d)[2] ) +{ + + real64 proj3d[3]; + computationalGeometry::LinePlaneIntersection( normal, + coord3d, + normal, + origin, + proj3d ); + + // find coordinates in the local plane + real64 tmp[3]; + + if( std::fabs( normal[0] ) < 0.9 ) + { + tmp[0] = 1.0; + tmp[1] = 0.0; + tmp[2] = 0.0; + } + else + { + tmp[0] = 0.0; + tmp[1] = 1.0; + tmp[2] = 0.0; + } + real64 u[3]; + real64 v[3]; + + LvArray::tensorOps::crossProduct( u, normal, tmp ); + LvArray::tensorOps::crossProduct( v, u, normal ); + LvArray::tensorOps::normalize< 3 >( u ); + LvArray::tensorOps::normalize< 3 >( v ); + + if constexpr ( side == MortarSide::Slave) + { + // flip v to ensure the slave polygon is in CCW order in 2D. + LvArray::tensorOps::scale< 3 >( v, -1 ); + } + + LvArray::tensorOps::subtract< 3 >( proj3d, origin ); + + proj2d[0] = LvArray::tensorOps::AiBi< 3 >( u, proj3d ); + proj2d[1] = LvArray::tensorOps::AiBi< 3 >( v, proj3d ); +} + +template< localIndex sizePoly, localIndex sizeClipper > +void SolidMechanicsMortarContact::polygonClipping( array2d< real64 > & poly, + array2d< real64 > & clipPoly ) +{ + // input polygons must be in CCW order! + for( localIndex i = 0; i < sizeClipper; ++i ) + { + localIndex k = (i + 1) % sizeClipper; + clip( poly, clipPoly( i, 0 ), clipPoly( i, 1 ), clipPoly( k, 0 ), clipPoly( k, 1 )); + } + +} + + +// intersect lines [(x1,y1),(x2,y2)] with [(x3,y3),(x4,y4)] +void SolidMechanicsMortarContact::intersect( real64 x1, real64 y1, real64 x2, real64 y2, + real64 x3, real64 y3, real64 x4, real64 y4, + real64 & xInt, real64 & yInt ) +{ + real64 numX = (x1*y2 - y1*x2) * (x3-x4) - + (x1-x2) * (x3*y4 - y3*x4); + real64 denX = (x1-x2) * (y3-y4) - (y1-y2) * (x3-x4); + + real64 numY = (x1*y2 - y1*x2) * (y3-y4) - + (y1-y2) * (x3*y4 - y3*x4); + real64 denY = (x1-x2) * (y3-y4) - (y1-y2) * (x3-x4); + + if( std::fabs( denX ) < 1e-10 || std::fabs( denY ) < 1e-10 ) + { + // check for coinciding vertices + if( std::fabs( x1-x3 ) < 1e-8 && std::fabs( y1-y3 ) < 1e-8 ) + { + xInt = x1; yInt = y1; return; + } + if( std::fabs( x1-x4 ) < 1e-8 && std::fabs( y1-y4 ) < 1e-8 ) + { + xInt = x1; yInt = y1; return; + } + if( std::fabs( x2-x3 ) < 1e-8 && std::fabs( y2-y3 ) < 1e-8 ) + { + xInt = x2; yInt = y2; return; + } + if( std::fabs( x2-x4 ) < 1e-8 && std::fabs( y2-y4 ) < 1e-8 ) + { + xInt = x2; yInt = y2; return; + } + std::cout << "Line 1: (" << x1 << ", " << y1 << ") - (" << x2 << ", " << y2 << ")" << std::endl; + std::cout << "Line 2: (" << x3 << ", " << y3 << ") - (" << x4 << ", " << y4 << ")" << std::endl; + GEOS_ERROR( "Lines are parallel" ); + } + else + { + xInt = numX / denX; + yInt = numY / denY; + } + +} + +bool SolidMechanicsMortarContact::validateClip( array2d< real64 > & clipPoly ) +{ + // remove duplicated vertices and check if the clipping polygon is valid (non-degenerate) + // return true if the resulting polygon is valid + + if( clipPoly.size( 0 ) < 3 ) + return false; + + // remove duplicated + real64 tol = 1e-8; + + // create a temporary container for filtered vertices + std::vector< std::array< real64, 2 > > filtered; + filtered.reserve( clipPoly.size( 0 )); + + filtered.push_back( {clipPoly( 0, 0 ), clipPoly( 0, 1 )} ); + + // remove consecutive duplicates + for( localIndex i = 1; i < clipPoly.size( 0 ); ++i ) + { + real64 dx = clipPoly( i, 0 ) - filtered.back()[0]; + real64 dy = clipPoly( i, 1 ) - filtered.back()[1]; + + if( std::sqrt( dx*dx + dy*dy ) > tol ) + { + filtered.push_back( {clipPoly( i, 0 ), clipPoly( i, 1 )} ); + } + } + + // remove last vertex if it is duplicate of the first (closed polygon) + if( static_cast< int >(filtered.size()) > 1 ) + { + real64 dx = filtered.front()[0] - filtered.back()[0]; + real64 dy = filtered.front()[1] - filtered.back()[1]; + if( std::sqrt( dx*dx + dy*dy ) < tol ) + { + filtered.pop_back(); + } + } + + if( static_cast< int >(filtered.size()) < 3 ) + return false; // polygon is degenerate + + clipPoly.resize( filtered.size(), 2 ); + + for( size_t i = 0; i < filtered.size(); ++i ) + { + int ii = static_cast< int >(i); + clipPoly( ii, 0 ) = filtered[i][0]; + clipPoly( ii, 1 ) = filtered[i][1]; + } + + return true; + +} + +template< ElementShape shape > +void SolidMechanicsMortarContact::projectGP( real64 const (&coordsTri)[3][2], + arrayView2d< real64 const > const & coordsElem, + real64 (& xi)[nGPtri][2] ) +{ + + constexpr localIndex numNodes = (shape == ElementShape::Triangle) ? 3 : 4; + + //auto const & FE = getFE< shape >(); + + using FEType = typename selectFE< shape >::type; + + // newton params + localIndex itMax = 3; + real64 tol = 1e-9; + + for( localIndex i = 0; i < nGPtri; ++i ) + { + real64 xiProj[2] = {0.0, 0.0}; + + real64 Ntri[3]; + + feTriangleCell::calcN( i, Ntri ); + real64 coordGP[2]; + LvArray::tensorOps::Ri_eq_AjiBj< 2, 3 >( coordGP, coordsTri, Ntri ); + + real64 Nq[numNodes]; + FEType::calcN( xiProj, Nq ); + permuteN< numNodes >( Nq ); + + real64 rhs[2] = {0.0, 0.0}; + LvArray::tensorOps::Ri_eq_AjiBj< 2, numNodes >( rhs, coordsElem, Nq ); + LvArray::tensorOps::subtract< 2 >( rhs, coordGP ); + + localIndex iter = 0; + + while( LvArray::tensorOps::l2Norm< 2 >( rhs ) > tol && iter < itMax ) + { + iter = iter + 1; + real64 dN[numNodes][2]= {{}}; + FEType::calcGradN( xiProj, dN ); + + real64 Jt[2][2] = {{}}; + real64 J[2][2] = {{}}; + LvArray::tensorOps::Rij_eq_AkiBkj< 2, 2, numNodes >( Jt, dN, coordsElem ); + LvArray::tensorOps::transpose< 2, 2 >( J, Jt ); + + real64 dxi[2]; + + bool success = denseLinearAlgebra::details::solveTwoByTwoSystem( J, rhs, dxi ); + + if( success ) + { + xiProj[0] -= dxi[0]; + xiProj[1] -= dxi[1]; + } + else + { + GEOS_ERROR( "Failed to solve linear system in GP projection algorithm" ); + } + + FEType::calcN( xiProj, Nq ); + permuteN< numNodes >( Nq ); + + LvArray::tensorOps::Ri_eq_AjiBj< 2, numNodes >( rhs, coordsElem, Nq ); + LvArray::tensorOps::subtract< 2 >( rhs, coordGP ); + } + + if( iter == itMax ) + { + GEOS_ERROR( "SolidMechanicsMortarContact::projectGP - Newton raphson not converged" ); + } + + if( !checkInFE< shape >( xiProj[0], xiProj[1] )) + { + GEOS_ERROR( "Projected GP is outside the reference element" ); + } + + xi[i][0] = xiProj[0]; + xi[i][1] = xiProj[1]; + + } + +} + +template< ElementShape shape > +bool SolidMechanicsMortarContact::checkInFE( real64 xi0, real64 xi1 ) +{ + real64 tol = 1e-7; + bool isInRange = false; + + if constexpr (shape == ElementShape::Triangle) + { + // barycentric coords must be >= 0 and xi1+xi2 <= 1 + isInRange = (xi0 >= -tol && xi1 >= -tol && + xi0 <= 1.0 + tol && xi1 <= 1.0 + tol && + xi0 + xi1 <= 1.0 + tol); + } + else if constexpr (shape == ElementShape::Quadrilateral) + { + // in square [-1,1] x [-1,1] + isInRange = (xi0 >= -1.0 - tol && xi0 <= 1.0 + tol && + xi1 >= -1.0 - tol && xi1 <= 1.0 + tol); + } + + return isInRange; +} + + +// clip polygon against clipping line [(xc1,yc1),(xc2,yc2)] +void SolidMechanicsMortarContact::clip( array2d< real64 > & poly, + real64 xc1, real64 yc1, real64 xc2, real64 yc2 ) +{ + constexpr localIndex maxPoints = 8; + + real64 clippedPoly[maxPoints][2]; + localIndex clippedSize = 0; + + localIndex polySize = poly.size( 0 ); + + for( localIndex i = 0; i < polySize; ++i ) + { + localIndex k = (i+1) % polySize; + real64 x1 = poly( i, 0 ); + real64 y1 = poly( i, 1 ); + real64 x2 = poly( k, 0 ); + real64 y2 = poly( k, 1 ); + + // positions of points w.r.t clipper line + real64 i_pos = (xc2 - xc1) * (y1 - yc1) - (yc2 - yc1) * (x1 - xc1); + real64 k_pos = (xc2 - xc1) * (y2 - yc1) - (yc2 - yc1) * (x2 - xc1); + + real64 xInt = 0.0; + real64 yInt = 0.0; + + // case 1: both points inside + if( i_pos >= 0 && k_pos >= 0 ) + { + // add only second point + clippedPoly[clippedSize][0] = x2; + clippedPoly[clippedSize][1] = y2; + ++clippedSize; + } + // case 2: only first point is outside + else if( i_pos < 0 && k_pos >= 0 ) + { + // compute intersection point + intersect( x1, y1, x2, y2, xc1, yc1, xc2, yc2, xInt, yInt ); + + // add intersection point and second point + clippedPoly[clippedSize][0] = xInt; + clippedPoly[clippedSize][1] = yInt; + ++clippedSize; + clippedPoly[clippedSize][0] = x2; + clippedPoly[clippedSize][1] = y2; + ++clippedSize; + } + // case 3: only second point is outside + else if( i_pos >= 0 && k_pos < 0 ) + { + // add intersection point + intersect( x1, y1, x2, y2, xc1, yc1, xc2, yc2, xInt, yInt ); + clippedPoly[clippedSize][0] = xInt; + clippedPoly[clippedSize][1] = yInt; + ++clippedSize; + } + // case 4: both points are outside, no points are added + } + + // Update the original polygon with the clipped points + poly.resize( clippedSize, 2 ); + for( localIndex i = 0; i < clippedSize; ++i ) + { + poly( i, 0 ) = clippedPoly[i][0]; + poly( i, 1 ) = clippedPoly[i][1]; + } +} + + +template< localIndex numNodes > +void SolidMechanicsMortarContact::permuteN( real64 (& N)[numNodes] ) +{ + if constexpr (numNodes == 3) + { + // permutation not need for triangles + return; + } + else if constexpr (numNodes == 4) + { + real64 Ntmp[numNodes]; + LvArray::tensorOps::copy< numNodes >( Ntmp, N ); + localIndex permutation[4]; + permutation[0] = 0; + permutation[1] = 1; + permutation[2] = 3; + permutation[3] = 2; + + for( localIndex i=0; i connections; + getMortarConnections( slaveShape, masterShape, connections ); + connectivityMap.get_inserted( {slaveShape, masterShape} ) = connections; + } + } +} + + +void SolidMechanicsMortarContact::getMortarConnections( ElementShape slaveShape, + ElementShape masterShape, + ArrayOfArrays< localIndex > & connections ) +{ + + // using smart pointers for the trees + std::unique_ptr< TreeNodeMortar > treeMaster = std::make_unique< TreeNodeMortar >(); + std::unique_ptr< TreeNodeMortar > treeSlave = std::make_unique< TreeNodeMortar >(); + + // get list of surfaces + array1d< localIndex > surfMaster = m_faceTypeToElementList.at( MortarSide::Master ).at( masterShape ); + array1d< localIndex > surfSlave = m_faceTypeToElementList.at( MortarSide::Slave ).at( slaveShape ); + + if( surfMaster.size() == 0 || surfSlave.size() == 0 ) + { + return; + } + + + array1d< localIndex > surfRootMaster( surfMaster.size()); + array1d< localIndex > surfRootSlave( surfSlave.size()); + for( localIndex i=0; icreateNode( *m_mortarSide.at( MortarSide::Master ).mesh, + m_mortarSide.at( MortarSide::Master ).surface->getUniqueSubRegion< FaceElementSubRegion >(), + surfMaster, + surfRootMaster ); + treeSlave->createNode( *m_mortarSide.at( MortarSide::Slave ).mesh, + m_mortarSide.at( MortarSide::Slave ).surface->getUniqueSubRegion< FaceElementSubRegion >(), + surfSlave, + surfRootSlave ); + + // initialize connectivity map + connections.resize( surfRootMaster.size()); + + // perform contact search and populate connectivity map + contactSearch( treeSlave, treeMaster, connections ); + + //localIndex numConnections = 0; + + for( localIndex i=0; i const & nodeSlave, + std::unique_ptr< TreeNodeMortar > const & nodeMaster, + ArrayOfArrays< localIndex > & connections ) +{ + if( !checkIntersection( nodeSlave, nodeMaster )) + return; + + if((nodeSlave->isLeaf) && (nodeMaster->isLeaf)) + { + connections.emplaceBack( nodeMaster->leafId, nodeSlave->leafId ); + return; + } + + if( !nodeSlave->isLeaf ) + { + contactSearch( nodeSlave->left, nodeMaster, connections ); + contactSearch( nodeSlave->right, nodeMaster, connections ); + } + else + { + contactSearch( nodeSlave, nodeMaster->left, connections ); + contactSearch( nodeSlave, nodeMaster->right, connections ); + } +} + +bool SolidMechanicsMortarContact::checkIntersection( std::unique_ptr< TreeNodeMortar > const & nodeSlave, std::unique_ptr< TreeNodeMortar > const & nodeMaster ) +{ + real64 * pS = nodeSlave->polytop; + real64 * pM = nodeMaster->polytop; + + for( localIndex k=0; k<9; ++k ) + { + if( pS[2*k]>=pM[2*k+1] ) + return false; + if( pM[2*k]>=pS[2*k+1] ) + return false; + } + + // if pass all checks, the two nodes intersect + return true; +} + + + +void TreeNodeMortar::createNode( MeshLevel const & mesh, + FaceElementSubRegion const & surf, + arrayView1d< localIndex > & surfId, + array1d< localIndex > & surfList ) +{ + + if( surfList.size() == 0 ) + { + GEOS_ERROR( "Mortar:createNode - empty list of surfaces" ); + } + + //GEOS_UNUSED_VAR(surfList); + FaceManager const & faceManager = mesh.getFaceManager(); + NodeManager const & nodeManager = mesh.getNodeManager(); + arrayView2d< localIndex const > const elemsToFaces = surf.faceList().toViewConst(); + ArrayOfArraysView< localIndex const > const & faceToNodeMap = faceManager.nodeList().toViewConst(); + arrayView2d< real64 const > const surfCenter = faceManager.faceCenter().toViewConst(); + localIndex nSurf = surfList.size(); + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const coords = nodeManager.referencePosition(); + + // compute primitive of polytopal bounding box for input list of surfaces + for( localIndex i=0; ipolytop[2*k] ) + this->polytop[2*k] = P; + // store maximum polytop primitive + if( P > this->polytop[2*k+1] ) + this->polytop[2*k+1] = P; + } + } + } + + // get direction of polytop maximum size + localIndex dir = 0; + real64 diff = -1; + for( localIndex k = 0; k<9; ++k ) + { + // add a small tolerance to avoid degenerate boxes + real64 d = polytop[2*k+1] - polytop[2*k] + 1e-3; + + // expand the polytop for safety + polytop[2*k] -= boundingBoxExpansion*d; + polytop[2*k+1] += boundingBoxExpansion*d; + + if( d > diff ) + { + dir = k; + diff = d; + } + } + + if( surfList.size() == 1 ) // leaf node + { + this->isLeaf = true; + this->leafId = surfList[0]; + return; + } + + // split the surface list into left and right child + std::vector< real64 > surfacePrimitive( nSurf ); + for( localIndex i=0; i surfSort = surfacePrimitive; + std::sort( surfSort.begin(), surfSort.end()); + real64 m = (nSurf % 2 == 1) ? surfSort[nSurf / 2] : (surfSort[nSurf / 2 - 1] + surfSort[nSurf / 2]) / 2.0; + + array1d< localIndex > surfLeft; + array1d< localIndex > surfRight; + + // divide left and right set of surface based on primitive value w.r.t median + for( localIndex i=0; ileft = std::make_unique< TreeNodeMortar >(); + this->left->createNode( mesh, surf, surfId, surfLeft ); + this->right = std::make_unique< TreeNodeMortar >(); + this->right->createNode( mesh, surf, surfId, surfRight ); +} + +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, SolidMechanicsMortarContact, string const &, Group * const ) +} +/* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.hpp new file mode 100644 index 00000000000..2e5e2ddb5cc --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.hpp @@ -0,0 +1,322 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +/* + * SolidMechanicsMortarContact.hpp + * + */ + +#ifndef GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSMORTARCONTACT_HPP_ +#define GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSMORTARCONTACT_HPP_ + +#include "physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp" +#include "mesh/utilities/ComputationalGeometry.hpp" + +namespace geos +{ + +class TreeNodeMortar; + +enum class ElementShape { Triangle, Quadrilateral }; +enum class MortarSide { Slave, Master }; + +using feTriangleCell = finiteElement::H1_TriangleFace_Lagrange1_Gauss4; + +struct MortarSurface +{ + MeshLevel * mesh = nullptr; + SurfaceElementRegion * surface = nullptr; +}; + +class SolidMechanicsMortarContact : public ContactSolverBase +{ +public: + SolidMechanicsMortarContact( const string & name, + Group * const parent ); + + ~SolidMechanicsMortarContact() override; + + /** + * @brief name of the node manager in the object catalog + * @return string that contains the catalog name to generate a new NodeManager object through the object catalog. + */ + static string catalogName() + { + return "SolidMechanicsMortarContact"; + } + /** + * @copydoc PhysicsSolverBase::getCatalogName() + */ + string getCatalogName() const override {return catalogName(); } + + using connectivityMapType = stdMap< std::pair< ElementShape, ElementShape >, ArrayOfArrays< localIndex > >; + + + virtual void registerDataOnMesh( dataRepository::Group & meshBodies ) override final; + + void registerMortarDataOnMesh( ); + + void setupDofs( DomainPartition const & domain, + DofManager & dofManager ) const override; + + virtual void setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity = true ) override final; + + virtual void implicitStepSetup( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) override final; + + virtual void implicitStepComplete( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) override final; + + virtual void assembleSystem( real64 const time, + real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override; + + virtual real64 calculateResidualNorm( real64 const & time_n, + real64 const & dt, + DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localRhs ) override; + + real64 calculateContactResidualNorm( DomainPartition const & domain, + DofManager const & dofManager, + arrayView1d< real64 const > const & localRhs ); + + virtual void applySystemSolution( DofManager const & dofManager, + arrayView1d< real64 const > const & localSolution, + real64 const scalingFactor, + real64 const dt, + DomainPartition & domain ) override; + + void updateState( DomainPartition & domain ) override final; + + void createFaceTypeListMortar( MortarSide side ); + + void createBubbleCellList( ) const; + + void setMortarSurfaces( DomainPartition & domain ); + + template< ElementShape shape > + void assembleMortarBubbles( DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + + +private: + + string m_meshSlaveName; + + stdMap< MortarSide, stdMap< ElementShape, array1d< localIndex > > > m_faceTypeToElementList; + + // holds pointers to master and slave MeshLevel and SurfaceElementRegion + stdMap< MortarSide, MortarSurface > m_mortarSide; + + // list of pairs of slave/master element ids for each mortar subtriangle + stdMap< MortarSide, stdMap< std::pair< ElementShape, ElementShape >, array1d< localIndex > > > m_triCells; + + // list of determinant of each mortar subtriangle + stdMap< std::pair< ElementShape, ElementShape >, array1d< real64 > > m_triCellsDet; + + // list local coordinates of gauss points of subtriangles on each mortar side + stdMap< MortarSide, stdMap< std::pair< ElementShape, ElementShape >, array3d< real64 > > > m_gpLocalCoords; + + // coordinates of triangle gauss points (they are private in the triangle class) + constexpr static localIndex nGPtri = feTriangleCell::numQuadraturePoints; + + + void addBubbleCouplingNumNonzeros( DofManager & dofManager, + arrayView1d< localIndex > const & rowLengths ) const; + + void addMortarCouplingNumNonzeros( DofManager & dofManager, + ElementShape const & slaveShape, + ElementShape const & masterShape, + connectivityMapType const & connectivityMap, + arrayView1d< localIndex > const & rowLengths ) const; + + void addBubbleCouplingSparsityPattern( DofManager const & dofManager, + SparsityPatternView< globalIndex > const & pattern ) const; + + void addMortarCouplingSparsityPattern( DofManager & dofManager, + ElementShape const & slaveShape, + ElementShape const & masterShape, + SparsityPatternView< globalIndex > const & pattern ) const; + + void assembleBubbles( real64 const dt, + DomainPartition & domain, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + + void assembleMortar( real64 const dt, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + template< ElementShape shpS, ElementShape shpM > + void assembleMortar( real64 const dt, + DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + void computeRotationMatrices( ); + + string m_slaveName; + string m_masterName; + + stdMap< ElementShape, std::unique_ptr< geos::finiteElement::FiniteElementBase > > m_faceTypeToMortarFiniteElements; + + struct viewKeyStruct : ContactSolverBase::viewKeyStruct + { + + constexpr static char const * masterString() { return "master"; } + + constexpr static char const * slaveString() { return "slave"; } + + }; + + void computeMortarInterpolation ( connectivityMapType & connectivityMap ); + + template< ElementShape slaveShape, ElementShape masterShape > + void computeMortarInterpolation( ArrayOfArrays< localIndex > const & connections ); + + template< ElementShape slaveShape, ElementShape masterShape > + localIndex processMortarPair( localIndex const slaveFaceId, + localIndex const masterFaceId, + arraySlice1d< localIndex const > const & nodesSlave, + arraySlice1d< localIndex const > const & nodesMaster, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & coordsSlave, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & coordsMaster, + arrayView2d< localIndex > & cellPairs, + arrayView1d< real64 > & subTriDeterminants, + arrayView3d< real64 > & localCoordsSlave, + arrayView3d< real64 > & localCoordsMaster, + localIndex const & kPair ); + + + template< MortarSide side > + void projectPointInPlane( real64 const (&coord3d)[3], + real64 const (&normal)[3], + real64 const (&origin)[3], + real64 ( &proj2d )[2] ); + + template< localIndex sizePoly, localIndex sizeClipper > + void polygonClipping( array2d< real64 > & poly, array2d< real64 > & clipPoly ); + + void intersect( real64 x1, real64 y1, real64 x2, real64 y2, + real64 x3, real64 y3, real64 x4, real64 y4, + real64 & xInt, real64 & yInt ); + + void clip( array2d< real64 > & poly, + real64 xc1, real64 yc1, real64 xc2, real64 yc2 ); + + bool validateClip( array2d< real64 > & clipPoly ); + + template< ElementShape shape > + void projectGP( real64 const (&coordsTri)[3][2], + arrayView2d< real64 const > const & coordsElem, + real64 ( &xi )[nGPtri][2] ); + + template< ElementShape shape > + bool checkInFE( real64 xi0, real64 xi1 ); + + template< localIndex numNodes > + void permuteN( real64 ( &N )[numNodes] ); + + void getConnectivityMap( connectivityMapType & connectivityMap ); + + void getMortarConnections( ElementShape slaveShape, ElementShape masterShape, ArrayOfArrays< localIndex > & connections ); + + // Tandem traversal contact search between master and slave nodes + void contactSearch( std::unique_ptr< TreeNodeMortar > const & nodeMaster, + std::unique_ptr< TreeNodeMortar > const & nodeSlave, + ArrayOfArrays< localIndex > & connections ); + + // Check intersection between two bounding boxes using polytops primitives + bool checkIntersection( std::unique_ptr< TreeNodeMortar > const & nodeMaster, + std::unique_ptr< TreeNodeMortar > const & nodeSlave ); + + template< ElementShape S > + struct selectFE + { + static auto getTypeImpl() + { + if constexpr (S == ElementShape::Triangle) + return finiteElement::H1_TriangleFace_Lagrange1_Gauss4{}; + else if constexpr (S == ElementShape::Quadrilateral) + return finiteElement::H1_QuadrilateralFace_Lagrange1_GaussLegendre2{}; + else + static_assert( S == ElementShape::Triangle || S == ElementShape::Quadrilateral, + "Unsupported ElementShape" ); + } + + using type = decltype(getTypeImpl()); + }; + +}; + +// binary tree for efficient contact search +class TreeNodeMortar +{ +public: + std::unique_ptr< TreeNodeMortar > left = nullptr; + std::unique_ptr< TreeNodeMortar > right = nullptr; + bool isLeaf = false; + + // primitives of polytopal bounding box + real64 polytop[18] = + { + 1e100, -1e100, 1e100, -1e100, 1e100, + -1e100, 1e100, -1e100, 1e100, -1e100, + 1e100, -1e100, 1e100, -1e100, 1e100, + -1e100, 1e100, -1e100 + }; + + static constexpr real64 boundingBoxExpansion = 0.025; + + static constexpr localIndex polytopPrimitives[3][9] = + { + { 1, 0, 0, 1, 1, 0, 1, 1, 0 }, + { 0, 1, 0, 1, 0, 1, -1, 0, 1 }, + { 0, 0, 1, 0, 1, 1, 0, -1, -1 } + }; + + // id of face corresponding to leaf nodes + localIndex leafId; + + // populate a tree node (use recursion) + void createNode( MeshLevel const & mesh, + FaceElementSubRegion const & surf, + arrayView1d< localIndex > & surfId, + array1d< localIndex > & surfList ); + +}; + + +} /* namespace geos */ + + + +#endif /* GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSMORTARCONTACT_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsMortarContactKernels.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsMortarContactKernels.hpp new file mode 100644 index 00000000000..63eeaebdd24 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/kernels/SolidMechanicsMortarContactKernels.hpp @@ -0,0 +1,435 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * 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 SolidMechanicsALMKernelsBase.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_CONTACT_KERNELS_SOLIDMECHANICSMORTARCONTACTKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_CONTACT_KERNELS_SOLIDMECHANICSMORTARCONTACTKERNELS_HPP_ + +#include "finiteElement/kernelInterface/InterfaceKernelBase.hpp" +#include "codingUtilities/Utilities.hpp" + +namespace geos +{ + +namespace solidMechanicsMortarContactKernels +{ + +/** + * @brief Implements kernels for ALM. + * @copydoc geos::finiteElement::InterfaceKernelBase + * + */ +template< typename CONSTITUTIVE_TYPE, + typename FE_TYPE > +class MortarContactKernels : + public finiteElement::InterfaceKernelBase< CONSTITUTIVE_TYPE, + FE_TYPE, + 3, 3 > +{ +public: + /// Alias for the base class; + using Base = finiteElement::InterfaceKernelBase< CONSTITUTIVE_TYPE, + FE_TYPE, + 3, 3 >; + + /// Alias for the sub triangle element type + using feTriangleCell = finiteElement::H1_TriangleFace_Lagrange1_Gauss4; + + /// Number of nodes per element...which is equal to the + /// numTestSupportPointPerElem and numTrialSupportPointPerElem by definition. + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + + /// Compile time value for the number of quadrature points per sub triangle. + static constexpr int numQuadraturePointsPerElem = feTriangleCell::numQuadraturePoints; + + /// The number of displacement dofs per element. + static constexpr int numUdofs = numNodesPerElem * 3; + + /// The number of lagrange multiplier dofs per element. + static constexpr int numTdofs = 3; + + using Base::m_dofNumber; + using Base::m_dofRankOffset; + using Base::m_finiteElementSpace; + using Base::m_matrix; + using Base::m_rhs; + + /** + * @brief Constructor + * @copydoc geos::finiteElement::InterfaceKernelBase::InterfaceKernelBase + */ + MortarContactKernels( NodeManager const & nodeManager, + EdgeManager const & edgeManager, + FaceManager const & faceManager, + localIndex const targetRegionIndex, + FaceElementSubRegion & subRegion, + FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, + arrayView1d< globalIndex const > const uDofNumber, + globalIndex const rankOffset, + CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > const inputRhs, + real64 const inputDt, + FaceElementSubRegion const & multipliersRegion, + FaceElementSubRegion const & primaryVarRegion, + arrayView1d< localIndex const > const faceElementList1, + arrayView1d< localIndex const > const faceElementList2, + arrayView1d< real64 const > const subTriangleDeterminants, + arrayView3d< real64 const > const gpLocalCoords, + string const tractionDofKey ): + Base( nodeManager, + edgeManager, + faceManager, + targetRegionIndex, + subRegion, + finiteElementSpace, + inputConstitutiveType, + uDofNumber, + rankOffset, + inputMatrix, + inputRhs, + inputDt ), + //m_X( nodeManager.referencePosition()), + m_faceToNodes( faceManager.nodeList().toViewConst()), + m_elemsToFaces( primaryVarRegion.faceList().toViewConst()), + m_faceElementList1( faceElementList1 ), + m_faceElementList2( faceElementList2 ), + m_subTriangleDeterminants( subTriangleDeterminants ), + m_gpLocalCoords( gpLocalCoords ), + m_traction( multipliersRegion.getField< fields::contact::traction >().toViewConst() ), + m_tDofNumber( multipliersRegion.getReference< globalIndex_array >( tractionDofKey ).toViewConst() ), + m_displacement( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ), + m_rotationMatrix( multipliersRegion.getField< fields::contact::rotationMatrix >().toViewConst()) + {} + + //*************************************************************************** + /** + * @copydoc finiteElement::InterfaceKernelBase::StackVariables + */ + struct StackVariables + { + +public: + + GEOS_HOST_DEVICE + StackVariables(): + dispEqnRowIndices{}, + dispColIndices{}, + tEqnRowIndices{}, + tColIndices{}, + localRu{}, + localRt{}, + localAtu{ {} }, + localAut{ {} }, + uLocal{}, + det{}, + localRotationMatrix{ {} } + {} + + /// C-array storage for the element local row degrees of freedom. + localIndex dispEqnRowIndices[numUdofs]; + + /// C-array storage for the element local column degrees of freedom. + globalIndex dispColIndices[numUdofs]; + + /// C-array storage for the traction local row degrees of freedom. + localIndex tEqnRowIndices[numTdofs]; + + /// C-array storage for the element local column degrees of freedom. + globalIndex tColIndices[numTdofs]; + + /// C-array storage for the element local Ru residual vector. + real64 localRu[numUdofs]; + + /// C-array storage for the element local Rt residual vector. + real64 localRt[numTdofs]; + + /// C-array storage for the element local Atu matrix. + real64 localAtu[numTdofs][numUdofs]; + + /// C-array storage for the element local Aut matrix. + real64 localAut[numUdofs][numTdofs]; + + /// C-array storage for the element local displacement vector + real64 uLocal[numUdofs]; + + /// C-array storage for the subtriangle cell determinants + real64 det[numQuadraturePointsPerElem]; + + /// C-array storage for rotation matrix + real64 localRotationMatrix[3][3]; + + + }; + + //*************************************************************************** + + /** + * @copydoc ::geos::finiteElement::InterfaceKernelBase::kernelLaunch + * + */ + //START_kernelLauncher + template< typename POLICY, + typename KERNEL_TYPE > + static + real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ) + { + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( numElems ); + + // Define a RAJA reduction variable to get the maximum residual contribution. + RAJA::ReduceMax< ReducePolicy< POLICY >, real64 > maxResidual( 0 ); + + // Loop over all existing triangular integration subcells + forAll< POLICY >( kernelComponent.m_faceElementList1.size(), + [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + typename KERNEL_TYPE::StackVariables stack; + + kernelComponent.setup( k, stack ); + for( integer q=0; q( matRRtAtu, stack.localRotationMatrix, stack.localAtu ); + LvArray::tensorOps::copy< numTdofs, numUdofs >( stack.localAtu, matRRtAtu ); + LvArray::tensorOps::transpose< numUdofs, numTdofs >( stack.localAut, stack.localAtu ); + + // Compute the traction contribute of the local residuals + LvArray::tensorOps::Ri_eq_AijBj< numUdofs, numTdofs >( rhsU, stack.localAut, m_traction[ke] ); + LvArray::tensorOps::scaledAdd< numUdofs >( stack.localRu, rhsU, 1.0 ); + + // Compute the displacemement contribute of the local residual + LvArray::tensorOps::Ri_eq_AijBj< numTdofs, numUdofs >( rhsT, stack.localAtu, stack.uLocal ); + LvArray::tensorOps::scaledAdd< numTdofs >( stack.localRt, rhsT, 1.0 ); + + fillGlobalMatrix( stack ); + + return 0.0; + } + + + /** + * @brief Fill global matrix and residual vector + * + * @param stack stack variables + */ + GEOS_HOST_DEVICE + void fillGlobalMatrix( StackVariables & stack ) const + { + + + for( localIndex i=0; i < numTdofs; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.tEqnRowIndices[ i ] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) continue; + + // TODO: May not need to be an atomic operation + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localRt[i] ); + + // Fill in matrix block Atu + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof, + stack.dispColIndices, + stack.localAtu[i], + numUdofs ); + + } + + for( localIndex i=0; i < numUdofs; ++i ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.dispEqnRowIndices[ i ] ); + + if( dof < 0 || dof >= m_matrix.numRows() ) continue; + + // Is it necessary? Each row should be indepenedent + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localRu[i] ); + + // Fill in matrix + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof, + stack.tColIndices, + stack.localAut[i], + numTdofs ); + + } + + } + + +protected: + + /// The array of array containing the face to node map. + ArrayOfArraysView< localIndex const > const m_faceToNodes; + + /// The array of array containing the element to face map. + arrayView2d< localIndex const > const m_elemsToFaces; + + /// The array containing the list of face element of the same type for the traction side + /// for each mortar subtriangle. + arrayView1d< localIndex const > const m_faceElementList1; + + /// The array containing the list of face element of the same type for the displacement side + /// for each mortar subtriangle. + arrayView1d< localIndex const > const m_faceElementList2; + + /// The array containing the determinants of the subtriangles for mortar integration. + arrayView1d< real64 const > const m_subTriangleDeterminants; + + /// The array containing the local gauss point coordinates projected on the displacement side + arrayView3d< real64 const > const m_gpLocalCoords; + + /// The array containing the traction field + arrayView2d< real64 const > const m_traction; + + /// The array holding the traction degrees of freedom + arrayView1d< globalIndex const > const m_tDofNumber; + + /// The array containing the displacement field + arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const m_displacement; + + + /// The array containing the rotation matrix for each element. + arrayView3d< real64 const > const m_rotationMatrix; + +}; + +/// The factory used to construct a Mortar kernel. +using MortarFactory = finiteElement::InterfaceKernelFactory< MortarContactKernels, + arrayView1d< globalIndex const > const, + globalIndex const, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + FaceElementSubRegion const &, + FaceElementSubRegion const &, + arrayView1d< localIndex const > const &, + arrayView1d< localIndex const > const &, + arrayView1d< real64 const > const &, + arrayView3d< real64 const > const &, + string const >; + + +} // namespace solidMechanicsMortarContactKernels + +} // namespace geos + + +#endif /* GEOS_PHYSICSSOLVERS_CONTACT_KERNELS_SOLIDMECHANICSMORTARCONTACTKERNELS_HPP_ */