diff --git a/inputFiles/mortar/CantileverBeam.xml b/inputFiles/mortar/CantileverBeam.xml new file mode 100644 index 00000000000..89348c71a66 --- /dev/null +++ b/inputFiles/mortar/CantileverBeam.xml @@ -0,0 +1,253 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt b/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt index ced15c47ba5..48ebd7ad1da 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 @@ -68,7 +69,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) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/CMakeLists.txt b/src/coreComponents/physicsSolvers/solidMechanics/contact/CMakeLists.txt index 6d6457f7597..9014b6551d4 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 @@ -30,4 +31,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 3b1b8aed17f..196ccc20355 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 { @@ -67,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 ); @@ -122,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, diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index f0fcbb3e604..e934ff50887 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -1,3 +1,4 @@ + /* * ------------------------------------------------------------------------------------------------------------ * SPDX-License-Identifier: LGPL-2.1-only diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.cpp new file mode 100644 index 00000000000..4c0193485f9 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.cpp @@ -0,0 +1,395 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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/LogLevelsInfo.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" +#include "physicsSolvers/LogLevelsInfo.hpp" + +#include "constitutive/contact/FrictionSelector.hpp" + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; +using namespace fields; + +SolidMechanicsMortarContact::SolidMechanicsMortarContact( const string & name, + Group * const parent ): + ContactSolverBase( name, parent ) +{} + +SolidMechanicsMortarContact::~SolidMechanicsMortarContact() +{} + +void SolidMechanicsMortarContact::registerDataOnMesh( dataRepository::Group & meshBodies ) +{ + GEOS_MARK_FUNCTION; + ContactSolverBase::registerDataOnMesh( meshBodies ); +} + +void SolidMechanicsMortarContact::setupDofs( DomainPartition const & domain, + DofManager & dofManager ) const +{ + GEOS_MARK_FUNCTION; + SolidMechanicsLagrangianFEM::setupDofs( domain, dofManager ); +} + +void SolidMechanicsMortarContact::setupSystem( DomainPartition & domain, + DofManager & dofManager, + CRSMatrix< real64, globalIndex > & localMatrix, + ParallelVector & rhs, + ParallelVector & solution, + bool const setSparsity ) +{ + + GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( domain, dofManager, localMatrix, rhs, solution, setSparsity ); + +} + +void SolidMechanicsMortarContact::implicitStepSetup( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) +{ + + GEOS_MARK_FUNCTION; + SolidMechanicsLagrangianFEM::implicitStepSetup( time_n, dt, domain ); + + ////////////////////////////////////////////////////////////////////////////////// + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshName, + MeshLevel & mesh, + string_array const & ) + { + // naive way to assign master and slave mesh levels + if (meshName == "meshMaster"){ + this->m_meshMaster = &mesh; + } + + if (meshName == "meshSlave"){ + this->m_meshSlave = &mesh; + } + + ElementRegionManager const & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< FaceElementSubRegion >([&](const FaceElementSubRegion & subRegion ) + { + // assign surface to master or slave depending on object path (again, naive) + string surfacePath = subRegion.getPath(); + if (surfacePath.find("surfaceSlave") != std::string::npos){ + GEOS_ERROR_IF(m_surfaceSlave != nullptr,"Slave surface has already been assigned."); + this->m_surfaceSlave = &subRegion; + } + + if (surfacePath.find("surfaceMaster") != std::string::npos){ + GEOS_ERROR_IF(m_surfaceMaster != nullptr,"Master surface has already been assigned."); + this->m_surfaceMaster = &subRegion; + } + }); + + } ); + + + std::cout << "Processing master/slave connectivity" << std::endl; + + this->getConnectivityMap(); + + ////////////////////////////////////////////////////////////////////////////////// + GEOS_ERROR("Mortar solver is not implemented yet. "); + +} + +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 ); + + //ParallelMatrix parallel_matrix; + //parallel_matrix.create( localMatrix.toViewConst(), dofManager.numLocalDofs(), MPI_COMM_GEOS ); + //parallel_matrix.write("mech.mtx"); + +} + +void SolidMechanicsMortarContact::implicitStepComplete( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) +{ + + SolidMechanicsLagrangianFEM::implicitStepComplete( time_n, dt, domain ); + +} + +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 solidResidualNorm = SolidMechanicsLagrangianFEM::calculateResidualNorm( time, dt, domain, dofManager, localRhs ); + + return solidResidualNorm; + +} + +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 ); + +} + +void SolidMechanicsMortarContact::updateState( DomainPartition & domain ) +{ + GEOS_UNUSED_VAR( domain ); +} + +void SolidMechanicsMortarContact::getConnectivityMap() +{ + // using smart pointers for the trees + std::unique_ptr treeMaster = std::make_unique(); + std::unique_ptr treeSlave = std::make_unique(); + + // progressive list of surfaces to create tree roots + array1d surfRootMaster(m_surfaceMaster->size()) ; + array1d surfRootSlave(m_surfaceSlave->size()) ; + + for (localIndex i=0; i < m_surfaceMaster->size(); ++i) + { + surfRootMaster(i) = i; + } + + for (localIndex i=0; i < m_surfaceSlave->size(); ++i) + { + surfRootSlave(i) = i; + } + + // create binary trees recursively + std::cout << "Populating master tree" << std::endl; + treeMaster->createNode(*m_meshMaster,*m_surfaceMaster,surfRootMaster); + std::cout << "Populating slave tree" << std::endl; + treeSlave->createNode(*m_meshSlave,*m_surfaceSlave,surfRootSlave); + + // initialize connectivity map + localIndex nSurfSlave = m_surfaceSlave->size(); + m_connectivityMap.resize(nSurfSlave); + + // perform contact search and populate connectivity map + contactSearch(treeSlave,treeMaster); + + for (localIndex i=0; i const & nodeSlave, std::unique_ptr const & nodeMaster) +{ + if (!checkIntersection(nodeSlave,nodeMaster)) return; + + if ((nodeSlave->isLeaf) && (nodeMaster->isLeaf)) + { + //std::cout << "found intersection between slave " << nodeSlave->leafId << " and master " << nodeMaster->leafId << std::endl; + m_connectivityMap.emplaceBack(nodeSlave->leafId,nodeMaster->leafId); + return; + } + + if (!nodeSlave->isLeaf) + { + contactSearch(nodeSlave->left,nodeMaster); + contactSearch(nodeSlave->right,nodeMaster); + } + else + { + contactSearch(nodeSlave,nodeMaster->left); + contactSearch(nodeSlave,nodeMaster->right); + } +} + +bool SolidMechanicsMortarContact::checkIntersection(std::unique_ptr const & nodeSlave, std::unique_ptr const & nodeMaster) +{ + double* pS = nodeSlave->polytop; + double* 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, + array1d & surfList) +{ + + //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< double const > const surfCenter = faceManager.faceCenter().toViewConst(); + localIndex nSurf = surfList.size(); + arrayView2d const coords = nodeManager.referencePosition(); + + double const pP[3][9] = POLYTOP_PRIMITIVES; + + // 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] -= BOUNDING_BOX_EXPANSION*d; + polytop[2*k+1] += BOUNDING_BOX_EXPANSION*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 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 surfLeft; + array1d surfRight; + + // divide left and right set of surface based on primitive value w.r.t median + for (localIndex i=0; ileft = std::make_unique(); + this->left->createNode(mesh,surf,surfLeft); + this->right = std::make_unique(); + this->right->createNode(mesh,surf,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..5d26c2eda46 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsMortarContact.hpp @@ -0,0 +1,157 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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_ + +#define POLYTOP { \ + 1e100, -1e100, 1e100, -1e100, 1e100, -1e100, 1e100, -1e100, 1e100, \ + -1e100, 1e100, -1e100, 1e100, -1e100, 1e100, -1e100, 1e100, -1e100 \ +} + +#define POLYTOP_PRIMITIVES { \ + {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} \ +} + +#define BOUNDING_BOX_EXPANSION 0.025 + + +#include "physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp" + +namespace geos +{ + +class TreeNodeMortar; + +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(); } + + virtual void registerDataOnMesh( dataRepository::Group & meshBodies ) override final; + + virtual 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; + + 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; + +private: + + // store pointers to relevant mesh objects of master and slave side + MeshLevel const* m_meshSlave = nullptr; + MeshLevel const* m_meshMaster = nullptr; + FaceElementSubRegion const* m_surfaceMaster = nullptr; + FaceElementSubRegion const* m_surfaceSlave = nullptr; + + // map id of slave elements to id of connected master elements + ArrayOfArrays< localIndex > m_connectivityMap; + + // tandem traversal contact search + void contactSearch(std::unique_ptr const & nodeMaster, + std::unique_ptr const & nodeSlave); + + // check intersection between two bounding boxes using polytops primitives + bool checkIntersection(std::unique_ptr const & nodeMaster, + std::unique_ptr const & nodeSlave); + + void getConnectivityMap(); + +}; + +// binary tree for efficient contact search +class TreeNodeMortar + { + public: + std::unique_ptr left = nullptr; + std::unique_ptr right = nullptr; + bool isLeaf = false; + + // primitives of polytopal bounding box + double polytop[18] = POLYTOP; + + // id of face corresponding to leaf nodes + localIndex leafId; + + // populate a tree node (use recursion) + void createNode(MeshLevel const & mesh, + FaceElementSubRegion const & surf, + array1d & surfList); + + + + }; + +}; /* namespace geos */ + + + +#endif /* GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSMORTARCONTACT_HPP_ */ diff --git a/src/main/GEOS.code-workspace b/src/main/GEOS.code-workspace new file mode 100644 index 00000000000..24463b69648 --- /dev/null +++ b/src/main/GEOS.code-workspace @@ -0,0 +1,12 @@ +{ + "folders": [ + { + "path": "../.." + } + ], + "settings": { + "files.associations": { + "stdexcept": "cpp" + } + } +} \ No newline at end of file