Improvements to EFEM solver.

host-configs/LLNL/quartz-base.cmake

@@ -31,3 +31,4 @@ set(ENABLE_PAMELA ON CACHE BOOL "")
 set(ENABLE_PVTPackage ON CACHE BOOL "")
 set(ENABLE_GEOSX_PTP ON CACHE BOOL "" FORCE)
 set(ENABLE_PETSC OFF CACHE BOOL "Enables PETSc." FORCE)
+

src/coreComponents/constitutive/solid/LinearElasticIsotropic.hpp

@@ -106,30 +106,6 @@ class LinearElasticIsotropicUpdates : public SolidBaseUpdates
     c[5][5] = G;
   }
 
-  void GetStiffness( localIndex const k, array2d< real64 > & c ) const
-  {
-    real64 const G = m_shearModulus[k];
-    real64 const Lame = m_bulkModulus[k] - 2.0/3.0 * G;
-
-    c[0][0] = Lame + 2 * G;
-    c[0][1] = Lame;
-    c[0][2] = Lame;
-
-    c[1][0] = Lame;
-    c[1][1] = Lame + 2 * G;
-    c[1][2] = Lame;
-
-    c[2][0] = Lame;
-    c[2][1] = Lame;
-    c[2][2] = Lame + 2 * G;
-
-    c[3][3] = G;
-
-    c[4][4] = G;
-
-    c[5][5] = G;
-  }
-
   GEOSX_FORCE_INLINE
   GEOSX_HOST_DEVICE
   void setDiscretizationOps( localIndex const k,

src/coreComponents/fileIO/schema/docs/SolidMechanicsEmbeddedFractures.rst

@@ -5,6 +5,7 @@ Name                      Type         Default  Description
 ========================= ============ ======== ====================================================================================================================================================================================================================================================================================================================== 
 cflFactor                 real64       0.5      Factor to apply to the `CFL condition <http://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition>`_ when calculating the maximum allowable time step. Values should be in the interval (0,1]                                                                                                                      
 contactRelationName       string       required Name of contact relation to enforce constraints on fracture boundary.                                                                                                                                                                                                                                                  
+fractureRegionName        string       required Name of the fracture region.                                                                                                                                                                                                                                                                                           
 initialDt                 real64       1e+99    Initial time-step value required by the solver to the event manager.                                                                                                                                                                                                                                                   
 logLevel                  integer      0        Log level                                                                                                                                                                                                                                                                                                              
 name                      string       required A name is required for any non-unique nodes                                                                                                                                                                                                                                                                            

src/coreComponents/fileIO/schema/schema.xsd

@@ -1359,6 +1359,8 @@ the relative residual norm satisfies:
 		<xsd:attribute name="cflFactor" type="real64" default="0.5" />
 		<!--contactRelationName => Name of contact relation to enforce constraints on fracture boundary.-->
 		<xsd:attribute name="contactRelationName" type="string" use="required" />
+		<!--fractureRegionName => Name of the fracture region.-->
+		<xsd:attribute name="fractureRegionName" type="string" use="required" />
 		<!--initialDt => Initial time-step value required by the solver to the event manager.-->
 		<xsd:attribute name="initialDt" type="real64" default="1e+99" />
 		<!--logLevel => Log level-->

src/coreComponents/linearAlgebra/DofManager.cpp

@@ -1790,9 +1790,10 @@ void DofManager::makeRestrictor( std::vector< SubComponent > const & selection,
   GEOSX_ERROR_IF( !m_reordered, "Cannot make restrictors before reorderByRank() has been called." );
 
   // 1. Populate selected fields and compute some basic dimensions
-  array1d< FieldDescription > fieldsSelected( selection.size() );
+  // array1d< FieldDescription > fieldsSelected( selection.size() );
+  std::vector< FieldDescription > fieldsSelected( selection.size() );
 
-  for( localIndex k = 0; k < fieldsSelected.size(); ++k )
+  for( std::size_t k = 0; k < fieldsSelected.size(); ++k )
   {
     SubComponent const & dof = selection[k];
     FieldDescription const & fieldOld = m_fields[getFieldIndex( dof.fieldName )];
@@ -1812,7 +1813,7 @@ void DofManager::makeRestrictor( std::vector< SubComponent > const & selection,
   // 2. Compute remaining offsets (we mostly just need globalOffset, but it depends on others)
 
   globalIndex blockOffset = 0;
-  for( localIndex k = 0; k < fieldsSelected.size(); ++k )
+  for( std::size_t k = 0; k < fieldsSelected.size(); ++k )
   {
     fieldsSelected[k].blockOffset = blockOffset;
     blockOffset += fieldsSelected[k].numGlobalDof;
@@ -1822,7 +1823,7 @@ void DofManager::makeRestrictor( std::vector< SubComponent > const & selection,
                                               []( localIndex const n, FieldDescription const & f )
   { return n + f.rankOffset; } );
 
-  for( localIndex k = 0; k < fieldsSelected.size(); ++k )
+  for( std::size_t k = 0; k < fieldsSelected.size(); ++k )
   {
     fieldsSelected[k].globalOffset = globalOffset;
     globalOffset += fieldsSelected[k].numLocalDof;
@@ -1840,7 +1841,7 @@ void DofManager::makeRestrictor( std::vector< SubComponent > const & selection,
   restrictor.createWithLocalSize( rowSize, colSize, 1, comm );
   restrictor.open();
 
-  for( localIndex k = 0; k < fieldsSelected.size(); ++k )
+  for( std::size_t k = 0; k < fieldsSelected.size(); ++k )
   {
     FieldDescription const & fieldNew = fieldsSelected[k];
     FieldDescription const & fieldOld = m_fields[getFieldIndex( fieldNew.name )];

src/coreComponents/mesh/CellElementSubRegion.cpp

@@ -33,6 +33,8 @@ CellElementSubRegion::CellElementSubRegion( string const & name, Group * const p
   registerWrapper( viewKeyStruct::dNdXString, &m_dNdX )->setSizedFromParent( 1 )->reference().resizeDimension< 3 >( 3 );
 
   registerWrapper( viewKeyStruct::detJString, &m_detJ )->setSizedFromParent( 1 )->reference();
+
+  registerWrapper( viewKeyStruct::toEmbSurfString, &m_toEmbeddedSurfaces )->setSizedFromParent( 1 );
 }
 
 CellElementSubRegion::~CellElementSubRegion()
@@ -82,6 +84,16 @@ void CellElementSubRegion::ConstructSubRegionFromFaceSet( FaceManager const * co
   this->resize( targetSet.size() );
 }
 
+
+void CellElementSubRegion::addFracturedElement( localIndex const cellElemIndex,
+                                                localIndex const embSurfIndex )
+{
+  // add the connection between the element and the embedded surface to the map
+  m_toEmbeddedSurfaces.emplaceBack( cellElemIndex, embSurfIndex );
+  // add the element to the fractured elements list
+  m_fracturedCells.insert( cellElemIndex );
+}
+
 void CellElementSubRegion::ViewPackingExclusionList( SortedArray< localIndex > & exclusionList ) const
 {
   ObjectManagerBase::ViewPackingExclusionList( exclusionList );

src/coreComponents/mesh/CellElementSubRegion.hpp

@@ -30,6 +30,10 @@ namespace geosx
 class CellElementSubRegion : public CellBlock
 {
 public:
+
+  /// Type of map between cell blocks and embedded elements
+  using EmbSurfMapType = InterObjectRelation< ArrayOfArrays< localIndex > >;
+
   /**
    * @name Constructor / Destructor
    */
@@ -70,6 +74,14 @@ class CellElementSubRegion : public CellBlock
 
   ///@}
 
+  /**
+   * @brief Add fractured element to list and relative entries to the map.
+   * @param cellElemIndex cell element index
+   * @param embSurfIndex embedded surface element index
+   */
+  void addFracturedElement( localIndex const cellElemIndex,
+                            localIndex const embSurfIndex );
+
   /**
    * @name Overriding packing / Unpacking functions
    */
@@ -129,6 +141,8 @@ class CellElementSubRegion : public CellBlock
     static constexpr auto constitutiveGroupingString = "ConstitutiveGrouping";
     /// String key for the constitutive map
     static constexpr auto constitutiveMapString = "ConstitutiveMap";
+    /// String key to embSurfMap
+    static constexpr auto toEmbSurfString = "ToEmbeddedSurfaces";
 
     /// ViewKey for the constitutive grouping
     dataRepository::ViewKey constitutiveGrouping  = { constitutiveGroupingString };
@@ -165,6 +179,28 @@ class CellElementSubRegion : public CellBlock
   arrayView2d< real64 const > detJ() const
   { return m_detJ; }
 
+  /**
+   * @brief @return The sorted array of fractured elements.
+   */
+  SortedArray< localIndex > & fracturedElementsList()
+  { return m_fracturedCells; }
+
+  /**
+   * @brief @return The sorted array view of fractured elements.
+   */
+  SortedArrayView< localIndex const > const fracturedElementsList() const
+  { return m_fracturedCells.toViewConst(); }
+
+  /**
+   * @brief @return The map to the embedded surfaces
+   */
+  EmbSurfMapType & embeddedSurfacesList() { return m_toEmbeddedSurfaces; }
+
+  /**
+   * @brief @return The map to the embedded surfaces
+   */
+  EmbSurfMapType const & embeddedSurfacesList() const { return m_toEmbeddedSurfaces; }
+
   /// Map used for constitutive grouping
   map< string, localIndex_array > m_constitutiveGrouping;
 
@@ -185,6 +221,12 @@ class CellElementSubRegion : public CellBlock
   /// Map of unmapped global indices in the element-to-face map
   map< localIndex, array1d< globalIndex > > m_unmappedGlobalIndicesInFacelist;
 
+  /// List of fractured elements
+  SortedArray< localIndex > m_fracturedCells;
+
+  /// Map from Cell Elements to Embedded Surfaces
+  EmbSurfMapType m_toEmbeddedSurfaces;
+
   /**
    * @brief Pack element-to-node and element-to-face maps
    * @tparam the flag for the bufferOps::Pack function

src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp

@@ -241,16 +241,4 @@ void EmbeddedSurfaceSubRegion::setupRelatedObjectsInRelations( MeshLevel const *
   this->m_toNodesRelation.SetRelatedObject( mesh->getNodeManager() );
 }
 
-real64 EmbeddedSurfaceSubRegion::ComputeHeavisideFunction( ArraySlice< real64 const, 1, nodes::REFERENCE_POSITION_USD - 1 > const nodeCoord,
-                                                           localIndex const k ) const
-{
-  real64 heaviside, distanceVector[3];
-  LvArray::tensorOps::copy< 3 >( distanceVector, nodeCoord );
-  LvArray::tensorOps::subtract< 3 >( distanceVector, m_elementCenter[k] );
-
-  heaviside = LvArray::tensorOps::AiBi< 3 >( distanceVector, m_normalVector[k] ) > 0 ? 1 : 0;
-
-  return heaviside;
-}
-
 } /* namespace geosx */

src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp

@@ -154,7 +154,22 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion
 
     /// Connectivity index string
     static constexpr auto connectivityIndexString     = "connectivityIndex";
-  };
+
+    /// Displacement jump string
+    static constexpr auto dispJumpString             = "displacementJump";
+
+    /// Delta displacement jump string
+    static constexpr auto deltaDispJumpString        = "deltaDisplacementJump";
+
+    /// Displacement jump key
+    dataRepository::ViewKey dispJump                  = {dispJumpString};
+
+    /// Delta displacement jump key
+    dataRepository::ViewKey deltaDispJump             = {deltaDispJumpString};
+
+  }
+  /// viewKey struct for the EmbeddedSurfaceSubRegion class
+  viewKeys;
 
   virtual void setupRelatedObjectsInRelations( MeshLevel const * const mesh ) override;
 
@@ -259,8 +274,40 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion
    */
   array1d< real64 > const & getConnectivityIndex() const { return m_connectivityIndex;}
 
-  ///@}
 
+  /**
+   * @brief Get a mutable total displacement array.
+   * @return the total displacement array if it exists, or an error is thrown if it does not exist
+   * @note An error is thrown if the displacement jump does not exist
+   */
+  array2d< real64 > & displacementJump()
+  { return getReference< array2d< real64 > >( viewKeys.dispJump ); }
+
+  /**
+   * @brief Provide an immutable arrayView to the total displacement array.
+   * @return immutable arrayView of the total displacement array if it exists, or an error is thrown if it does not exist
+   * @note An error is thrown if the displacement jump does not exist
+   */
+  arrayView2d< real64 const > displacementJump() const
+  {return getReference< array2d< real64 > >( viewKeys.dispJump ); }
+
+  /**
+   * @brief Get a mutable incremental displacement array.
+   * @return the incremental displacement array if it exists, or an error is thrown if it does not exist
+   * @note An error is thrown if the incremental displacement jump does not exist
+   */
+  array2d< real64 > & incrementalDisplacementJump()
+  { return getReference< array2d< real64 > >( viewKeys.deltaDispJump ); }
+
+  /**
+   * @brief Provide an immutable arrayView to the incremental displacement jump array.
+   * @return immutable arrayView of the incremental displacement array if it exists, or an error is thrown if it does not exist
+   * @note An error is thrown if the incremental displacement jump does not exist
+   */
+  arrayView2d< real64 const > incrementalDisplacementJump() const
+  { return getReference< array2d< real64 > >( viewKeys.deltaDispJump ); }
+
+  ///@}
 
 private:
 

src/coreComponents/physicsSolvers/CMakeLists.txt

@@ -45,6 +45,8 @@ set( physicsSolvers_headers
      solidMechanics/SolidMechanicsLagrangianFEM.hpp
      solidMechanics/SolidMechanicsLagrangianSSLE.hpp
      solidMechanics/SolidMechanicsLagrangianFEMKernels.hpp
+     solidMechanics/SolidMechanicsEFEMKernels.hpp
+     solidMechanics/SolidMechanicsEFEMKernelsHelper.hpp
      solidMechanics/SolidMechanicsPoroElasticKernel.hpp
      solidMechanics/SolidMechanicsSmallStrainQuasiStaticKernel.hpp
      solidMechanics/SolidMechanicsSmallStrainImplicitNewmarkKernel.hpp