Skip to content

Commit

Permalink
move creation of EmbeddedFractures upon input.
Browse files Browse the repository at this point in the history
  • Loading branch information
@CusiniM
CusiniM committed Jan 15, 2025
1 parent 1bf96de commit b1fa4c4
Show file tree
Hide file tree
Showing 4 changed files with 128 additions and 121 deletions.
44 changes: 44 additions & 0 deletions src/coreComponents/mainInterface/ProblemManager.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -685,6 +685,9 @@ void ProblemManager::generateMesh()
domain.setupCommunications( useNonblockingMPI );
domain.outputPartitionInformation();

// Create Embedded fractures here.
generateEmbeddedFractures();

domain.forMeshBodies( [&]( MeshBody & meshBody )
{
if( meshBody.hasGroup( keys::particleManager ) )
Expand Down Expand Up @@ -728,6 +731,47 @@ void ProblemManager::generateMesh()

}

void ProblemManager::generateEmbeddedFractures()
{
DomainPartition & domain = getDomainPartition();
MeshManager & meshManager = this->getGroup< MeshManager >( groupKeys.meshManager );

meshManager.generateEmbeddedFractures( domain );

MeshBody & meshBody = domain.getMeshBody( 0 );
CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager();
Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks();
string const & faceBlockName = embeddedSurfaceRegion.getFaceBlockName();

if( embSurfBlocks.hasGroup( faceBlockName ))
{
EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( faceBlockName );

elemManager.forElementSubRegionsComplete< CellElementSubRegion >( [&]( localIndex const er,
localIndex const esr,
ElementRegionBase &,
CellElementSubRegion & subRegion )
{
FixedOneToManyRelation const & cellToEdges = subRegion.edgeList();

embeddedSurfaceSubRegion.copyFromEmbeddedSurfaceBlock( er,
esr,
nodeManager,
embSurfNodeManager,
edgeManager,
cellToEdges,
embSurf );

// Add all the fracture information to the CellElementSubRegion
for( localIndex edfmIndex=0; edfmIndex < embSurf.numEmbeddedSurfElem(); ++edfmIndex )
{
localIndex cellIndex = embSurf.getEmbeddedSurfElemTo3dElem().toCellIndex[edfmIndex][0];
subRegion.addFracturedElement( cellIndex, edfmIndex );
newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( edfmIndex );
}
} ); // end loop over subregions
}
}

void ProblemManager::importFields()
{
Expand Down
23 changes: 7 additions & 16 deletions src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,11 +27,6 @@
#include "BufferOps.hpp"
#include "mesh/MeshFields.hpp"


#include "mainInterface/GeosxState.hpp"
#include "mainInterface/ProblemManager.hpp"
#include "mesh/DomainPartition.hpp"

namespace geos
{
using namespace dataRepository;
Expand Down Expand Up @@ -345,14 +340,13 @@ array1d< localIndex > EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex( Arra
return parentEdgeIndex;
}

bool EmbeddedSurfaceSubRegion::addAllEmbeddedSurfaces(
localIndex const regionIndex,
localIndex const subRegionIndex,
NodeManager const & nodeManager,
EmbeddedSurfaceNodeManager & embSurfNodeManager,
EdgeManager const & edgeManager,
FixedOneToManyRelation const & cellToEdges,
EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock )
void EmbeddedSurfaceSubRegion::copyFromEmbeddedSurfaceBlock( localIndex const regionIndex,
localIndex const subRegionIndex,
NodeManager const & nodeManager,
EmbeddedSurfaceNodeManager & embSurfNodeManager,
EdgeManager const & edgeManager,
FixedOneToManyRelation const & cellToEdges,
EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock )
{

arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoord = nodeManager.referencePosition();
Expand Down Expand Up @@ -419,9 +413,6 @@ bool EmbeddedSurfaceSubRegion::addAllEmbeddedSurfaces(
LvArray::tensorOps::copy< 3 >( m_tangentVector2[i], elemTangentialVectors2[i] );
this->calculateElementGeometricQuantities( elemNodeCoords.toViewConst(), i );
}

return true;

}

void EmbeddedSurfaceSubRegion::inheritGhostRank( array1d< array1d< arrayView1d< integer const > > > const & cellGhostRank )
Expand Down
12 changes: 12 additions & 0 deletions src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,6 +290,18 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion

///@}

/**
* @brief Fill the EmbeddedSurfaceSubRegion by copying those of the source CellBlock
* @param embeddedSurfaceBlock the CellBlEmbeddedSurfaceBlock which properties (connectivity info) will be copied.
*/
void copyFromEmbeddedSurfaceBlock( localIndex const regionIndex,
localIndex const subRegionIndex,
NodeManager const & nodeManager,
EmbeddedSurfaceNodeManager & embSurfNodeManager,
EdgeManager const & edgeManager,
FixedOneToManyRelation const & cellToEdges,
EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock );

private:

/**
Expand Down
170 changes: 65 additions & 105 deletions src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ EmbeddedSurfaceGenerator::EmbeddedSurfaceGenerator( const string & name,

registerWrapper( viewKeyStruct::targetObjectsNameString(), &m_targetObjectsName ).
setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ).
setInputFlag( dataRepository::InputFlags::OPTIONAL ).
setInputFlag( dataRepository::InputFlags::REQUIRED ).
setDescription( "List of geometric objects that will be used to initialized the embedded surfaces/fractures." );

registerWrapper( viewKeyStruct::mpiCommOrderString(), &m_mpiCommOrder ).
Expand Down Expand Up @@ -130,121 +130,81 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups()

NewObjectLists newObjects;

if( m_targetObjectsName.empty() ) // when targetObjects were not provided - try to import from vtk
// Loop over all the fracture planes
geometricObjManager.forGeometricObject< PlanarGeometricObject >( m_targetObjectsName, [&]( localIndex const,
PlanarGeometricObject & fracture )
{
MeshBody & meshBody = domain.getMeshBody( 0 );
CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager();
Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks();
string const & faceBlockName = embeddedSurfaceRegion.getFaceBlockName();
if( embSurfBlocks.hasGroup( faceBlockName ))
/* 1. Find out if an element is cut by the fracture or not.
* Loop over all the elements and for each one of them loop over the nodes and compute the
* dot product between the distance between the plane center and the node and the normal
* vector defining the plane. If two scalar products have different signs the plane cuts the
* cell. If a nodes gives a 0 dot product it has to be neglected or the method won't work.
*/
real64 const planeCenter[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( fracture.getCenter() );
real64 const normalVector[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( fracture.getNormal() );
// Initialize variables
globalIndex nodeIndex;
integer isPositive, isNegative;
real64 distVec[ 3 ];

elemManager.forElementSubRegionsComplete< CellElementSubRegion >(
[&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion )
{
EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( faceBlockName );

elemManager.forElementSubRegionsComplete< CellElementSubRegion >(
[&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion )
{
FixedOneToManyRelation const & cellToEdges = subRegion.edgeList();
arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList();
FixedOneToManyRelation const & cellToEdges = subRegion.edgeList();

bool added = embeddedSurfaceSubRegion.addAllEmbeddedSurfaces( er,
esr,
nodeManager,
embSurfNodeManager,
edgeManager,
cellToEdges,
embSurf );
arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

if( added )
{
// Add all the fracture information to the CellElementSubRegion
for( localIndex edfmIndex=0; edfmIndex < embSurf.numEmbeddedSurfElem(); ++edfmIndex )
{
localIndex cellIndex = embSurf.getEmbeddedSurfElemTo3dElem().toCellIndex[edfmIndex][0];
subRegion.addFracturedElement( cellIndex, edfmIndex );
newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( edfmIndex );
}
}
} ); // end loop over subregions
}
}
else
{
// Loop over all the fracture planes
geometricObjManager.forGeometricObject< PlanarGeometricObject >( m_targetObjectsName, [&]( localIndex const,
PlanarGeometricObject & fracture )
{
/* 1. Find out if an element is cut by the fracture or not.
* Loop over all the elements and for each one of them loop over the nodes and compute the
* dot product between the distance between the plane center and the node and the normal
* vector defining the plane. If two scalar products have different signs the plane cuts the
* cell. If a nodes gives a 0 dot product it has to be neglected or the method won't work.
*/
real64 const planeCenter[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( fracture.getCenter() );
real64 const normalVector[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3( fracture.getNormal() );
// Initialize variables
globalIndex nodeIndex;
integer isPositive, isNegative;
real64 distVec[ 3 ];

elemManager.forElementSubRegionsComplete< CellElementSubRegion >(
[&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion )
forAll< serialPolicy >( subRegion.size(), [ &, ghostRank,
cellToNodes,
nodesCoord ] ( localIndex const cellIndex )
{
arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList();
FixedOneToManyRelation const & cellToEdges = subRegion.edgeList();

arrayView1d< integer const > const ghostRank = subRegion.ghostRank();

forAll< serialPolicy >( subRegion.size(), [ &, ghostRank,
cellToNodes,
nodesCoord ] ( localIndex const cellIndex )
if( ghostRank[cellIndex] < 0 )
{
if( ghostRank[cellIndex] < 0 )
isPositive = 0;
isNegative = 0;
for( localIndex kn = 0; kn < subRegion.numNodesPerElement(); kn++ )
{
isPositive = 0;
isNegative = 0;
for( localIndex kn = 0; kn < subRegion.numNodesPerElement(); kn++ )
nodeIndex = cellToNodes[cellIndex][kn];
LvArray::tensorOps::copy< 3 >( distVec, nodesCoord[nodeIndex] );
LvArray::tensorOps::subtract< 3 >( distVec, planeCenter );
// check if the dot product is zero
if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) > 0 )
{
nodeIndex = cellToNodes[cellIndex][kn];
LvArray::tensorOps::copy< 3 >( distVec, nodesCoord[nodeIndex] );
LvArray::tensorOps::subtract< 3 >( distVec, planeCenter );
// check if the dot product is zero
if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) > 0 )
{
isPositive = 1;
}
else if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) < 0 )
{
isNegative = 1;
}
} // end loop over nodes

if( isPositive * isNegative == 1 )
isPositive = 1;
}
else if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) < 0 )
{
bool added = embeddedSurfaceSubRegion.addNewEmbeddedSurface( cellIndex,
er,
esr,
nodeManager,
embSurfNodeManager,
edgeManager,
cellToEdges,
&fracture );

if( added )
{
GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, "Element " << cellIndex << " is fractured" );

// Add the information to the CellElementSubRegion
subRegion.addFracturedElement( cellIndex, localNumberOfSurfaceElems );

newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( localNumberOfSurfaceElems );

localNumberOfSurfaceElems++;
}
isNegative = 1;
}
} // end loop over nodes
if( isPositive * isNegative == 1 )
{
bool added = embeddedSurfaceSubRegion.addNewEmbeddedSurface( cellIndex,
er,
esr,
nodeManager,
embSurfNodeManager,
edgeManager,
cellToEdges,
&fracture );

if( added )
{
GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, "Element " << cellIndex << " is fractured" );

// Add the information to the CellElementSubRegion
subRegion.addFracturedElement( cellIndex, localNumberOfSurfaceElems );

newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( localNumberOfSurfaceElems );

localNumberOfSurfaceElems++;
}
}
} );// end loop over cells
} );// end loop over subregions
} );// end loop over planes
}
}
} );// end loop over cells
} );// end loop over subregions
} );// end loop over planes

// Launch kernel to compute connectivity index of each fractured element.
elemManager.forElementSubRegionsComplete< CellElementSubRegion >(
Expand Down

0 comments on commit b1fa4c4

Please sign in to comment.