From ac45ee2a112dca04a2ecddacfc8252d2a486c29e Mon Sep 17 00:00:00 2001
From: Matteo Cusini <49037133+CusiniM@users.noreply.github.com>
Date: Tue, 28 Jan 2025 22:10:29 -0800
Subject: [PATCH 1/5] feat: MpiWrapper::allReduce overload for arrays. (#3446)

* feat: MpiWrapper::allReduce overload for arrays.
---
 host-configs/apple/macOS_base.cmake           |   2 +-
 src/coreComponents/common/CMakeLists.txt      |   1 +
 src/coreComponents/common/MpiWrapper.hpp      |  88 +++++++++++--
 src/coreComponents/common/TypesHelpers.hpp    | 124 ++++++++++++++++++
 .../common/initializeEnvironment.cpp          |   2 +-
 .../fileIO/timeHistory/HDFHistoryIO.cpp       |   5 +-
 .../TwoPointFluxApproximation.cpp             |   7 +-
 src/coreComponents/mesh/ParticleManager.cpp   |   8 +-
 .../mesh/generators/InternalMeshGenerator.cpp |   7 +-
 .../mesh/generators/VTKUtilities.cpp          |   3 +-
 .../PhysicsSolverBaseKernels.hpp              |  24 ++--
 .../contact/ContactSolverBase.cpp             |   7 +-
 .../SolidMechanicsEmbeddedFractures.cpp       |   9 +-
 .../fluidFlow/FlowSolverBase.cpp              |  28 ++--
 .../multiphysics/HydrofractureSolver.cpp      |   8 +-
 .../SolidMechanicsLagrangianFEM.cpp           |   8 +-
 .../solidMechanics/SolidMechanicsMPM.cpp      |  18 +--
 .../SolidMechanicsStatistics.cpp              |  14 +-
 .../surfaceGeneration/SurfaceGenerator.cpp    |   9 +-
 .../isotropic/ElasticWaveEquationSEM.cpp      |   6 +-
 20 files changed, 270 insertions(+), 108 deletions(-)
 create mode 100644 src/coreComponents/common/TypesHelpers.hpp

diff --git a/host-configs/apple/macOS_base.cmake b/host-configs/apple/macOS_base.cmake
index ec74d15b9c1..6d56f70755c 100644
--- a/host-configs/apple/macOS_base.cmake
+++ b/host-configs/apple/macOS_base.cmake
@@ -25,7 +25,7 @@ set(ENABLE_CALIPER "OFF" CACHE PATH "" FORCE )
 set( BLAS_LIBRARIES ${HOMEBREW_DIR}/opt/lapack/lib/libblas.dylib CACHE PATH "" FORCE )
 set( LAPACK_LIBRARIES ${HOMEBREW_DIR}/opt/lapack/lib/liblapack.dylib CACHE PATH "" FORCE )
 
-set(ENABLE_DOXYGEN ON CACHE BOOL "" FORCE)
+set(ENABLE_DOXYGEN OFF CACHE BOOL "" FORCE)
 set(ENABLE_SPHINX ON CACHE BOOL "" FORCE)
 set(ENABLE_MATHPRESSO ON CACHE BOOL "" FORCE )
 
diff --git a/src/coreComponents/common/CMakeLists.txt b/src/coreComponents/common/CMakeLists.txt
index 830be7b7b1f..9a4de42a717 100644
--- a/src/coreComponents/common/CMakeLists.txt
+++ b/src/coreComponents/common/CMakeLists.txt
@@ -48,6 +48,7 @@ set( common_headers
      Tensor.hpp
      TimingMacros.hpp
      TypeDispatch.hpp
+     TypesHelpers.hpp
      initializeEnvironment.hpp
      LifoStorage.hpp
      LifoStorageCommon.hpp
diff --git a/src/coreComponents/common/MpiWrapper.hpp b/src/coreComponents/common/MpiWrapper.hpp
index 4d017ef8393..de529215ea0 100644
--- a/src/coreComponents/common/MpiWrapper.hpp
+++ b/src/coreComponents/common/MpiWrapper.hpp
@@ -22,6 +22,7 @@
 
 #include "common/DataTypes.hpp"
 #include "common/Span.hpp"
+#include "common/TypesHelpers.hpp"
 
 #if defined(GEOS_USE_MPI)
   #include <mpi.h>
@@ -128,6 +129,8 @@ struct MpiWrapper
     Min,  //!< Min
     Sum,  //!< Sum
     Prod, //!< Prod
+    LogicalAnd, //!< Logical and
+    LogicalOr, //!< Logical or
   };
 
   MpiWrapper() = delete;
@@ -351,18 +354,6 @@ struct MpiWrapper
                         array1d< T > & recvbuf,
                         MPI_Comm comm = MPI_COMM_GEOS );
 
-  /**
-   * @brief Strongly typed wrapper around MPI_Allreduce.
-   * @param[in] sendbuf The pointer to the sending buffer.
-   * @param[out] recvbuf The pointer to the receive buffer.
-   * @param[in] count The number of values to send/receive.
-   * @param[in] op The MPI_Op to perform.
-   * @param[in] comm The MPI_Comm over which the gather operates.
-   * @return The return value of the underlying call to MPI_Allreduce().
-   */
-  template< typename T >
-  static int allReduce( T const * sendbuf, T * recvbuf, int count, MPI_Op op, MPI_Comm comm = MPI_COMM_GEOS );
-
   /**
    * @brief Convenience wrapper for the MPI_Allreduce function.
    * @tparam T type of data to reduce. Must correspond to a valid MPI_Datatype.
@@ -385,6 +376,29 @@ struct MpiWrapper
   template< typename T >
   static void allReduce( Span< T const > src, Span< T > dst, Reduction const op, MPI_Comm comm = MPI_COMM_GEOS );
 
+  /**
+   * @brief Convenience wrapper for the MPI_Allreduce function. Version for arrays.
+   * @tparam T type of data to reduce. Must correspond to a valid MPI_Datatype.
+   * @param src[in] The values to send to the reduction.
+   * @param dst[out] The resulting values.
+   * @param op The Reduction enum to perform.
+   * @param comm The communicator.
+   */
+  template< typename SRC_CONTAINER_TYPE, typename DST_CONTAINER_TYPE >
+  static void allReduce( SRC_CONTAINER_TYPE const & src, DST_CONTAINER_TYPE & dst, Reduction const op, MPI_Comm const comm = MPI_COMM_GEOS );
+
+  /**
+   * @brief Convenience wrapper for the MPI_Allreduce function. Version for arrays.
+   * @tparam T type of data to reduce. Must correspond to a valid MPI_Datatype.
+   * @param src[in] The values to send to the reduction.
+   * @param dst[out] The resulting values.
+   * @param count The number of contiguos elements of the arrays to perform the reduction on (must be leq than the size).
+   * @param op The Reduction enum to perform.
+   * @param comm The communicator.
+   */
+  template< typename SRC_CONTAINER_TYPE, typename DST_CONTAINER_TYPE >
+  static void allReduce( SRC_CONTAINER_TYPE const & src, DST_CONTAINER_TYPE & dst, int const count, Reduction const op, MPI_Comm const comm );
+
 
   /**
    * @brief Strongly typed wrapper around MPI_Reduce.
@@ -639,6 +653,19 @@ struct MpiWrapper
    */
   template< typename T > static T maxValLoc( T localValueLocation, MPI_Comm comm = MPI_COMM_GEOS );
 
+private:
+
+  /**
+   * @brief Strongly typed wrapper around MPI_Allreduce.
+   * @param[in] sendbuf The pointer to the sending buffer.
+   * @param[out] recvbuf The pointer to the receive buffer.
+   * @param[in] count The number of values to send/receive.
+   * @param[in] op The MPI_Op to perform.
+   * @param[in] comm The MPI_Comm over which the gather operates.
+   * @return The return value of the underlying call to MPI_Allreduce().
+   */
+  template< typename T >
+  static int allReduce( T const * sendbuf, T * recvbuf, int count, MPI_Op op, MPI_Comm comm = MPI_COMM_GEOS );
 };
 
 namespace internal
@@ -701,6 +728,14 @@ inline MPI_Op MpiWrapper::getMpiOp( Reduction const op )
     {
       return MPI_PROD;
     }
+    case Reduction::LogicalAnd:
+    {
+      return MPI_LAND;
+    }
+    case Reduction::LogicalOr:
+    {
+      return MPI_LOR;
+    }
     default:
       GEOS_ERROR( "Unsupported reduction operation" );
       return MPI_NO_OP;
@@ -1113,6 +1148,35 @@ void MpiWrapper::allReduce( Span< T const > const src, Span< T > const dst, Redu
   allReduce( src.data(), dst.data(), LvArray::integerConversion< int >( src.size() ), getMpiOp( op ), comm );
 }
 
+template< typename SRC_CONTAINER_TYPE, typename DST_CONTAINER_TYPE >
+void MpiWrapper::allReduce( SRC_CONTAINER_TYPE const & src, DST_CONTAINER_TYPE & dst, int const count, Reduction const op, MPI_Comm const comm )
+{
+  static_assert( std::is_trivially_copyable< typename get_value_type< SRC_CONTAINER_TYPE >::type >::value,
+                 "The type in the source container must be trivially copyable." );
+  static_assert( std::is_trivially_copyable< typename get_value_type< DST_CONTAINER_TYPE >::type >::value,
+                 "The type in the destination container must be trivially copyable." );
+  static_assert( std::is_same< typename get_value_type< SRC_CONTAINER_TYPE >::type,
+                               typename get_value_type< DST_CONTAINER_TYPE >::type >::value,
+                 "Source and destination containers must have the same value type." );
+  GEOS_ASSERT_GE( src.size(), count );
+  GEOS_ASSERT_GE( dst.size(), count );
+  allReduce( src.data(), dst.data(), count, getMpiOp( op ), comm );
+}
+
+template< typename SRC_CONTAINER_TYPE, typename DST_CONTAINER_TYPE >
+void MpiWrapper::allReduce( SRC_CONTAINER_TYPE const & src, DST_CONTAINER_TYPE & dst, Reduction const op, MPI_Comm const comm )
+{
+  static_assert( std::is_trivially_copyable< typename get_value_type< SRC_CONTAINER_TYPE >::type >::value,
+                 "The type in the source container must be trivially copyable." );
+  static_assert( std::is_trivially_copyable< typename get_value_type< DST_CONTAINER_TYPE >::type >::value,
+                 "The type in the destination container must be trivially copyable." );
+  static_assert( std::is_same< typename get_value_type< SRC_CONTAINER_TYPE >::type,
+                               typename get_value_type< DST_CONTAINER_TYPE >::type >::value,
+                 "Source and destination containers must have the same value type." );
+  GEOS_ASSERT_EQ( src.size(), dst.size() );
+  allReduce( src.data(), dst.data(), LvArray::integerConversion< int >( src.size() ), getMpiOp( op ), comm );
+}
+
 template< typename T >
 T MpiWrapper::sum( T const & value, MPI_Comm comm )
 {
diff --git a/src/coreComponents/common/TypesHelpers.hpp b/src/coreComponents/common/TypesHelpers.hpp
new file mode 100644
index 00000000000..f0e514e35a1
--- /dev/null
+++ b/src/coreComponents/common/TypesHelpers.hpp
@@ -0,0 +1,124 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 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 TypesHelpers.hpp
+ *
+ */
+
+#ifndef TYPES_HELPERS_HPP
+#define TYPES_HELPERS_HPP
+
+#include <type_traits>
+
+namespace geos
+{
+
+namespace internal
+{
+/**
+ * @brief Trait to determine if a type defines a `value_type` member.
+ *
+ * This primary template defaults to `std::false_type`, indicating that
+ * the type `T` does not define a `value_type` member.
+ *
+ * @tparam T The type to check.
+ * @tparam void A SFINAE parameter used to specialize the trait.
+ */
+template< typename T, typename = void >
+struct has_value_type : std::false_type {};
+
+/**
+ * @brief Specialization of `has_value_type` for types with a `value_type` member.
+ *
+ * If the type `T` defines a `value_type` member, this specialization
+ * is used, which inherits from `std::true_type`.
+ *
+ * @tparam T The type to check.
+ */
+template< typename T >
+struct has_value_type< T, std::void_t< typename T::value_type > > : std::true_type {};
+
+/**
+ * @brief Trait to determine if a type defines a `ValueType` member.
+ *
+ * This primary template defaults to `std::false_type`, indicating that
+ * the type `T` does not define a `ValueType` member.
+ *
+ * @tparam T The type to check.
+ * @tparam void A SFINAE parameter used to specialize the trait.
+ */
+template< typename T, typename = void >
+struct has_ValueType : std::false_type {};
+
+/**
+ * @brief Specialization of `has_ValueType` for types with a `ValueType` member.
+ *
+ * If the type `T` defines a `ValueType` member, this specialization
+ * is used, which inherits from `std::true_type`.
+ *
+ * @tparam T The type to check.
+ */
+template< typename T >
+struct has_ValueType< T, std::void_t< typename T::ValueType > > : std::true_type {};
+
+}   // namespace internal
+
+/**
+ * @brief Trait to retrieve the `value_type` or `ValueType` of a type `T`.
+ *
+ * This primary template provides a static assertion error if `T` does not
+ * define either `value_type` or `ValueType`.
+ *
+ * @tparam T The type from which to extract the type alias.
+ * @tparam Enable A SFINAE parameter used for specialization.
+ */
+template< typename T, typename Enable = void >
+struct get_value_type
+{
+  static_assert( sizeof(T) == 0, "T must define either value_type or ValueType." );
+};
+
+/**
+ * @brief Specialization of `get_value_type` for types with a `value_type` member.
+ *
+ * If the type `T` defines a `value_type` member, this specialization
+ * retrieves it as the alias `type`.
+ *
+ * @tparam T The type from which to extract `value_type`.
+ */
+template< typename T >
+struct get_value_type< T, std::enable_if_t< internal::has_value_type< T >::value > >
+{
+  using type = typename T::value_type;
+};
+
+/**
+ * @brief Specialization of `get_value_type` for types with a `ValueType` member.
+ *
+ * If the type `T` does not define a `value_type` but defines a `ValueType`,
+ * this specialization retrieves it as the alias `type`.
+ *
+ * @tparam T The type from which to extract `ValueType`.
+ */
+template< typename T >
+struct get_value_type< T, std::enable_if_t< !internal::has_value_type< T >::value && internal::has_ValueType< T >::value > >
+{
+  using type = typename T::ValueType;
+};
+
+} // namespace geos
+
+#endif /* TYPES_HELPERS_HPP */
diff --git a/src/coreComponents/common/initializeEnvironment.cpp b/src/coreComponents/common/initializeEnvironment.cpp
index 832f4fa542f..64387060583 100644
--- a/src/coreComponents/common/initializeEnvironment.cpp
+++ b/src/coreComponents/common/initializeEnvironment.cpp
@@ -290,7 +290,7 @@ static void addUmpireHighWaterMarks()
     string allocatorNameMinChars = string( MAX_NAME_LENGTH, '\0' );
 
     // Make sure that each rank is looking at the same allocator.
-    MpiWrapper::allReduce( allocatorNameFixedSize.c_str(), &allocatorNameMinChars.front(), MAX_NAME_LENGTH, MPI_MIN, MPI_COMM_GEOS );
+    MpiWrapper::allReduce( allocatorNameFixedSize, allocatorNameMinChars, MpiWrapper::Reduction::Min, MPI_COMM_GEOS );
     if( allocatorNameFixedSize != allocatorNameMinChars )
     {
       GEOS_WARNING( "Not all ranks have an allocator named " << allocatorNameFixedSize << ", cannot compute high water mark." );
diff --git a/src/coreComponents/fileIO/timeHistory/HDFHistoryIO.cpp b/src/coreComponents/fileIO/timeHistory/HDFHistoryIO.cpp
index 24ae0ddddf6..fdfbf95c61d 100644
--- a/src/coreComponents/fileIO/timeHistory/HDFHistoryIO.cpp
+++ b/src/coreComponents/fileIO/timeHistory/HDFHistoryIO.cpp
@@ -258,8 +258,9 @@ void HDFHistoryIO::init( bool existsOkay )
 void HDFHistoryIO::write()
 {
   // check if the size has changed on any process in the primary comm
-  int anyChanged = false;
-  MpiWrapper::allReduce( &m_sizeChanged, &anyChanged, 1, MPI_LOR, m_comm );
+  int const anyChanged = MpiWrapper::allReduce( m_sizeChanged,
+                                                MpiWrapper::Reduction::LogicalOr,
+                                                m_comm );
   m_sizeChanged = anyChanged;
 
   // this will set the first dim large enough to hold all the rows we're about to write
diff --git a/src/coreComponents/finiteVolume/TwoPointFluxApproximation.cpp b/src/coreComponents/finiteVolume/TwoPointFluxApproximation.cpp
index 25c128770fa..665e1eed140 100644
--- a/src/coreComponents/finiteVolume/TwoPointFluxApproximation.cpp
+++ b/src/coreComponents/finiteVolume/TwoPointFluxApproximation.cpp
@@ -951,10 +951,9 @@ void TwoPointFluxApproximation::computeAquiferStencil( DomainPartition & domain,
     localSumFaceAreasView[aquiferIndex] += targetSetSumFaceAreas.get();
   } );
 
-  MpiWrapper::allReduce( localSumFaceAreas.data(),
-                         globalSumFaceAreas.data(),
-                         localSumFaceAreas.size(),
-                         MpiWrapper::getMpiOp( MpiWrapper::Reduction::Sum ),
+  MpiWrapper::allReduce( localSumFaceAreas,
+                         globalSumFaceAreas,
+                         MpiWrapper::Reduction::Sum,
                          MPI_COMM_GEOS );
 
   // Step 3: compute the face area fraction for each connection, and insert into boundary stencil
diff --git a/src/coreComponents/mesh/ParticleManager.cpp b/src/coreComponents/mesh/ParticleManager.cpp
index 057533e4f7e..57bd1687935 100644
--- a/src/coreComponents/mesh/ParticleManager.cpp
+++ b/src/coreComponents/mesh/ParticleManager.cpp
@@ -77,11 +77,9 @@ void ParticleManager::setMaxGlobalIndex()
     m_localMaxGlobalIndex = std::max( m_localMaxGlobalIndex, subRegion.maxGlobalIndex() );
   } );
 
-  MpiWrapper::allReduce( &m_localMaxGlobalIndex,
-                         &m_maxGlobalIndex,
-                         1,
-                         MPI_MAX,
-                         MPI_COMM_GEOS );
+  m_maxGlobalIndex = MpiWrapper::allReduce( m_localMaxGlobalIndex,
+                                            MpiWrapper::Reduction::Max,
+                                            MPI_COMM_GEOS );
 }
 
 Group * ParticleManager::createChild( string const & childKey, string const & childName )
