diff --git a/inputFiles/singlePhaseFlowFractures/edfm_vtk_fractures.vtu b/inputFiles/singlePhaseFlowFractures/edfm_vtk_fractures.vtu new file mode 100644 index 00000000000..abc49433c04 --- /dev/null +++ b/inputFiles/singlePhaseFlowFractures/edfm_vtk_fractures.vtu @@ -0,0 +1,78 @@ + + + + + + + 0 + 1 + 2 + 3 + 4 + 5 + + + + + 18 19 + + + 18 4 + 19 13 + + + 0.1 + 0.1 + + + 0.0001 + 0.0001 + + + 0 0 1 + 0 0 1 + + + 0 1 0 + 0 1 0 + + + 1 0 0 + 1 0 0 + + + + + 0.5 0.3333 0 + 0.5 0.6666 0 + 0.5 0.3333 0.5 + 0.5 0.6666 0.5 + 0.5 0.3333 1 + 0.5 0.6666 1 + + + 0 + + + 1 + + + + + + + 1 0 2 3 + 3 2 4 5 + + + 4 + 8 + + + 9 + 9 + + + + + diff --git a/inputFiles/singlePhaseFlowFractures/edfm_vtk_matrix.vtu b/inputFiles/singlePhaseFlowFractures/edfm_vtk_matrix.vtu new file mode 100644 index 00000000000..9d43724aea3 --- /dev/null +++ b/inputFiles/singlePhaseFlowFractures/edfm_vtk_matrix.vtu @@ -0,0 +1,122 @@ + + + + + + 0 1 2 3 4 5 + 6 7 8 9 10 11 + 12 13 14 15 16 17 + 18 19 20 21 22 23 + 24 25 26 27 28 29 + 30 31 32 33 34 35 + 36 37 38 39 40 41 + 42 43 44 45 46 47 + + + + + 0 0 0 0 0 0 + 0 0 0 0 0 0 + 0 0 0 0 0 0 + + + 0 1 2 3 4 5 + 6 7 8 9 10 11 + 12 13 14 15 16 17 + + + + + 0 0 0 + 0.333 0 0 + 0.666 0 0 + 1 0 0 + 0 0.333 0 + 0.333 0.333 0 + 0.666 0.333 0 + 1 0.333 0 + 0 0.666 0 + 0.333 0.666 0 + 0.666 0.666 0 + 1 0.666 0 + 0 1 0 + 0.333 1 0 + 0.666 1 0 + 1 1 0 + 0 0 0.5 + 0.333 0 0.5 + 0.666 0 0.5 + 1 0 0.5 + 0 0.333 0.5 + 0.333 0.333 0.5 + 0.666 0.333 0.5 + 1 0.333 0.5 + 0 0.666 0.5 + 0.333 0.666 0.5 + 0.666 0.666 0.5 + 1 0.666 0.5 + 0 1 0.5 + 0.333 1 0.5 + 0.666 1 0.5 + 1 1 0.5 + 0 0 1 + 0.333 0 1 + 0.666 0 1 + 1 0 1 + 0 0.333 1 + 0.333 0.333 1 + 0.666 0.333 1 + 1 0.333 1 + 0 0.666 1 + 0.333 0.666 1 + 0.666 0.666 1 + 1 0.666 1 + 0 1 1 + 0.333 1 1 + 0.666 1 1 + 1 1 1 + + + 0 + + + 1.7320508076 + + + + + + + 0 1 5 4 16 17 21 20 + 1 2 6 5 17 18 22 21 + 2 3 7 6 18 19 23 22 + 4 5 9 8 20 21 25 24 + 5 6 10 9 21 22 26 25 + 6 7 11 10 22 23 27 26 + 8 9 13 12 24 25 29 28 + 9 10 14 13 25 26 30 29 + 10 11 15 14 26 27 31 30 + 16 17 21 20 32 33 37 36 + 17 18 22 21 33 34 38 37 + 18 19 23 22 34 35 39 38 + 20 21 25 24 36 37 41 40 + 21 22 26 25 37 38 42 41 + 22 23 27 26 38 39 43 42 + 24 25 29 28 40 41 45 44 + 25 26 30 29 41 42 46 45 + 26 27 31 30 42 43 47 46 + + + 8 16 24 32 40 48 + 56 64 72 80 88 96 + 104 112 120 128 136 144 + + + 12 12 12 12 12 12 + 12 12 12 12 12 12 + 12 12 12 12 12 12 + + + + + diff --git a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk.xml new file mode 100644 index 00000000000..e147b352089 --- /dev/null +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk.xml @@ -0,0 +1,137 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk_smoke.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk_smoke.xml new file mode 100644 index 00000000000..bc1a60fa14e --- /dev/null +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk_smoke.xml @@ -0,0 +1,61 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlowFractures/main_edfm_vtk.vtm b/inputFiles/singlePhaseFlowFractures/main_edfm_vtk.vtm new file mode 100644 index 00000000000..1e5eda4c330 --- /dev/null +++ b/inputFiles/singlePhaseFlowFractures/main_edfm_vtk.vtm @@ -0,0 +1,6 @@ + + + + + + diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp index 652c1d3193d..7cc0ecfc9d0 100644 --- a/src/coreComponents/mainInterface/ProblemManager.cpp +++ b/src/coreComponents/mainInterface/ProblemManager.cpp @@ -694,6 +694,8 @@ void ProblemManager::generateMesh() else if( meshBody.hasGroup( keys::cellManager ) ) { // meshBody.deregisterGroup( keys::cellManager ); + // Cell block manager is needed for EDFM loading from VTK + // Can't be deregistred here meshBody.deregisterCellBlockManager(); } diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt index 05f09e22dc8..080d86dd8c8 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -79,6 +79,8 @@ set( mesh_headers generators/ExternalMeshGeneratorBase.hpp generators/FaceBlock.hpp generators/FaceBlockABC.hpp + generators/EmbeddedSurfaceBlock.hpp + generators/EmbeddedSurfaceBlockABC.hpp generators/InternalMeshGenerator.hpp generators/InternalWellGenerator.hpp generators/InternalWellboreGenerator.hpp @@ -158,6 +160,7 @@ set( mesh_sources generators/LineBlock.cpp generators/ExternalMeshGeneratorBase.cpp generators/FaceBlock.cpp + generators/EmbeddedSurfaceBlock.cpp generators/InternalMeshGenerator.cpp generators/InternalWellGenerator.cpp generators/InternalWellboreGenerator.cpp @@ -191,6 +194,7 @@ if( ENABLE_VTK ) set( mesh_headers ${mesh_headers} generators/CollocatedNodes.hpp generators/VTKFaceBlockUtilities.hpp + generators/VTKEmbeddedSurfaceBlockUtilities.hpp generators/VTKHierarchicalDataSource.hpp generators/VTKMeshGenerator.hpp generators/VTKMeshGeneratorTools.hpp @@ -200,6 +204,7 @@ if( ENABLE_VTK ) set( mesh_sources ${mesh_sources} generators/CollocatedNodes.cpp generators/VTKFaceBlockUtilities.cpp + generators/VTKEmbeddedSurfaceBlockUtilities.cpp generators/VTKHierarchicalDataSource.cpp generators/VTKMeshGenerator.cpp generators/VTKMeshGeneratorTools.cpp @@ -241,4 +246,3 @@ endif() if( GEOS_ENABLE_TESTS ) add_subdirectory( unitTests ) endif( ) - diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index 18da04ec066..cb7a2d2d04d 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -122,7 +122,6 @@ void ElementRegionManager::setSchemaDeviations( xmlWrapper::xmlNode schemaRoot, } } - void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockManager ) { { // cellBlocks loading @@ -141,7 +140,16 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa this->forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion & elemRegion ) { - elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); + + if( elemRegion.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::faceElement ) + { + elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); + } + else if( elemRegion.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::embeddedElement ) + { + elemRegion.generateMesh( cellBlockManager.getEmbeddedSurfaceBlocks() ); + } + } ); // Some mappings of the surfaces subregions point to elements in other subregions and regions. diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index 875d23b289f..c4a1993dbba 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -27,6 +27,11 @@ #include "BufferOps.hpp" #include "mesh/MeshFields.hpp" + +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/DomainPartition.hpp" + namespace geos { using namespace dataRepository; @@ -102,12 +107,11 @@ void EmbeddedSurfaceSubRegion::calculateElementGeometricQuantities( arrayView2d< { LvArray::tensorOps::add< 3 >( m_elementCenter[k], intersectionPoints[ p ] ); } + LvArray::tensorOps::scale< 3 >( m_elementCenter[ k ], 1.0 / intersectionPoints.size( 0 ) ); // update area m_elementArea[ k ] = computationalGeometry::ComputeSurfaceArea( intersectionPoints, m_normalVector[k] ); - LvArray::tensorOps::scale< 3 >( m_elementCenter[ k ], 1.0 / intersectionPoints.size( 0 ) ); - // update volume m_elementVolume[k] = m_elementAperture[k] * m_elementArea[k]; } @@ -277,6 +281,148 @@ bool EmbeddedSurfaceSubRegion::addNewEmbeddedSurface( localIndex const cellIndex } return addEmbeddedElem; } +array1d< localIndex > EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex( ArrayOfArraysView< real64 > const & elemNodesLocations, + ArrayOfArraysView< localIndex > const & elemToNodes, + ToCellRelation< ArrayOfArrays< localIndex > > const & elemTo3dElem, + FixedOneToManyRelation const & cellToEdges, + arrayView2d< localIndex const > const & edgeToNodes, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoord ) +{ + // first we build a node index to parent cell index map (we only need one) + array1d< int > nodeToElem( elemNodesLocations.size()); + for( int i = 0; i parentEdgeIndex( elemNodesLocations.size()); + for( localIndex i = 0; i < elemNodesLocations.size(); ++i ) + { + + auto elementIndex =nodeToElem[i]; + auto cell3dIndex = elemTo3dElem.toCellIndex[elementIndex][0]; + double elemMinDistance = std::numeric_limits< double >::max(); + localIndex targetEdgeIndex = -1; + + R1Tensor nodexyz; + LvArray::tensorOps::copy< 3 >( nodexyz, elemNodesLocations[i] ); + for( localIndex ke = 0; ke < cellToEdges.size( 1 ); ke++ ) + { + R1Tensor v1; + R1Tensor v2; + R1Tensor crossProd; + localIndex edgeIndex = cellToEdges[cell3dIndex][ke]; + LvArray::tensorOps::copy< 3 >( v1, nodesCoord[edgeToNodes[edgeIndex][0]] ); + LvArray::tensorOps::subtract< 3 >( v1, nodesCoord[edgeToNodes[edgeIndex][1]] ); + + LvArray::tensorOps::copy< 3 >( v2, nodesCoord[edgeToNodes[edgeIndex][1]] ); + LvArray::tensorOps::subtract< 3 >( v2, nodexyz ); + + LvArray::tensorOps::crossProduct( crossProd, v1, v2 ); + auto numeratorNorm = LvArray::tensorOps::l2Norm< 3 >( crossProd ); + auto denomNorm = LvArray::tensorOps::l2Norm< 3 >( v1 ); + auto dist = numeratorNorm / denomNorm; + if( dist < elemMinDistance ) + { + elemMinDistance = dist; + targetEdgeIndex = edgeIndex; + } + } + + parentEdgeIndex[i] = targetEdgeIndex; + } + 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 ) +{ + + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoord = nodeManager.referencePosition(); + arrayView2d< localIndex const > const edgeToNodes = edgeManager.nodeList(); + arrayView1d< integer const > const edgeGhostRank = edgeManager.ghostRank(); + arrayView1d< globalIndex const > const & edgeLocalToGlobal = edgeManager.localToGlobalMap(); + + localIndex const numElems = embeddedSurfaceBlock.numEmbeddedSurfElem(); + + resize( embeddedSurfaceBlock.numEmbeddedSurfElem()); + + ToCellRelation< ArrayOfArrays< localIndex > > elemTo3dElem = embeddedSurfaceBlock.getEmbeddedSurfElemTo3dElem(); + ArrayOfArrays< localIndex > elemToNodes = embeddedSurfaceBlock.getEmbeddedSurfElemToNodes(); + ArrayOfArrays< real64 > elemNodesLocations = embeddedSurfaceBlock.getEmbeddedSurfElemNodesCoords(); + ArrayOfArrays< real64 > elemNormalVectors = embeddedSurfaceBlock.getEmbeddedSurfElemNormalVectors(); + ArrayOfArrays< real64 > elemTangentialVectors1 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialWidthVectors(); + ArrayOfArrays< real64 > elemTangentialVectors2 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialLengthVectors(); + + localIndex numEdfmNodes = elemNodesLocations.size(); + array1d< localIndex > allPointsParentIndex = getEdfmNodeParentEdgeIndex( elemNodesLocations.toView(), elemToNodes.toView(), elemTo3dElem, cellToEdges, edgeToNodes, nodesCoord ); + + // EDFM nodes get their gost rank and global parent index + // from their parent 3d element local edge index. + array1d< integer > allPointsGhostRank( numEdfmNodes ); + array1d< integer > allPointsGlobalParenIndex( numEdfmNodes ); + for( auto i = 0; i < numEdfmNodes; ++i ) + { + integer pointGhostRank = edgeGhostRank[allPointsParentIndex[i]]; + localIndex pointLocalParentEdgeIndex = allPointsParentIndex[i]; + globalIndex pointGlobalParentEdgeIndex = edgeLocalToGlobal[allPointsParentIndex[i]]; + + embSurfNodeManager.appendNode( elemNodesLocations[i].toSliceConst(), pointGhostRank ); + + arrayView1d< localIndex > const & parentIndex = embSurfNodeManager.getField< fields::parentEdgeIndex >(); + + parentIndex[i] = pointLocalParentEdgeIndex; + + array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); + parentEdgeGlobalIndex[i] = pointGlobalParentEdgeIndex; + } + + for( auto i = 0; i < numElems; ++i ) + { + m_toNodesRelation.resizeArray( i, 4 ); + localIndex elem3dIndex = elemTo3dElem.toCellIndex[i][0]; + //localIndex elemRegionIndex = elemTo3dElem.toBlockIndex[i][0]; + // region and subregion will be filled later. + m_2dElemToElems.m_toElementIndex.emplaceBack( i, elem3dIndex ); + m_2dElemToElems.m_toElementSubRegion.emplaceBack( i, subRegionIndex ); + m_2dElemToElems.m_toElementRegion.emplaceBack( i, regionIndex ); + + // we only accept quads + array2d< real64 > elemNodeCoords( 4, 3 ); + for( localIndex inode = 0; inode < 4; inode++ ) + { + auto nodeGlobalIndex = elemToNodes[i][inode]; + m_toNodesRelation( i, inode ) = nodeGlobalIndex; + LvArray::tensorOps::copy< 3 >( elemNodeCoords[inode], elemNodesLocations[nodeGlobalIndex] ); + } + + m_parentPlaneName[i] = "elem_" + std::to_string( i ); + LvArray::tensorOps::copy< 3 >( m_normalVector[i], elemNormalVectors[i] ); + LvArray::tensorOps::copy< 3 >( m_tangentVector1[i], elemTangentialVectors1[i] ); + 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 ) { diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp index 5baa73ec274..249862f4835 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp @@ -26,6 +26,7 @@ #include "EdgeManager.hpp" #include "EmbeddedSurfaceNodeManager.hpp" #include "CellElementSubRegion.hpp" +#include "mesh/generators/EmbeddedSurfaceBlockABC.hpp" //Do we really need this include Rectangle? #include "simpleGeometricObjects/Rectangle.hpp" @@ -146,6 +147,44 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion FixedOneToManyRelation const & cellToEdges, PlanarGeometricObject const * fracture ); + /** + * @brief Function that maps each node of an embedded surface element to its parent edge on the 3D element + * @param elemNodesLocations geometric coordinates of each embedded surface element node (no duplicates for shared nodes) + * @param elemToNodes embedded surface node index to the array of its node indices map + * @param elemTo3dElem embedded element to its parent 3D cell element map + * @param cellToEdges 3D cells to its edges map + * @param edgeToNodes edge to its node indices map + * @param nodesCoord geometric coordinates of 3D cells + * @return embedded node to 3D cell edge mapping + */ + array1d< localIndex > getEdfmNodeParentEdgeIndex( ArrayOfArraysView< real64 > const & elemNodesLocations, + ArrayOfArraysView< localIndex > const & elemToNodes, + ToCellRelation< ArrayOfArrays< localIndex > > const & elemTo3dElem, + FixedOneToManyRelation const & cellToEdges, + arrayView2d< localIndex const > const & edgeToNodes, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & nodesCoord ); + + /** + * @brief Function to add a all the embedded surface elements (loaded from VTK). + * @param regionIndex cell element region index + * @param subRegionIndex cell element subregion index + * @param nodeManager the nodemanager group + * @param embSurfNodeManager the embSurfNodeManager group + * @param edgeManager the edgemanager group + * @param cellToEdges cellElement to edges map + * @param embeddedSurfaceBlock embeddedSurfaceBloc holds the information for all embedded element network + * @return boolean defining whether all the embedded elements were added or not + */ + + bool addAllEmbeddedSurfaces( localIndex const regionIndex, + localIndex const subRegionIndex, + NodeManager const & nodeManager, + EmbeddedSurfaceNodeManager & embSurfNodeManager, + EdgeManager const & edgeManager, + FixedOneToManyRelation const & cellToEdges, + EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock ); + + /** * @brief inherit ghost rank from cell elements. * @param cellGhostRank cell element ghost ranks diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index e6b47dbae4d..dd4c7b7dbc3 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -22,7 +22,6 @@ #include "SurfaceElementRegion.hpp" #include "common/MpiWrapper.hpp" - namespace geos { using namespace dataRepository; @@ -55,8 +54,10 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) { Group & elementSubRegions = this->getGroup( viewKeyStruct::elementSubRegions() ); + if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) { + // We just register the subregion copying of data is done at the EmbeddedSurfaceGenerator elementSubRegions.registerGroup< EmbeddedSurfaceSubRegion >( m_faceBlockName ); } else if( m_subRegionType == SurfaceSubRegionType::faceElement ) diff --git a/src/coreComponents/mesh/SurfaceElementRegion.hpp b/src/coreComponents/mesh/SurfaceElementRegion.hpp index f15de2d63f6..705938958bf 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.hpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.hpp @@ -134,6 +134,12 @@ class SurfaceElementRegion : public ElementRegionBase */ SurfaceSubRegionType subRegionType() const { return m_subRegionType; } + /** + * @brief Get face block name. + * @return face block name + */ + string const & getFaceBlockName() const { return m_faceBlockName; } + /** * @brief Returns the unique sub-region of type @p SUBREGION_TYPE for the current @p SurfaceElementRegion. * @tparam SUBREGION_TYPE The type of the sub region we're looking for. diff --git a/src/coreComponents/mesh/generators/CellBlock.hpp b/src/coreComponents/mesh/generators/CellBlock.hpp index 6fa0de501a9..e51f2fed324 100644 --- a/src/coreComponents/mesh/generators/CellBlock.hpp +++ b/src/coreComponents/mesh/generators/CellBlock.hpp @@ -196,6 +196,20 @@ class CellBlock : public CellBlockABC arrayView1d< globalIndex const > localToGlobalMapConstView() const { return m_localToGlobalMap.toViewConst(); } + /** + * @brief Get global to local map. + * @return The mapping relationship. + */ + unordered_map< globalIndex, localIndex > & globalToLocalMap() + { return m_globalToLocalMap; } + + /** + * @brief Get global to local map, const version. + * @return The mapping relationship. + */ + unordered_map< globalIndex, localIndex > const & globalToLocalMap() const + { return m_globalToLocalMap; } + /** * @brief Resize the cell block to hold @p numElements * @param numElements The new number of elements. @@ -242,6 +256,9 @@ class CellBlock : public CellBlockABC /// Contains the global index of each object. array1d< globalIndex > m_localToGlobalMap; + /// Map from object global index to the local index. + unordered_map< globalIndex, localIndex > m_globalToLocalMap; + /// Name of the properties registered from an external mesh string_array m_externalPropertyNames; diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index 2dd41c0987a..284c27bdbfe 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.cpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp @@ -32,6 +32,7 @@ CellBlockManager::CellBlockManager( string const & name, Group * const parent ): this->registerGroup< Group >( viewKeyStruct::cellBlocks() ); this->registerGroup< Group >( viewKeyStruct::faceBlocks() ); this->registerGroup< Group >( viewKeyStruct::lineBlocks() ); + this->registerGroup< Group >( viewKeyStruct::embeddedSurfaceBlocks() ); } void CellBlockManager::resize( integer_array const & numElements, @@ -707,6 +708,15 @@ Group & CellBlockManager::getFaceBlocks() return this->getGroup( viewKeyStruct::faceBlocks() ); } +Group const & CellBlockManager::getEmbeddedSurfaceBlocks() const +{ + return this->getGroup( viewKeyStruct::embeddedSurfaceBlocks() ); +} + +Group & CellBlockManager::getEmbeddedSurfaceBlocks() +{ + return this->getGroup( viewKeyStruct::embeddedSurfaceBlocks() ); +} Group & CellBlockManager::getLineBlocks() { return this->getGroup( viewKeyStruct::lineBlocks() ); @@ -787,6 +797,12 @@ LineBlock & CellBlockManager::registerLineBlock( string const & name ) return this->getLineBlocks().registerGroup< LineBlock >( name ); } +EmbeddedSurfaceBlock & CellBlockManager::registerEmbeddedSurfaceBlock( string const & name ) +{ + + return this->getEmbeddedSurfaceBlocks().registerGroup< EmbeddedSurfaceBlock >( name ); +} + array2d< real64, nodes::REFERENCE_POSITION_PERM > CellBlockManager::getNodePositions() const { return m_nodesPositions; diff --git a/src/coreComponents/mesh/generators/CellBlockManager.hpp b/src/coreComponents/mesh/generators/CellBlockManager.hpp index 4e8d2688564..80caf28087b 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.hpp @@ -25,6 +25,7 @@ #include "mesh/generators/InternalWellGenerator.hpp" #include "mesh/generators/LineBlock.hpp" #include "mesh/generators/LineBlockABC.hpp" +#include "mesh/generators/EmbeddedSurfaceBlock.hpp" #include "mesh/generators/CellBlockManagerABC.hpp" #include "mesh/generators/PartitionDescriptor.hpp" @@ -165,6 +166,10 @@ class CellBlockManager : public CellBlockManagerABC Group & getFaceBlocks() override; + Group const & getEmbeddedSurfaceBlocks() const override; + + Group & getEmbeddedSurfaceBlocks() override; + LineBlockABC const & getLineBlock( string name ) const override; /** @@ -195,6 +200,12 @@ class CellBlockManager : public CellBlockManagerABC * @return A reference to the new line block. The CellBlockManager owns this new instance. */ LineBlock & registerLineBlock( string const & name ); + /** + * @brief Registers and returns an embedded surface block of name @p name. + * @param name The name of the created embedded surface block. + * @return A reference to the new embedded surface block. The CellBlockManager owns this new instance. + */ + EmbeddedSurfaceBlock & registerEmbeddedSurfaceBlock( string const & name ); /** * @brief Launch kernel function over all the sub-regions * @tparam LAMBDA type of the user-provided function @@ -214,6 +225,15 @@ class CellBlockManager : public CellBlockManagerABC */ void setGlobalLength( real64 globalLength ) { m_globalLength = globalLength; } + /** + * @brief Get cell block at index @p blockIndex. + * @param[in] blockIndex The cell block index. + * @return Const reference to the instance. + * + * @note Mainly useful for iteration purposes. + */ + CellBlock const & getCellBlock( localIndex const blockIndex ) const; + private: struct viewKeyStruct @@ -229,6 +249,9 @@ class CellBlockManager : public CellBlockManagerABC /// Line blocks key static constexpr char const * lineBlocks() { return "lineBlocks"; } + /// Embedded Surface blocks key + static constexpr char const * embeddedSurfaceBlocks() + { return "embeddedSurfaceBlocks"; } }; /** @@ -239,15 +262,6 @@ class CellBlockManager : public CellBlockManagerABC */ Group & getLineBlocks(); - /** - * @brief Get cell block at index @p blockIndex. - * @param[in] blockIndex The cell block index. - * @return Const reference to the instance. - * - * @note Mainly useful for iteration purposes. - */ - CellBlock const & getCellBlock( localIndex const blockIndex ) const; - /** * @brief Returns the number of cells blocks * @return Number of cell blocks diff --git a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp index c68c2a1bf13..583071a28f2 100644 --- a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp @@ -101,6 +101,14 @@ class CellBlockManagerABC : public dataRepository::Group */ virtual Group & getFaceBlocks() = 0; + /** + * @brief Returns a group containing the embedded surface block as @p EmbeddedSurfaceBlockABC instances. + * @return Mutable reference to the embedded surface blocks group. + * + * @note It should probably be better not to expose a non-const accessor here. + */ + virtual Group & getEmbeddedSurfaceBlocks() = 0; + /** * @brief Returns LineBlockABC corresponding to the given identifier * @param name the name of the required LineBlockABC @@ -126,6 +134,12 @@ class CellBlockManagerABC : public dataRepository::Group */ virtual std::map< integer, std::set< string > > const & getRegionAttributesCellBlocks() const = 0; + /** + * @brief Returns a group containing the embedded surfaces blocks as EmbeddedSurfaceBlockABC instances + * @return Const reference to the Group instance. + */ + virtual Group const & getEmbeddedSurfaceBlocks() const = 0; + /** * @brief Total number of nodes across all the cell blocks. * @return The total number of nodes. diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp new file mode 100644 index 00000000000..46f7d713376 --- /dev/null +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp @@ -0,0 +1,114 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2020- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + + +#include "EmbeddedSurfaceBlock.hpp" + +namespace geos +{ + +localIndex EmbeddedSurfaceBlock::numEmbeddedSurfElem() const +{ + return m_numEmbeddedSurfaces; +} + +ArrayOfArrays< localIndex > EmbeddedSurfaceBlock::getEmbeddedSurfElemToNodes() const +{ + return m_embeddedSurfElemToNodes; +} + +ToCellRelation< ArrayOfArrays< localIndex > > EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const +{ + return m_embeddedSurfElemTo3dElem; +} + +ArrayOfArrays< real64 > EmbeddedSurfaceBlock::getEmbeddedSurfElemNodesCoords() const +{ + return m_embeddedSurfElemNodesCoords; +} + +ArrayOfArrays< real64 > EmbeddedSurfaceBlock::getEmbeddedSurfElemNormalVectors() const +{ + return m_embeddedSurfElemNormals; +} +ArrayOfArrays< real64 > EmbeddedSurfaceBlock::getEmbeddedSurfElemTangentialLengthVectors() const +{ + return m_embeddedSurfElemLengthVectors; +} +ArrayOfArrays< real64 > EmbeddedSurfaceBlock::getEmbeddedSurfElemTangentialWidthVectors() const +{ + return m_embeddedSurfElemWidthVectors; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemNormalVectors( ArrayOfArrays< real64 > && _normals ) +{ + m_embeddedSurfElemNormals = _normals; +} +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialLengthVectors( ArrayOfArrays< real64 > && _lengthVectors ) +{ + m_embeddedSurfElemLengthVectors = _lengthVectors; +} +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialWidthVectors( ArrayOfArrays< real64 > && _widthVectors ) +{ + m_embeddedSurfElemWidthVectors= _widthVectors; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemAperture( array1d< real64 > && _apertures ) +{ + m_embeddedSurfElemApertures= _apertures; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemPermeability( array1d< real64 > && _perms ) +{ + m_embeddedSurfElemPermeability= _perms; +} + +void EmbeddedSurfaceBlock::setNumEmbeddedSurfElem( localIndex _numEmbeddedSurfaces ) +{ + + m_numEmbeddedSurfaces = _numEmbeddedSurfaces; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemToNodes( ArrayOfArrays< localIndex > && _embeddedSurfElemToNodes ) +{ + m_embeddedSurfElemToNodes = _embeddedSurfElemToNodes; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem( ToCellRelation< ArrayOfArrays< localIndex > > && _embeddedSurfElemTo3dElem ) +{ + m_embeddedSurfElemTo3dElem = _embeddedSurfElemTo3dElem; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemNodes( ArrayOfArrays< real64 > && _embeddedSurfElemNodes ) +{ + m_embeddedSurfElemNodesCoords = _embeddedSurfElemNodes; +} + +unordered_map< globalIndex, localIndex > & EmbeddedSurfaceBlock::globalToLocalMap() +{ + return m_globalToLocalMap; +} + +unordered_map< globalIndex, localIndex > const & EmbeddedSurfaceBlock::globalToLocalMap() const +{ + return m_globalToLocalMap; +} + +void EmbeddedSurfaceBlock::setGlobalToLocalMap( unordered_map< globalIndex, localIndex > && g2l ) +{ + m_globalToLocalMap = g2l; +} + +} diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp new file mode 100644 index 00000000000..d1048113522 --- /dev/null +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -0,0 +1,142 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2020- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_EMBEDDEDSURFACEBLOCK_HPP +#define GEOS_EMBEDDEDSURFACEBLOCK_HPP + +#include "EmbeddedSurfaceBlockABC.hpp" + +namespace geos +{ + +/** + * @brief Simple implementation of the @p EmbeddedSurfaceBlockABC contract + * + */ + +class EmbeddedSurfaceBlock : public EmbeddedSurfaceBlockABC +{ +public: + /** + * @brief Constructor. + * @param[in] name Name of this EmbeddedSurfaceBlock. + * @param[in] parent Parent group. + */ + EmbeddedSurfaceBlock( string const & name, + Group * const parent ) + : + EmbeddedSurfaceBlockABC( name, parent ) + { } + + localIndex numEmbeddedSurfElem() const override; + ArrayOfArrays< localIndex > getEmbeddedSurfElemToNodes() const override; + ToCellRelation< ArrayOfArrays< localIndex > > getEmbeddedSurfElemTo3dElem() const override; + ArrayOfArrays< real64 > getEmbeddedSurfElemNodesCoords() const override; + ArrayOfArrays< real64 > getEmbeddedSurfElemNormalVectors() const override; + ArrayOfArrays< real64 > getEmbeddedSurfElemTangentialLengthVectors() const override; + ArrayOfArrays< real64 > getEmbeddedSurfElemTangentialWidthVectors() const override; + + + /** + * @brief Sets the number of embedded elements + * @param _numEmbeddedSurfElem the input value + */ + void setNumEmbeddedSurfElem( localIndex _numEmbeddedSurfElem ); + + /** + * @brief Sets the embedded elements to nodes mapping + * @param _embeddedSurfElemToNodes the input mapping. + */ + void setEmbeddedSurfElemToNodes( ArrayOfArrays< localIndex > && _embeddedSurfElemToNodes ); + + /** + * @brief Sets the embedded elements to 1 3d element mapping + * @param _embeddedSurfElemTo3dElem the input mapping. + */ + void setEmbeddedSurfElemTo3dElem( ToCellRelation< ArrayOfArrays< localIndex > > && _embeddedSurfElemTo3dElem ); + + /** + * @brief Sets the embedded elements nodes coordinates + * @param _embeddedSurfElemNodes the coordinates of embedded elements nodes. + */ + void setEmbeddedSurfElemNodes( ArrayOfArrays< real64 > && _embeddedSurfElemNodes ); + + /** + * @brief Sets the embedded elements normals + * @param _normals the coordinates of embedded elements normals. + */ + void setEmbeddedSurfElemNormalVectors( ArrayOfArrays< real64 > && _normals ); + + /** + * @brief Sets the embedded elements tangential length vectors (horizontal plan) + * @param _lengthVectors the coordinates of embedded elements tangential length vectors. + */ + void setEmbeddedSurfElemTangentialLengthVectors( ArrayOfArrays< real64 > && _lengthVectors ); + + /** + * @brief Sets the embedded elements tangential width vectors (vertical plan) + * @param _widthVectors the coordinates of embedded elements tangential width vectors. + */ + void setEmbeddedSurfElemTangentialWidthVectors( ArrayOfArrays< real64 > && _widthVectors ); + + /** + * @brief Sets the embedded elements Apertures + * @param _apertures the aperture values of embedded elements. + */ + void setEmbeddedSurfElemAperture( array1d< real64 > && _apertures ); + + /** + * @brief Sets the embedded elements Permeabilities + * @param _perms the permeability values of embedded elements. + */ + void setEmbeddedSurfElemPermeability( array1d< real64 > && _perms ); + + /** + * @brief Get global to local map. + * @return The mapping relationship. + */ + unordered_map< globalIndex, localIndex > & globalToLocalMap(); + + /** + * @brief Get global to local map, const version. + * @return The mapping relationship. + */ + unordered_map< globalIndex, localIndex > const & globalToLocalMap() const; + + /** + * @brief Set global to local map from a provided mapping. + * @param g2l provided mapping. + */ + void setGlobalToLocalMap( unordered_map< globalIndex, localIndex > && g2l ); + +private: + + localIndex m_numEmbeddedSurfaces; + ArrayOfArrays< localIndex > m_embeddedSurfElemToNodes; + ToCellRelation< ArrayOfArrays< localIndex > > m_embeddedSurfElemTo3dElem; + ArrayOfArrays< real64 > m_embeddedSurfElemNodesCoords; + ArrayOfArrays< real64 > m_embeddedSurfElemNormals; + ArrayOfArrays< real64 > m_embeddedSurfElemLengthVectors; + ArrayOfArrays< real64 > m_embeddedSurfElemWidthVectors; + array1d< real64 > m_embeddedSurfElemApertures; + array1d< real64 > m_embeddedSurfElemPermeability; + + /// Map from object global index to the local index. + unordered_map< globalIndex, localIndex > m_globalToLocalMap; + +}; + +} + +#endif //inlcude guard diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp new file mode 100644 index 00000000000..7ced1c7c089 --- /dev/null +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -0,0 +1,124 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 Total, S.A + * Copyright (c) 2020- GEOSX Contributors + * All right reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_EMBEDDEDSURFACEBLOCKABC_HPP +#define GEOS_EMBEDDEDSURFACEBLOCKABC_HPP + + +#include "CellBlockUtilities.hpp" +#include "dataRepository/Group.hpp" +#include "common/DataTypes.hpp" + + +namespace geos +{ + +/** + * @brief An embedded 2d element within a 3d element + * + * @details The @p EmbeddedSurfaceBlockABC represents an array of 2d elements (@e surfacic) embedded within + * a 3d elements grid (2d element can intersect one 3d element). The 2d element is assumed to be a quad. + * and 3d elements are assumed to be regular hexa. + * @details In this class, we'll use the term @e 2d @e element for the elements of the @p EmbeddedSurfaceBlockABC, + * which are geometrical quad surfaces (in 3d). + * In the same way, we'll use the wording @e 2d @e face + * to refer to the 2d boundaries of the @e 2d @e elements. + * The @e 2d @e face are geometrical segments (in 3d). + * + */ +class EmbeddedSurfaceBlockABC : public dataRepository::Group +{ +public: + /** + * @brief Constructor + * @param name The name of this Group + * @param parent The parent Group + */ + + EmbeddedSurfaceBlockABC( string const & name, + Group * const parent ): + Group( name, parent ) + { } + + /** + * @brief Get the number of embedded surface elements + * @return Number of embedded surface elements + * @details Return the number of embedded surface elements, each surface element + * can intersect 1 3d elements + */ + virtual localIndex numEmbeddedSurfElem() const = 0; + + /** + * @brief Get the indices of the nodes of all embedded surface elements + * @return The mapping of first dimension @p numEmbeddedSurfElem. + * Second dimension is 4 (quad), and represents the indices of each embedded surface element. + * + * @details each embedded surface element is supposed to have 4 nodes. Node indices for all the embedded + * surface elements are given by getEmbeddedSurfaceElementsNodes. This method returns the mapping between + * an embedded surface element and the 4 indices of its nodes. + */ + virtual ArrayOfArrays< localIndex > getEmbeddedSurfElemToNodes() const = 0; + + /** + * @brief Get the indices of the parent 3d element of each embedded surface element (1 elem) + * @return The mapping of first dimension @p numEmbeddedSurfElem. + * Second dimension 1 to 1 3d element that it intersects. + * 2d element numbering is local to the @p EmbeddedSurfaceBlockABC + * 3d element numbering is local to the @p CellBlockABC + * + * @details each embedded surface element intersects 1 3d element only. Element and block indices + * of this 3d element (ToCellRelation) are returned for each embedded surface element. + */ + virtual ToCellRelation< ArrayOfArrays< localIndex > > getEmbeddedSurfElemTo3dElem() const = 0; + + + /** + * @brief Get the x, y, and z coordinates of the embedded surface elements nodes. + * @return An array of x, y, z coordinates of all the embedded surface elements nodes. + * first dimension is a maximum of @p numEmbeddedSurfaceElements*4 as nodes can be shared between elements. + * Second dimension is 3, and represents the x, y and z. + */ + virtual ArrayOfArrays< real64 > getEmbeddedSurfElemNodesCoords() const = 0; + + /** + * @brief Get all the normal vectors for all the embedded surface elements. + * @return An array of x, y, z coordinates of all the embedded surface elements normals. + * first dimension is @p numEmbeddedSurfaceElements. + * Second dimension is 3, and represents the x, y and z coordinates of the normal vector. + */ + virtual ArrayOfArrays< real64 > getEmbeddedSurfElemNormalVectors() const = 0; + + + /** + * @brief Get all the tangential length vectors for all the embedded surface elements. + * @return An array of x, y, z coordinates of all the embedded surface elements tangential length vectors + * in the orthonormal basis with the normal vectors. For 2.5D case, this is along the horizontal direction. + * first dimension is @p numEmbeddedSurfaceElements. + * Second dimension is 3, and represents the x, y and z coordinates of the length vector. + */ + virtual ArrayOfArrays< real64 > getEmbeddedSurfElemTangentialLengthVectors() const = 0; + + /** + * @brief Get all the tangential width vectors for all the embedded surface elements. + * @return An array of x, y, z coordinates of all the embedded surface elements tangential width vectors + * in the orthonormal basis with the normal vectors. For 2.5D case, this is along vertical direction. + * first dimension is @p numEmbeddedSurfaceElements. + * Second dimension is 3, and represents the x, y and z coordinates of the width vector. + */ + virtual ArrayOfArrays< real64 > getEmbeddedSurfElemTangentialWidthVectors() const = 0; +}; + +} + +#endif // include guard diff --git a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp new file mode 100644 index 00000000000..78cfd84115e --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp @@ -0,0 +1,168 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2020- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +#include "VTKEmbeddedSurfaceBlockUtilities.hpp" +#include +#include +#include +#include "CellBlockManager.hpp" + +namespace geos::vtk +{ + +unordered_map< globalIndex, localIndex > buildGlobalToLocal( vtkIdTypeArray const * edfmMeshCellGlobalIds ) +{ + unordered_map< globalIndex, localIndex > g2l; + + vtkIdType const numEdfms = edfmMeshCellGlobalIds ? edfmMeshCellGlobalIds->GetNumberOfTuples() : 0; + + for( auto i = 0; i < numEdfms; ++i ) + { + // Note that `numEdfms` is zero if `edfmMeshCellGlobalIds` is 0 too. + // This prevents from calling a member on a null instance. + g2l[edfmMeshCellGlobalIds->GetValue( i )] = i; + } + + return g2l; +} + +void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName, + vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, + CellBlockManager & cellBlockManager ) +{ + EmbeddedSurfaceBlock & embeddedSurfBlock = cellBlockManager.registerEmbeddedSurfaceBlock( embeddedSurfaceBlockName ); + vtkIdType const numEdfmFracs = embeddedSurfaceMesh->GetNumberOfCells(); + + vtkIdTypeArray const * edfmMeshCellGlobalIds = vtkIdTypeArray::FastDownCast( embeddedSurfaceMesh->GetCellData()->GetGlobalIds() ); + + GEOS_ERROR_IF( numEdfmFracs > 0 && !edfmMeshCellGlobalIds, "GlobalIds needs to be provided for " << embeddedSurfaceBlockName ); + + // build and set global to local mapping to use later for 'fracture_to_parent_matrix_cell_mapping' + embeddedSurfBlock.setGlobalToLocalMap( buildGlobalToLocal( edfmMeshCellGlobalIds ) ); + + //1.Get EDFM vertices + vtkUnstructuredGrid * const grid = vtkUnstructuredGrid::SafeDownCast( embeddedSurfaceMesh ); + vtkPoints * const nodes = grid->GetPoints(); + vtkIdType const numNodes = nodes ? nodes->GetNumberOfPoints() : 0; + + ArrayOfArrays< localIndex > elem2dToNodes( numEdfmFracs ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto const val = grid->GetCell( i )->GetPointIds(); + elem2dToNodes.emplaceBack( i, val->GetId( 0 )); + elem2dToNodes.emplaceBack( i, val->GetId( 1 )); + elem2dToNodes.emplaceBack( i, val->GetId( 2 )); + elem2dToNodes.emplaceBack( i, val->GetId( 3 )); + } + + ArrayOfArrays< real64 > embeddedSurfaceElemNodes( LvArray::integerConversion< localIndex >( numNodes ) ); + for( vtkIdType i = 0; i < numNodes; ++i ) + { + double point[3]; + nodes->GetPoint( i, point ); + embeddedSurfaceElemNodes.emplaceBack( i, point[0] ); + embeddedSurfaceElemNodes.emplaceBack( i, point[1] ); + embeddedSurfaceElemNodes.emplaceBack( i, point[2] ); + } + + //2. Get EDFM normal vectors + ArrayOfArrays< real64 > embeddedSurfaceNormalVectors( numEdfmFracs ); + const char * const normal_str = "normal_vectors"; + vtkDataArray * const normals = grid->GetCellData()->GetArray( normal_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto const val = normals->GetTuple3( i ); + embeddedSurfaceNormalVectors.emplaceBack( i, val[0] ); + embeddedSurfaceNormalVectors.emplaceBack( i, val[1] ); + embeddedSurfaceNormalVectors.emplaceBack( i, val[2] ); + } + + //3. Get EDFM tangential length vectors (horizontal) + ArrayOfArrays< real64 > embeddedSurfaceTangentialLengthVectors( numEdfmFracs ); + const char * const tl_str = "tangential_length_vectors"; + vtkDataArray * const tangential_len = grid->GetCellData()->GetArray( tl_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto const val = tangential_len->GetTuple3( i ); + embeddedSurfaceTangentialLengthVectors.emplaceBack( i, val[0] ); + embeddedSurfaceTangentialLengthVectors.emplaceBack( i, val[1] ); + embeddedSurfaceTangentialLengthVectors.emplaceBack( i, val[2] ); + } + + //4. Get EDFM tangential width vectors (vertical) + ArrayOfArrays< real64 > embeddedSurfaceTangentialWidthVectors( numEdfmFracs ); + const char * const tw_str = "tangential_width_vectors"; + vtkDataArray * const tangential_width = grid->GetCellData()->GetArray( tw_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto const val = tangential_width->GetTuple3( i ); + embeddedSurfaceTangentialWidthVectors.emplaceBack( i, val[0] ); + embeddedSurfaceTangentialWidthVectors.emplaceBack( i, val[1] ); + embeddedSurfaceTangentialWidthVectors.emplaceBack( i, val[2] ); + } + + //4. Get EDFM aperture + array1d< real64 > embeddedSurfaceAperture( numEdfmFracs ); + const char * const aperture_str = "aperture"; + vtkDataArray * const apertures = grid->GetCellData()->GetArray( aperture_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto const val = apertures->GetTuple1( i ); + embeddedSurfaceAperture[i] = val; + } + + //5. Get EDFM permeability + array1d< real64 > embeddedSurfacePermeability( numEdfmFracs ); + const char * const perm_str = "permeability"; + vtkDataArray * const perms = grid->GetCellData()->GetArray( perm_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto val =perms->GetTuple1( i ); + embeddedSurfaceAperture[i] = val; + } + + //6. Get edfm to matrix cell mapping + const char * const fid_to_mid_mapping = "fracture_to_parent_matrix_cell_mapping"; + vtkDataArray * const fid_mid = grid->GetCellData()->GetArray( fid_to_mid_mapping ); + localIndex const blockIndex = 0; // TODO what happens when there are multiple blocks? + CellBlock const & cellBlock = cellBlockManager.getCellBlock( blockIndex ); + auto const & globalToLocalCell = cellBlock.globalToLocalMap(); + auto const & globalToLocalEdfm = embeddedSurfBlock.globalToLocalMap(); + + ArrayOfArrays< localIndex > toBlockIndex( numEdfmFracs ); + ArrayOfArrays< localIndex > toCellIndex( numEdfmFracs ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + // these indexes are global and needs to be converted to local + auto const fidmid = fid_mid->GetTuple2( i ); + // convert global to local using the maps + localIndex const id0 = globalToLocalEdfm.at( fidmid[0] ); + localIndex const id1 = globalToLocalCell.at( fidmid[1] ); + toBlockIndex.emplaceBack( id0, blockIndex );// cell block is set to 0 for now + toCellIndex.emplaceBack( id0, id1 ); + } + ToCellRelation< ArrayOfArrays< localIndex > > elem2dTo3d( std::move( toBlockIndex ), std::move( toCellIndex ) ); + + embeddedSurfBlock.setNumEmbeddedSurfElem( numEdfmFracs ); + embeddedSurfBlock.setEmbeddedSurfElemNodes( std::move( embeddedSurfaceElemNodes )); + embeddedSurfBlock.setEmbeddedSurfElemNormalVectors( std::move( embeddedSurfaceNormalVectors )); + embeddedSurfBlock.setEmbeddedSurfElemTangentialLengthVectors( std::move( embeddedSurfaceTangentialLengthVectors )); + embeddedSurfBlock.setEmbeddedSurfElemTangentialWidthVectors( std::move( embeddedSurfaceTangentialWidthVectors )); + embeddedSurfBlock.setEmbeddedSurfElemAperture( std::move( embeddedSurfaceAperture )); + embeddedSurfBlock.setEmbeddedSurfElemPermeability( std::move( embeddedSurfacePermeability )); + embeddedSurfBlock.setEmbeddedSurfElemTo3dElem( std::move( elem2dTo3d )); + embeddedSurfBlock.setEmbeddedSurfElemToNodes( std::move( elem2dToNodes )); +} +} diff --git a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp new file mode 100644 index 00000000000..b7f20b0e33a --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp @@ -0,0 +1,39 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2020- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_VTKEMBEDDEDSURFACEBLOCKUTILITIES_HPP +#define GEOS_VTKEMBEDDEDSURFACEBLOCKUTILITIES_HPP + +#include "CellBlockManager.hpp" + +#include "common/DataTypes.hpp" + +#include +#include + +namespace geos::vtk +{ + +/** + * @brief Attach the embedded surface block information to the cell block manager. + * @param embeddedSurfaceBlockName[in] The name of the embedded surface block. + * @param embeddedSurfaceMesh[in] The vtk mesh for the embedded surface block. + * @param cellBlockManager[inout] The cell block manager that will receive the embedded surface block information. + */ +void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName, + vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, + CellBlockManager & cellBlockManager ); +} + +#endif // include guard diff --git a/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp index 2618d8138f9..a6ed881f9f1 100644 --- a/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKFaceBlockUtilities.cpp @@ -605,7 +605,6 @@ array1d< globalIndex > buildLocalToGlobal( vtkIdTypeArray const * faceMeshCellGl return l2g; } - void importFractureNetwork( string const & faceBlockName, vtkSmartPointer< vtkDataSet > faceMesh, vtkSmartPointer< vtkDataSet > mesh, diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp index 6d791068dfd..c6ce7daaff2 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp @@ -21,6 +21,7 @@ #include "mesh/ExternalDataSourceManager.hpp" #include "mesh/generators/VTKFaceBlockUtilities.hpp" +#include "mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp" #include "mesh/generators/VTKMeshGeneratorTools.hpp" #include "mesh/generators/CellBlockManager.hpp" #include "mesh/generators/Region.hpp" @@ -65,6 +66,11 @@ VTKMeshGenerator::VTKMeshGenerator( string const & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "For multi-block files, names of the face mesh block." ); + registerWrapper( viewKeyStruct::embeddedSurfaceBlockNamesString(), &m_embeddedSurfaceBlockNames ). + setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "For multi-block files, names of the EDFM mesh block." ); + registerWrapper( viewKeyStruct::partitionRefinementString(), &m_partitionRefinement ). setInputFlag( InputFlags::OPTIONAL ). setApplyDefaultValue( 1 ). @@ -130,7 +136,7 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager if( !m_filePath.empty()) { GEOS_LOG_RANK_0( GEOS_FMT( "{} '{}': reading mesh from {}", catalogName(), getName(), m_filePath ) ); - allMeshes = vtk::loadAllMeshes( m_filePath, m_mainBlockName, m_faceBlockNames ); + allMeshes = vtk::loadAllMeshes( m_filePath, m_mainBlockName, m_faceBlockNames, m_embeddedSurfaceBlockNames ); } else if( !m_dataSourceName.empty()) { @@ -185,9 +191,11 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{} '{}': redistributing mesh...", catalogName(), getName() ) ); vtk::AllMeshes redistributedMeshes = - vtk::redistributeMeshes( getLogLevel(), allMeshes.getMainMesh(), allMeshes.getFaceBlocks(), comm, m_partitionMethod, m_partitionRefinement, m_useGlobalIds ); + vtk::redistributeMeshes( getLogLevel(), allMeshes.getMainMesh(), allMeshes.getFaceBlocks(), allMeshes.getEmbeddedSurfaceBlocks(), comm, m_partitionMethod, m_partitionRefinement, + m_useGlobalIds ); m_vtkMesh = redistributedMeshes.getMainMesh(); m_faceBlockMeshes = redistributedMeshes.getFaceBlocks(); + m_embeddedSurfaceBlockMeshes = redistributedMeshes.getEmbeddedSurfaceBlocks(); GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{} '{}': finding neighbor ranks...", catalogName(), getName() ) ); std::vector< vtkBoundingBox > boxes = vtk::exchangeBoundingBoxes( *m_vtkMesh, comm ); std::vector< int > const neighbors = vtk::findNeighborRanks( std::move( boxes ) ); @@ -214,8 +222,14 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager for( auto const & [name, mesh]: m_faceBlockMeshes ) { + GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{} '{}': importing fracture network {}...", catalogName(), getName(), name ) ); vtk::importFractureNetwork( name, mesh, m_vtkMesh, cellBlockManager ); } + for( auto const & [name, mesh]: m_embeddedSurfaceBlockMeshes ) + { + GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{} '{}': importing embedded fracture network {}...", catalogName(), getName(), name ) ); + vtk::importEmbeddedFractureNetwork( name, mesh, cellBlockManager ); + } GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{} '{}': done!", catalogName(), getName() ) ); vtk::printMeshStatistics( *m_vtkMesh, m_cellMap, comm ); diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp index 20747c848c2..ed12d06d6e6 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp @@ -109,6 +109,7 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase constexpr static char const * regionAttributeString() { return "regionAttribute"; } constexpr static char const * mainBlockNameString() { return "mainBlockName"; } constexpr static char const * faceBlockNamesString() { return "faceBlocks"; } + constexpr static char const * embeddedSurfaceBlockNamesString() { return "embeddedSurfaceBlocks"; } constexpr static char const * nodesetNamesString() { return "nodesetNames"; } constexpr static char const * partitionRefinementString() { return "partitionRefinement"; } constexpr static char const * partitionMethodString() { return "partitionMethod"; } @@ -148,9 +149,16 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase /// Name of the face blocks to be imported (for multi-block files). array1d< string > m_faceBlockNames; + /// Name of the edfm blocks to be imported (for multi-block files). + array1d< string > m_embeddedSurfaceBlockNames; + /// Maps the face block name to its vtk mesh instance. std::map< string, vtkSmartPointer< vtkDataSet > > m_faceBlockMeshes; + // Maps the edfm surface block name to its vtk mesh instance. + + std::map< string, vtkSmartPointer< vtkDataSet > > m_embeddedSurfaceBlockMeshes; + /// Names of VTK nodesets to import string_array m_nodesetNames; diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index a18b3e6767a..f719303a143 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -565,18 +565,24 @@ loadMesh( Path const & filePath, AllMeshes loadAllMeshes( Path const & filePath, string const & mainBlockName, - array1d< string > const & faceBlockNames ) + array1d< string > const & faceBlockNames, + array1d< string > const & edfmSurfBlockNames ) { int const lastRank = MpiWrapper::commSize() - 1; vtkSmartPointer< vtkDataSet > main = loadMesh( filePath, mainBlockName ); std::map< string, vtkSmartPointer< vtkDataSet > > faces; + std::map< string, vtkSmartPointer< vtkDataSet > > edfmSurfaces; for( string const & faceBlockName: faceBlockNames ) { faces[faceBlockName] = loadMesh( filePath, faceBlockName, lastRank ); } - return AllMeshes( main, faces ); + for( string const & edfmSurfBlockName: edfmSurfBlockNames ) + { + edfmSurfaces[edfmSurfBlockName] = loadMesh( filePath, edfmSurfBlockName, lastRank ); + } + return AllMeshes( main, faces, edfmSurfaces ); } @@ -688,7 +694,9 @@ AllMeshes redistributeByCellGraph( AllMeshes & input, finalFractures[fractureName] = finalFracMesh; } - return AllMeshes( finalMesh, finalFractures ); + // Ouassim: just add the edfm mesh at the moment and see. + auto edfmMesh = input.getEmbeddedSurfaceBlocks(); + return AllMeshes( finalMesh, finalFractures, edfmMesh ); } /** @@ -746,6 +754,7 @@ vtkSmartPointer< vtkDataSet > manageGlobalIds( vtkSmartPointer< vtkDataSet > mes { // Add global ids on the fly if needed int const me = hasGlobalIds( mesh ); + GEOS_LOG_RANK( "me hasGlobalIds="< manageGlobalIds( vtkSmartPointer< vtkDataSet > mes { mesh->GetPointData()->SetGlobalIds( vtkIdTypeArray::New() ); mesh->GetCellData()->SetGlobalIds( vtkIdTypeArray::New() ); + GEOS_LOG_RANK( "SetGlobalIds" ); } } @@ -900,6 +910,7 @@ AllMeshes redistributeMeshes( integer const logLevel, vtkSmartPointer< vtkDataSet > loadedMesh, std::map< string, vtkSmartPointer< vtkDataSet > > & namesToFractures, + std::map< string, vtkSmartPointer< vtkDataSet > > & namesToEdfmFractures, MPI_Comm const comm, PartitionMethod const method, int const partitionRefinement, @@ -913,8 +924,15 @@ redistributeMeshes( integer const logLevel, fractures.push_back( nameToFracture.second ); } + std::vector< vtkSmartPointer< vtkDataSet > > edfms; + for( auto & nameToEdfm: namesToEdfmFractures ) + { + edfms.push_back( nameToEdfm.second ); + } + // Generate global IDs for vertices and cells, if needed vtkSmartPointer< vtkDataSet > mesh = manageGlobalIds( loadedMesh, useGlobalIds, !std::empty( fractures ) ); + //vtkSmartPointer< vtkDataSet > mesh = manageGlobalIds( loadedMesh, useGlobalIds, !std::empty( fractures ) || !std::empty( edfms ) ); if( MpiWrapper::commRank( comm ) != ( MpiWrapper::commSize( comm ) - 1 ) ) { @@ -922,6 +940,11 @@ redistributeMeshes( integer const logLevel, { GEOS_ASSERT_EQ( nameToFracture.second->GetNumberOfCells(), 0 ); } + + for( auto nameToEdfm: namesToEdfmFractures ) + { + GEOS_ASSERT_EQ( nameToEdfm.second->GetNumberOfCells(), 0 ); + } } // Determine if redistribution is required @@ -944,13 +967,14 @@ redistributeMeshes( integer const logLevel, // Redistribute the mesh again using higher-quality graph partitioner if( partitionRefinement > 0 ) { - AllMeshes input( mesh, namesToFractures ); + AllMeshes input( mesh, namesToFractures, namesToEdfmFractures ); result = redistributeByCellGraph( input, method, comm, partitionRefinement - 1 ); } else { result.setMainMesh( mesh ); result.setFaceBlocks( namesToFractures ); + result.setEmbeddedSurfaceBlocks( namesToEdfmFractures ); } // Logging some information about the redistribution. @@ -962,6 +986,10 @@ redistributeMeshes( integer const logLevel, { messages.push_back( GEOS_FMT( pattern, faceName, faceMesh->GetNumberOfCells() ) ); } + for( auto const & [edfmName, edfmMesh]: result.getEmbeddedSurfaceBlocks() ) + { + messages.push_back( GEOS_FMT( pattern, edfmName, edfmMesh->GetNumberOfCells() ) ); + } if( logLevel >= 5 ) { GEOS_LOG_RANK( stringutilities::join( messages, ", " ) ); @@ -1770,6 +1798,7 @@ void fillCellBlock( vtkDataSet & mesh, localIndex const numNodesPerElement = cellBlock.numNodesPerElement(); arrayView2d< localIndex, cells::NODE_MAP_USD > const cellToVertex = cellBlock.getElemToNode(); arrayView1d< globalIndex > const & localToGlobal = cellBlock.localToGlobalMap(); + auto & globalToLocal = cellBlock.globalToLocalMap(); vtkIdTypeArray const * const globalCellId = vtkIdTypeArray::FastDownCast( mesh.GetCellData()->GetGlobalIds() ); GEOS_ERROR_IF( !cellIds.empty() && globalCellId == nullptr, "Global cell IDs have not been generated" ); @@ -1780,7 +1809,10 @@ void fillCellBlock( vtkDataSet & mesh, { cellToVertex[cellCount][v] = cell->GetPointId( nodeOrder[v] ); } - localToGlobal[cellCount++] = globalCellId->GetValue( c ); + globalIndex g = globalCellId->GetValue( c ); + localToGlobal[cellCount] = g; + globalToLocal[g] = cellCount; + cellCount++; }; // Writing connectivity and Local to Global diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp index 0a2f5f3903a..eefbf8eebaf 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.hpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.hpp @@ -76,11 +76,14 @@ class AllMeshes * @brief Builds the compound from values. * @param main The main 3d mesh (the matrix). * @param faceBlocks The fractures meshes. + * @param edfmSurfBlocks The edfm fractures meshes. */ AllMeshes( vtkSmartPointer< vtkDataSet > const & main, - std::map< string, vtkSmartPointer< vtkDataSet > > const & faceBlocks ) + std::map< string, vtkSmartPointer< vtkDataSet > > const & faceBlocks, + std::map< string, vtkSmartPointer< vtkDataSet > > const & edfmSurfBlocks ) : m_main( main ), - m_faceBlocks( faceBlocks ) + m_faceBlocks( faceBlocks ), + m_embeddedSurfaceBlocks( edfmSurfBlocks ) { } /** @@ -99,6 +102,14 @@ class AllMeshes return m_faceBlocks; } + /** + * @return a mapping linking the name of each edfm block to its mesh. + */ + std::map< string, vtkSmartPointer< vtkDataSet > > & getEmbeddedSurfaceBlocks() + { + return m_embeddedSurfaceBlocks; + } + /** * @brief Defines the main 3d mesh for the simulation. * @param main The new 3d mesh. @@ -117,12 +128,23 @@ class AllMeshes m_faceBlocks = faceBlocks; } + /** + * @brief Defines the face blocks/fractures. + * @param edfmSurfBlocks A map which connects each name of the face block to its mesh. + */ + void setEmbeddedSurfaceBlocks( std::map< string, vtkSmartPointer< vtkDataSet > > const & edfmSurfBlocks ) + { + m_embeddedSurfaceBlocks = edfmSurfBlocks; + } private: /// The main 3d mesh (namely the matrix). vtkSmartPointer< vtkDataSet > m_main; /// The face meshes (namely the fractures). std::map< string, vtkSmartPointer< vtkDataSet > > m_faceBlocks; + + /// The EDFM meshes (namely the EDFM fractures). + std::map< string, vtkSmartPointer< vtkDataSet > > m_embeddedSurfaceBlocks; }; /** @@ -130,11 +152,13 @@ class AllMeshes * @param[in] filePath the Path of the file to load * @param[in] mainBlockName The name of the block to import (will be considered for multi-block files only). * @param[in] faceBlockNames The names of the face blocks to import (will be considered for multi-block files only). + * @param[in] edfmSurfBlockNames The names of the EDFM blocks to import (will be considered for multi-block files only). * @return The compound of the main mesh and the face block meshes. */ AllMeshes loadAllMeshes( Path const & filePath, string const & mainBlockName, - array1d< string > const & faceBlockNames ); + array1d< string > const & faceBlockNames, + array1d< string > const & edfmSurfBlockNames ); /** * @brief Compute the rank neighbor candidate list. @@ -149,6 +173,7 @@ findNeighborRanks( std::vector< vtkBoundingBox > boundingBoxes ); * @param[in] logLevel the log level * @param[in] loadedMesh the mesh that was loaded on one or several MPI ranks * @param[in] namesToFractures the fracture meshes + * @param[in] namesToEdfmFractures the EDFM meshes * @param[in] comm the MPI communicator * @param[in] method the partitionning method * @param[in] partitionRefinement number of graph partitioning refinement cycles @@ -159,6 +184,7 @@ AllMeshes redistributeMeshes( integer const logLevel, vtkSmartPointer< vtkDataSet > loadedMesh, std::map< string, vtkSmartPointer< vtkDataSet > > & namesToFractures, + std::map< string, vtkSmartPointer< vtkDataSet > > & namesToEdfmFractures, MPI_Comm const comm, PartitionMethod const method, int const partitionRefinement, diff --git a/src/coreComponents/mesh/generators/VTKWellGenerator.cpp b/src/coreComponents/mesh/generators/VTKWellGenerator.cpp index b52409d3f68..fc570fdb7ac 100644 --- a/src/coreComponents/mesh/generators/VTKWellGenerator.cpp +++ b/src/coreComponents/mesh/generators/VTKWellGenerator.cpp @@ -50,7 +50,7 @@ void VTKWellGenerator::fillPolylineDataStructure( ) GEOS_LOG_RANK_0( GEOS_FMT( "{} '{}': reading well from {}", catalogName(), getName(), m_filePath ) ); { GEOS_LOG_LEVEL_RANK_0( 2, " reading the dataset..." ); - vtk::AllMeshes allMeshes = vtk::loadAllMeshes( m_filePath, "main", array1d< string >()); + vtk::AllMeshes allMeshes = vtk::loadAllMeshes( m_filePath, "main", array1d< string >(), array1d< string >()); vtkSmartPointer< vtkDataSet > loadedMesh = allMeshes.getMainMesh(); controller->Broadcast( loadedMesh, 0 ); diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index da4e2856368..f34c1a2218d 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -75,7 +75,7 @@ EmbeddedSurfaceGenerator::EmbeddedSurfaceGenerator( const string & name, registerWrapper( viewKeyStruct::targetObjectsNameString(), &m_targetObjectsName ). setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ). - setInputFlag( dataRepository::InputFlags::REQUIRED ). + setInputFlag( dataRepository::InputFlags::OPTIONAL ). setDescription( "List of geometric objects that will be used to initialized the embedded surfaces/fractures." ); registerWrapper( viewKeyStruct::mpiCommOrderString(), &m_mpiCommOrder ). @@ -130,81 +130,121 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() NewObjectLists newObjects; - // Loop over all the fracture planes - geometricObjManager.forGeometricObject< PlanarGeometricObject >( m_targetObjectsName, [&]( localIndex const, - PlanarGeometricObject & fracture ) + if( m_targetObjectsName.empty() ) // when targetObjects were not provided - try to import from vtk { - /* 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 ) + MeshBody & meshBody = domain.getMeshBody( 0 ); + CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager(); + Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks(); + string const & faceBlockName = embeddedSurfaceRegion.getFaceBlockName(); + if( embSurfBlocks.hasGroup( faceBlockName )) { - arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList(); - FixedOneToManyRelation const & cellToEdges = subRegion.edgeList(); + EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( faceBlockName ); + + elemManager.forElementSubRegionsComplete< CellElementSubRegion >( + [&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion ) + { + FixedOneToManyRelation const & cellToEdges = subRegion.edgeList(); - arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + bool added = embeddedSurfaceSubRegion.addAllEmbeddedSurfaces( er, + esr, + nodeManager, + embSurfNodeManager, + edgeManager, + cellToEdges, + embSurf ); - forAll< serialPolicy >( subRegion.size(), [ &, ghostRank, - cellToNodes, - nodesCoord ] ( localIndex const cellIndex ) + 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 ) { - if( ghostRank[cellIndex] < 0 ) + 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 ) { - isPositive = 0; - isNegative = 0; - for( localIndex kn = 0; kn < subRegion.numNodesPerElement(); kn++ ) + if( ghostRank[cellIndex] < 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 ) + isPositive = 0; + isNegative = 0; + for( localIndex kn = 0; kn < subRegion.numNodesPerElement(); kn++ ) { - isNegative = 1; - } - } // end loop over nodes - if( isPositive * isNegative == 1 ) - { - bool added = embeddedSurfaceSubRegion.addNewEmbeddedSurface( cellIndex, - er, - esr, - nodeManager, - embSurfNodeManager, - edgeManager, - cellToEdges, - &fracture ); - - if( added ) + 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 ) { - 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++; + 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 >( @@ -267,9 +307,6 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() embSurfNodeManager.compressRelationMaps(); } -void EmbeddedSurfaceGenerator::initializePostInitialConditionsPreSubGroups() -{} - real64 EmbeddedSurfaceGenerator::solverStep( real64 const & GEOS_UNUSED_PARAM( time_n ), real64 const & GEOS_UNUSED_PARAM( dt ), const int GEOS_UNUSED_PARAM( cycleNumber ), diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.hpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.hpp index 5f9219a4ea0..5bfbe9dfa5a 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.hpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.hpp @@ -97,8 +97,6 @@ class EmbeddedSurfaceGenerator : public PhysicsSolverBase virtual void initializePostSubGroups() override final; - virtual void initializePostInitialConditionsPreSubGroups() override final; - virtual void postRestartInitialization() override final { GEOS_ERROR( "Restarting is not supported for cases involving EmbeddedSurfaceGenerator" ); diff --git a/src/coreComponents/schema/docs/VTKMesh.rst b/src/coreComponents/schema/docs/VTKMesh.rst new file mode 100644 index 00000000000..458ef31976e --- /dev/null +++ b/src/coreComponents/schema/docs/VTKMesh.rst @@ -0,0 +1,27 @@ + + +====================== ======================== ========= ============================================================================================================================================================================================================================================================================================================================================================================================================================================================================ +Name Type Default Description +====================== ======================== ========= ============================================================================================================================================================================================================================================================================================================================================================================================================================================================================ +embeddedSurfaceBlocks groupNameRef_array {} For multi-block files, names of the EDFM mesh block. +faceBlocks groupNameRef_array {} For multi-block files, names of the face mesh block. +fieldNamesInGEOSX groupNameRef_array {} Names of the volumic fields in GEOSX to import into +fieldsToImport groupNameRef_array {} Volumic fields to be imported from the external mesh file +file path required Path to the mesh file +logLevel integer 0 Log level +mainBlockName groupNameRef main For multi-block files, name of the 3d mesh block. +name groupName required A name is required for any non-unique nodes +nodesetNames groupNameRef_array {} Names of the VTK nodesets to import +partitionMethod geos_vtk_PartitionMethod parmetis Method (library) used to partition the mesh +partitionRefinement integer 1 Number of partitioning refinement iterations (defaults to 1, recommended value).A value of 0 disables graph partitioning and keeps simple kd-tree partitions (not recommended). Values higher than 1 may lead to slightly improved partitioning, but yield diminishing returns. +regionAttribute groupNameRef attribute Name of the VTK cell attribute to use as region marker +scale R1Tensor {1,1,1} Scale the coordinates of the vertices by given scale factors (after translation) +surfacicFieldsInGEOSX groupNameRef_array {} Names of the surfacic fields in GEOSX to import into +surfacicFieldsToImport groupNameRef_array {} Surfacic fields to be imported from the external mesh file +translate R1Tensor {0,0,0} Translate the coordinates of the vertices by a given vector (prior to scaling) +useGlobalIds integer 0 Controls the use of global IDs in the input file for cells and points. If set to 0 (default value), the GlobalId arrays in the input mesh are used if available, and generated otherwise. If set to a negative value, the GlobalId arrays in the input mesh are not used, and generated global Ids are automatically generated. If set to a positive value, the GlobalId arrays in the input mesh are used and required, and the simulation aborts if they are not available +InternalWell node :ref:`XML_InternalWell` +VTKWell node :ref:`XML_VTKWell` +====================== ======================== ========= ============================================================================================================================================================================================================================================================================================================================================================================================================================================================================ + + diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index e130b1e56d4..0c8988ee9b8 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -1765,6 +1765,8 @@ stress - traction is applied to the faces as specified by the inner product of i + +