feat: Mortar Method

inputFiles/mortar/CantileverBeam.xml

@@ -0,0 +1,253 @@
+<?xml version="1.0" ?>
+
+<Problem>
+  <Events
+    maxTime="1.0">
+
+    <SoloEvent
+      name="genSurfaceMaster"
+      target="/Solvers/SurfaceGenMaster"/>
+
+    <SoloEvent
+      name="genSurfaceSlave"
+      target="/Solvers/SurfaceGenSlave"/>
+
+    <PeriodicEvent
+      name="solverApplications"
+      beginTime="0.0"
+      forceDt="1.0"
+      target="/Solvers/Mortar"/>
+
+  </Events>
+
+  <ElementRegions>
+    <CellElementRegion
+      name="regionMaster"
+      meshBody="meshMaster"
+      cellBlocks="{ cb1 }"
+      materialList="{ steel }"/>
+
+    <CellElementRegion
+      name="regionSlave"
+      meshBody="meshSlave"
+      cellBlocks="{ cb2 }"
+      materialList="{ steel }"/>
+
+    <SurfaceElementRegion
+      name="surfaceMaster"
+      meshBody="meshMaster"
+      faceBlock="faceElementSubRegion"
+      defaultAperture="1.0e-6"
+      materialList="{ frictionLaw }"/>
+
+    <SurfaceElementRegion
+      name="surfaceSlave"
+      meshBody="meshSlave"
+      faceBlock="faceElementSubRegion"
+      defaultAperture="1.0e-6"
+      materialList="{ frictionLaw }"/>
+  </ElementRegions>
+
+  <Constitutive>
+    <ElasticIsotropic
+      name="steel"
+      defaultDensity="0.0"
+      defaultBulkModulus="175.0"
+      defaultShearModulus="80.77"/>
+
+    <Coulomb
+      name="frictionLaw"
+      cohesion="0.0"
+      frictionCoefficient="0.577350269"/>
+  </Constitutive>
+
+  <FieldSpecifications>
+    <FieldSpecification
+      name="surf"
+      initialCondition="1"
+      setNames="{ surface }"
+      objectPath="faceManager"
+      fieldName="ruptureState"
+      scale="1"/>
+
+    <FieldSpecification
+      name="x_fix_left"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="0"
+      scale="0.0"
+      setNames="{ fixed_left}"/>
+
+    <FieldSpecification
+      name="y_fix_left"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="1"
+      scale="0.0"
+      setNames="{ fixed_left}"/>
+
+    <FieldSpecification
+      name="z_fix_left"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="2"
+      scale="0.0"
+      setNames="{ fixed_left }"/>
+
+    <FieldSpecification
+      name="x_fix_right"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="0"
+      scale="0.0"
+      setNames="{ fixed_right}"/>
+
+    <FieldSpecification
+      name="y_fix_right"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="1"
+      scale="0.0"
+      setNames="{ fixed_right}"/>
+
+    <FieldSpecification
+      name="z_fix_right"
+      objectPath="nodeManager"
+      fieldName="totalDisplacement"
+      component="2"
+      scale="0.0"
+      setNames="{ fixed_right}"/>
+
+    <Traction
+      name="distributed_load_left"
+      objectPath="faceManager"
+      scale="1.0e-2"
+      direction="{ 0, 1, 0 }"
+      functionName="loadFunction"
+      setNames="{ top_surface_left}"/>
+
+    <Traction
+      name="distributed_load_right"
+      objectPath="faceManager"
+      scale="1.0e-2"
+      direction="{ 0, 1, 0 }"
+      functionName="loadFunction"
+      setNames="{ top_surface_right}"/>
+
+  </FieldSpecifications>
+
+  <Functions>
+    <TableFunction
+      name="loadFunction"
+      inputVarNames="{ time }"
+      coordinates="{ 0.0, 1.0 }"
+      values="{ -1.0, -1.0 }"/>
+  </Functions>
+
+  <Outputs>
+    <VTK
+      name="vtkOutput_2beams"
+      plotFileRoot="VTKplot"/>
+
+  </Outputs>
+
+  <Solvers
+    gravityVector="{ 0.0, 0.0, 0.0 }">
+    <SurfaceGenerator
+      name="SurfaceGenMaster"
+      logLevel="0"
+      fractureRegion="surfaceMaster"
+      targetRegions="{ meshMaster/regionMaster}"
+      initialRockToughness="1.0e6"
+      mpiCommOrder="1"/>
+
+    <SurfaceGenerator
+      name="SurfaceGenSlave"
+      logLevel="0"
+      fractureRegion="surfaceSlave"
+      targetRegions="{ meshSlave/regionSlave}"
+      initialRockToughness="1.0e6"
+      mpiCommOrder="1"/>
+
+    <SolidMechanicsMortarContact
+      name="Mortar"
+      timeIntegrationOption="QuasiStatic"
+      logLevel="2"
+      discretization="FE1"
+      targetRegions="{ meshMaster/regionMaster, meshMaster/surfaceMaster, meshSlave/regionSlave, meshSlave/surfaceSlave  }">
+      <NonlinearSolverParameters
+        newtonTol="5.0e-6"
+        logLevel="2"
+        newtonMaxIter="20"
+        lineSearchAction="None"
+        lineSearchMaxCuts="4"/>
+      <LinearSolverParameters
+        solverType="direct"
+       logLevel="1"/> 
+    </SolidMechanicsMortarContact>
+
+  </Solvers>
+
+  <NumericalMethods>
+    <FiniteElements>
+      <FiniteElementSpace
+        name="FE1"
+        order="1"/>
+    </FiniteElements>
+  </NumericalMethods>
+
+  <Mesh>
+    <InternalMesh
+      name="meshMaster"
+      elementTypes="{ C3D8 }"
+      xCoords="{ 0, 8 }"
+      yCoords="{ 0, 1 }"
+      zCoords="{ 0, 1 }"
+      nx="{ 32 }"
+      ny="{ 4 }"
+      nz="{ 4 }"
+      cellBlockNames="{ cb1 }"/>
+
+    <InternalMesh
+      name="meshSlave"
+      elementTypes="{ C3D8 }"
+      xCoords="{ 8, 16 }"
+      yCoords="{ 0, 1 }"
+      zCoords="{ 0, 1 }"
+      nx="{ 48 }"
+      ny="{ 6 }"
+      nz="{ 6 }"
+      cellBlockNames="{ cb2 }"/>
+
+  </Mesh>
+
+  <Geometry>
+    <Box
+      name = "fixed_left"
+      xMin = "{-0.01,-0.01,-0.01}"
+      xMax = "{0.01,1.01,1.01}"/>
+    <Box
+      name = "fixed_right"
+      xMin = "{15.99,-0.01,-0.01}"
+      xMax = "{16.01,1.01,1.01}"/>
+
+    <Box
+      name = "top_surface_left"
+      xMin = "{-0.01,0.99,-0.01}"
+      xMax = "{8.01,1.01,1.01}"/>
+    <Box
+      name = "top_surface_right"
+      xMin = "{7.99,0.99,-0.01}"
+      xMax = "{16.01,1.01,1.01}"/>
+
+    <Rectangle
+      name="surface"
+      normal="{1.0, 0.0, 0.0}"
+      origin="{8.0, 0.5, 0.5}"
+      lengthVector="{0.0, 1.0, 0.0}"
+      widthVector="{0.0, 0.0, 1.0}"
+      dimensions="{ 1.0, 1.0 }"/>
+
+  </Geometry>
+
+ </Problem>

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)
 

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 )
+     contact/SolidMechanicsMortarContact.cpp
+     PARENT_SCOPE )

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,

src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp

@@ -1,3 +1,4 @@
+
 /*
  * ------------------------------------------------------------------------------------------------------------
  * SPDX-License-Identifier: LGPL-2.1-only