diff --git a/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp b/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp
index ace6dd8846e..99fe38fb059 100644
--- a/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp
+++ b/src/coreComponents/mesh/generators/InternalMeshGenerator.cpp
@@ -603,10 +603,9 @@ void InternalMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockMa
     {
       elemCenterCoordsLocal[k] = m_min[dim] + ( m_max[dim] - m_min[dim] ) * ( k + 0.5 ) / m_numElemsTotal[dim];
     }
-    MpiWrapper::allReduce( elemCenterCoordsLocal.data(),
-                           elemCenterCoords[dim].data(),
-                           m_numElemsTotal[dim],
-                           MPI_MAX,
+    MpiWrapper::allReduce( elemCenterCoordsLocal,
+                           elemCenterCoords[dim],
+                           MpiWrapper::Reduction::Max,
                            MPI_COMM_GEOS );
   }
 
diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp
index a18b3e6767a..1fb198d70ce 100644
--- a/src/coreComponents/mesh/generators/VTKUtilities.cpp
+++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp
@@ -746,8 +746,7 @@ vtkSmartPointer< vtkDataSet > manageGlobalIds( vtkSmartPointer< vtkDataSet > mes
   {
     // Add global ids on the fly if needed
     int const me = hasGlobalIds( mesh );
-    int everyone;
-    MpiWrapper::allReduce( &me, &everyone, 1, MPI_MAX, MPI_COMM_GEOS );
+    int const everyone = MpiWrapper::allReduce( me, MpiWrapper::Reduction::Max, MPI_COMM_GEOS );
 
     if( everyone and not me )
     {
diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBaseKernels.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBaseKernels.hpp
index ac8433b208b..45e2ade8570 100644
--- a/src/coreComponents/physicsSolvers/PhysicsSolverBaseKernels.hpp
+++ b/src/coreComponents/physicsSolvers/PhysicsSolverBaseKernels.hpp
@@ -267,10 +267,9 @@ class LinfResidualNormHelper
   static void computeGlobalNorm( array1d< real64 > const & localResidualNorm,
                                  array1d< real64 > & globalResidualNorm )
   {
-    MpiWrapper::allReduce( localResidualNorm.data(),
-                           globalResidualNorm.data(),
-                           localResidualNorm.size(),
-                           MpiWrapper::getMpiOp( MpiWrapper::Reduction::Max ),
+    MpiWrapper::allReduce( localResidualNorm,
+                           globalResidualNorm,
+                           MpiWrapper::Reduction::Max,
                            MPI_COMM_GEOS );
   }
 };
@@ -309,16 +308,17 @@ class L2ResidualNormHelper
   {
     array1d< real64 > sumLocalResidualNorm( localResidualNorm.size() );
     array1d< real64 > sumLocalResidualNormalizer( localResidualNormalizer.size() );
-    MpiWrapper::allReduce( localResidualNorm.data(),
-                           sumLocalResidualNorm.data(),
-                           localResidualNorm.size(),
-                           MpiWrapper::getMpiOp( MpiWrapper::Reduction::Sum ),
+
+    MpiWrapper::allReduce( localResidualNorm,
+                           sumLocalResidualNorm,
+                           MpiWrapper::Reduction::Sum,
                            MPI_COMM_GEOS );
-    MpiWrapper::allReduce( localResidualNormalizer.data(),
-                           sumLocalResidualNormalizer.data(),
-                           localResidualNormalizer.size(),
-                           MpiWrapper::getMpiOp( MpiWrapper::Reduction::Sum ),
+
+    MpiWrapper::allReduce( localResidualNormalizer,
+                           sumLocalResidualNormalizer,
+                           MpiWrapper::Reduction::Sum,
                            MPI_COMM_GEOS );
+
     for( integer i = 0; i < localResidualNorm.size(); ++i )
     {
       globalResidualNorm[i] = sqrt( sumLocalResidualNorm[i] ) / sqrt( sumLocalResidualNormalizer[i] );
diff --git a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp
index edd3bd4ec24..9555658db65 100644
--- a/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp
+++ b/src/coreComponents/physicsSolvers/contact/ContactSolverBase.cpp
@@ -184,10 +184,9 @@ void ContactSolverBase::computeFractureStateStatistics( MeshLevel const & mesh,
 
   array1d< globalIndex > totalCounter( 4 );
 
-  MpiWrapper::allReduce( localCounter.data(),
-                         totalCounter.data(),
-                         4,
-                         MPI_SUM,
+  MpiWrapper::allReduce( localCounter,
+                         totalCounter,
+                         MpiWrapper::Reduction::Sum,
                          MPI_COMM_GEOS );
 
   numStick    = totalCounter[0];
diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp
index b78afce691a..4e31c7c4eaa 100644
--- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp
+++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp
@@ -841,12 +841,9 @@ bool SolidMechanicsEmbeddedFractures::updateConfiguration( DomainPartition & dom
   synchronizeFractureState( domain );
 
   // Compute if globally the fracture state has changed
-  int hasConfigurationConvergedGlobally;
-  MpiWrapper::allReduce( &hasConfigurationConverged,
-                         &hasConfigurationConvergedGlobally,
-                         1,
-                         MPI_LAND,
-                         MPI_COMM_GEOS );
+  int const hasConfigurationConvergedGlobally = MpiWrapper::allReduce( hasConfigurationConverged,
+                                                                       MpiWrapper::Reduction::LogicalAnd,
+                                                                       MPI_COMM_GEOS );
 
   // for this solver it makes sense to reset the state.
   // if( !hasConfigurationConvergedGlobally )
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
index e9a55b86c46..35c1a9a6b40 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
@@ -671,15 +671,13 @@ void FlowSolverBase::findMinMaxElevationInEquilibriumTarget( DomainPartition & d
 
   } );
 
-  MpiWrapper::allReduce( localMaxElevation.data(),
-                         maxElevation.data(),
-                         localMaxElevation.size(),
-                         MpiWrapper::getMpiOp( MpiWrapper::Reduction::Max ),
+  MpiWrapper::allReduce( localMaxElevation.toView(),
+                         maxElevation,
+                         MpiWrapper::Reduction::Max,
                          MPI_COMM_GEOS );
-  MpiWrapper::allReduce( localMinElevation.data(),
-                         minElevation.data(),
-                         localMinElevation.size(),
-                         MpiWrapper::getMpiOp( MpiWrapper::Reduction::Min ),
+  MpiWrapper::allReduce( localMinElevation.toView(),
+                         minElevation,
+                         MpiWrapper::Reduction::Min,
                          MPI_COMM_GEOS );
 }
 
@@ -726,10 +724,9 @@ void FlowSolverBase::computeSourceFluxSizeScalingFactor( real64 const & time,
   } );
 
   // synchronize the set size over all the MPI ranks
-  MpiWrapper::allReduce( bcAllSetsSize.data(),
-                         bcAllSetsSize.data(),
-                         bcAllSetsSize.size(),
-                         MpiWrapper::getMpiOp( MpiWrapper::Reduction::Sum ),
+  MpiWrapper::allReduce( bcAllSetsSize,
+                         bcAllSetsSize,
+                         MpiWrapper::Reduction::Sum,
                          MPI_COMM_GEOS );
 }
 
@@ -809,10 +806,9 @@ void FlowSolverBase::saveAquiferConvergedState( real64 const & time,
     localSumFluxes[aquiferIndex] += targetSetSumFluxes;
   } );
 
-  MpiWrapper::allReduce( localSumFluxes.data(),
-                         globalSumFluxes.data(),
-                         localSumFluxes.size(),
-                         MpiWrapper::getMpiOp( MpiWrapper::Reduction::Sum ),
+  MpiWrapper::allReduce( localSumFluxes,
+                         globalSumFluxes,
+                         MpiWrapper::Reduction::Sum,
                          MPI_COMM_GEOS );
 
   // Step 3: we are ready to save the summed fluxes for each individual aquifer
diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp
index dbf00bbb3d7..3175f518666 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp
@@ -226,11 +226,9 @@ real64 HydrofractureSolver< POROMECHANICS_SOLVER >::fullyCoupledSolverStep( real
     {
       locallyFractured = 1;
     }
-    MpiWrapper::allReduce( &locallyFractured,
-                           &globallyFractured,
-                           1,
-                           MPI_MAX,
-                           MPI_COMM_GEOS );
+    globallyFractured = MpiWrapper::allReduce( locallyFractured,
+                                               MpiWrapper::Reduction::Max,
+                                               MPI_COMM_GEOS );
 
     if( globallyFractured == 0 )
     {
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
index fa5c1685db3..dd19b20d9b3 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
@@ -513,11 +513,9 @@ real64 SolidMechanicsLagrangianFEM::solverStep( real64 const & time_n,
         {
           locallyFractured = 1;
         }
-        MpiWrapper::allReduce( &locallyFractured,
-                               &globallyFractured,
-                               1,
-                               MPI_MAX,
-                               MPI_COMM_GEOS );
+        globallyFractured = MpiWrapper::allReduce( locallyFractured,
+                                                   MpiWrapper::Reduction::Max,
+                                                   MPI_COMM_GEOS );
       }
       if( globallyFractured == 0 )
       {
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
index a2b1d3ec2ff..c576f6513e2 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
@@ -3313,23 +3313,17 @@ void SolidMechanicsMPM::printProfilingResults()
 
   // Get total CPU time for the entire time step
   real64 totalStepTimeThisRank = m_profilingTimes[numIntervals] - m_profilingTimes[0];
-  real64 totalStepTimeAllRanks;
-  MpiWrapper::allReduce< real64 >( &totalStepTimeThisRank,
-                                   &totalStepTimeAllRanks,
-                                   1,
-                                   MPI_SUM,
-                                   MPI_COMM_GEOS );
+  real64 const totalStepTimeAllRanks = MpiWrapper::allReduce( totalStepTimeThisRank,
+                                                              MpiWrapper::Reduction::Sum,
+                                                              MPI_COMM_GEOS );
 
   // Get total CPU times for each queried time interval
   for( unsigned int i = 0; i < numIntervals; i++ )
   {
     real64 timeIntervalThisRank = ( m_profilingTimes[i+1] - m_profilingTimes[i] );
-    real64 timeIntervalAllRanks;
-    MpiWrapper::allReduce< real64 >( &timeIntervalThisRank,
-                                     &timeIntervalAllRanks,
-                                     1,
-                                     MPI_SUM,
-                                     MPI_COMM_GEOS );
+    real64 const timeIntervalAllRanks = MpiWrapper::allReduce( timeIntervalThisRank,
+                                                               MpiWrapper::Reduction::Sum,
+                                                               MPI_COMM_GEOS );
     if( rank == 0 )
     {
       timeIntervalsAllRanks[i] = timeIntervalAllRanks;
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStatistics.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStatistics.cpp
index d16dd1ea519..12f47f2b3e3 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStatistics.cpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStatistics.cpp
@@ -146,16 +146,14 @@ void SolidMechanicsStatistics::computeNodeStatistics( MeshLevel & mesh, real64 c
   nodeStatistics.minDisplacement[1] = minDispY.get();
   nodeStatistics.minDisplacement[2] = minDispZ.get();
 
-  MpiWrapper::allReduce( nodeStatistics.maxDisplacement.data(),
-                         nodeStatistics.maxDisplacement.data(),
-                         3,
-                         MpiWrapper::getMpiOp( MpiWrapper::Reduction::Max ),
+  MpiWrapper::allReduce( nodeStatistics.maxDisplacement,
+                         nodeStatistics.maxDisplacement,
+                         MpiWrapper::Reduction::Max,
                          MPI_COMM_GEOS );
 
-  MpiWrapper::allReduce( nodeStatistics.minDisplacement.data(),
-                         nodeStatistics.minDisplacement.data(),
-                         3,
-                         MpiWrapper::getMpiOp( MpiWrapper::Reduction::Min ),
+  MpiWrapper::allReduce( nodeStatistics.minDisplacement,
+                         nodeStatistics.minDisplacement,
+                         MpiWrapper::Reduction::Min,
                          MPI_COMM_GEOS );
 
   TableData mechanicsData;
diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
index 788c078ada3..e41674c9f82 100644
--- a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
+++ b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
@@ -4580,12 +4580,9 @@ SurfaceGenerator::calculateRuptureRate( SurfaceElementRegion & faceElementRegion
     maxRuptureRate = std::max( maxRuptureRate, ruptureRate( faceElemIndex ) );
   }
 
-  real64 globalMaxRuptureRate;
-  MpiWrapper::allReduce( &maxRuptureRate,
-                         &globalMaxRuptureRate,
-                         1,
-                         MPI_MAX,
-                         MPI_COMM_GEOS );
+  real64 const globalMaxRuptureRate = MpiWrapper::allReduce( maxRuptureRate,
+                                                             MpiWrapper::Reduction::Max,
+                                                             MPI_COMM_GEOS );
 
   return globalMaxRuptureRate;
 }
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp
index acf9572b216..84277473919 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/elastic/secondOrderEqn/isotropic/ElasticWaveEquationSEM.cpp
@@ -1013,10 +1013,10 @@ void ElasticWaveEquationSEM::cleanup( real64 const time_n,
       computeAllSeismoTraces( time_n, 0.0, uy_np1, uy_n, dasReceivers, m_linearDASVectorY.toView(), true );
       computeAllSeismoTraces( time_n, 0.0, uz_np1, uz_n, dasReceivers, m_linearDASVectorZ.toView(), true );
       // sum contributions from all MPI ranks, since some receivers might be split among multiple ranks
-      MpiWrapper::allReduce( dasReceivers.data(),
-                             dasReceivers.data(),
+      MpiWrapper::allReduce( dasReceivers,
+                             dasReceivers,
                              m_linearDASGeometry.size( 0 ),
-                             MpiWrapper::getMpiOp( MpiWrapper::Reduction::Sum ),
+                             MpiWrapper::Reduction::Sum,
                              MPI_COMM_GEOS );
       WaveSolverUtils::writeSeismoTrace( "dasTraceReceiver", getName(), m_outputSeismoTrace, m_linearDASGeometry.size( 0 ),
                                          m_receiverIsLocal, m_nsamplesSeismoTrace, dasReceivers );

From bd709034a07eb9cb4192cc7bbc420da30e147a24 Mon Sep 17 00:00:00 2001
From: Pavel Tomin <paveltomin@users.noreply.github.com>
Date: Wed, 29 Jan 2025 14:24:58 -0600
Subject: [PATCH 2/5] refactor: unify flow solver discretization name check
 (#3514)

---
 .../fluidFlow/CompositionalMultiphaseFVM.cpp         |  9 +--------
 .../physicsSolvers/fluidFlow/FlowSolverBase.cpp      | 12 ++++++++++++
 .../physicsSolvers/fluidFlow/FlowSolverBase.hpp      |  2 ++
 .../fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp |  9 +--------
 .../physicsSolvers/fluidFlow/SinglePhaseFVM.cpp      |  9 +--------
 5 files changed, 17 insertions(+), 24 deletions(-)

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
index aeb88ef68e3..67ee8bf2f9a 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp
@@ -165,14 +165,7 @@ void CompositionalMultiphaseFVM::initializePreSubGroups()
                                                 ? LinearSolverParameters::MGR::StrategyType::thermalCompositionalMultiphaseFVM
                                                 : LinearSolverParameters::MGR::StrategyType::compositionalMultiphaseFVM;
 
-  DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
-  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
-  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
-  if( !fvManager.hasGroup< FluxApproximationBase >( m_discretizationName ) )
-  {
-    GEOS_ERROR( "A discretization deriving from FluxApproximationBase must be selected with CompositionalMultiphaseFlow" );
-  }
-
+  checkDiscretizationName();
 }
 
 void CompositionalMultiphaseFVM::setupDofs( DomainPartition const & domain,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
index 35c1a9a6b40..fe3bad87c88 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
@@ -368,6 +368,18 @@ void FlowSolverBase::initializePreSubGroups()
   }
 }
 
+void FlowSolverBase::checkDiscretizationName() const
+{
+  DomainPartition const & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
+  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
+  if( !fvManager.hasGroup< FluxApproximationBase >( m_discretizationName ) )
+  {
+    GEOS_ERROR( GEOS_FMT( "{}: can not find discretization named '{}' (a discretization deriving from FluxApproximationBase must be selected for {} solver '{}' )",
+                          getDataContext(), m_discretizationName, getCatalogName(), getName()));
+  }
+}
+
 void FlowSolverBase::validatePoreVolumes( DomainPartition const & domain ) const
 {
   real64 minPoreVolume = LvArray::NumericLimits< real64 >::max;
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
index b97716d0b00..19e79d90ab2 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
@@ -230,6 +230,8 @@ class FlowSolverBase : public PhysicsSolverBase
 
   virtual void initializePreSubGroups() override;
 
+  void checkDiscretizationName() const;
+
   virtual void initializePostInitialConditionsPreSubGroups() override;
 
   void initializeState( DomainPartition & domain );
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp
index 07c38cd2726..4670d6425a4 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/ReactiveCompositionalMultiphaseOBL.cpp
@@ -125,14 +125,7 @@ void ReactiveCompositionalMultiphaseOBL::initializePreSubGroups()
 {
   FlowSolverBase::initializePreSubGroups();
 
-  DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
-  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
-  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
-  if( !fvManager.hasGroup< FluxApproximationBase >( m_discretizationName ) )
-  {
-    GEOS_ERROR( "A discretization deriving from FluxApproximationBase must be selected with ReactiveCompositionalMultiphaseOBL" );
-  }
-
+  checkDiscretizationName();
 }
 
 void ReactiveCompositionalMultiphaseOBL::setupDofs( DomainPartition const & domain,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp
index cc9cdc6aaf6..5da5846273c 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseFVM.cpp
@@ -68,14 +68,7 @@ void SinglePhaseFVM< BASE >::initializePreSubGroups()
 {
   BASE::initializePreSubGroups();
 
-  DomainPartition & domain = this->template getGroupByPath< DomainPartition >( "/Problem/domain" );
-  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
-  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
-
-  if( !fvManager.hasGroup< FluxApproximationBase >( m_discretizationName ) )
-  {
-    GEOS_ERROR( "A discretization deriving from FluxApproximationBase must be selected with SinglePhaseFVM" );
-  }
+  this->checkDiscretizationName();
 
   if( m_isThermal )
   {

From 9cfb92a9618f2e6f7224a43d07be16afce5e8767 Mon Sep 17 00:00:00 2001
From: Pavel Tomin <paveltomin@users.noreply.github.com>
Date: Tue, 4 Feb 2025 13:22:10 -0600
Subject: [PATCH 3/5] fix: issue 3499 (#3500)

---
 src/coreComponents/common/Units.hpp                           | 2 ++
 .../capillaryPressure/JFunctionCapillaryPressure.cpp          | 4 +++-
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/coreComponents/common/Units.hpp b/src/coreComponents/common/Units.hpp
index 8bfdc04ca80..693c3a13f75 100644
--- a/src/coreComponents/common/Units.hpp
+++ b/src/coreComponents/common/Units.hpp
@@ -30,6 +30,8 @@ namespace geos
 namespace units
 {
 
+/// Darcy to m^2 conversion factor
+static constexpr double DarcyToSqM = 9.869233e-13;
 
 /**
  * @return the input Kelvin degrees converted in Celsius
diff --git a/src/coreComponents/constitutive/capillaryPressure/JFunctionCapillaryPressure.cpp b/src/coreComponents/constitutive/capillaryPressure/JFunctionCapillaryPressure.cpp
index 908c42604b4..3e67db4b27b 100644
--- a/src/coreComponents/constitutive/capillaryPressure/JFunctionCapillaryPressure.cpp
+++ b/src/coreComponents/constitutive/capillaryPressure/JFunctionCapillaryPressure.cpp
@@ -22,6 +22,7 @@
 #include "constitutive/capillaryPressure/CapillaryPressureFields.hpp"
 #include "constitutive/capillaryPressure/TableCapillaryPressureHelpers.hpp"
 #include "functions/FunctionManager.hpp"
+#include "common/Units.hpp"
 
 namespace geos
 {
@@ -250,7 +251,8 @@ void JFunctionCapillaryPressure::saveConvergedRockState( arrayView2d< real64 con
     {
       permeability = convergedPermeability[ei][0][2];
     }
-    GEOS_ERROR_IF( permeability < LvArray::NumericLimits< real64 >::epsilon, "Zero permeability in J-function capillary pressure" );
+    // multiply epsilon by Darcy to sq m factor
+    GEOS_ERROR_IF( permeability < LvArray::NumericLimits< real64 >::epsilon * units::DarcyToSqM, "Zero permeability in J-function capillary pressure" );
 
     // here we compute an average of the porosity over quadrature points
     // this average is exact for tets, regular pyramids/wedges/hexes, or for VEM

From 0b8a7e25d54afd1140801db276a31f50e98730a8 Mon Sep 17 00:00:00 2001
From: acitrain <60715545+acitrain@users.noreply.github.com>
Date: Wed, 5 Feb 2025 11:01:35 +0100
Subject: [PATCH 4/5] fix: Add source value array for fwi in wave solvers
 (#3502)

This PR is adding an array to store the values of the sources in time to the acoustic wave solver.
This array is useful for FWI and pygeos.

---------

Co-authored-by: j0405284 <stefano.frambati@total.com>
Co-authored-by: BESSET Julien <julien.besset@inria.fr>
Co-authored-by: Stefano Frambati <stefano.frambati@totalenergies.com>
---
 .integrated_tests.yaml                        |   4 +-
 BASELINE_NOTES.md                             |   6 +-
 .../isotropic/AcousticWaveEquationSEM.cpp     |  85 ++++++++---
 .../isotropic/AcousticWaveEquationSEM.hpp     |   4 +-
 .../AcousticElasticWaveEquationSEM.cpp        |   4 +-
 .../PrecomputeSourcesAndReceiversKernel.hpp   |  23 ++-
 .../wavePropagation/shared/WaveSolverBase.cpp |   7 +
 .../wavePropagation/shared/WaveSolverBase.hpp |   4 +
 src/coreComponents/schema/schema.xsd          | 141 ++++++++++++------
 src/coreComponents/schema/schema.xsd.other    |  43 ++++--
 .../testWavePropagationAdjoint1.cpp           |   4 +-
 11 files changed, 244 insertions(+), 81 deletions(-)

diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
index ea7e4b5b516..35bfac096d1 100644
--- a/.integrated_tests.yaml
+++ b/.integrated_tests.yaml
@@ -1,7 +1,7 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3395-9832-31f145e
-  
+  baseline: integratedTests/baseline_integratedTests-pr3502-10013-5b4e015
+
 allow_fail:
   all: ''
   streak: ''
diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md
index fe597c2c6c8..cb85ce3d56a 100644
--- a/BASELINE_NOTES.md
+++ b/BASELINE_NOTES.md
@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3502 (2025-02-04)
+=====================
+Add array to store the source values in time inside wave solvers
+
 PR  #3395 (2024-01-22)
 =====================
 Add new fields and change the default input for some tests.
@@ -128,7 +132,7 @@ Add routine for automatic time steps in waveSolvers with new attributes
 
 PR #3156 (2024-10-29)
 ====================
-Restart check errors due to 1) schema node added to enable thermal option in well model and 2) arrays removed/added for option.  Max difference errors due treatment of shutin wells.  Previously non-zero rate value reported for shutin well, new code will set rate arrays to zero. 
+Restart check errors due to 1) schema node added to enable thermal option in well model and 2) arrays removed/added for option.  Max difference errors due treatment of shutin wells.  Previously non-zero rate value reported for shutin well, new code will set rate arrays to zero.
 
 PR #2878 (2024-10-17)
 =====================
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp
index c47e6e83fc1..1289cbd95b0 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.cpp
@@ -146,6 +146,35 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseM
   receiverConstants.setValues< EXEC_POLICY >( -1 );
   receiverIsLocal.zero();
 
+  arrayView1d< TableFunction::KernelWrapper const > const sourceWaveletTableWrappers = m_sourceWaveletTableWrappers.toViewConst();
+  bool useSourceWaveletTables = m_useSourceWaveletTables;
+
+  //Correct size for sourceValue
+  EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
+  real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() );
+  real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() );
+  real64 dt = 0;
+  for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent )
+  {
+    EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] );
+    if( subEvent->getEventName() == "/Solvers/" + this->getName() )
+    {
+      dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() );
+    }
+  }
+
+  real64 dtCompute;
+
+  localIndex nSubSteps = (int) ceil( dt/m_timeStep );
+  dtCompute = dt/nSubSteps;
+
+  localIndex const nsamples = int( (maxTime - minTime) / dtCompute) + 1;
+
+  localIndex const numSourcesGlobal = m_sourceCoordinates.size( 0 );
+  m_sourceValue.resize( nsamples, numSourcesGlobal );
+
+  arrayView2d< real32 > const sourceValue = m_sourceValue.toView();
+
   mesh.getElemManager().forElementSubRegionsComplete< CellElementSubRegion >( regionNames, [&]( localIndex const,
                                                                                                 localIndex const er,
                                                                                                 localIndex const esr,
@@ -192,7 +221,14 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseM
           receiverCoordinates,
           receiverIsLocal,
           receiverNodeIds,
-          receiverConstants );
+          receiverConstants,
+          sourceValue,
+          dtCompute,
+          m_timeSourceFrequency,
+          m_timeSourceDelay,
+          m_rickerOrder,
+          sourceWaveletTableWrappers,
+          useSourceWaveletTables );
       }
     } );
     elementSubRegion.faceList().freeOnDevice();
@@ -209,25 +245,24 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & baseM
   nodesToElements.freeOnDevice();
 }
 
-void AcousticWaveEquationSEM::addSourceToRightHandSide( real64 const & time_n, arrayView1d< real32 > const rhs )
+void AcousticWaveEquationSEM::addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs )
 {
   arrayView2d< localIndex const > const sourceNodeIds = m_sourceNodeIds.toViewConst();
   arrayView2d< real64 const > const sourceConstants   = m_sourceConstants.toViewConst();
   arrayView1d< localIndex const > const sourceIsAccessible = m_sourceIsAccessible.toViewConst();
-  real32 const timeSourceFrequency = m_timeSourceFrequency;
-  real32 const timeSourceDelay = m_timeSourceDelay;
-  localIndex const rickerOrder = m_rickerOrder;
-  bool useSourceWaveletTables = m_useSourceWaveletTables;
-  arrayView1d< TableFunction::KernelWrapper const > const sourceWaveletTableWrappers = m_sourceWaveletTableWrappers.toViewConst();
+  arrayView2d< real32 const > const sourceValue   = m_sourceValue.toViewConst();
+
+  GEOS_THROW_IF( cycleNumber > sourceValue.size( 0 ),
+                 getDataContext() << ": Too many steps compared to array size",
+                 std::runtime_error );
   forAll< EXEC_POLICY >( sourceConstants.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const isrc )
   {
     if( sourceIsAccessible[isrc] == 1 )
     {
-      real64 const srcValue =
-        useSourceWaveletTables ? sourceWaveletTableWrappers[ isrc ].compute( &time_n ) : WaveSolverUtils::evaluateRicker( time_n, timeSourceFrequency, timeSourceDelay, rickerOrder );
+
       for( localIndex inode = 0; inode < sourceConstants.size( 1 ); ++inode )
       {
-        real32 const localIncrement = sourceConstants[isrc][inode] * srcValue;
+        real32 const localIncrement = sourceConstants[isrc][inode] * sourceValue[cycleNumber][isrc];
         RAJA::atomicAdd< ATOMIC_POLICY >( &rhs[sourceNodeIds[isrc][inode]], localIncrement );
       }
     }
@@ -252,12 +287,10 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
 
 
 
-  forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName,
+  forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
                                                                 MeshLevel & mesh,
                                                                 arrayView1d< string const > const & regionNames )
   {
-    MeshLevel & baseMesh = domain.getMeshBodies().getGroup< MeshBody >( meshBodyName ).getBaseDiscretization();
-    precomputeSourceAndReceiverTerm( baseMesh, mesh, regionNames );
 
     NodeManager & nodeManager = mesh.getNodeManager();
     FaceManager & faceManager = mesh.getFaceManager();
@@ -338,6 +371,16 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
     m_timeStep=dtOut;
   }
 
+  forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName,
+                                                                MeshLevel & mesh,
+                                                                arrayView1d< string const > const & regionNames )
+  {
+    GEOS_UNUSED_VAR( meshBodyName );
+    MeshLevel & baseMesh = domain.getMeshBodies().getGroup< MeshBody >( meshBodyName ).getBaseDiscretization();
+    precomputeSourceAndReceiverTerm( baseMesh, mesh, regionNames );
+
+  } );
+
 
   WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ), m_receiverIsLocal );
 }
@@ -875,7 +918,7 @@ real64 AcousticWaveEquationSEM::explicitStepForward( real64 const & time_n,
                                                      DomainPartition & domain,
                                                      bool computeGradient )
 {
-  real64 dtCompute = explicitStepInternal( time_n, dt, domain );
+  real64 dtCompute = explicitStepInternal( time_n, dt, cycleNumber, domain );
 
   forDiscretizationOnMeshTargets( domain.getMeshBodies(),
                                   [&] ( string const &,
@@ -946,7 +989,7 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
                                                       DomainPartition & domain,
                                                       bool computeGradient )
 {
-  real64 dtCompute = explicitStepInternal( time_n, dt, domain );
+  real64 dtCompute = explicitStepInternal( time_n, dt, cycleNumber, domain );
   forDiscretizationOnMeshTargets( domain.getMeshBodies(),
                                   [&] ( string const &,
                                         MeshLevel & mesh,
@@ -1076,6 +1119,7 @@ void AcousticWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh )
 
 void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n,
                                                real64 const & dt,
+                                               integer const & cycleNumber,
                                                DomainPartition & domain,
                                                MeshLevel & mesh,
                                                arrayView1d< string const > const & regionNames )
@@ -1103,8 +1147,14 @@ void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n,
                                                           getDiscretizationName(),
                                                           "",
                                                           kernelFactory );
+
   //Modification of cycleNember useful when minTime < 0
-  addSourceToRightHandSide( time_n, rhs );
+  EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
+  real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() );
+  //localIndex const cycleNumber = time_n/dt;
+  integer const cycleForSource = int(round( -minTime / dt + cycleNumber ));
+
+  addSourceToRightHandSide( cycleForSource, rhs );
 
   /// calculate your time integrators
   real64 const dt2 = pow( dt, 2 );
@@ -1228,6 +1278,7 @@ void AcousticWaveEquationSEM::synchronizeUnknowns( real64 const & time_n,
 
 real64 AcousticWaveEquationSEM::explicitStepInternal( real64 const & time_n,
                                                       real64 const & dt,
+                                                      integer const & cycleNumber,
                                                       DomainPartition & domain )
 {
   GEOS_MARK_FUNCTION;
@@ -1241,7 +1292,7 @@ real64 AcousticWaveEquationSEM::explicitStepInternal( real64 const & time_n,
   {
     localIndex nSubSteps = (int) ceil( dt/m_timeStep );
     dtCompute = dt/nSubSteps;
-    computeUnknowns( time_n, dtCompute, domain, mesh, regionNames );
+    computeUnknowns( time_n, dtCompute, cycleNumber, domain, mesh, regionNames );
     synchronizeUnknowns( time_n, dtCompute, domain, mesh, regionNames );
   } );
 
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp
index 5f3f1c70885..973c0856c03 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustic/secondOrderEqn/isotropic/AcousticWaveEquationSEM.hpp
@@ -87,7 +87,7 @@ class AcousticWaveEquationSEM : public WaveSolverBase
    * @param cycleNumber the cycle number/step number of evaluation of the source
    * @param rhs the right hand side vector to be computed
    */
-  virtual void addSourceToRightHandSide( real64 const & time_n, arrayView1d< real32 > const rhs );
+  virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs );
 
 
   /**
@@ -123,10 +123,12 @@ class AcousticWaveEquationSEM : public WaveSolverBase
    */
   real64 explicitStepInternal( real64 const & time_n,
                                real64 const & dt,
+                               integer const & cycleNumber,
                                DomainPartition & domain );
 
   void computeUnknowns( real64 const & time_n,
                         real64 const & dt,
+                        integer const & cycleNumber,
                         DomainPartition & domain,
                         MeshLevel & mesh,
                         arrayView1d< string const > const & regionNames );
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEM.cpp
index ce244edf2ec..ef965933529 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/sem/acoustoelastic/secondOrderEqn/isotropic/AcousticElasticWaveEquationSEM.cpp
@@ -129,7 +129,7 @@ void AcousticElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups
 
 real64 AcousticElasticWaveEquationSEM::solverStep( real64 const & time_n,
                                                    real64 const & dt,
-                                                   int const,
+                                                   integer const cycleNumber,
                                                    DomainPartition & domain )
 {
   GEOS_MARK_FUNCTION;
@@ -174,7 +174,7 @@ real64 AcousticElasticWaveEquationSEM::solverStep( real64 const & time_n,
 
     elasSolver->synchronizeUnknowns( time_n, dt, domain, mesh, m_elasRegions );
 
-    acousSolver->computeUnknowns( time_n, dt, domain, mesh, m_acousRegions );
+    acousSolver->computeUnknowns( time_n, dt, cycleNumber, domain, mesh, m_acousRegions );
 
     forAll< EXEC_POLICY >( interfaceNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const in )
     {
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp
index 6fd876ff6ef..cae03815697 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/PrecomputeSourcesAndReceiversKernel.hpp
@@ -72,7 +72,14 @@ struct PreComputeSourcesAndReceivers
                                        arrayView2d< real64 const > const receiverCoordinates,
                                        arrayView1d< localIndex > const receiverIsLocal,
                                        arrayView2d< localIndex > const receiverNodeIds,
-                                       arrayView2d< real64 > const receiverConstants )
+                                       arrayView2d< real64 > const receiverConstants,
+                                       arrayView2d< real32 > const sourceValue,
+                                       real64 const dt,
+                                       real32 const timeSourceFrequency,
+                                       real32 const timeSourceDelay,
+                                       localIndex const rickerOrder,
+                                       arrayView1d< TableFunction::KernelWrapper const > const sourceWaveletTableWrappers,
+                                       bool useSourceWaveletTables )
   {
     constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
 
@@ -121,6 +128,20 @@ struct PreComputeSourcesAndReceivers
               sourceNodeIds[isrc][a] = elemsToNodes( k, a );
               sourceConstants[isrc][a] = Ntest[a];
             }
+
+            for( localIndex cycle = 0; cycle < sourceValue.size( 0 ); ++cycle )
+            {
+              real64 const time_n = cycle * dt;
+              if( useSourceWaveletTables )
+              {
+                sourceValue[cycle][isrc]= sourceWaveletTableWrappers[ isrc ].compute( &time_n );
+              }
+              else
+              {
+                sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( cycle * dt, timeSourceFrequency, timeSourceDelay, rickerOrder );
+              }
+            }
+
           }
         }
       } // end loop over all sources
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp
index b55ca1f38d7..342033418f0 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.cpp
@@ -54,6 +54,13 @@ WaveSolverBase::WaveSolverBase( const std::string & name,
     setSizedFromParent( 0 ).
     setDescription( "Coordinates (x,y,z) of the receivers" );
 
+  registerWrapper( viewKeyStruct::sourceValueString(), &m_sourceValue ).
+    setInputFlag( InputFlags::FALSE ).
+    setRestartFlags( RestartFlags::NO_WRITE ).
+    setSizedFromParent( 0 ).
+    setDescription( "Array which contains the value of the Ricker wavelets at each time-steps" );
+
+
   registerWrapper( viewKeyStruct::timeSourceDelayString(), &m_timeSourceDelay ).
     setInputFlag( InputFlags::OPTIONAL ).
     setApplyDefaultValue( -1 ).
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp
index 2c319ca4f38..06c4dfe752f 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/shared/WaveSolverBase.hpp
@@ -83,6 +83,7 @@ class WaveSolverBase : public PhysicsSolverBase
   struct viewKeyStruct : PhysicsSolverBase::viewKeyStruct
   {
     static constexpr char const * sourceCoordinatesString() { return "sourceCoordinates"; }
+    static constexpr char const * sourceValueString() { return "sourceValue"; }
 
     static constexpr char const * timeSourceFrequencyString() { return "timeSourceFrequency"; }
     static constexpr char const * timeSourceDelayString() { return "timeSourceDelay"; }
@@ -254,6 +255,9 @@ class WaveSolverBase : public PhysicsSolverBase
 
   localIndex getNumNodesPerElem();
 
+  /// Precomputed value of the source terms
+  array2d< real32 > m_sourceValue;
+
   /// Coordinates of the sources in the mesh
   array2d< real64 > m_sourceCoordinates;
 
diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index 87f70867a1e..9c0418e4f5f 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -222,18 +222,10 @@
 				</xsd:unique>
 			</xsd:element>
 			<xsd:element name="Functions" type="FunctionsType" maxOccurs="1">
-				<xsd:unique name="FunctionsCompositeFunctionUniqueName">
-					<xsd:selector xpath="CompositeFunction" />
-					<xsd:field xpath="@name" />
-				</xsd:unique>
 				<xsd:unique name="FunctionsMultivariableTableFunctionUniqueName">
 					<xsd:selector xpath="MultivariableTableFunction" />
 					<xsd:field xpath="@name" />
 				</xsd:unique>
-				<xsd:unique name="FunctionsSymbolicFunctionUniqueName">
-					<xsd:selector xpath="SymbolicFunction" />
-					<xsd:field xpath="@name" />
-				</xsd:unique>
 				<xsd:unique name="FunctionsTableFunctionUniqueName">
 					<xsd:selector xpath="TableFunction" />
 					<xsd:field xpath="@name" />
@@ -395,6 +387,10 @@
 					<xsd:selector xpath="MultiphasePoromechanics" />
 					<xsd:field xpath="@name" />
 				</xsd:unique>
+				<xsd:unique name="SolversMultiphasePoromechanicsConformingFracturesUniqueName">
+					<xsd:selector xpath="MultiphasePoromechanicsConformingFractures" />
+					<xsd:field xpath="@name" />
+				</xsd:unique>
 				<xsd:unique name="SolversMultiphasePoromechanicsReservoirUniqueName">
 					<xsd:selector xpath="MultiphasePoromechanicsReservoir" />
 					<xsd:field xpath="@name" />
@@ -521,6 +517,10 @@
 					<xsd:selector xpath="HydrofractureInitialization" />
 					<xsd:field xpath="@name" />
 				</xsd:unique>
+				<xsd:unique name="TasksMultiphasePoromechanicsConformingFracturesInitializationUniqueName">
+					<xsd:selector xpath="MultiphasePoromechanicsConformingFracturesInitialization" />
+					<xsd:field xpath="@name" />
+				</xsd:unique>
 				<xsd:unique name="TasksMultiphasePoromechanicsInitializationUniqueName">
 					<xsd:selector xpath="MultiphasePoromechanicsInitialization" />
 					<xsd:field xpath="@name" />
@@ -1387,40 +1387,16 @@ stress - traction is applied to the faces as specified by the inner product of i
 	</xsd:simpleType>
 	<xsd:complexType name="FunctionsType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
-			<xsd:element name="CompositeFunction" type="CompositeFunctionType" />
 			<xsd:element name="MultivariableTableFunction" type="MultivariableTableFunctionType" />
-			<xsd:element name="SymbolicFunction" type="SymbolicFunctionType" />
 			<xsd:element name="TableFunction" type="TableFunctionType" />
 		</xsd:choice>
 	</xsd:complexType>
-	<xsd:complexType name="CompositeFunctionType">
-		<!--expression => Composite math expression-->
-		<xsd:attribute name="expression" type="string" default="" />
-		<!--functionNames => List of source functions. The order must match the variableNames argument.-->
-		<xsd:attribute name="functionNames" type="string_array" default="{}" />
-		<!--inputVarNames => Name of fields are input to function.-->
-		<xsd:attribute name="inputVarNames" type="groupNameRef_array" default="{}" />
-		<!--variableNames => List of variables in expression-->
-		<xsd:attribute name="variableNames" type="groupNameRef_array" default="{}" />
-		<!--name => A name is required for any non-unique nodes-->
-		<xsd:attribute name="name" type="groupName" use="required" />
-	</xsd:complexType>
 	<xsd:complexType name="MultivariableTableFunctionType">
 		<!--inputVarNames => Name of fields are input to function.-->
 		<xsd:attribute name="inputVarNames" type="groupNameRef_array" default="{}" />
 		<!--name => A name is required for any non-unique nodes-->
 		<xsd:attribute name="name" type="groupName" use="required" />
 	</xsd:complexType>
-	<xsd:complexType name="SymbolicFunctionType">
-		<!--expression => Symbolic math expression-->
-		<xsd:attribute name="expression" type="string" use="required" />
-		<!--inputVarNames => Name of fields are input to function.-->
-		<xsd:attribute name="inputVarNames" type="groupNameRef_array" default="{}" />
-		<!--variableNames => List of variables in expression.  The order must match the evaluate argument-->
-		<xsd:attribute name="variableNames" type="groupNameRef_array" use="required" />
-		<!--name => A name is required for any non-unique nodes-->
-		<xsd:attribute name="name" type="groupName" use="required" />
-	</xsd:complexType>
 	<xsd:complexType name="TableFunctionType">
 		<!--coordinateFiles => List of coordinate file names for ND Table-->
 		<xsd:attribute name="coordinateFiles" type="path_array" default="{}" />
@@ -1849,6 +1825,8 @@ stress - traction is applied to the faces as specified by the inner product of i
 		<xsd:attribute name="formulation" type="geos_FiniteElementDiscretization_Formulation" default="default" />
 		<!--order => The order of the finite element basis.-->
 		<xsd:attribute name="order" type="integer" use="required" />
+		<!--useHighOrderQuadratureRule => Specifier to indicate whether to use a high order quadrature rule-->
+		<xsd:attribute name="useHighOrderQuadratureRule" type="integer" default="0" />
 		<!--useVirtualElements => Specifier to indicate whether to force the use of VEM-->
 		<xsd:attribute name="useVirtualElements" type="integer" default="0" />
 		<!--name => A name is required for any non-unique nodes-->
@@ -2302,6 +2280,7 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 			<xsd:element name="ImplicitSpringSlider" type="ImplicitSpringSliderType" />
 			<xsd:element name="LaplaceFEM" type="LaplaceFEMType" />
 			<xsd:element name="MultiphasePoromechanics" type="MultiphasePoromechanicsType" />
+			<xsd:element name="MultiphasePoromechanicsConformingFractures" type="MultiphasePoromechanicsConformingFracturesType" />
 			<xsd:element name="MultiphasePoromechanicsReservoir" type="MultiphasePoromechanicsReservoirType" />
 			<xsd:element name="OneWayCoupledFractureFlowContactMechanics" type="OneWayCoupledFractureFlowContactMechanicsType" />
 			<xsd:element name="PhaseFieldDamageFEM" type="PhaseFieldDamageFEMType" />
@@ -2625,6 +2604,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<xsd:attribute name="continuationDBC" type="integer" default="1" />
 		<!--discretization => Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified.-->
 		<xsd:attribute name="discretization" type="groupNameRef" use="required" />
+		<!--gravityDensityScheme => Scheme for density treatment in gravity-->
+		<xsd:attribute name="gravityDensityScheme" type="geos_GravityDensityScheme" default="ArithmeticAverage" />
 		<!--initialDt => Initial time-step value required by the solver to the event manager.-->
 		<xsd:attribute name="initialDt" type="real64" default="1e+99" />
 		<!--isThermal => Flag indicating whether the problem is thermal or not.-->
@@ -2671,7 +2652,7 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<!--scalingType => Solution scaling type.Valid options:
 * Global
 * Local-->
-		<xsd:attribute name="scalingType" type="geos_CompositionalMultiphaseFVM_ScalingType" default="Global" />
+		<xsd:attribute name="scalingType" type="geos_compositionalMultiphaseUtilities_ScalingType" default="Global" />
 		<!--solutionChangeScalingFactor => Damping factor for solution change targets-->
 		<xsd:attribute name="solutionChangeScalingFactor" type="real64" default="0.5" />
 		<!--targetFlowCFL => Target CFL condition `CFL condition <http://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition>`_ when computing the next timestep.-->
@@ -2692,8 +2673,6 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<xsd:attribute name="useDBC" type="integer" default="0" />
 		<!--useMass => Use mass formulation instead of molar. Warning : Affects SourceFlux rates units.-->
 		<xsd:attribute name="useMass" type="integer" default="0" />
-		<!--useNewGravity => Flag indicating whether new gravity treatment is used-->
-		<xsd:attribute name="useNewGravity" type="integer" default="0" />
 		<!--useSimpleAccumulation => Flag indicating whether simple accumulation form is used-->
 		<xsd:attribute name="useSimpleAccumulation" type="integer" default="1" />
 		<!--useTotalMassEquation => Flag indicating whether total mass equation is used-->
@@ -2708,7 +2687,12 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 			<xsd:pattern value=".*[\[\]`$].*|Linfinity|L2" />
 		</xsd:restriction>
 	</xsd:simpleType>
-	<xsd:simpleType name="geos_CompositionalMultiphaseFVM_ScalingType">
+	<xsd:simpleType name="geos_GravityDensityScheme">
+		<xsd:restriction base="xsd:string">
+			<xsd:pattern value=".*[\[\]`$].*|ArithmeticAverage|PhasePresence" />
+		</xsd:restriction>
+	</xsd:simpleType>
+	<xsd:simpleType name="geos_compositionalMultiphaseUtilities_ScalingType">
 		<xsd:restriction base="xsd:string">
 			<xsd:pattern value=".*[\[\]`$].*|Global|Local" />
 		</xsd:restriction>
@@ -2765,8 +2749,6 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<xsd:attribute name="minScalingFactor" type="real64" default="0.01" />
 		<!--solutionChangeScalingFactor => Damping factor for solution change targets-->
 		<xsd:attribute name="solutionChangeScalingFactor" type="real64" default="0.5" />
-		<!--targetFlowCFL => Target CFL condition `CFL condition <http://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition>`_ when computing the next timestep.-->
-		<xsd:attribute name="targetFlowCFL" type="real64" default="-1" />
 		<!--targetPhaseVolFractionChangeInTimeStep => Target (absolute) change in phase volume fraction in a time step-->
 		<xsd:attribute name="targetPhaseVolFractionChangeInTimeStep" type="real64" default="0.2" />
 		<!--targetRegions => Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.-->
@@ -2781,8 +2763,6 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<xsd:attribute name="temperature" type="real64" use="required" />
 		<!--useMass => Use mass formulation instead of molar. Warning : Affects SourceFlux rates units.-->
 		<xsd:attribute name="useMass" type="integer" default="0" />
-		<!--useNewGravity => Flag indicating whether new gravity treatment is used-->
-		<xsd:attribute name="useNewGravity" type="integer" default="0" />
 		<!--useSimpleAccumulation => Flag indicating whether simple accumulation form is used-->
 		<xsd:attribute name="useSimpleAccumulation" type="integer" default="1" />
 		<!--useTotalMassEquation => Flag indicating whether total mass equation is used-->
@@ -3506,6 +3486,51 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<!--stabilizationType => StabilizationType. Options are:
 None- Add no stabilization to mass equation 
 Global- Add jump stabilization to all faces 
+Local- Add jump stabilization on interior of macro elements-->
+		<xsd:attribute name="stabilizationType" type="geos_stabilization_StabilizationType" default="None" />
+		<!--targetRegions => Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.-->
+		<xsd:attribute name="targetRegions" type="groupNameRef_array" use="required" />
+		<!--writeLinearSystem => Write matrix, rhs, solution to screen ( = 1) or file ( = 2).-->
+		<xsd:attribute name="writeLinearSystem" type="integer" default="0" />
+		<!--name => A name is required for any non-unique nodes-->
+		<xsd:attribute name="name" type="groupName" use="required" />
+	</xsd:complexType>
+	<xsd:complexType name="MultiphasePoromechanicsConformingFracturesType">
+		<xsd:choice minOccurs="0" maxOccurs="unbounded">
+			<xsd:element name="LinearSolverParameters" type="LinearSolverParametersType" maxOccurs="1" />
+			<xsd:element name="NonlinearSolverParameters" type="NonlinearSolverParametersType" maxOccurs="1" />
+		</xsd:choice>
+		<!--cflFactor => 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] -->
+		<xsd:attribute name="cflFactor" type="real64" default="0.5" />
+		<!--flowSolverName => Name of the flow solver used by the coupled solver-->
+		<xsd:attribute name="flowSolverName" type="groupNameRef" use="required" />
+		<!--initialDt => Initial time-step value required by the solver to the event manager.-->
+		<xsd:attribute name="initialDt" type="real64" default="1e+99" />
+		<!--isThermal => Flag indicating whether the problem is thermal or not. Set isThermal="1" to enable the thermal coupling-->
+		<xsd:attribute name="isThermal" type="integer" default="0" />
+		<!--logLevel => Sets the level of information to write in the standard output (the console typically).
+Level 0 outputs no specific information for this solver. Higher levels require more outputs.
+1
+ - Line search information
+ - Solution information (scaling, maximum changes, quality check)
+ - Convergence information
+ - Time step information
+ - Linear solver information
+ - Nonlinear solver information
+ - Solver timers information
+ - Coupling information
+2
+ - The summary of declared fields and coupling-->
+		<xsd:attribute name="logLevel" type="integer" default="0" />
+		<!--solidSolverName => Name of the solid solver used by the coupled solver-->
+		<xsd:attribute name="solidSolverName" type="groupNameRef" use="required" />
+		<!--stabilizationMultiplier => Constant multiplier of stabilization strength-->
+		<xsd:attribute name="stabilizationMultiplier" type="real64" default="1" />
+		<!--stabilizationRegionNames => Regions where stabilization is applied.-->
+		<xsd:attribute name="stabilizationRegionNames" type="groupNameRef_array" default="{}" />
+		<!--stabilizationType => StabilizationType. Options are:
+None- Add no stabilization to mass equation 
+Global- Add jump stabilization to all faces 
 Local- Add jump stabilization on interior of macro elements-->
 		<xsd:attribute name="stabilizationType" type="geos_stabilization_StabilizationType" default="None" />
 		<!--targetRegions => Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.-->
@@ -4343,6 +4368,10 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<xsd:attribute name="discretization" type="groupNameRef" use="required" />
 		<!--initialDt => Initial time-step value required by the solver to the event manager.-->
 		<xsd:attribute name="initialDt" type="real64" default="1e+99" />
+		<!--iterPenaltyN => Factor for tuning the iterative penalty coefficient for normal traction-->
+		<xsd:attribute name="iterPenaltyN" type="real64" default="10" />
+		<!--iterPenaltyT => Factor for tuning the iterative penalty coefficient for tangential traction-->
+		<xsd:attribute name="iterPenaltyT" type="real64" default="0.1" />
 		<!--logLevel => Sets the level of information to write in the standard output (the console typically).
 Level 0 outputs no specific information for this solver. Higher levels require more outputs.
 1
@@ -4365,12 +4394,16 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<xsd:attribute name="newmarkBeta" type="real64" default="0.25" />
 		<!--newmarkGamma => Value of :math:`\gamma` in the Newmark Method for Implicit Dynamic time integration option-->
 		<xsd:attribute name="newmarkGamma" type="real64" default="0.5" />
+		<!--simultaneous => Flag to update the Lagrange multiplier at each Newton iteration (true), or only after the Newton loop has converged (false)-->
+		<xsd:attribute name="simultaneous" type="integer" default="1" />
 		<!--stiffnessDamping => Value of stiffness based damping coefficient. -->
 		<xsd:attribute name="stiffnessDamping" type="real64" default="0" />
 		<!--strainTheory => Indicates whether or not to use `Infinitesimal Strain Theory <https://en.wikipedia.org/wiki/Infinitesimal_strain_theory>`_, or `Finite Strain Theory <https://en.wikipedia.org/wiki/Finite_strain_theory>`_. Valid Inputs are:
  0 - Infinitesimal Strain 
  1 - Finite Strain-->
 		<xsd:attribute name="strainTheory" type="integer" default="0" />
+		<!--symmetric => Flag to neglect the non-symmetric contribution in the tangential matrix-->
+		<xsd:attribute name="symmetric" type="integer" default="1" />
 		<!--targetRegions => Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.-->
 		<xsd:attribute name="targetRegions" type="groupNameRef_array" use="required" />
 		<!--timeIntegrationOption => Time integration method. Options are:
@@ -4378,6 +4411,14 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 * ImplicitDynamic
 * ExplicitDynamic-->
 		<xsd:attribute name="timeIntegrationOption" type="geos_SolidMechanicsLagrangianFEM_TimeIntegrationOption" default="ExplicitDynamic" />
+		<!--tolJumpN => Factor to adjust the tolerance for normal jump-->
+		<xsd:attribute name="tolJumpN" type="real64" default="1e-07" />
+		<!--tolJumpT => Factor to adjust the tolerance for tangential jump-->
+		<xsd:attribute name="tolJumpT" type="real64" default="1e-05" />
+		<!--tolNormalTrac => Factor to adjust the tolerance for normal traction-->
+		<xsd:attribute name="tolNormalTrac" type="real64" default="0.5" />
+		<!--tolTauLimit => Tolerance for the sliding check-->
+		<xsd:attribute name="tolTauLimit" type="real64" default="0.05" />
 		<!--writeLinearSystem => Write matrix, rhs, solution to screen ( = 1) or file ( = 2).-->
 		<xsd:attribute name="writeLinearSystem" type="integer" default="0" />
 		<!--name => A name is required for any non-unique nodes-->
@@ -4797,6 +4838,7 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 			<xsd:element name="CompositionalMultiphaseReservoirPoromechanicsInitialization" type="CompositionalMultiphaseReservoirPoromechanicsInitializationType" />
 			<xsd:element name="CompositionalMultiphaseStatistics" type="CompositionalMultiphaseStatisticsType" />
 			<xsd:element name="HydrofractureInitialization" type="HydrofractureInitializationType" />
+			<xsd:element name="MultiphasePoromechanicsConformingFracturesInitialization" type="MultiphasePoromechanicsConformingFracturesInitializationType" />
 			<xsd:element name="MultiphasePoromechanicsInitialization" type="MultiphasePoromechanicsInitializationType" />
 			<xsd:element name="PVTDriver" type="PVTDriverType" />
 			<xsd:element name="PackCollection" type="PackCollectionType" />
@@ -4860,6 +4902,19 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 	<xsd:complexType name="HydrofractureInitializationType">
 		<!--logLevel => Sets the level of information to write in the standard output (the console typically).
 Level 0 outputs no specific information for this solver. Higher levels require more outputs.
+1
+ - Initialization information-->
+		<xsd:attribute name="logLevel" type="integer" default="0" />
+		<!--poromechanicsSolverName => Name of the poromechanics solver-->
+		<xsd:attribute name="poromechanicsSolverName" type="groupNameRef" use="required" />
+		<!--solidMechanicsStatisticsName => Name of the solid mechanics statistics-->
+		<xsd:attribute name="solidMechanicsStatisticsName" type="groupNameRef" default="" />
+		<!--name => A name is required for any non-unique nodes-->
+		<xsd:attribute name="name" type="groupName" use="required" />
+	</xsd:complexType>
+	<xsd:complexType name="MultiphasePoromechanicsConformingFracturesInitializationType">
+		<!--logLevel => Sets the level of information to write in the standard output (the console typically).
+Level 0 outputs no specific information for this solver. Higher levels require more outputs.
 1
  - Initialization information-->
 		<xsd:attribute name="logLevel" type="integer" default="0" />
@@ -5071,11 +5126,11 @@ Level 0 outputs no specific information for this solver. Higher levels require m
 		<!--logLevel => Sets the level of information to write in the standard output (the console typically).
 Level 0 outputs no specific information for this solver. Higher levels require more outputs.
 1
- - Print source flux statistics resulting from all flux in a mesh
+ - Print aggregated statistics of all source fluxes in a mesh
 2
- - Print source flux statistics for each flux in a mesh
+ - Print statistics for each source flux in a mesh
 3
- - Print source flux statistics for each regions -->
+ - Print statistics for each source flux in each regions-->
 		<xsd:attribute name="logLevel" type="integer" default="0" />
 		<!--writeCSV => Write statistics into a CSV file-->
 		<xsd:attribute name="writeCSV" type="integer" default="0" />
diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other
index bdfbfbec836..91ec1f05a2c 100644
--- a/src/coreComponents/schema/schema.xsd.other
+++ b/src/coreComponents/schema/schema.xsd.other
@@ -341,15 +341,11 @@
 	</xsd:complexType>
 	<xsd:complexType name="FunctionsType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
-			<xsd:element name="CompositeFunction" type="CompositeFunctionType" />
 			<xsd:element name="MultivariableTableFunction" type="MultivariableTableFunctionType" />
-			<xsd:element name="SymbolicFunction" type="SymbolicFunctionType" />
 			<xsd:element name="TableFunction" type="TableFunctionType" />
 		</xsd:choice>
 	</xsd:complexType>
-	<xsd:complexType name="CompositeFunctionType" />
 	<xsd:complexType name="MultivariableTableFunctionType" />
-	<xsd:complexType name="SymbolicFunctionType" />
 	<xsd:complexType name="TableFunctionType" />
 	<xsd:complexType name="GeometryType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -537,6 +533,7 @@
 			<xsd:element name="ImplicitSpringSlider" type="ImplicitSpringSliderType" />
 			<xsd:element name="LaplaceFEM" type="LaplaceFEMType" />
 			<xsd:element name="MultiphasePoromechanics" type="MultiphasePoromechanicsType" />
+			<xsd:element name="MultiphasePoromechanicsConformingFractures" type="MultiphasePoromechanicsConformingFracturesType" />
 			<xsd:element name="MultiphasePoromechanicsReservoir" type="MultiphasePoromechanicsReservoirType" />
 			<xsd:element name="OneWayCoupledFractureFlowContactMechanics" type="OneWayCoupledFractureFlowContactMechanicsType" />
 			<xsd:element name="PhaseFieldDamageFEM" type="PhaseFieldDamageFEMType" />
@@ -635,6 +632,8 @@
 		<xsd:attribute name="sourceNodeIds" type="integer_array2d" />
 		<!--sourceRegion => Region containing the sources-->
 		<xsd:attribute name="sourceRegion" type="integer_array" />
+		<!--sourceValue => Source Value of the sources-->
+		<xsd:attribute name="sourceValue" type="real32_array2d" />
 		<!--timeStep => TimeStep computed with the power iteration method (if we don't want to compute it, it is initialized with the XML value-->
 		<xsd:attribute name="timeStep" type="real64" />
 		<!--usePML => Flag to apply PML-->
@@ -682,6 +681,8 @@
 		<xsd:attribute name="sourceIsAccessible" type="integer_array" />
 		<!--sourceNodeIds => Indices of the nodes (in the right order) for each source point-->
 		<xsd:attribute name="sourceNodeIds" type="integer_array2d" />
+		<!--sourceValue => Source Value of the sources-->
+		<xsd:attribute name="sourceValue" type="real32_array2d" />
 		<!--timeStep => TimeStep computed with the power iteration method (if we don't want to compute it, it is initialized with the XML value-->
 		<xsd:attribute name="timeStep" type="real64" />
 		<!--usePML => Flag to apply PML-->
@@ -723,6 +724,8 @@
 		<xsd:attribute name="sourceIsAccessible" type="integer_array" />
 		<!--sourceNodeIds => Indices of the nodes (in the right order) for each source point-->
 		<xsd:attribute name="sourceNodeIds" type="integer_array2d" />
+		<!--sourceValue => Source Value of the sources-->
+		<xsd:attribute name="sourceValue" type="real32_array2d" />
 		<!--timeStep => TimeStep computed with the power iteration method (if we don't want to compute it, it is initialized with the XML value-->
 		<xsd:attribute name="timeStep" type="real64" />
 		<!--usePML => Flag to apply PML-->
@@ -857,6 +860,8 @@
 		<xsd:attribute name="sourceNodeIds" type="integer_array2d" />
 		<!--sourceRegion => Region containing the sources-->
 		<xsd:attribute name="sourceRegion" type="integer_array" />
+		<!--sourceValue => Source Value of the sources-->
+		<xsd:attribute name="sourceValue" type="real32_array2d" />
 		<!--timeStep => TimeStep computed with the power iteration method (if we don't want to compute it, it is initialized with the XML value-->
 		<xsd:attribute name="timeStep" type="real64" />
 		<!--usePML => Flag to apply PML-->
@@ -904,6 +909,8 @@
 		<xsd:attribute name="sourceIsAccessible" type="integer_array" />
 		<!--sourceNodeIds => Indices of the nodes (in the right order) for each source point-->
 		<xsd:attribute name="sourceNodeIds" type="integer_array2d" />
+		<!--sourceValue => Source Value of the sources-->
+		<xsd:attribute name="sourceValue" type="real32_array2d" />
 		<!--timeStep => TimeStep computed with the power iteration method (if we don't want to compute it, it is initialized with the XML value-->
 		<xsd:attribute name="timeStep" type="real64" />
 		<!--usePML => Flag to apply PML-->
@@ -929,7 +936,7 @@
 		<!--maxStableDt => Value of the Maximum Stable Timestep for this solver.-->
 		<xsd:attribute name="maxStableDt" type="real64" />
 		<!--meshTargets => MeshBody/Region combinations that the solver will be applied to.-->
-		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std___1_pair_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__gt__cm_-LvArray_Array_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std___1_integral_constant_lt_bool_cm_-true_gt__gt_" />
+		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std_pair_lt_string_cm_-string-_gt__cm_-LvArray_Array_lt_string_cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std_integral_constant_lt_bool_cm_-true_gt_-_gt_" />
 	</xsd:complexType>
 	<xsd:complexType name="ExplicitSpringSliderType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -940,7 +947,7 @@
 		<!--maxStableDt => Value of the Maximum Stable Timestep for this solver.-->
 		<xsd:attribute name="maxStableDt" type="real64" />
 		<!--meshTargets => MeshBody/Region combinations that the solver will be applied to.-->
-		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std___1_pair_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__gt__cm_-LvArray_Array_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std___1_integral_constant_lt_bool_cm_-true_gt__gt_" />
+		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std_pair_lt_string_cm_-string-_gt__cm_-LvArray_Array_lt_string_cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std_integral_constant_lt_bool_cm_-true_gt_-_gt_" />
 	</xsd:complexType>
 	<xsd:complexType name="FlowProppantTransportType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -990,11 +997,7 @@
 		<!--maxStableDt => Value of the Maximum Stable Timestep for this solver.-->
 		<xsd:attribute name="maxStableDt" type="real64" />
 		<!--meshTargets => MeshBody/Region combinations that the solver will be applied to.-->
-		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std___1_pair_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__gt__cm_-LvArray_Array_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std___1_integral_constant_lt_bool_cm_-true_gt__gt_" />
-
 		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std_pair_lt_string_cm_-string-_gt__cm_-LvArray_Array_lt_string_cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std_integral_constant_lt_bool_cm_-true_gt_-_gt_" />
-		<!--performStressInitialization => Flag to indicate that the solver is going to perform stress initialization-->
-		<xsd:attribute name="performStressInitialization" type="integer" />
 	</xsd:complexType>
 	<xsd:complexType name="LaplaceFEMType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -1005,7 +1008,7 @@
 		<!--maxStableDt => Value of the Maximum Stable Timestep for this solver.-->
 		<xsd:attribute name="maxStableDt" type="real64" />
 		<!--meshTargets => MeshBody/Region combinations that the solver will be applied to.-->
-		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std___1_pair_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__gt__cm_-LvArray_Array_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std___1_integral_constant_lt_bool_cm_-true_gt__gt_" />
+		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std_pair_lt_string_cm_-string-_gt__cm_-LvArray_Array_lt_string_cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std_integral_constant_lt_bool_cm_-true_gt_-_gt_" />
 	</xsd:complexType>
 	<xsd:complexType name="MultiphasePoromechanicsType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -1018,10 +1021,24 @@
 		<!--maxStableDt => Value of the Maximum Stable Timestep for this solver.-->
 		<xsd:attribute name="maxStableDt" type="real64" />
 		<!--meshTargets => MeshBody/Region combinations that the solver will be applied to.-->
-		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std___1_pair_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__gt__cm_-LvArray_Array_lt_std___1_basic_string_lt_char_cm_-std___1_char_traits_lt_char_gt__cm_-std___1_allocator_lt_char_gt__gt__cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std___1_integral_constant_lt_bool_cm_-true_gt__gt_" />
+		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std_pair_lt_string_cm_-string-_gt__cm_-LvArray_Array_lt_string_cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std_integral_constant_lt_bool_cm_-true_gt_-_gt_" />
 		<!--performStressInitialization => Flag to indicate that the solver is going to perform stress initialization-->
 		<xsd:attribute name="performStressInitialization" type="integer" />
+	</xsd:complexType>
+	<xsd:complexType name="MultiphasePoromechanicsConformingFracturesType">
+		<xsd:choice minOccurs="0" maxOccurs="unbounded">
+			<xsd:element name="LinearSolverParameters" type="LinearSolverParametersType" maxOccurs="1" />
+			<xsd:element name="NonlinearSolverParameters" type="NonlinearSolverParametersType" maxOccurs="1" />
+			<xsd:element name="SolverStatistics" type="SolverStatisticsType" />
+		</xsd:choice>
+		<!--discretization => Name of discretization object (defined in the :ref:`NumericalMethodsManager`) to use for this solver. For instance, if this is a Finite Element Solver, the name of a :ref:`FiniteElement` should be specified. If this is a Finite Volume Method, the name of a :ref:`FiniteVolume` discretization should be specified.-->
+		<xsd:attribute name="discretization" type="groupNameRef" />
+		<!--maxStableDt => Value of the Maximum Stable Timestep for this solver.-->
+		<xsd:attribute name="maxStableDt" type="real64" />
+		<!--meshTargets => MeshBody/Region combinations that the solver will be applied to.-->
 		<xsd:attribute name="meshTargets" type="geos_mapBase_lt_std_pair_lt_string_cm_-string-_gt__cm_-LvArray_Array_lt_string_cm_-1_cm_-camp_int_seq_lt_long_cm_-0l_gt__cm_-int_cm_-LvArray_ChaiBuffer_gt__cm_-std_integral_constant_lt_bool_cm_-true_gt_-_gt_" />
+		<!--performStressInitialization => Flag to indicate that the solver is going to perform stress initialization-->
+		<xsd:attribute name="performStressInitialization" type="integer" />
 	</xsd:complexType>
 	<xsd:complexType name="MultiphasePoromechanicsReservoirType">
 		<xsd:choice minOccurs="0" maxOccurs="unbounded">
@@ -1445,6 +1462,7 @@
 			<xsd:element name="CompositionalMultiphaseReservoirPoromechanicsInitialization" type="CompositionalMultiphaseReservoirPoromechanicsInitializationType" />
 			<xsd:element name="CompositionalMultiphaseStatistics" type="CompositionalMultiphaseStatisticsType" />
 			<xsd:element name="HydrofractureInitialization" type="HydrofractureInitializationType" />
+			<xsd:element name="MultiphasePoromechanicsConformingFracturesInitialization" type="MultiphasePoromechanicsConformingFracturesInitializationType" />
 			<xsd:element name="MultiphasePoromechanicsInitialization" type="MultiphasePoromechanicsInitializationType" />
 			<xsd:element name="PVTDriver" type="PVTDriverType" />
 			<xsd:element name="PackCollection" type="PackCollectionType" />
@@ -1475,6 +1493,7 @@
 	<xsd:complexType name="CompositionalMultiphaseReservoirPoromechanicsInitializationType" />
 	<xsd:complexType name="CompositionalMultiphaseStatisticsType" />
 	<xsd:complexType name="HydrofractureInitializationType" />
+	<xsd:complexType name="MultiphasePoromechanicsConformingFracturesInitializationType" />
 	<xsd:complexType name="MultiphasePoromechanicsInitializationType" />
 	<xsd:complexType name="PVTDriverType" />
 	<xsd:complexType name="PackCollectionType" />
diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAdjoint1.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAdjoint1.cpp
index fb4bcd943e3..1837b7e793b 100644
--- a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAdjoint1.cpp
+++ b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagationAdjoint1.cpp
@@ -390,8 +390,8 @@ TEST_F( AcousticWaveEquationSEMTest, SeismoTrace )
   std::cout << " / ||f'||.||u||=" << std::sqrt( sum_fb2*sum_u2 ) << " / ||f||.||f'||=" << std::sqrt( sum_ff2*sum_fb2 ) << std::endl;
   real32 diffToCheck;
   diffToCheck=std::abs( sum_ufb-sum_qff ) / std::max( std::sqrt( sum_fb2*sum_u2 ), std::sqrt( sum_q2*sum_ff2 ));
-  std::cout << " Diff to compare with 2.e-4: " << diffToCheck << std::endl;
-  ASSERT_TRUE( diffToCheck < 2.e-4 );
+  std::cout << " Diff to compare with 9e-3: " << diffToCheck << std::endl;
+  ASSERT_TRUE( diffToCheck < 9e-3 );
 }
 
 int main( int argc, char * * argv )

From 3cddb9f5afd714e8fefd2ad79d70d92bc45d91be Mon Sep 17 00:00:00 2001
From: Pavel Tomin <paveltomin@users.noreply.github.com>
Date: Thu, 6 Feb 2025 10:32:04 -0600
Subject: [PATCH 5/5] refactor: Preparatory for well to frac connection (#3227)

---
 .integrated_tests.yaml                        |   2 +-
 BASELINE_NOTES.md                             |   4 +
 ...echanics_FaultModel_well_fim_new_smoke.xml | 116 +++++++
 ...Poromechanics_FaultModel_well_new_base.xml | 235 +++++++++++++
 ...echanics_FaultModel_well_seq_new_smoke.xml | 120 +++++++
 .../finiteVolume/FiniteVolumeManager.cpp      |   1 -
 .../mesh/MeshForLoopInterface.hpp             |  66 +++-
 src/coreComponents/mesh/Perforation.cpp       |   5 +
 src/coreComponents/mesh/Perforation.hpp       |  17 +-
 .../mesh/WellElementSubRegion.cpp             | 312 ++++++++++--------
 .../mesh/generators/LineBlock.hpp             |  10 +
 .../mesh/generators/LineBlockABC.hpp          |   6 +
 .../mesh/generators/MeshGeneratorBase.cpp     |   1 +
 .../mesh/generators/WellGeneratorBase.cpp     |   2 +
 .../mesh/generators/WellGeneratorBase.hpp     |   9 +
 .../contact/SolidMechanicsLagrangeContact.cpp |   1 -
 .../multiphysics/PhaseFieldFractureSolver.cpp |   1 -
 .../SolidMechanicsLagrangianFEM.cpp           |   1 +
 .../SolidMechanicsLagrangianFEM.hpp           |   1 -
 .../solidMechanics/SolidMechanicsMPM.cpp      |   1 +
 .../solidMechanics/SolidMechanicsMPM.hpp      |   1 -
 .../EmbeddedSurfaceGenerator.cpp              |   2 +-
 .../surfaceGeneration/SurfaceGenerator.cpp    |   2 +-
 23 files changed, 770 insertions(+), 146 deletions(-)
 create mode 100644 inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_new_smoke.xml
 create mode 100644 inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml
 create mode 100644 inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml

diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
index 35bfac096d1..af490bff930 100644
--- a/.integrated_tests.yaml
+++ b/.integrated_tests.yaml
@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3502-10013-5b4e015
+  baseline: integratedTests/baseline_integratedTests-pr3227-10037-26bd879
 
 allow_fail:
   all: ''
diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md
index cb85ce3d56a..1da00fbf042 100644
--- a/BASELINE_NOTES.md
+++ b/BASELINE_NOTES.md
@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR  #3227 (2024-02-01)
+=====================
+Add targetRegion for perforations (optional).
+
 PR #3502 (2025-02-04)
 =====================
 Add array to store the source values in time inside wave solvers
diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_new_smoke.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_new_smoke.xml
new file mode 100644
index 00000000000..05679f48257
--- /dev/null
+++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_fim_new_smoke.xml
@@ -0,0 +1,116 @@
+<Problem>
+  <Included>
+    <File name="./singlePhasePoromechanics_FaultModel_well_new_base.xml"/>
+  </Included>
+
+  <Solvers
+    gravityVector="{ 0.0, 0.0, 0.0 }">
+
+    <SinglePhasePoromechanicsConformingFracturesReservoir
+      name="reservoirSolver"
+      targetRegions="{ Region, Fault, wellRegion1, wellRegion2 }" 
+      initialDt="1e-3"
+      poromechanicsConformingFracturesSolverName="poroFractureSolver"
+      wellSolverName="singlePhaseWellSolver"
+      logLevel="1">
+      <NonlinearSolverParameters
+        lineSearchAction="None"
+        maxTimeStepCuts="10"
+        newtonTol="1e-3"
+        newtonMaxIter="20"
+        maxAllowedResidualNorm="1e+25"/>
+      <LinearSolverParameters
+        solverType="direct"/> 
+    </SinglePhasePoromechanicsConformingFracturesReservoir>
+
+    <SinglePhasePoromechanicsConformingFractures
+      name="poroFractureSolver"
+      targetRegions="{ Region, Fault }" 
+      flowSolverName="singlePhaseFlowSolver"
+      solidSolverName="fractureMechSolver"
+      logLevel="1">
+    </SinglePhasePoromechanicsConformingFractures>
+
+    <SolidMechanicsLagrangeContact
+      name="fractureMechSolver"
+      stabilizationName="singlePhaseTPFA_stab"
+      logLevel="1"
+      targetRegions="{ Region, Fault }" 
+      timeIntegrationOption="QuasiStatic"
+      discretization="FE1"/>  
+
+    <SinglePhaseFVM 
+      name="singlePhaseFlowSolver"
+      logLevel="1"
+      allowNegativePressure="1"
+      discretization="singlePhaseTPFA"
+      targetRegions="{ Region, Fault}" />
+
+    <SinglePhaseWell
+      name="singlePhaseWellSolver"
+      logLevel="1"
+      targetRegions="{ wellRegion1, wellRegion2 }">
+      <WellControls
+        name="wellControls1"
+        logLevel="2"
+        type="producer"
+        control="BHP"
+        referenceElevation="-2500"
+        targetBHP="1e7"
+        targetTotalRate="1e10" />
+      <WellControls
+        name="wellControls2"
+        logLevel="2"
+        type="injector"
+        control="BHP"
+        referenceElevation="-2500"
+        targetBHP="10e7"
+        targetTotalRate="1e10" />
+    </SinglePhaseWell>
+  </Solvers>
+
+  <Tasks>
+      <SinglePhasePoromechanicsConformingFracturesInitialization 
+        logLevel="1" 
+        name="SinglePhasePoroelasticityPreEquilibrationStep"  
+        poromechanicsSolverName="poroFractureSolver"/>
+  </Tasks>
+
+  <Events
+    minTime="-1e11"
+    maxTime="3.0e7">
+
+    <!-- Initialization -->
+    <SoloEvent 
+      name="SinglePhasePoroelasticityPreEquilibrationStep" 
+      targetTime="-1e11" 
+      beginTime="-1e11" 
+      target="/Tasks/SinglePhasePoroelasticityPreEquilibrationStep"
+    />
+    <!-- Injection -->
+    <SoloEvent
+      name="outputPostEquilibrationStep1"
+      targetTime="0"
+      targetExactTimestep="1"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="outputs"
+      beginTime="0.0"
+      endTime="3.0e7"
+      forceDt="3.0e7"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="solverApplication"
+      endTime="3.0e7"
+      forceDt="3.0e7"
+      target="/Solvers/reservoirSolver"/>
+
+    <PeriodicEvent
+      name="restarts"
+      beginTime="0.0"
+      endTime="3.0e7"
+      target="/Outputs/restartOutput"/>
+  </Events>
+</Problem>
diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml
new file mode 100644
index 00000000000..f91d3f1e8f1
--- /dev/null
+++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_new_base.xml
@@ -0,0 +1,235 @@
+<Problem>
+  <Constitutive>
+    <PorousElasticIsotropic
+      name="porousRock"
+      permeabilityModelName="rockPerm"
+      porosityModelName="rockPorosity"
+      solidModelName="rock" />
+
+    <ElasticIsotropic
+      name="rock"
+      defaultDensity="2927"
+      defaultPoissonRatio="0.15"
+      defaultYoungModulus="14950000000" />
+
+    <CompressibleSinglePhaseFluid
+      name="water"
+      compressibility="2.09028227021e-10"
+      defaultDensity="1325"
+      defaultViscosity="0.001"
+      referenceDensity="1325"
+      referencePressure="0e6"
+      referenceViscosity="0.001"
+      viscosibility="0.0" />
+
+    <BiotPorosity
+      name="rockPorosity"
+      defaultReferencePorosity="0.3"
+      defaultGrainBulkModulus="7.12e30" />
+
+    <ConstantPermeability
+      name="rockPerm"
+      permeabilityComponents="{ 1.0e-18, 1.0e-18, 1.0e-18 }" />
+      <!-- Material inside the fault -->
+    <CompressibleSolidParallelPlatesPermeability
+      name="faultFilling"
+      solidModelName="nullSolid"
+      porosityModelName="fracturePorosity"
+      permeabilityModelName="fracturePerm"/>
+
+    <PressurePorosity
+      name="fracturePorosity"
+      defaultReferencePorosity="1.00"
+      referencePressure="0e6"
+      compressibility="0.0"/>
+
+    <ParallelPlatesPermeability
+      name="fracturePerm"/>
+
+    <Coulomb
+      name="fractureContact"
+      cohesion="40.0e6"
+      frictionCoefficient="0.01"/>
+
+    <HydraulicApertureTable
+    name="hApertureModel"
+    apertureTableName="apertureTable"/>   
+
+    <NullModel
+      name="nullSolid"/>
+
+  </Constitutive>
+
+  <FieldSpecifications>
+
+    <FieldSpecification
+      name="fracPorosity"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions/Fault"
+      fieldName="fracturePorosity_porosity"
+      scale="1.00"/>
+
+    <FieldSpecification
+      name="xconstraint"
+      component="0"
+      fieldName="totalDisplacement"
+      objectPath="nodeManager"
+      scale="0.0"
+      setNames="{ xnegFace, xposFace }" />
+
+    <FieldSpecification
+      name="yconstraint"
+      component="1"
+      fieldName="totalDisplacement"
+      objectPath="nodeManager"
+      scale="0.0"
+      setNames="{ ynegFace, yposFace }" />
+
+    <FieldSpecification
+      name="zconstraint"
+      component="2"
+      fieldName="totalDisplacement"
+      objectPath="nodeManager"
+      scale="0.0"
+      setNames="{ znegFace }" /> 
+
+    <Traction
+      name="tractionTop"
+      objectPath="faceManager"
+      tractionType="normal"
+      scale="-4e7"      
+      setNames="{ zposFace }"/>
+
+    <FieldSpecification
+      name="init_pressure"
+      initialCondition="1"
+      setNames="{ all }"
+      objectPath="ElementRegions"
+      fieldName="pressure"
+      scale="4e7"/>
+  </FieldSpecifications>
+
+  <ElementRegions>
+    <CellElementRegion
+      name="Region"
+      cellBlocks="{ 2_hexahedra, 1_hexahedra, 0_hexahedra }"
+      materialList="{ water, porousRock }"/>  
+
+    <SurfaceElementRegion
+      name="Fault"
+      faceBlock="fracture"
+      materialList="{ water, faultFilling, fractureContact, hApertureModel }"
+      defaultAperture="1e-3"/>
+
+    <WellElementRegion
+      name="wellRegion1"
+      materialList="{ water }"/>
+
+    <WellElementRegion
+      name="wellRegion2"
+      materialList="{ water }"/>
+  </ElementRegions>
+
+  <NumericalMethods>
+    <FiniteElements>
+      <FiniteElementSpace
+        name="FE1"
+        order="1" />
+    </FiniteElements>
+
+    <FiniteVolume>
+      <TwoPointFluxApproximation
+        name="singlePhaseTPFA"/>
+      <TwoPointFluxApproximation
+        name="singlePhaseTPFA_stab"/>  
+    </FiniteVolume>
+  </NumericalMethods>
+
+  <Outputs>
+    <VTK
+      name="vtkOutput" />
+
+    <Restart
+      name="restartOutput"/>
+  </Outputs>
+
+  <Functions>
+    <TableFunction
+      name="apertureTable"
+      coordinates="{ -1.0e-3, 0.0 }"
+      values="{ 1.0e-3, 1.0e-3 }"
+    />
+  </Functions>
+
+  <Mesh>
+    <VTKMesh
+      name="mesh1"
+      useGlobalIds="1"
+      faceBlocks="{ fracture }"
+      file="verticalFault_ExternalMesh.vtm">
+      <InternalWell
+        name="well_producer1"
+        logLevel="1"
+        wellRegionName="wellRegion1"
+        wellControlsName="wellControls1"
+        polylineNodeCoords="{ { 5400, 5400, 0 },
+                              { 5400, 5400, -2500 } }"
+        polylineSegmentConn="{ { 0, 1 } }"
+        radius="0.1"
+        numElementsPerSegment="1">
+        <Perforation
+          name="producer1_perf1"
+          targetRegion="Region"
+          distanceFromHead="2500"/>
+      </InternalWell>
+      <InternalWell
+        name="well_injector1"
+        logLevel="1"
+        wellRegionName="wellRegion2"
+        wellControlsName="wellControls2"
+        polylineNodeCoords="{ { -5400, -5400, 0 },
+                              { -5400, -5400, -2500 } }"
+        polylineSegmentConn="{ { 0, 1 } }"
+        radius="0.1"
+        numElementsPerSegment="1">
+        <Perforation
+          name="injector1_perf1"
+          distanceFromHead="2500"/>
+      </InternalWell>
+    </VTKMesh>
+  </Mesh>
+
+  <Geometry>
+    <Box
+      name="xposFace"
+      xMin="{ 5999.99,  -6000.01,  -4000.01 }"
+      xMax="{ 6000.01,   6000.01,      0.01 }"/>
+
+    <Box
+      name="xnegFace"
+      xMin="{-6000.01,  -6000.01,  -4000.01 }"
+      xMax="{-5999.99,   6000.01,      0.01 }"/>
+
+    <Box
+      name="yposFace"
+      xMin="{-6000.01,   5999.99,  -4000.01 }"
+      xMax="{ 6000.01,   6000.01,      0.01 }"/>
+
+    <Box
+      name="ynegFace"
+      xMin="{-6000.01,  -6000.01,  -4000.01 }"
+      xMax="{ 6000.01,  -5999.99,      0.01 }"/>
+
+    <Box
+      name="zposFace"
+      xMin="{-6000.01,  -6000.01,     -0.01 }"
+      xMax="{ 6000.01,   6000.01,      0.01 }"/>
+
+    <Box
+      name="znegFace"
+      xMin="{-6000.01,  -6000.01,  -4000.01 }"
+      xMax="{ 6000.01,   6000.01,  -3999.99 }"/>
+  </Geometry>
+
+</Problem>
diff --git a/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml
new file mode 100644
index 00000000000..2f56fcd4625
--- /dev/null
+++ b/inputFiles/poromechanicsFractures/singlePhasePoromechanics_FaultModel_well_seq_new_smoke.xml
@@ -0,0 +1,120 @@
+<Problem>
+  <Included>
+    <File name="./singlePhasePoromechanics_FaultModel_well_new_base.xml"/>
+  </Included>
+
+  <Solvers
+    gravityVector="{ 0.0, 0.0, 0.0 }">
+
+    <SinglePhaseReservoirPoromechanicsConformingFractures
+      name="reservoirSolver"
+      targetRegions="{ Region, Fault }" 
+      initialDt="1e-3"
+      reservoirAndWellsSolverName="reservoirAndWellsSolver"
+      solidSolverName="fractureMechSolver"
+      logLevel="1">
+      <NonlinearSolverParameters
+        couplingType="Sequential"
+        lineSearchAction="None"
+        maxTimeStepCuts="10"
+        newtonTol="1e-3"
+        newtonMaxIter="20"/>
+      <LinearSolverParameters
+        solverType="direct"/> 
+    </SinglePhaseReservoirPoromechanicsConformingFractures>
+
+    <SolidMechanicsLagrangeContact
+      name="fractureMechSolver"
+      stabilizationName="singlePhaseTPFA_stab"
+      logLevel="1"
+      targetRegions="{ Region, Fault }" 
+      timeIntegrationOption="QuasiStatic"
+      discretization="FE1"/>  
+
+    <SinglePhaseReservoir
+      name="reservoirAndWellsSolver"
+      flowSolverName="singlePhaseFlowSolver"
+      wellSolverName="singlePhaseWellSolver"
+      logLevel="1"
+      targetRegions="{ Region, Fault, wellRegion1, wellRegion2 }">
+      <NonlinearSolverParameters
+        newtonMaxIter="40"
+        maxAllowedResidualNorm="1e+25"/>
+      <LinearSolverParameters
+        directParallel="0"/>
+    </SinglePhaseReservoir>
+
+    <SinglePhaseFVM 
+      name="singlePhaseFlowSolver"
+      logLevel="1"
+      discretization="singlePhaseTPFA"
+      targetRegions="{ Region, Fault}" />
+
+    <SinglePhaseWell
+      name="singlePhaseWellSolver"
+      logLevel="1"
+      targetRegions="{ wellRegion1, wellRegion2 }">
+      <WellControls
+        name="wellControls1"
+        logLevel="2"
+        type="producer"
+        control="BHP"
+        referenceElevation="-2500"
+        targetBHP="1e7"
+        targetTotalRate="1e10" />
+      <WellControls
+        name="wellControls2"
+        logLevel="2"
+        type="injector"
+        control="BHP"
+        referenceElevation="-2500"
+        targetBHP="10e7"
+        targetTotalRate="1e10" />
+    </SinglePhaseWell>
+  </Solvers>
+
+  <Tasks>
+      <SinglePhaseReservoirPoromechanicsConformingFracturesInitialization 
+        logLevel="1" 
+        name="SinglePhasePoroelasticityPreEquilibrationStep"  
+        poromechanicsSolverName="reservoirSolver"/>
+  </Tasks>
+
+  <Events
+    minTime="-1e11"
+    maxTime="3.0e7">
+
+    <!-- Initialization -->
+    <SoloEvent 
+      name="SinglePhasePoroelasticityPreEquilibrationStep" 
+      targetTime="-1e11" 
+      beginTime="-1e11" 
+      target="/Tasks/SinglePhasePoroelasticityPreEquilibrationStep"
+    />
+    <!-- Injection -->
+    <SoloEvent
+      name="outputPostEquilibrationStep1"
+      targetTime="0"
+      targetExactTimestep="1"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="outputs"
+      beginTime="0.0"
+      endTime="3.0e7"
+      forceDt="3.0e7"
+      target="/Outputs/vtkOutput"/>
+
+    <PeriodicEvent
+      name="solverApplication"
+      endTime="3.0e7"
+      forceDt="3.0e7"
+      target="/Solvers/reservoirSolver"/>
+
+    <PeriodicEvent
+      name="restarts"
+      beginTime="0.0"
+      endTime="3.0e7"
+      target="/Outputs/restartOutput"/>
+  </Events>
+</Problem>
diff --git a/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp b/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp
index b4329194c67..ef6c94b5fb5 100644
--- a/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp
+++ b/src/coreComponents/finiteVolume/FiniteVolumeManager.cpp
@@ -22,7 +22,6 @@
 
 #include "finiteVolume/FluxApproximationBase.hpp"
 #include "finiteVolume/HybridMimeticDiscretization.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
 #include "common/GEOS_RAJA_Interface.hpp"
 
 namespace geos
diff --git a/src/coreComponents/mesh/MeshForLoopInterface.hpp b/src/coreComponents/mesh/MeshForLoopInterface.hpp
index b7c126b1e23..eb1e6040d89 100644
--- a/src/coreComponents/mesh/MeshForLoopInterface.hpp
+++ b/src/coreComponents/mesh/MeshForLoopInterface.hpp
@@ -69,11 +69,11 @@ minLocOverElemsInMesh( MeshLevel const & mesh, LAMBDA && lambda )
 
   ElementRegionManager const & elemManager = mesh.getElemManager();
 
-  for( localIndex er=0; er<elemManager.numRegions(); ++er )
+  for( localIndex er = 0; er < elemManager.numRegions(); ++er )
   {
     ElementRegionBase const & elemRegion = elemManager.getRegion( er );
 
-    elemRegion.forElementSubRegionsIndex< CellElementSubRegion >( [&]( localIndex const esr, CellElementSubRegion const & subRegion )
+    elemRegion.forElementSubRegionsIndex< ElementSubRegionBase >( [&]( localIndex const esr, ElementSubRegionBase const & subRegion )
     {
       localIndex const size = subRegion.size();
       for( localIndex k = 0; k < size; ++k )
@@ -93,6 +93,68 @@ minLocOverElemsInMesh( MeshLevel const & mesh, LAMBDA && lambda )
   return std::make_pair( minVal, std::make_tuple( minReg, minSubreg, minIndex ));
 }
 
+/**
+ * @brief @return Return the minimum location/indices for a value condition specified by @p lambda.
+ * @tparam LAMBDA The type of the lambda function to be used to specify the minimum condition.
+ * @param region The region that will have all of its elements processed by this function.
+ * @param lambda  A lambda function that returns as value that will be used in the minimum comparison.
+ */
+template< typename LAMBDA >
+auto
+minLocOverElemsInRegion( ElementRegionBase const & region, LAMBDA && lambda )
+{
+  using NUMBER = decltype( lambda( 0, 0 ) );
+
+  NUMBER minVal = std::numeric_limits< NUMBER >::max();
+  localIndex minSubreg = -1, minIndex = -1;
+
+  region.forElementSubRegionsIndex< ElementSubRegionBase >( [&]( localIndex const esr, ElementSubRegionBase const & subRegion )
+  {
+    localIndex const size = subRegion.size();
+    for( localIndex k = 0; k < size; ++k )
+    {
+      NUMBER const val = lambda( esr, k );
+      if( val < minVal )
+      {
+        minVal = val;
+        minSubreg = esr;
+        minIndex = k;
+      }
+    }
+  } );
+
+  return std::make_pair( minVal, std::make_tuple( minSubreg, minIndex ));
+}
+
+/**
+ * @brief @return Return the minimum location/indices for a value condition specified by @p lambda.
+ * @tparam LAMBDA The type of the lambda function to be used to specify the minimum condition.
+ * @param subRegion The subregion that will have all of its elements processed by this function.
+ * @param lambda  A lambda function that returns as value that will be used in the minimum comparison.
+ */
+template< typename LAMBDA >
+auto
+minLocOverElemsInSubRegion( ElementSubRegionBase const & subRegion, LAMBDA && lambda )
+{
+  using NUMBER = decltype( lambda( 0 ) );
+
+  NUMBER minVal = std::numeric_limits< NUMBER >::max();
+  localIndex minIndex = -1;
+
+  localIndex const size = subRegion.size();
+  for( localIndex k = 0; k < size; ++k )
+  {
+    NUMBER const val = lambda( k );
+    if( val < minVal )
+    {
+      minVal = val;
+      minIndex = k;
+    }
+  }
+
+  return std::make_pair( minVal, minIndex );
+}
+
 } // namespace geos
 
 #endif // GEOS_MESH_MESHFORLOOPINTERFACE_HPP
diff --git a/src/coreComponents/mesh/Perforation.cpp b/src/coreComponents/mesh/Perforation.cpp
index 74ade6d4269..77c5ea57a81 100644
--- a/src/coreComponents/mesh/Perforation.cpp
+++ b/src/coreComponents/mesh/Perforation.cpp
@@ -46,6 +46,11 @@ Perforation::Perforation( string const & name, Group * const parent )
     setApplyDefaultValue( 0.0 ).
     setInputFlag( InputFlags::OPTIONAL ).
     setDescription( "Perforation skin factor" );
+
+  registerWrapper( viewKeyStruct::targetRegionString(), &m_targetRegionName ).
+    setRTTypeName( rtTypes::CustomTypes::groupNameRef ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "Target region to connect the perforation" );
 }
 
 
diff --git a/src/coreComponents/mesh/Perforation.hpp b/src/coreComponents/mesh/Perforation.hpp
index 5a3a3ae8721..fb3636ef5dc 100644
--- a/src/coreComponents/mesh/Perforation.hpp
+++ b/src/coreComponents/mesh/Perforation.hpp
@@ -102,6 +102,12 @@ class Perforation : public dataRepository::Group
    */
   real64 getWellSkinFactor() const { return m_wellSkinFactor; }
 
+  /**
+   * @brief Get the target region for the perforation.
+   * @return the list of target regions
+   */
+  string const & getTargetRegion() const { return m_targetRegionName; }
+
   ///@}
 
   /**
@@ -116,12 +122,8 @@ class Perforation : public dataRepository::Group
     static constexpr char const *wellTransmissibilityString() { return "transmissibility"; }
     /// @return String key for the well skin factor at this perforation
     static constexpr char const *wellSkinFactorString() { return "skinFactor"; }
-    /// ViewKey for the linear distance from well head
-    dataRepository::ViewKey distanceFromHead = {distanceFromHeadString()};
-    /// ViewKey for the well transmissibility at this perforation
-    dataRepository::ViewKey wellTransmissibility = {wellTransmissibilityString()};
-    /// ViewKey for the well transmissibility at this perforation
-    dataRepository::ViewKey wellSkinFactor = { wellSkinFactorString() };
+    /// @return Target region for this perforation
+    static constexpr char const *targetRegionString() { return "targetRegion"; }
   }
   /// ViewKey struct for the Perforation class
   viewKeysPerforation;
@@ -138,6 +140,9 @@ class Perforation : public dataRepository::Group
 
   /// Well skin factor at this perforation
   real64 m_wellSkinFactor;
+
+  /// Name of region the perforation will be connected to
+  string m_targetRegionName;
 };
 
 } // namespace geos
diff --git a/src/coreComponents/mesh/WellElementSubRegion.cpp b/src/coreComponents/mesh/WellElementSubRegion.cpp
index cd8d5cde24f..c16ea1c6d00 100644
--- a/src/coreComponents/mesh/WellElementSubRegion.cpp
+++ b/src/coreComponents/mesh/WellElementSubRegion.cpp
@@ -120,7 +120,8 @@ void collectLocalAndBoundaryNodes( LineBlockABC const & lineBlock,
  * @param[in] ei the index of the reservoir element
  * @param[inout] nodes the nodes that have already been visited
  */
-void collectElementNodes( CellElementSubRegion const & subRegion,
+template< typename SUBREGION_TYPE >
+void collectElementNodes( SUBREGION_TYPE const & subRegion,
                           localIndex ei,
                           SortedArray< localIndex > & nodes )
 {
@@ -137,6 +138,33 @@ void collectElementNodes( CellElementSubRegion const & subRegion,
   }
 }
 
+template< typename SUBREGION_TYPE >
+bool isPointInsideElement( SUBREGION_TYPE const & GEOS_UNUSED_PARAM( subRegion ),
+                           arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & GEOS_UNUSED_PARAM( referencePosition ),
+                           localIndex const & GEOS_UNUSED_PARAM( eiLocal ),
+                           ArrayOfArraysView< localIndex const > const & GEOS_UNUSED_PARAM( facesToNodes ),
+                           real64 const (&GEOS_UNUSED_PARAM( elemCenter ))[3],
+                           real64 const (&GEOS_UNUSED_PARAM( location ))[3] )
+{
+  // only CellElementSubRegion is currently supported
+  return false;
+}
+
+bool isPointInsideElement( CellElementSubRegion const & subRegion,
+                           arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & referencePosition,
+                           localIndex const & eiLocal,
+                           ArrayOfArraysView< localIndex const > const & facesToNodes,
+                           real64 const (&elemCenter)[3],
+                           real64 const (&location)[3] )
+{
+  arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList();
+  return computationalGeometry::isPointInsidePolyhedron( referencePosition,
+                                                         elemsToFaces[eiLocal],
+                                                         facesToNodes,
+                                                         elemCenter,
+                                                         location );
+}
+
 /**
  * @brief Search the reservoir elements that can be accessed from the set "nodes".
           Stop if a reservoir element containing the perforation is found.
@@ -145,17 +173,18 @@ void collectElementNodes( CellElementSubRegion const & subRegion,
  * @param[in] location the location of that we are trying to match with a reservoir element
  * @param[inout] nodes the nodes that have already been visited
  * @param[inout] elements the reservoir elements that have already been visited
- * @param[inout] erMatched the region index of the reservoir element that contains "location", if any
- * @param[inout] esrMatched the subregion index of the reservoir element that contains "location", if any
+ * @param[in] targetRegionIndex the target region index for the reservoir element
+ * @param[in] targetSubRegionIndex the target subregion index for the reservoir element
  * @param[inout] eiMatched the element index of the reservoir element that contains "location", if any
  * @param[inout] giMatched the element global index of the reservoir element that contains "location", if any
  */
+template< typename SUBREGION_TYPE >
 bool visitNeighborElements( MeshLevel const & mesh,
                             real64 const (&location)[3],
                             SortedArray< localIndex > & nodes,
                             SortedArray< globalIndex > & elements,
-                            localIndex & erMatched,
-                            localIndex & esrMatched,
+                            localIndex const & targetRegionIndex,
+                            localIndex const & targetSubRegionIndex,
                             localIndex & eiMatched,
                             globalIndex & giMatched )
 {
@@ -166,6 +195,7 @@ bool visitNeighborElements( MeshLevel const & mesh,
   ArrayOfArraysView< localIndex const > const & toElementRegionList    = nodeManager.elementRegionList();
   ArrayOfArraysView< localIndex const > const & toElementSubRegionList = nodeManager.elementSubRegionList();
   ArrayOfArraysView< localIndex const > const & toElementList          = nodeManager.elementList();
+
   arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const referencePosition =
     nodeManager.referencePosition().toViewConst();
 
@@ -189,14 +219,17 @@ bool visitNeighborElements( MeshLevel const & mesh,
   for( localIndex currNode : currNodes )
   {
     // collect the elements that have not been visited yet
-    for( localIndex b=0; b<toElementRegionList.sizeOfArray( currNode ); ++b )
+    for( localIndex b = 0; b < toElementRegionList.sizeOfArray( currNode ); ++b )
     {
       localIndex const er      = toElementRegionList[currNode][b];
       localIndex const esr     = toElementSubRegionList[currNode][b];
       localIndex const eiLocal = toElementList[currNode][b];
-      CellElementRegion const & region = elemManager.getRegion< CellElementRegion >( er );
-      CellElementSubRegion const & subRegion = region.getSubRegion< CellElementSubRegion >( esr );
-      arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList();
+
+      if( er != targetRegionIndex || esr != targetSubRegionIndex )
+        continue;
+
+      ElementRegionBase const & region = elemManager.getRegion< ElementRegionBase >( er );
+      SUBREGION_TYPE const & subRegion = region.getSubRegion< SUBREGION_TYPE >( esr );
       arrayView2d< real64 const > const elemCenters = subRegion.getElementCenter();
 
       globalIndex const eiGlobal = subRegion.localToGlobalMap()[eiLocal];
@@ -212,17 +245,11 @@ bool visitNeighborElements( MeshLevel const & mesh,
 
         // perform the test to see if the point is in this reservoir element
         // if the point is in the resevoir element, save the indices and stop the search
-        if( computationalGeometry::isPointInsidePolyhedron( referencePosition,
-                                                            elemsToFaces[eiLocal],
-                                                            facesToNodes,
-                                                            elemCenter,
-                                                            location ) )
+        if( isPointInsideElement( subRegion, referencePosition, eiLocal, facesToNodes, elemCenter, location ) )
         {
-          erMatched  = er;
-          esrMatched = esr;
-          eiMatched  = eiLocal;
+          eiMatched = eiLocal;
           giMatched = eiGlobal;
-          matched    = true;
+          matched   = true;
           break;
         }
         // otherwise add the nodes of this element to the set of new nodes to visit
@@ -250,36 +277,33 @@ bool visitNeighborElements( MeshLevel const & mesh,
           contains the well element.
  * @param[in] meshLevel the mesh object (single level only)
  * @param[in] location the location of that we are trying to match with a reservoir element
- * @param[inout] erInit the region index of the reservoir element from which we start the search
- * @param[inout] esrInit the subregion index of the reservoir element from which we start the search
+ * @param[in] targetRegionIndex the region index of the reservoir element from which we start the search
+ * @param[in] targetSubRegionIndex the subregion index of the reservoir element from which we start the search
  * @param[inout] eiInit the element index of the reservoir element from which we start the search
  */
 void initializeLocalSearch( MeshLevel const & mesh,
                             real64 const (&location)[3],
-                            localIndex & erInit,
-                            localIndex & esrInit,
+                            localIndex const & targetRegionIndex,
+                            localIndex const & targetSubRegionIndex,
                             localIndex & eiInit )
 {
+  ElementSubRegionBase const & subRegion = mesh.getElemManager().getRegion( targetRegionIndex ).getSubRegion( targetSubRegionIndex );
   ElementRegionManager::ElementViewAccessor< arrayView2d< real64 const > >
   resElemCenter = mesh.getElemManager().constructViewAccessor< array2d< real64 >,
                                                                arrayView2d< real64 const > >( ElementSubRegionBase::viewKeyStruct::elementCenterString() );
   // to initialize the local search for the reservoir element that contains "location",
   // we find the reservoir element that minimizes the distance from "location" to the reservoir element center
-  auto ret = minLocOverElemsInMesh( mesh, [&] ( localIndex const er,
-                                                localIndex const esr,
-                                                localIndex const ei )
+  auto ret = minLocOverElemsInSubRegion( subRegion, [&] ( localIndex const ei )
   {
     real64 v[3] = { location[0], location[1], location[2] };
-    LvArray::tensorOps::subtract< 3 >( v, resElemCenter[er][esr][ei] );
-    return LvArray::tensorOps::l2Norm< 3 >( v );
+    LvArray::tensorOps::subtract< 3 >( v, resElemCenter[targetRegionIndex][targetSubRegionIndex][ei] );
+    auto dist = LvArray::tensorOps::l2Norm< 3 >( v );
+    return dist;
   } );
 
-  // save the region, subregion and index of the reservoir element
+  // save the index of the reservoir element
   // note that this reservoir element does not necessarily contains "location"
-  erInit  = std::get< 0 >( ret.second );
-  esrInit = std::get< 1 >( ret.second );
-  eiInit  = std::get< 2 >( ret.second );
-
+  eiInit  = ret.second;
 }
 
 /**
@@ -287,10 +311,7 @@ void initializeLocalSearch( MeshLevel const & mesh,
           To do that, loop over the reservoir elements that are in the neighborhood of (erInit,esrInit,eiInit)
  * @param[in] meshLevel the mesh object (single level only)
  * @param[in] location the location of that we are trying to match with a reservoir element
- * @param[in] erInit the region index of the reservoir element from which we start the search
- * @param[in] esrInit the subregion index of the reservoir element from which we start the search
- * @param[in] eiInit the element index of the reservoir element from which we start the search
- * @param[inout] erMatched the region index of the reservoir element that contains "location", if any
+ * @param[in] targetRegionIndex the target region index for the reservoir element
  * @param[inout] esrMatched the subregion index of the reservoir element that contains "location", if any
  * @param[inout] eiMatched the element index of the reservoir element that contains "location", if any
  * @param[inout] giMatched the element global index of the reservoir element that contains "location", if any
@@ -298,50 +319,82 @@ void initializeLocalSearch( MeshLevel const & mesh,
 bool searchLocalElements( MeshLevel const & mesh,
                           real64 const (&location)[3],
                           localIndex const & searchDepth,
-                          localIndex const & erInit,
-                          localIndex const & esrInit,
-                          localIndex const & eiInit,
-                          localIndex & erMatched,
+                          localIndex const & targetRegionIndex,
                           localIndex & esrMatched,
                           localIndex & eiMatched,
                           globalIndex & giMatched )
 {
-  // search locally, starting from the location of the previous perforation
-  // the assumption here is that perforations have been entered in order of depth
-  bool resElemFound = false;
+  ElementRegionBase const & region = mesh.getElemManager().getRegion< ElementRegionBase >( targetRegionIndex );
 
-  CellElementRegion const & region = mesh.getElemManager().getRegion< CellElementRegion >( erInit );
-  CellElementSubRegion const & subRegion = region.getSubRegion< CellElementSubRegion >( esrInit );
-
-  SortedArray< localIndex >  nodes;
-  SortedArray< globalIndex > elements;
-
-  // here is how the search is done:
-  //   1 - We check if "location" is within the "init" reservoir element defined by (erInit,esrMatched,eiMatched)
-  //   2 - If yes, stop
-  //     - If not, a) collect the nodes of the reservoir element defined by (erInit,esrMatched,eiMatched)
-  //               b) use these nodes to grab the neighbors of (erInit,esrMatched,eiMatched)
-  //               c) check if "location" is within the neighbors. If not, grab the neighbors of the neighbors, and so
-  // on...
-
-  // collect the nodes of the current element
-  // they will be used to access the neighbors and check if they contain the perforation
-  collectElementNodes( subRegion, eiInit, nodes );
-  // if no match is found, enlarge the neighborhood m_searchDepth'th times
-  for( localIndex d = 0; d < searchDepth; ++d )
+  bool resElemFound = false;
+  for( localIndex esr = 0; esr < region.numSubRegions(); ++esr )
   {
-    localIndex nNodes = nodes.size();
-
-    // search the reservoir elements that can be accessed from the set "nodes"
-    // stop if a reservoir element containing the perforation is found
-    // if not, enlarge the set "nodes"
-    resElemFound = visitNeighborElements( mesh, location, nodes, elements,
-                                          erMatched, esrMatched, eiMatched, giMatched );
-    if( resElemFound || nNodes == nodes.size())
+    ElementSubRegionBase const & subRegionBase = region.getSubRegion( esr );
+    region.applyLambdaToContainer< CellElementSubRegion, SurfaceElementSubRegion >( subRegionBase, [&]( auto const & subRegion )
+    {
+      GEOS_LOG_RANK_0( GEOS_FMT( "  searching well connections with region/subregion: {}/{}", region.getName(), subRegion.getName() ) );
+
+      // first, we search for the reservoir element that is the *closest* from the center of well element
+      // note that this reservoir element does not necessarily contain the center of the well element
+      // this "init" reservoir element will be used later to find the reservoir element that
+      // contains the well element
+      localIndex eiInit = -1;
+
+      initializeLocalSearch( mesh, location, targetRegionIndex, esr, eiInit );
+
+      if( eiInit < 0 ) // nothing found, skip the rest
+        return;
+
+      // loop over the reservoir elements that are in the neighborhood of (esrInit,eiInit)
+      // search locally, starting from the location of the previous perforation
+      // the assumption here is that perforations have been entered in order of depth
+
+      SortedArray< localIndex >  nodes;
+      SortedArray< globalIndex > elements;
+
+      // here is how the search is done:
+      //   1 - We check if "location" is within the "init" reservoir element defined by (erInit,esrMatched,eiMatched)
+      //   2 - If yes, stop
+      //     - If not, a) collect the nodes of the reservoir element defined by (erInit,esrMatched,eiMatched)
+      //               b) use these nodes to grab the neighbors of (erInit,esrMatched,eiMatched)
+      //               c) check if "location" is within the neighbors. If not, grab the neighbors of the neighbors, and so
+      // on...
+
+      // collect the nodes of the current element
+      // they will be used to access the neighbors and check if they contain the perforation
+      collectElementNodes( subRegion, eiInit, nodes );
+
+      // if no match is found, enlarge the neighborhood m_searchDepth'th times
+      for( localIndex d = 0; d < searchDepth; ++d )
+      {
+        localIndex nNodes = nodes.size();
+
+        // search the reservoir elements that can be accessed from the set "nodes"
+        // stop if a reservoir element containing the perforation is found
+        // if not, enlarge the set "nodes"
+
+        resElemFound =
+          visitNeighborElements< TYPEOFREF( subRegion ) >( mesh, location, nodes, elements, targetRegionIndex, esr, eiMatched, giMatched );
+
+        if( resElemFound || nNodes == nodes.size())
+        {
+          if( resElemFound )
+          {
+            esrMatched = esr;
+            GEOS_LOG( GEOS_FMT( "    found {}/{}/{}", region.getName(), subRegion.getName(), giMatched ) );
+          }
+          // TODO learn how to exit forElementSubRegionsIndex
+          break;
+        }
+      }
+    } );
+
+    if( resElemFound )
     {
       break;
     }
   }
+
   return resElemFound;
 }
 
@@ -448,10 +501,11 @@ void WellElementSubRegion::generate( MeshLevel & mesh,
 
 void WellElementSubRegion::assignUnownedElementsInReservoir( MeshLevel & mesh,
                                                              LineBlockABC const & lineBlock,
-                                                             SortedArray< globalIndex >      const & unownedElems,
+                                                             SortedArray< globalIndex > const & unownedElems,
                                                              SortedArray< globalIndex > & localElems,
                                                              arrayView1d< integer > & elemStatusGlobal ) const
 {
+  ElementRegionManager const & elemManager = mesh.getElemManager();
   // get the well and reservoir element coordinates
   arrayView2d< real64 const > const & wellElemCoordsGlobal = lineBlock.getElemCoords();
 
@@ -464,31 +518,21 @@ void WellElementSubRegion::assignUnownedElementsInReservoir( MeshLevel & mesh,
                                  wellElemCoordsGlobal[currGlobal][1],
                                  wellElemCoordsGlobal[currGlobal][2] };
 
-    // this will contain the indices of the reservoir element
-    // in which the center of the well element is located
-    localIndex erMatched  = -1;
-    localIndex esrMatched = -1;
-    localIndex eiMatched  = -1;
-
-    // this will contain the indices of the reservoir element
-    // from which we are going to start the search
-    localIndex erInit     = -1;
-    localIndex esrInit    = -1;
-    localIndex eiInit     = -1;
-    globalIndex giMatched   = -1;
-
-    // Step 1: first, we search for the reservoir element that is the *closest* from the center of well element
-    //         note that this reservoir element does not necessarily contain the center of the well element
-    //         this "init" reservoir element will be used in SearchLocalElements to find the reservoir element that
-    //         contains the well element
-    initializeLocalSearch( mesh, location,
-                           erInit, esrInit, eiInit );
-
-    // Step 2: then, search for the reservoir element that contains the well element
-    //         to do that, we loop over the reservoir elements that are in the neighborhood of (erInit,esrInit,eiInit)
-    bool resElemFound = searchLocalElements( mesh, location, m_searchDepth,
-                                             erInit, esrInit, eiInit,
-                                             erMatched, esrMatched, eiMatched, giMatched );
+    // for each perforation, we have to find the reservoir element that contains the perforation
+    bool resElemFound = false;
+    for( localIndex er = 0; er < elemManager.numRegions(); er++ )
+    {
+      // search for the reservoir element that contains the well element
+      localIndex esrMatched = -1;
+      localIndex eiMatched  = -1;
+      globalIndex giMatched = -1;
+      resElemFound = searchLocalElements( mesh, location, m_searchDepth, er, esrMatched, eiMatched, giMatched );
+
+      if( resElemFound )
+      {
+        break;
+      }
+    }
 
     // if the element was found
     if( resElemFound )
@@ -770,59 +814,67 @@ void WellElementSubRegion::connectPerforationsToMeshElements( MeshLevel & mesh,
   arrayView2d< real64 const > const perfCoordsGlobal = lineBlock.getPerfCoords();
   arrayView1d< real64 const > const perfWellTransmissibilityGlobal = lineBlock.getPerfTransmissibility();
   arrayView1d< real64 const > const perfWellSkinFactorGlobal = lineBlock.getPerfSkinFactor();
+  arrayView1d< string const > const perfTargetRegionGlobal = lineBlock.getPerfTargetRegion();
 
   m_perforationData.resize( perfCoordsGlobal.size( 0 ) );
   localIndex iperfLocal = 0;
 
   arrayView2d< real64 > const perfLocation = m_perforationData.getLocation();
 
+  ElementRegionManager const & elemManager = mesh.getElemManager();
+
   // loop over all the perforations
   for( globalIndex iperfGlobal = 0; iperfGlobal < perfCoordsGlobal.size( 0 ); ++iperfGlobal )
   {
     real64 const location[3] = { perfCoordsGlobal[iperfGlobal][0],
                                  perfCoordsGlobal[iperfGlobal][1],
                                  perfCoordsGlobal[iperfGlobal][2] };
+    GEOS_LOG_RANK_0( GEOS_FMT( "{}: perforation {} location = ({}, {}, {})",
+                               lineBlock.getName(), iperfGlobal,
+                               location[0], location[1], location[2] ) );
 
-    localIndex erMatched  = -1;
-    localIndex esrMatched = -1;
-    localIndex eiMatched  = -1;
+    localIndex erStart = -1, erEnd = -1;
 
-    localIndex erInit     = -1;
-    localIndex esrInit    = -1;
-    localIndex eiInit     = -1;
-    globalIndex giMatched  = -1;
+    localIndex const targetRegionIndex = elemManager.getRegions().getIndex( perfTargetRegionGlobal[iperfGlobal] );
+    if( targetRegionIndex >= 0 )
+    {
+      erStart = targetRegionIndex;
+      erEnd = erStart + 1;
+    }
+    else // default is all regions
+    {
+      erStart = 0;
+      erEnd = elemManager.numRegions();
+    }
 
     // for each perforation, we have to find the reservoir element that contains the perforation
+    for( localIndex er = erStart; er < erEnd; er++ )
+    {
+      // search for the reservoir element that contains the well element
+      localIndex esrMatched = -1;
+      localIndex eiMatched  = -1;
+      globalIndex giMatched  = -1;
+      bool const resElemFound = searchLocalElements( mesh, location, m_searchDepth, er, esrMatched, eiMatched, giMatched );
+
+      // if the element was found
+      if( resElemFound )
+      {
+        // set the indices for the matched reservoir element
+        m_perforationData.getMeshElements().m_toElementRegion[iperfLocal] = er;
+        m_perforationData.getMeshElements().m_toElementSubRegion[iperfLocal] = esrMatched;
+        m_perforationData.getMeshElements().m_toElementIndex[iperfLocal] = eiMatched;
+        m_perforationData.getReservoirElementGlobalIndex()[iperfLocal] = giMatched;
 
-    // Step 1: first, we search for the reservoir element that is the *closest* from the center of well element
-    //         note that this reservoir element does not necessarily contain the center of the well element
-    //         this "init" reservoir element will be used in SearchLocalElements to find the reservoir element that
-    //         contains the well element
-    initializeLocalSearch( mesh, location,
-                           erInit, esrInit, eiInit );
+        // construct the local wellTransmissibility and location maps
+        m_perforationData.getWellTransmissibility()[iperfLocal] = perfWellTransmissibilityGlobal[iperfGlobal];
+        m_perforationData.getWellSkinFactor()[iperfLocal] = perfWellSkinFactorGlobal[iperfGlobal];
+        LvArray::tensorOps::copy< 3 >( perfLocation[iperfLocal], location );
 
-    // Step 2: then, search for the reservoir element that contains the well element
-    //         to do that, we loop over the reservoir elements that are in the neighborhood of (erInit,esrInit,eiInit)
-    bool resElemFound = searchLocalElements( mesh, location, m_searchDepth,
-                                             erInit, esrInit, eiInit,
-                                             erMatched, esrMatched, eiMatched, giMatched );
+        // increment the local to global map
+        m_perforationData.localToGlobalMap()[iperfLocal++] = iperfGlobal;
 
-    // if the element was found
-    if( resElemFound )
-    {
-      // set the indices for the matched reservoir element
-      m_perforationData.getMeshElements().m_toElementRegion[iperfLocal] = erMatched;
-      m_perforationData.getMeshElements().m_toElementSubRegion[iperfLocal] = esrMatched;
-      m_perforationData.getMeshElements().m_toElementIndex[iperfLocal] = eiMatched;
-      m_perforationData.getReservoirElementGlobalIndex()[iperfLocal] = giMatched;
-
-      // construct the local wellTransmissibility and location maps
-      m_perforationData.getWellTransmissibility()[iperfLocal] = perfWellTransmissibilityGlobal[iperfGlobal];
-      m_perforationData.getWellSkinFactor()[iperfLocal] = perfWellSkinFactorGlobal[iperfGlobal];
-      LvArray::tensorOps::copy< 3 >( perfLocation[iperfLocal], location );
-
-      // increment the local to global map
-      m_perforationData.localToGlobalMap()[iperfLocal++] = iperfGlobal;
+        break;
+      }
     }
   }
 
diff --git a/src/coreComponents/mesh/generators/LineBlock.hpp b/src/coreComponents/mesh/generators/LineBlock.hpp
index e2af5724ec5..6d660510816 100644
--- a/src/coreComponents/mesh/generators/LineBlock.hpp
+++ b/src/coreComponents/mesh/generators/LineBlock.hpp
@@ -154,6 +154,7 @@ class LineBlock : public LineBlockABC
 
   arrayView1d< real64 const > getPerfSkinFactor() const override final { return m_perfSkinFactor; }
 
+  arrayView1d< string const > getPerfTargetRegion() const override final { return m_perfTargetRegion; }
 
   /**
    * @brief Set the well skin factor at the perforations.
@@ -161,6 +162,12 @@ class LineBlock : public LineBlockABC
    */
   void setPerfSkinFactor( arrayView1d< real64 const > perfSkinFactor ) { m_perfSkinFactor = perfSkinFactor; }
 
+  /**
+   * @brief Set the target region for the perforations.
+   * @param perfTargetRegion list of target regions for all the perforations on the well
+   */
+  void setPerfTargetRegion( arrayView1d< string const > perfTargetRegion ) { m_perfTargetRegion = perfTargetRegion; }
+
   arrayView1d< globalIndex const > getPerfElemIndex() const override final { return m_perfElemId; }
 
   /**
@@ -242,6 +249,9 @@ class LineBlock : public LineBlockABC
   /// Well skin factor at the perforation
   array1d< real64 > m_perfSkinFactor;
 
+  /// Target region for the perforation
+  array1d< string > m_perfTargetRegion;
+
   /// Global index of the well element
   array1d< globalIndex > m_perfElemId;
 
diff --git a/src/coreComponents/mesh/generators/LineBlockABC.hpp b/src/coreComponents/mesh/generators/LineBlockABC.hpp
index 2855789ed66..ad56b5b1caa 100644
--- a/src/coreComponents/mesh/generators/LineBlockABC.hpp
+++ b/src/coreComponents/mesh/generators/LineBlockABC.hpp
@@ -145,6 +145,12 @@ class LineBlockABC : public dataRepository::Group
    */
   virtual arrayView1d< real64 const > getPerfSkinFactor() const = 0;
 
+  /**
+   * @brief Get the target region for the perforations.
+   * @return list of target regions for all the perforations on the well
+   */
+  virtual arrayView1d< string const > getPerfTargetRegion() const = 0;
+
   /**
    * @brief Get the global indices of the well elements connected to each perforation.
    * @return list providing the global index of the connected well element for each perforation
diff --git a/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp b/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp
index 86eb7f0b849..da8e7ad355b 100644
--- a/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp
+++ b/src/coreComponents/mesh/generators/MeshGeneratorBase.cpp
@@ -90,6 +90,7 @@ void MeshGeneratorBase::attachWellInfo( CellBlockManager & cellBlockManager )
     lb.setPerfCoords( wellGen.getPerfCoords() );
     lb.setPerfTransmissibility( wellGen.getPerfTransmissibility() );
     lb.setPerfSkinFactor( wellGen.getPerfSkinFactor() );
+    lb.setPerfTargetRegion( wellGen.getPerfTargetRegion() );
     lb.setPerfElemIndex( wellGen.getPerfElemIndex() );
     lb.setWellControlsName( wellGen.getWellControlsName() );
     lb.setWellGeneratorName( wellGen.getName() );
diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp
index 5debc76589f..e9b9e5387b0 100644
--- a/src/coreComponents/mesh/generators/WellGeneratorBase.cpp
+++ b/src/coreComponents/mesh/generators/WellGeneratorBase.cpp
@@ -113,6 +113,7 @@ void WellGeneratorBase::generateWellGeometry( )
   m_perfDistFromHead.resize( m_numPerforations );
   m_perfTransmissibility.resize( m_numPerforations );
   m_perfSkinFactor.resize( m_numPerforations );
+  m_perfTargetRegion.resize( m_numPerforations );
   m_perfElemId.resize( m_numPerforations );
 
   // construct a reverse map from the polyline nodes to the segments
@@ -343,6 +344,7 @@ void WellGeneratorBase::connectPerforationsToWellElements()
     m_perfDistFromHead[iperf]  = perf.getDistanceFromWellHead();
     m_perfTransmissibility[iperf] = perf.getWellTransmissibility();
     m_perfSkinFactor[iperf] = perf.getWellSkinFactor();
+    m_perfTargetRegion[iperf] = perf.getTargetRegion();
 
     // search in all the elements of this well between head and bottom
     globalIndex iwelemTop    = 0;
diff --git a/src/coreComponents/mesh/generators/WellGeneratorBase.hpp b/src/coreComponents/mesh/generators/WellGeneratorBase.hpp
index fbe05157d6f..db7505d8f02 100644
--- a/src/coreComponents/mesh/generators/WellGeneratorBase.hpp
+++ b/src/coreComponents/mesh/generators/WellGeneratorBase.hpp
@@ -185,6 +185,12 @@ class WellGeneratorBase : public MeshComponentBase
    */
   arrayView1d< real64 const > getPerfSkinFactor() const { return m_perfSkinFactor; };
 
+  /**
+   * @brief Get the target region for a perforation.
+   * @return the target regions for a perforation
+   */
+  arrayView1d< string const > getPerfTargetRegion() const { return m_perfTargetRegion; };
+
   /**
    * @brief Get the global indices of the well elements connected to each perforation.
    * @return list providing the global index of the connected well element for each perforation
@@ -393,6 +399,9 @@ class WellGeneratorBase : public MeshComponentBase
   /// Physical location of the perforation wrt to well head
   array1d< real64 > m_perfDistFromHead;
 
+  /// Target region for the perforation
+  array1d< string > m_perfTargetRegion;
+
 };
 }
 #endif /* GEOS_MESH_GENERATORS_WELLGENERATORBASE_HPP_ */
diff --git a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
index 41f39281a2d..af569535465 100644
--- a/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
+++ b/src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp
@@ -30,7 +30,6 @@
 #include "discretizationMethods/NumericalMethodsManager.hpp"
 #include "mainInterface/ProblemManager.hpp"
 #include "mesh/SurfaceElementRegion.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
 #include "mesh/mpiCommunications/NeighborCommunicator.hpp"
 #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" // needed to register pressure(_n)
 #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp"
diff --git a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp
index c9fe35f6145..a9589a4424a 100644
--- a/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp
+++ b/src/coreComponents/physicsSolvers/multiphysics/PhaseFieldFractureSolver.cpp
@@ -25,7 +25,6 @@
 #include "finiteElement/Kinematics.h"
 #include "finiteElement/FiniteElementDispatch.hpp"
 #include "mesh/DomainPartition.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
 #include "mesh/utilities/ComputationalGeometry.hpp"
 
 namespace geos
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
index dd19b20d9b3..3622a33879b 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp
@@ -33,6 +33,7 @@
 #include "fieldSpecification/FieldSpecificationManager.hpp"
 #include "fieldSpecification/TractionBoundaryCondition.hpp"
 #include "finiteElement/FiniteElementDiscretizationManager.hpp"
+#include "finiteElement/FiniteElementDiscretization.hpp"
 #include "LvArray/src/output.hpp"
 #include "mainInterface/ProblemManager.hpp"
 #include "mesh/DomainPartition.hpp"
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp
index 76d98cb5143..ee5d8049f0a 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp
@@ -24,7 +24,6 @@
 #include "common/TimingMacros.hpp"
 #include "kernels/SolidMechanicsLagrangianFEMKernels.hpp"
 #include "kernels/StrainHelper.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
 #include "mesh/mpiCommunications/CommunicationTools.hpp"
 #include "mesh/mpiCommunications/MPI_iCommData.hpp"
 #include "physicsSolvers/PhysicsSolverBase.hpp"
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
index c576f6513e2..ac93eec252e 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
@@ -31,6 +31,7 @@
 #include "constitutive/ConstitutiveManager.hpp"
 #include "constitutive/contact/FrictionBase.hpp"
 #include "finiteElement/FiniteElementDiscretizationManager.hpp"
+#include "finiteElement/FiniteElementDiscretization.hpp"
 #include "finiteElement/Kinematics.h"
 #include "LvArray/src/output.hpp"
 #include "mesh/DomainPartition.hpp"
diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp
index 571e365f14a..672ce06c89a 100644
--- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp
+++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp
@@ -24,7 +24,6 @@
 #include "common/TimingMacros.hpp"
 #include "kernels/SolidMechanicsLagrangianFEMKernels.hpp"
 #include "kernels/ExplicitMPM.hpp"
-#include "mesh/MeshForLoopInterface.hpp"
 #include "mesh/mpiCommunications/CommunicationTools.hpp"
 #include "mesh/mpiCommunications/MPI_iCommData.hpp"
 #include "physicsSolvers/PhysicsSolverBase.hpp"
diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp
index da4e2856368..dd8fe2d42ce 100644
--- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp
+++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp
@@ -191,7 +191,7 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups()
 
             if( added )
             {
-              GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, "Element " << cellIndex << " is fractured" );
+              GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, GEOS_FMT( "{}: element {} is fractured", subRegion.getName(), cellIndex ) );
 
               // Add the information to the CellElementSubRegion
               subRegion.addFracturedElement( cellIndex, localNumberOfSurfaceElems );
diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
index e41674c9f82..76d6dd85a9f 100644
--- a/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
+++ b/src/coreComponents/physicsSolvers/surfaceGeneration/SurfaceGenerator.cpp
@@ -22,7 +22,7 @@
 #include "mesh/mpiCommunications/CommunicationTools.hpp"
 #include "mesh/mpiCommunications/NeighborCommunicator.hpp"
 #include "mesh/mpiCommunications/SpatialPartition.hpp"
-#include "finiteElement/FiniteElementDiscretizationManager.hpp"
+#include "finiteElement/FiniteElementDiscretization.hpp"
 #include "finiteVolume/FiniteVolumeManager.hpp"
 #include "finiteVolume/FluxApproximationBase.hpp"
 #include "mesh/SurfaceElementRegion.hpp"