From d77d0825126423835a7e82d4e1101c02bee891e8 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Wed, 6 Sep 2023 16:16:30 -0500 Subject: [PATCH 01/32] initial embedded surface block abc --- .../generators/EmbeddedSurfaceBlockABC.hpp | 91 +++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp new file mode 100644 index 00000000000..558fb96dea4 --- /dev/null +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -0,0 +1,91 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" + + +/** + * @brief An embedded 2d element within a 3d element + * + * @details The @p EmbeddedSurfaceBlockABC represents a 2d element (@e surfacic) embedded within + * a 3d element. The 2d element is assumed to be a quad. + * @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 or N 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 getEmbeddedSurfElemToNodes() const = 0; + + /** + * @brief Get the indices of the parent 3d elements of eacch embedded surface element (1 or more) + * @return The mapping of first dimension @p numEmbeddedSurfElem. + * Second dimension is heterogeneous, and depends on how many 3d element does the embedded surface element + * intersect. + * + * @details each embedded surface element is intersects 1 or more 3d elements. Indices of these 3d elements + * are returned for each embedded surface element + */ + virtual ArrayOfArrays 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 @p numEmbeddedSurfaceElements*4. + * Second dimension is 3, and represents the x, y and z. + */ + virtual ArrayOfArrays getEmbeddedSurfaceElementsNodes() const = 0; + + + +} \ No newline at end of file From 35196c8bb3faf2b7e5189acb2054e5e2e076d807 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Thu, 7 Sep 2023 09:11:41 -0500 Subject: [PATCH 02/32] minor change --- .../mesh/generators/EmbeddedSurfaceBlockABC.hpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp index 558fb96dea4..1b1fb2ffab3 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -67,12 +67,12 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group virtual ArrayOfArrays getEmbeddedSurfElemToNodes() const = 0; /** - * @brief Get the indices of the parent 3d elements of eacch embedded surface element (1 or more) + * @brief Get the indices of the parent 3d elements of each embedded surface element (1 or more) * @return The mapping of first dimension @p numEmbeddedSurfElem. - * Second dimension is heterogeneous, and depends on how many 3d element does the embedded surface element + * Second dimension 1 to number of 3d elements, and depends on how many 3d elements does the embedded surface element * intersect. * - * @details each embedded surface element is intersects 1 or more 3d elements. Indices of these 3d elements + * @details each embedded surface element intersects 1 or more 3d elements. Indices of these 3d elements * are returned for each embedded surface element */ virtual ArrayOfArrays getEmbeddedSurfElemTo3dElem() const = 0; @@ -84,8 +84,6 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group * first dimension is @p numEmbeddedSurfaceElements*4. * Second dimension is 3, and represents the x, y and z. */ - virtual ArrayOfArrays getEmbeddedSurfaceElementsNodes() const = 0; - - + virtual ArrayOfArrays getEmbeddedSurfElemNodes() const = 0; } \ No newline at end of file From 2135a211874903938f83988e909cfbd43dabdd44 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Wed, 13 Sep 2023 13:02:00 -0500 Subject: [PATCH 03/32] start concrete implementation --- .../mesh/generators/EmbeddedSurfaceBlock.hpp | 85 +++++++++++++++++++ .../generators/EmbeddedSurfaceBlockABC.hpp | 12 ++- 2 files changed, 94 insertions(+), 3 deletions(-) create mode 100644 src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp new file mode 100644 index 00000000000..7fd9bf3d8e9 --- /dev/null +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -0,0 +1,85 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 ) + : + EmbeddedSurfaceBlock( name, parent ) + { } + + localIndex numEmbeddedSurfElem() const override; + ArrayOfArrays getEmbeddedSurfElemToNodes() const override; + ArrayOfArrays getEmbeddedSurfElemTo3dElem() const override; + ArrayOfArrays getEmbeddedSurfElemNodes() 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 && _embeddedSurfElemToNodes); + + /** + * @brief Sets the embedded elements to 3d elements mapping + * @param _embeddedSurfElemTo3dElem the input mapping. + */ + void setEmbeddedSurfElemTo3dElem(ArrayOfArrays && _embeddedSurfElemTo3dElem); + + /** + * @brief Sets the embedded elements nodes coordinates + * @param _embeddedSurfElemNodes the coordinates of embedded elements nodes. + */ + void setEmbeddedSurfElemNodes(ArrayOfArrays && _embeddedSurfElemNodes); + +private: + + localIndex m_numEmbeddedSurfElem; + + ArrayOfArrays< localIndex > m_embeddedSurfElemToNodes; + + ArrayOfArrays< localIndex > m_embeddedSurfElemTo3dElem; + + ArrayOfArrays< real64 > m_embeddedSurfElemNodes; +}; + +} + +#endif \ No newline at end of file diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp index 1b1fb2ffab3..c5a1536cc3d 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -21,11 +21,14 @@ #include "common/DataTypes.hpp" +namespace geos +{ + /** * @brief An embedded 2d element within a 3d element * - * @details The @p EmbeddedSurfaceBlockABC represents a 2d element (@e surfacic) embedded within - * a 3d element. The 2d element is assumed to be a quad. + * @details The @p EmbeddedSurfaceBlockABC represents an array of 2d elements (@e surfacic) embedded within + * a 3d elements grid (2d element can intersect one or more 3d element). The 2d element is assumed to be a quad. * @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 @@ -85,5 +88,8 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group * Second dimension is 3, and represents the x, y and z. */ virtual ArrayOfArrays getEmbeddedSurfElemNodes() const = 0; +}; -} \ No newline at end of file +} + +#endif // include guard \ No newline at end of file From 7eb79a56f4c4c71971ef2afa085c4fa719b76773 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Wed, 13 Sep 2023 16:36:42 -0500 Subject: [PATCH 04/32] add the cpp --- .../mesh/generators/EmbeddedSurfaceBlock.cpp | 59 +++++++++++++++++++ .../mesh/generators/EmbeddedSurfaceBlock.hpp | 5 +- 2 files changed, 60 insertions(+), 4 deletions(-) create mode 100644 src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp new file mode 100644 index 00000000000..9addec5ef2b --- /dev/null +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp @@ -0,0 +1,59 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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_numEmbeddedSurfElem; +} + +ArrayOfArrays EmbeddedSurfaceBLock::getEmbeddedSurfElemToNodes() const { + return m_embeddedSurfElemToNodes; +} + +ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const { + return m_embeddedSurfElemTo3dElem; +}; + +ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNodes() const { + return m_embeddedSurfElemNodes; +} + + +void EmbeddedSurfaceBlock::setNumEmbeddedSurfElem(localIndex _numEmbeddedSurfElem){ + + m_numEmbeddedSurfElem = _numEmbeddedSurfElem; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemToNodes(ArrayOfArrays && _embeddedSurfElemToNodes){ + m_embeddedSurfElemToNodes = _embeddedSurfElemToNodes; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem(ArrayOfArrays && _embeddedSurfElemTo3dElem){ + m_embeddedSurfElemTo3dElem = _embeddedSurfElemTo3dElem; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemNodes(ArrayOfArrays && _embeddedSurfElemNodes){ + m_embeddedSurfElemNodes = _embeddedSurfElemNodes; +} + +} \ No newline at end of file diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp index 7fd9bf3d8e9..d49e9b8633b 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -72,14 +72,11 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC private: localIndex m_numEmbeddedSurfElem; - ArrayOfArrays< localIndex > m_embeddedSurfElemToNodes; - ArrayOfArrays< localIndex > m_embeddedSurfElemTo3dElem; - ArrayOfArrays< real64 > m_embeddedSurfElemNodes; }; } -#endif \ No newline at end of file +#endif //inlcude guard \ No newline at end of file From b20dd27b1028c94accae50dee0792264f1776b03 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Fri, 15 Sep 2023 10:49:51 -0500 Subject: [PATCH 05/32] add embedded surfaces to CellBlockManager --- .../mesh/generators/CellBlockManager.cpp | 1 + .../mesh/generators/CellBlockManager.hpp | 21 +++++++++++++++++++ .../mesh/generators/CellBlockManagerABC.hpp | 19 +++++++++++++++++ src/docs/doxygen/GeosxConfig.hpp | 2 +- 4 files changed, 42 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index 2db3b3cfff2..4fa118c241f 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, diff --git a/src/coreComponents/mesh/generators/CellBlockManager.hpp b/src/coreComponents/mesh/generators/CellBlockManager.hpp index dcb15629251..901379a96f0 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" @@ -81,6 +82,8 @@ class CellBlockManager : public CellBlockManagerABC ToCellRelation< array2d< localIndex > > getFaceToElements() const override; + ToCellRelation< ArrayOfArrays< localIndex > > getEmbeddedSurfaceToElements() const override; + array1d< globalIndex > getNodeLocalToGlobal() const override; /** @@ -125,6 +128,8 @@ class CellBlockManager : public CellBlockManagerABC localIndex numFaces() const override; + localIndex numEmbeddedSurfaces() const override; + using Group::resize; /** @@ -165,6 +170,10 @@ class CellBlockManager : public CellBlockManagerABC Group & getFaceBlocks() override; + Group & const getEmbeddedSurfaceBlocks() override; + + Group & getEmbeddedSurfaceBlocks() override; + LineBlockABC const & getLineBlock( string name ) const override; /** @@ -195,6 +204,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 @@ -229,6 +244,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"; } }; /** @@ -279,6 +297,8 @@ class CellBlockManager : public CellBlockManagerABC ArrayOfArrays< localIndex > m_faceToEdges; ToCellRelation< array2d< localIndex > > m_faceToCells; + ToCellRelation< ArrayOfArrays< localIndex > > m_embeddedSurfToCells; + array1d< globalIndex > m_nodeLocalToGlobal; std::map< string, SortedArray< localIndex > > m_nodeSets; @@ -290,6 +310,7 @@ class CellBlockManager : public CellBlockManagerABC localIndex m_numNodes; localIndex m_numFaces; localIndex m_numEdges; + localIndex m_numEmbeddedSurfElem; }; } diff --git a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp index 3c474520389..b99da778dfa 100644 --- a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp @@ -126,6 +126,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 const Group & getEmbeddedSurfaceBlocks() const = 0; + /** * @brief Total number of nodes across all the cell blocks. * @return The total number of nodes. @@ -146,6 +152,12 @@ class CellBlockManagerABC : public dataRepository::Group */ virtual localIndex numFaces() const = 0; + /** + * @brief Total number of embedded surfaces across all the cell blocks. + * @return The total number of embedded surfaces. + */ + virtual localIndex numEmbeddedSurfaces() const = 0; + /** * @brief Returns the node coordinates in a (numNodes, 3) 2d array. * @return A const view to the array. @@ -201,6 +213,13 @@ class CellBlockManagerABC : public dataRepository::Group */ virtual ToCellRelation< array2d< localIndex > > getFaceToElements() const = 0; + /** + * @brief Returns the embedded surface to elements mapping. + * @return A 1 to many relationship. The result is jagged array mapping each embedded element + * To the elements it intersect. + * + */ + virtual ToCellRelation> getEmbeddedSurfaceToElements() const = 0; /** * @brief The node to global mapping for nodes. * @return The mapping as an array of size numNodes. diff --git a/src/docs/doxygen/GeosxConfig.hpp b/src/docs/doxygen/GeosxConfig.hpp index 44a7b9be831..4fcd8ccc5bd 100644 --- a/src/docs/doxygen/GeosxConfig.hpp +++ b/src/docs/doxygen/GeosxConfig.hpp @@ -171,7 +171,7 @@ #define fmt_VERSION 10.0.0 /// Version information for python -#define Python3_VERSION 3.10.6 +#define Python3_VERSION 3.10.12 /// Version information for CUDAToolkit /* #undef CUDAToolkit_VERSION */ From 738c14ac7572396aa8d48c2d30f712197a74cdbe Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 29 Sep 2023 11:24:41 -0500 Subject: [PATCH 06/32] this should not be commited --- src/docs/doxygen/GeosxConfig.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/docs/doxygen/GeosxConfig.hpp b/src/docs/doxygen/GeosxConfig.hpp index 4fcd8ccc5bd..f264cce1d54 100644 --- a/src/docs/doxygen/GeosxConfig.hpp +++ b/src/docs/doxygen/GeosxConfig.hpp @@ -171,7 +171,7 @@ #define fmt_VERSION 10.0.0 /// Version information for python -#define Python3_VERSION 3.10.12 +#define Python3_VERSION 3.10.8 /// Version information for CUDAToolkit /* #undef CUDAToolkit_VERSION */ From 2fcdeed861d337fc0e570d9de0e728980753f3ad Mon Sep 17 00:00:00 2001 From: Ouassim Date: Wed, 27 Sep 2023 11:11:42 -0500 Subject: [PATCH 07/32] upate cellBlock Manger and embedded Ruface block --- .../mesh/generators/CellBlockManager.cpp | 17 +++++++++++++++++ .../mesh/generators/CellBlockManager.hpp | 2 +- .../mesh/generators/CellBlockManagerABC.hpp | 2 +- .../mesh/generators/EmbeddedSurfaceBlock.cpp | 6 +++--- .../mesh/generators/EmbeddedSurfaceBlock.hpp | 2 +- 5 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index 4fa118c241f..92f0a3c1d80 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.cpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp @@ -669,6 +669,8 @@ void CellBlockManager::buildMaps() buildNodeToEdges(); fillElementToEdgesOfCellBlocks( m_faceToEdges.toViewConst(), this->getCellBlocks() ); + + buildEmbeddedSurfaceMaps(); } ArrayOfArrays< localIndex > CellBlockManager::getFaceToNodes() const @@ -688,6 +690,11 @@ ToCellRelation< array2d< localIndex > > CellBlockManager::getFaceToElements() co return m_faceToCells; } +ToCellRelation< ArrayOfArrays< localIndex > > getEmbeddedSurfaceToElements() const { + + return m_embeddedSurfToCells +} + const Group & CellBlockManager::getCellBlocks() const { return this->getGroup( viewKeyStruct::cellBlocks() ); @@ -738,6 +745,11 @@ localIndex CellBlockManager::numFaces() const return m_numFaces; } +localIndex CellBlockManager::numEmbeddedSurfaces() const +{ + return m_numEmbeddedSurfaces; +} + ArrayOfArrays< localIndex > CellBlockManager::getEdgeToFaces() const { return m_edgeToFaces; @@ -788,6 +800,11 @@ LineBlock & CellBlockManager::registerLineBlock( string const & name ) return this->getLineBlocks().registerGroup< LineBlock >( name ); } +EmbeddedSurfaceBlock & CellBlockManager::registerEmbeddedSurfaceBlock( string const & name ){ + + return this->getLineBlocks().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 901379a96f0..edcff4b217c 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.hpp @@ -310,7 +310,7 @@ class CellBlockManager : public CellBlockManagerABC localIndex m_numNodes; localIndex m_numFaces; localIndex m_numEdges; - localIndex m_numEmbeddedSurfElem; + localIndex m_numEmbeddedSurfaces; }; } diff --git a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp index b99da778dfa..45e9d80e595 100644 --- a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp @@ -216,7 +216,7 @@ class CellBlockManagerABC : public dataRepository::Group /** * @brief Returns the embedded surface to elements mapping. * @return A 1 to many relationship. The result is jagged array mapping each embedded element - * To the elements it intersect. + * To the 3d cell elements it intersects. * */ virtual ToCellRelation> getEmbeddedSurfaceToElements() const = 0; diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp index 9addec5ef2b..70613a3908b 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp @@ -23,7 +23,7 @@ namespace geos localIndex EmbeddedSurfaceBlock::numEmbeddedSurfElem() const{ - return m_numEmbeddedSurfElem; + return m_numEmbeddedSurfaces; } ArrayOfArrays EmbeddedSurfaceBLock::getEmbeddedSurfElemToNodes() const { @@ -39,9 +39,9 @@ ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNodes() const { } -void EmbeddedSurfaceBlock::setNumEmbeddedSurfElem(localIndex _numEmbeddedSurfElem){ +void EmbeddedSurfaceBlock::setNumEmbeddedSurfElem(localIndex _numEmbeddedSurfaces){ - m_numEmbeddedSurfElem = _numEmbeddedSurfElem; + m_numEmbeddedSurfaces = _numEmbeddedSurfaces; } void EmbeddedSurfaceBlock::setEmbeddedSurfElemToNodes(ArrayOfArrays && _embeddedSurfElemToNodes){ diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp index d49e9b8633b..be248c85c4f 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -71,7 +71,7 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC private: - localIndex m_numEmbeddedSurfElem; + localIndex m_numEmbeddedSurfaces; ArrayOfArrays< localIndex > m_embeddedSurfElemToNodes; ArrayOfArrays< localIndex > m_embeddedSurfElemTo3dElem; ArrayOfArrays< real64 > m_embeddedSurfElemNodes; From 186a222ba2e939554998601840112fd6c57d8b2a Mon Sep 17 00:00:00 2001 From: Ouassim Date: Thu, 26 Oct 2023 12:03:28 -0500 Subject: [PATCH 08/32] Correct edfm to cell type --- src/coreComponents/mesh/generators/CellBlockManagerABC.hpp | 7 ------- .../mesh/generators/EmbeddedSurfaceBlock.cpp | 4 ++-- .../mesh/generators/EmbeddedSurfaceBlock.hpp | 6 +++--- .../mesh/generators/EmbeddedSurfaceBlockABC.hpp | 2 +- 4 files changed, 6 insertions(+), 13 deletions(-) diff --git a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp index 45e9d80e595..0175293ccd5 100644 --- a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp @@ -213,13 +213,6 @@ class CellBlockManagerABC : public dataRepository::Group */ virtual ToCellRelation< array2d< localIndex > > getFaceToElements() const = 0; - /** - * @brief Returns the embedded surface to elements mapping. - * @return A 1 to many relationship. The result is jagged array mapping each embedded element - * To the 3d cell elements it intersects. - * - */ - virtual ToCellRelation> getEmbeddedSurfaceToElements() const = 0; /** * @brief The node to global mapping for nodes. * @return The mapping as an array of size numNodes. diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp index 70613a3908b..7931478b6d7 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp @@ -30,7 +30,7 @@ ArrayOfArrays EmbeddedSurfaceBLock::getEmbeddedSurfElemToNodes() con return m_embeddedSurfElemToNodes; } -ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const { +ToCellRelation> EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const { return m_embeddedSurfElemTo3dElem; }; @@ -48,7 +48,7 @@ void EmbeddedSurfaceBlock::setEmbeddedSurfElemToNodes(ArrayOfArrays m_embeddedSurfElemToNodes = _embeddedSurfElemToNodes; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem(ArrayOfArrays && _embeddedSurfElemTo3dElem){ +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem(ToCellRelation> && _embeddedSurfElemTo3dElem){ m_embeddedSurfElemTo3dElem = _embeddedSurfElemTo3dElem; } diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp index be248c85c4f..31d47315f64 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -41,7 +41,7 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC localIndex numEmbeddedSurfElem() const override; ArrayOfArrays getEmbeddedSurfElemToNodes() const override; - ArrayOfArrays getEmbeddedSurfElemTo3dElem() const override; + ToCellRelation> getEmbeddedSurfElemTo3dElem() const override; ArrayOfArrays getEmbeddedSurfElemNodes() const override; @@ -61,7 +61,7 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC * @brief Sets the embedded elements to 3d elements mapping * @param _embeddedSurfElemTo3dElem the input mapping. */ - void setEmbeddedSurfElemTo3dElem(ArrayOfArrays && _embeddedSurfElemTo3dElem); + void setEmbeddedSurfElemTo3dElem(ToCellRelation> && _embeddedSurfElemTo3dElem); /** * @brief Sets the embedded elements nodes coordinates @@ -73,7 +73,7 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC localIndex m_numEmbeddedSurfaces; ArrayOfArrays< localIndex > m_embeddedSurfElemToNodes; - ArrayOfArrays< localIndex > m_embeddedSurfElemTo3dElem; + ToCellRelation> m_embeddedSurfElemTo3dElem; ArrayOfArrays< real64 > m_embeddedSurfElemNodes; }; diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp index c5a1536cc3d..a1fb49e79a0 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -78,7 +78,7 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group * @details each embedded surface element intersects 1 or more 3d elements. Indices of these 3d elements * are returned for each embedded surface element */ - virtual ArrayOfArrays getEmbeddedSurfElemTo3dElem() const = 0; + virtual ToCellRelation> getEmbeddedSurfElemTo3dElem() const = 0; /** From 2401cc7ba2c4cb117f2dec0dfd6d3c2d3325a32e Mon Sep 17 00:00:00 2001 From: Ouassim Date: Tue, 16 Apr 2024 07:38:22 -0500 Subject: [PATCH 09/32] corrections --- src/coreComponents/mesh/generators/CellBlockManager.cpp | 2 +- src/coreComponents/mesh/generators/CellBlockManager.hpp | 2 +- src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp | 2 +- src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index 92f0a3c1d80..8c0f468b484 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.cpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp @@ -692,7 +692,7 @@ ToCellRelation< array2d< localIndex > > CellBlockManager::getFaceToElements() co ToCellRelation< ArrayOfArrays< localIndex > > getEmbeddedSurfaceToElements() const { - return m_embeddedSurfToCells + return m_embeddedSurfToCells; } const Group & CellBlockManager::getCellBlocks() const diff --git a/src/coreComponents/mesh/generators/CellBlockManager.hpp b/src/coreComponents/mesh/generators/CellBlockManager.hpp index edcff4b217c..5f3cb3dbae4 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.hpp @@ -170,7 +170,7 @@ class CellBlockManager : public CellBlockManagerABC Group & getFaceBlocks() override; - Group & const getEmbeddedSurfaceBlocks() override; + Group & const getEmbeddedSurfaceBlocks() const override; Group & getEmbeddedSurfaceBlocks() override; diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp index 31d47315f64..e51c4a36a97 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -36,7 +36,7 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC EmbeddedSurfaceBlock( string const & name, Group * const parent ) : - EmbeddedSurfaceBlock( name, parent ) + EmbeddedSurfaceBlockABC( name, parent ) { } localIndex numEmbeddedSurfElem() const override; diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp index a1fb49e79a0..72d0f9fb3a0 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -56,7 +56,7 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group * @details Return the number of embedded surface elements, each surface element * can intersect 1 or N 3d elements */ - virtual localIndex numEmbeddedSurfElem() const = 0 + virtual localIndex numEmbeddedSurfElem() const = 0; /** * @brief Get the indices of the nodes of all embedded surface elements From e2346e7a69399eafb72069f8e8f45f8868de7886 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Tue, 9 Apr 2024 15:51:41 -0500 Subject: [PATCH 10/32] create edfm patches in bulk --- .../fractureMatrixFlow_edfm_base.xml | 4 +- ...reMatrixFlow_edfm_horizontalFrac_smoke.xml | 8 ++-- src/coreComponents/mesh/CMakeLists.txt | 3 ++ .../mesh/ElementRegionManager.cpp | 12 ++++- .../mesh/EmbeddedSurfaceSubRegion.cpp | 45 +++++++++++++++++++ .../mesh/EmbeddedSurfaceSubRegion.hpp | 7 +++ .../mesh/SurfaceElementRegion.cpp | 11 ++++- .../mesh/SurfaceElementRegion.hpp | 1 + .../mesh/generators/CellBlockManager.cpp | 17 +++++-- .../mesh/generators/CellBlockManager.hpp | 8 ++-- .../mesh/generators/CellBlockManagerABC.hpp | 14 ++++++ .../mesh/generators/EmbeddedSurfaceBlock.cpp | 30 ++++++++++--- .../mesh/generators/EmbeddedSurfaceBlock.hpp | 32 +++++++++++-- .../generators/EmbeddedSurfaceBlockABC.hpp | 38 +++++++++++++--- 14 files changed, 196 insertions(+), 34 deletions(-) diff --git a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base.xml index be4523c1da2..f9c78d0843f 100644 --- a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base.xml +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base.xml @@ -15,14 +15,14 @@ directParallel="0"/> - + mpiCommOrder="1"/> --> diff --git a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml index 2088dfd7d56..45a837b1436 100644 --- a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml @@ -41,18 +41,18 @@ - + target="/Solvers/SurfaceGenerator"/> --> - + target="/Outputs/vtkOutput"/> --> forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion & elemRegion ) { - elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); + // testing only + // + CellBlockManagerABC & cellBlockManagerNoConst = const_cast(cellBlockManager); + CellBlockManager & cellBlockManagerConcrete = dynamic_cast(cellBlockManagerNoConst); + string name = "EmbeddedSurface"; + createDummyEmbeddedSurfaceBlock(name, cellBlockManagerConcrete); + + //elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); + + 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 83f603d9992..9937edbfbb9 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -262,6 +262,51 @@ bool EmbeddedSurfaceSubRegion::addNewEmbeddedSurface( localIndex const cellIndex return addEmbeddedElem; } + +bool EmbeddedSurfaceSubRegion::copyFromCellBlock(EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock){ + + localIndex const numElems= embeddedSurfaceBlock.numEmbeddedSurfElem(); + + resize( embeddedSurfaceBlock.numEmbeddedSurfElem() ); + + ArrayOfArrays< localIndex > elemTo3dElem = embeddedSurfaceBlock.getEmbeddedSurfElemTo3dElem(); + ArrayOfArrays< localIndex > elemToNodes = embeddedSurfaceBlock.getEmbeddedSurfElemToNodes(); + ArrayOfArrays elemNodesLocations = embeddedSurfaceBlock.getEmbeddedSurfElemNodes(); + ArrayOfArrays elemNormalVectors = embeddedSurfaceBlock.getEmbeddedSurfElemNormalVectors(); + ArrayOfArrays elemTangentialVectors1 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialWidthVectors(); + ArrayOfArrays elemTangentialVectors2 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialLengthVectors(); + + + + + for(auto i = 0; i elemNodeCoords(4,3); + for( localIndex inode = 0; inode < 4; inode++ ) + { + auto nodeGlobalIndex = elemToNodes[elemIndex][ inode ]; + m_toNodesRelation( elemIndex, inode ) = nodeGlobalIndex; + LvArray::tensorOps::copy<3>(elemNodeCoords[inode], elemNodesLocations[nodeGlobalIndex]); + } + + m_parentPlaneName[elemIndex] = "elem_" + std::to_string(elemIndex); + LvArray::tensorOps::copy<3>(m_normalVector[elemIndex], elemNormalVectors[elemIndex]); + LvArray::tensorOps::copy<3>(m_tangentVector1[elemIndex], elemTangentialVectors1[elemIndex]); + LvArray::tensorOps::copy<3>(m_tangentVector2[elemIndex], elemTangentialVectors2[elemIndex]); + this->calculateElementGeometricQuantities(elemNodeCoords.toViewConst(), elemIndex); + } + return true; +} + void EmbeddedSurfaceSubRegion::inheritGhostRank( array1d< array1d< arrayView1d< integer const > > > const & cellGhostRank ) { arrayView1d< integer > const & ghostRank = this->ghostRank(); diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp index 311b8aec610..60abe76f7f5 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,12 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion FixedOneToManyRelation const & cellToEdges, PlanarGeometricObject const * fracture ); + /** + * @brief Fill @p EmbeddedSurfaceSubRegion by copying the data from the source embedded surface block + * + * **/ + bool copyFromCellBlock(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 78233703bd1..36293fadb59 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -56,7 +56,16 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) { - elementSubRegions.registerGroup< EmbeddedSurfaceSubRegion >( m_faceBlockName ); + EmbeddedSurfaceSubRegion & subRegion = elementSubRegions.registerGroup< EmbeddedSurfaceSubRegion >( m_faceBlockName ); + if( faceBlocks.hasGroup( "EmbeddedSurface" ) ) + { + EmbeddedSurfaceBlockABC const & source = faceBlocks.getGroup< EmbeddedSurfaceBlockABC >( "EmbeddedSurface" ); + subRegion.copyFromCellBlock( source ); + } + else + { + GEOS_LOG_RANK_0( "No face block \"" << m_faceBlockName << "\" was found in the mesh. Empty surface region was created." ); + } } else if( m_subRegionType == SurfaceSubRegionType::faceElement ) { diff --git a/src/coreComponents/mesh/SurfaceElementRegion.hpp b/src/coreComponents/mesh/SurfaceElementRegion.hpp index 182d4f730b0..feb4d4633eb 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.hpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.hpp @@ -224,6 +224,7 @@ ENUM_STRINGS( SurfaceElementRegion::SurfaceSubRegionType, "faceElement", "embeddedElement" ); +EmbeddedSurfaceBlockABC & createDummyEmbeddedSurfaceBlock(string embeddedSurfFace); } /* namespace geos */ #endif /* CORECOMPONENTS_MESH_SurfaceElementRegion_HPP_ */ diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index 8c0f468b484..b24188abc09 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.cpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp @@ -670,7 +670,7 @@ void CellBlockManager::buildMaps() fillElementToEdgesOfCellBlocks( m_faceToEdges.toViewConst(), this->getCellBlocks() ); - buildEmbeddedSurfaceMaps(); + //buildEmbeddedSurfaceMaps(); } ArrayOfArrays< localIndex > CellBlockManager::getFaceToNodes() const @@ -690,7 +690,7 @@ ToCellRelation< array2d< localIndex > > CellBlockManager::getFaceToElements() co return m_faceToCells; } -ToCellRelation< ArrayOfArrays< localIndex > > getEmbeddedSurfaceToElements() const { +ToCellRelation< localIndex > CellBlockManager::getEmbeddedSurfaceToElements() const { return m_embeddedSurfToCells; } @@ -715,6 +715,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() ); @@ -747,7 +756,7 @@ localIndex CellBlockManager::numFaces() const localIndex CellBlockManager::numEmbeddedSurfaces() const { - return m_numEmbeddedSurfaces; + return m_numEmbeddedSurfElem; } ArrayOfArrays< localIndex > CellBlockManager::getEdgeToFaces() const @@ -802,7 +811,7 @@ LineBlock & CellBlockManager::registerLineBlock( string const & name ) EmbeddedSurfaceBlock & CellBlockManager::registerEmbeddedSurfaceBlock( string const & name ){ - return this->getLineBlocks().registerGroup< EmbeddedSurfaceBlock >( name ); + return this->getEmbeddedSurfaceBlocks().registerGroup< EmbeddedSurfaceBlock >( name ); } array2d< real64, nodes::REFERENCE_POSITION_PERM > CellBlockManager::getNodePositions() const diff --git a/src/coreComponents/mesh/generators/CellBlockManager.hpp b/src/coreComponents/mesh/generators/CellBlockManager.hpp index 5f3cb3dbae4..64a80422d26 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.hpp @@ -82,7 +82,7 @@ class CellBlockManager : public CellBlockManagerABC ToCellRelation< array2d< localIndex > > getFaceToElements() const override; - ToCellRelation< ArrayOfArrays< localIndex > > getEmbeddedSurfaceToElements() const override; + ToCellRelation< localIndex > getEmbeddedSurfaceToElements() const override; array1d< globalIndex > getNodeLocalToGlobal() const override; @@ -170,7 +170,7 @@ class CellBlockManager : public CellBlockManagerABC Group & getFaceBlocks() override; - Group & const getEmbeddedSurfaceBlocks() const override; + Group const & getEmbeddedSurfaceBlocks() const override; Group & getEmbeddedSurfaceBlocks() override; @@ -297,7 +297,7 @@ class CellBlockManager : public CellBlockManagerABC ArrayOfArrays< localIndex > m_faceToEdges; ToCellRelation< array2d< localIndex > > m_faceToCells; - ToCellRelation< ArrayOfArrays< localIndex > > m_embeddedSurfToCells; + ToCellRelation< localIndex > m_embeddedSurfToCells; array1d< globalIndex > m_nodeLocalToGlobal; @@ -310,7 +310,7 @@ class CellBlockManager : public CellBlockManagerABC localIndex m_numNodes; localIndex m_numFaces; localIndex m_numEdges; - localIndex m_numEmbeddedSurfaces; + localIndex m_numEmbeddedSurfElem; }; } diff --git a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp index 0175293ccd5..bb6044aa7cd 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 @@ -213,6 +221,12 @@ class CellBlockManagerABC : public dataRepository::Group */ virtual ToCellRelation< array2d< localIndex > > getFaceToElements() const = 0; + /** + * @brief Returns the embedded surface to elements mapping. + * @return A 1 to 1 relationship. The result is mapping from 1 fracture element to its parent matrix element. + * + */ + virtual ToCellRelation< localIndex > getEmbeddedSurfaceToElements() const = 0; /** * @brief The node to global mapping for nodes. * @return The mapping as an array of size numNodes. diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp index 7931478b6d7..8edd794d75f 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp @@ -18,26 +18,42 @@ namespace geos { - - - localIndex EmbeddedSurfaceBlock::numEmbeddedSurfElem() const{ return m_numEmbeddedSurfaces; } -ArrayOfArrays EmbeddedSurfaceBLock::getEmbeddedSurfElemToNodes() const { +ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemToNodes() const { return m_embeddedSurfElemToNodes; } -ToCellRelation> EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const { +ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const { return m_embeddedSurfElemTo3dElem; -}; +} ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNodes() const { return m_embeddedSurfElemNodes; } +ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNormalVectors() const { + return m_embeddedSurfElemNormals; +} +ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemTangentialLengthVectors() const{ + return m_embeddedSurfElemLengthVectors; +} +ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemTangentialWidthVectors() const{ + return m_embeddedSurfElemWidthVectors; +} + +void EmbeddedSurfaceBlock::setEmbeddedSurfElemNormalVectors(ArrayOfArrays&& _normals){ + m_embeddedSurfElemNormals = _normals; +} +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialLengthVectors(ArrayOfArrays && _lengthVectors){ + m_embeddedSurfElemLengthVectors = _lengthVectors; +} +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialWidthVectors(ArrayOfArrays && _widthVectors){ + m_embeddedSurfElemWidthVectors= _widthVectors; +} void EmbeddedSurfaceBlock::setNumEmbeddedSurfElem(localIndex _numEmbeddedSurfaces){ @@ -48,7 +64,7 @@ void EmbeddedSurfaceBlock::setEmbeddedSurfElemToNodes(ArrayOfArrays m_embeddedSurfElemToNodes = _embeddedSurfElemToNodes; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem(ToCellRelation> && _embeddedSurfElemTo3dElem){ +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem(ArrayOfArrays && _embeddedSurfElemTo3dElem){ m_embeddedSurfElemTo3dElem = _embeddedSurfElemTo3dElem; } diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp index e51c4a36a97..f2c8f6fa26e 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -41,8 +41,11 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC localIndex numEmbeddedSurfElem() const override; ArrayOfArrays getEmbeddedSurfElemToNodes() const override; - ToCellRelation> getEmbeddedSurfElemTo3dElem() const override; + ArrayOfArrays getEmbeddedSurfElemTo3dElem() const override; ArrayOfArrays getEmbeddedSurfElemNodes() const override; + ArrayOfArrays getEmbeddedSurfElemNormalVectors() const override; + ArrayOfArrays getEmbeddedSurfElemTangentialLengthVectors() const override; + ArrayOfArrays getEmbeddedSurfElemTangentialWidthVectors() const override; /** @@ -58,23 +61,44 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC void setEmbeddedSurfElemToNodes(ArrayOfArrays && _embeddedSurfElemToNodes); /** - * @brief Sets the embedded elements to 3d elements mapping + * @brief Sets the embedded elements to 1 3d element mapping * @param _embeddedSurfElemTo3dElem the input mapping. */ - void setEmbeddedSurfElemTo3dElem(ToCellRelation> && _embeddedSurfElemTo3dElem); + void setEmbeddedSurfElemTo3dElem(ArrayOfArrays && _embeddedSurfElemTo3dElem); /** * @brief Sets the embedded elements nodes coordinates * @param _embeddedSurfElemNodes the coordinates of embedded elements nodes. */ void setEmbeddedSurfElemNodes(ArrayOfArrays && _embeddedSurfElemNodes); + + /** + * @brief Sets the embedded elements normals + * @param _normals the coordinates of embedded elements normals. + */ + void setEmbeddedSurfElemNormalVectors(ArrayOfArrays && _normals); + + /** + * @brief Sets the embedded elements tangential length vectors (horizontal plan) + * @param _lengthVectors the coordinates of embedded elements tangential length vectors. + */ + void setEmbeddedSurfElemTangentialLengthVectors(ArrayOfArrays && _lengthVectors); + + /** + * @brief Sets the embedded elements tangential width vectors (vertical plan) + * @param _widthVectors the coordinates of embedded elements tangential width vectors. + */ + void setEmbeddedSurfElemTangentialWidthVectors(ArrayOfArrays && _widthVectors); private: localIndex m_numEmbeddedSurfaces; ArrayOfArrays< localIndex > m_embeddedSurfElemToNodes; - ToCellRelation> m_embeddedSurfElemTo3dElem; + ArrayOfArrays m_embeddedSurfElemTo3dElem; ArrayOfArrays< real64 > m_embeddedSurfElemNodes; + ArrayOfArrays< real64 > m_embeddedSurfElemNormals; + ArrayOfArrays< real64 > m_embeddedSurfElemLengthVectors; + ArrayOfArrays< real64 > m_embeddedSurfElemWidthVectors; }; } diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp index 72d0f9fb3a0..90b70fe1075 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -70,15 +70,14 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group virtual ArrayOfArrays getEmbeddedSurfElemToNodes() const = 0; /** - * @brief Get the indices of the parent 3d elements of each embedded surface element (1 or more) + * @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 number of 3d elements, and depends on how many 3d elements does the embedded surface element - * intersect. + * Second dimension 1 to 1 3d element that it intersects. * - * @details each embedded surface element intersects 1 or more 3d elements. Indices of these 3d elements - * are returned for each embedded surface element + * @details each embedded surface element intersects 1 3d element only. Index of this 3d element + * is returned for each embedded surface element. */ - virtual ToCellRelation> getEmbeddedSurfElemTo3dElem() const = 0; + virtual ArrayOfArrays getEmbeddedSurfElemTo3dElem() const = 0; /** @@ -88,6 +87,33 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group * Second dimension is 3, and represents the x, y and z. */ virtual ArrayOfArrays getEmbeddedSurfElemNodes() 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 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 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 getEmbeddedSurfElemTangentialWidthVectors() const = 0; }; } From c42b37f0e9bf3cedb7ea63788aa4eed5fcb4bb4e Mon Sep 17 00:00:00 2001 From: Ouassim Date: Wed, 17 Apr 2024 11:32:01 -0500 Subject: [PATCH 11/32] small changes --- .../fractureMatrixFlow_edfm_base.xml | 4 +- .../fractureMatrixFlow_edfm_base_modified.xml | 137 ++++++++++++++++++ ...reMatrixFlow_edfm_horizontalFrac_smoke.xml | 8 +- ...low_edfm_horizontalFrac_smoke_modified.xml | 73 ++++++++++ .../mesh/EmbeddedSurfaceSubRegion.cpp | 28 ++-- .../mesh/SurfaceElementRegion.cpp | 18 ++- .../mesh/generators/EmbeddedSurfaceBlock.cpp | 10 +- .../mesh/generators/EmbeddedSurfaceBlock.hpp | 10 +- .../generators/EmbeddedSurfaceBlockABC.hpp | 17 ++- 9 files changed, 263 insertions(+), 42 deletions(-) create mode 100644 inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base_modified.xml create mode 100644 inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke_modified.xml diff --git a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base.xml index f9c78d0843f..be4523c1da2 100644 --- a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base.xml +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base.xml @@ -15,14 +15,14 @@ directParallel="0"/> - + mpiCommOrder="1"/> diff --git a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base_modified.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base_modified.xml new file mode 100644 index 00000000000..d928a856304 --- /dev/null +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base_modified.xml @@ -0,0 +1,137 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml index 45a837b1436..2088dfd7d56 100644 --- a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml @@ -41,18 +41,18 @@ - + target="/Solvers/SurfaceGenerator"/> - + target="/Outputs/vtkOutput"/> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index 9937edbfbb9..e43fae47ae7 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -269,9 +269,9 @@ bool EmbeddedSurfaceSubRegion::copyFromCellBlock(EmbeddedSurfaceBlockABC const & resize( embeddedSurfaceBlock.numEmbeddedSurfElem() ); - ArrayOfArrays< localIndex > elemTo3dElem = embeddedSurfaceBlock.getEmbeddedSurfElemTo3dElem(); + ToCellRelation> elemTo3dElem = embeddedSurfaceBlock.getEmbeddedSurfElemTo3dElem(); ArrayOfArrays< localIndex > elemToNodes = embeddedSurfaceBlock.getEmbeddedSurfElemToNodes(); - ArrayOfArrays elemNodesLocations = embeddedSurfaceBlock.getEmbeddedSurfElemNodes(); + ArrayOfArrays elemNodesLocations = embeddedSurfaceBlock.getEmbeddedSurfElemNodesCoords(); ArrayOfArrays elemNormalVectors = embeddedSurfaceBlock.getEmbeddedSurfElemNormalVectors(); ArrayOfArrays elemTangentialVectors1 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialWidthVectors(); ArrayOfArrays elemTangentialVectors2 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialLengthVectors(); @@ -282,27 +282,27 @@ bool EmbeddedSurfaceSubRegion::copyFromCellBlock(EmbeddedSurfaceBlockABC const & for(auto i = 0; i elemNodeCoords(4,3); for( localIndex inode = 0; inode < 4; inode++ ) { - auto nodeGlobalIndex = elemToNodes[elemIndex][ inode ]; - m_toNodesRelation( elemIndex, inode ) = nodeGlobalIndex; + auto nodeGlobalIndex = elemToNodes[i][ inode ]; + m_toNodesRelation( i, inode ) = nodeGlobalIndex; LvArray::tensorOps::copy<3>(elemNodeCoords[inode], elemNodesLocations[nodeGlobalIndex]); } - m_parentPlaneName[elemIndex] = "elem_" + std::to_string(elemIndex); - LvArray::tensorOps::copy<3>(m_normalVector[elemIndex], elemNormalVectors[elemIndex]); - LvArray::tensorOps::copy<3>(m_tangentVector1[elemIndex], elemTangentialVectors1[elemIndex]); - LvArray::tensorOps::copy<3>(m_tangentVector2[elemIndex], elemTangentialVectors2[elemIndex]); - this->calculateElementGeometricQuantities(elemNodeCoords.toViewConst(), elemIndex); + 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; } diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index 36293fadb59..ffb90aff3d2 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -56,15 +56,23 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) { - EmbeddedSurfaceSubRegion & subRegion = elementSubRegions.registerGroup< EmbeddedSurfaceSubRegion >( m_faceBlockName ); - if( faceBlocks.hasGroup( "EmbeddedSurface" ) ) + bool experiment = false; + if (experiment) { - EmbeddedSurfaceBlockABC const & source = faceBlocks.getGroup< EmbeddedSurfaceBlockABC >( "EmbeddedSurface" ); - subRegion.copyFromCellBlock( source ); + EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); + if (faceBlocks.hasGroup("EmbeddedSurface")) + { + EmbeddedSurfaceBlockABC const &source = faceBlocks.getGroup("EmbeddedSurface"); + subRegion.copyFromCellBlock(source); + } + else + { + GEOS_LOG_RANK_0("No face block \"" << m_faceBlockName << "\" was found in the mesh. Empty surface region was created."); + } } else { - GEOS_LOG_RANK_0( "No face block \"" << m_faceBlockName << "\" was found in the mesh. Empty surface region was created." ); + elementSubRegions.registerGroup(m_faceBlockName); } } else if( m_subRegionType == SurfaceSubRegionType::faceElement ) diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp index 8edd794d75f..0b86848bd34 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp @@ -27,12 +27,12 @@ ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemToNodes() con return m_embeddedSurfElemToNodes; } -ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const { +ToCellRelation> EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const { return m_embeddedSurfElemTo3dElem; } -ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNodes() const { - return m_embeddedSurfElemNodes; +ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNodesCoords() const { + return m_embeddedSurfElemNodesCoords; } ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNormalVectors() const { @@ -64,12 +64,12 @@ void EmbeddedSurfaceBlock::setEmbeddedSurfElemToNodes(ArrayOfArrays m_embeddedSurfElemToNodes = _embeddedSurfElemToNodes; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem(ArrayOfArrays && _embeddedSurfElemTo3dElem){ +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem(ToCellRelation< ArrayOfArrays> && _embeddedSurfElemTo3dElem){ m_embeddedSurfElemTo3dElem = _embeddedSurfElemTo3dElem; } void EmbeddedSurfaceBlock::setEmbeddedSurfElemNodes(ArrayOfArrays && _embeddedSurfElemNodes){ - m_embeddedSurfElemNodes = _embeddedSurfElemNodes; + m_embeddedSurfElemNodesCoords = _embeddedSurfElemNodes; } } \ No newline at end of file diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp index f2c8f6fa26e..44f6b2bf632 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -41,8 +41,8 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC localIndex numEmbeddedSurfElem() const override; ArrayOfArrays getEmbeddedSurfElemToNodes() const override; - ArrayOfArrays getEmbeddedSurfElemTo3dElem() const override; - ArrayOfArrays getEmbeddedSurfElemNodes() const override; + ToCellRelation> getEmbeddedSurfElemTo3dElem() const override; + ArrayOfArrays getEmbeddedSurfElemNodesCoords() const override; ArrayOfArrays getEmbeddedSurfElemNormalVectors() const override; ArrayOfArrays getEmbeddedSurfElemTangentialLengthVectors() const override; ArrayOfArrays getEmbeddedSurfElemTangentialWidthVectors() const override; @@ -64,7 +64,7 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC * @brief Sets the embedded elements to 1 3d element mapping * @param _embeddedSurfElemTo3dElem the input mapping. */ - void setEmbeddedSurfElemTo3dElem(ArrayOfArrays && _embeddedSurfElemTo3dElem); + void setEmbeddedSurfElemTo3dElem(ToCellRelation> && _embeddedSurfElemTo3dElem); /** * @brief Sets the embedded elements nodes coordinates @@ -94,8 +94,8 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC localIndex m_numEmbeddedSurfaces; ArrayOfArrays< localIndex > m_embeddedSurfElemToNodes; - ArrayOfArrays m_embeddedSurfElemTo3dElem; - ArrayOfArrays< real64 > m_embeddedSurfElemNodes; + ToCellRelation> m_embeddedSurfElemTo3dElem; + ArrayOfArrays< real64 > m_embeddedSurfElemNodesCoords; ArrayOfArrays< real64 > m_embeddedSurfElemNormals; ArrayOfArrays< real64 > m_embeddedSurfElemLengthVectors; ArrayOfArrays< real64 > m_embeddedSurfElemWidthVectors; diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp index 90b70fe1075..57da75c2564 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -28,7 +28,8 @@ 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 or more 3d element). The 2d element is assumed to be a quad. + * 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 @@ -54,7 +55,7 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group * @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 or N 3d elements + * can intersect 1 3d elements */ virtual localIndex numEmbeddedSurfElem() const = 0; @@ -73,20 +74,22 @@ class EmbeddedSurfaceBlockABC : public dataRepository::Group * @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. Index of this 3d element - * is returned for each embedded surface element. + * @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 ArrayOfArrays getEmbeddedSurfElemTo3dElem() const = 0; + virtual ToCellRelation> 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 @p numEmbeddedSurfaceElements*4. + * 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 getEmbeddedSurfElemNodes() const = 0; + virtual ArrayOfArrays getEmbeddedSurfElemNodesCoords() const = 0; /** * @brief Get all the normal vectors for all the embedded surface elements. From 55239bded6e8a946e11c2758a712d5cebaf4c9b5 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Thu, 18 Apr 2024 12:46:13 -0500 Subject: [PATCH 12/32] add flag to enable experimentation --- src/coreComponents/mesh/ElementRegionManager.cpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index 7b8f40a10e3..459b4700eae 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -140,14 +140,22 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa { // testing only // + + bool experiment = true; + if (experiment) + { + CellBlockManagerABC & cellBlockManagerNoConst = const_cast(cellBlockManager); CellBlockManager & cellBlockManagerConcrete = dynamic_cast(cellBlockManagerNoConst); + string name = "EmbeddedSurface"; createDummyEmbeddedSurfaceBlock(name, cellBlockManagerConcrete); - - //elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); - elemRegion.generateMesh( cellBlockManager.getEmbeddedSurfaceBlocks() ); + } + else{ + elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); + } + } ); // Some mappings of the surfaces subregions point to elements in other subregions and regions. From 4b2363e71b5a907612f0cfdde4edad83ce1c4147 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Thu, 18 Apr 2024 12:46:55 -0500 Subject: [PATCH 13/32] missed files --- .../mesh/EmbeddedSurfaceSubRegion.cpp | 16 ++++++++++++++++ src/coreComponents/mesh/SurfaceElementRegion.cpp | 5 +++-- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index e43fae47ae7..395a03f0fb6 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -27,6 +27,16 @@ #include "BufferOps.hpp" #include "mesh/MeshFields.hpp" + +#include "mainInterface/GeosxState.hpp" +#include "common/DataTypes.hpp" +#include "linearAlgebra/DofManager.hpp" +#include "mainInterface/initialization.hpp" +#include "mainInterface/ProblemManager.hpp" +#include "mesh/DomainPartition.hpp" +#include "mesh/MeshManager.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" + namespace geos { using namespace dataRepository; @@ -265,6 +275,12 @@ bool EmbeddedSurfaceSubRegion::addNewEmbeddedSurface( localIndex const cellIndex bool EmbeddedSurfaceSubRegion::copyFromCellBlock(EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock){ + + DomainPartition & domain = getGlobalState().getProblemManager().getDomainPartition(); + MeshLevel & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); + EmbeddedSurfaceNodeManager & embSurfNodeManager = meshLevel.getEmbSurfNodeManager(); + EdgeManager & edgeManager = meshLevel.getEdgeManager(); + localIndex const numElems= embeddedSurfaceBlock.numEmbeddedSurfElem(); resize( embeddedSurfaceBlock.numEmbeddedSurfElem() ); diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index ffb90aff3d2..9c523adc133 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -21,7 +21,6 @@ #include "EdgeManager.hpp" #include "SurfaceElementRegion.hpp" - namespace geos { using namespace dataRepository; @@ -54,11 +53,13 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) { Group & elementSubRegions = this->getGroup( viewKeyStruct::elementSubRegions() ); + if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) { - bool experiment = false; + bool experiment = true; if (experiment) { + EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); if (faceBlocks.hasGroup("EmbeddedSurface")) { From 0fde8c6eedfc901a91aa010299aa89b0541299d6 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Fri, 26 Apr 2024 16:12:14 -0500 Subject: [PATCH 14/32] add bulk edfm fracs --- src/coreComponents/mesh/CMakeLists.txt | 3 + .../mesh/ElementRegionManager.cpp | 4 +- .../mesh/EmbeddedSurfaceSubRegion.cpp | 175 +++++++++++++++++- .../mesh/EmbeddedSurfaceSubRegion.hpp | 18 ++ .../mesh/SurfaceElementRegion.cpp | 40 ++-- .../EmbeddedSurfaceGenerator.cpp | 50 ++++- 6 files changed, 263 insertions(+), 27 deletions(-) diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt index 9963d8b1a7f..6375ef6fcd5 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -209,3 +209,6 @@ if( GEOS_ENABLE_TESTS ) add_subdirectory( unitTests ) endif( ) + +#add_dependencies(mesh mainInterface) +#add_dependencies(mesh linearAlgebra) \ No newline at end of file diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index 459b4700eae..e6d75e41e3f 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -141,7 +141,7 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa // testing only // - bool experiment = true; + bool experiment = false; if (experiment) { @@ -149,7 +149,7 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa CellBlockManager & cellBlockManagerConcrete = dynamic_cast(cellBlockManagerNoConst); string name = "EmbeddedSurface"; - createDummyEmbeddedSurfaceBlock(name, cellBlockManagerConcrete); + createDummyEmbeddedSurfaceBlockOrig(name, cellBlockManagerConcrete); elemRegion.generateMesh( cellBlockManager.getEmbeddedSurfaceBlocks() ); } else{ diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index 395a03f0fb6..0c4a4e1df67 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -29,13 +29,13 @@ #include "mainInterface/GeosxState.hpp" -#include "common/DataTypes.hpp" -#include "linearAlgebra/DofManager.hpp" -#include "mainInterface/initialization.hpp" +//#include "common/DataTypes.hpp" +//#include "linearAlgebra/DofManager.hpp" +//#include "mainInterface/initialization.hpp" #include "mainInterface/ProblemManager.hpp" #include "mesh/DomainPartition.hpp" -#include "mesh/MeshManager.hpp" -#include "mesh/mpiCommunications/CommunicationTools.hpp" +//#include "mesh/MeshManager.hpp" +//#include "mesh/mpiCommunications/CommunicationTools.hpp" namespace geos { @@ -271,7 +271,129 @@ bool EmbeddedSurfaceSubRegion::addNewEmbeddedSurface( localIndex const cellIndex } return addEmbeddedElem; } +array1d EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex(ArrayOfArraysView const &elemNodesLocations, + ToCellRelation> const &elemTo3dElem, + FixedOneToManyRelation const &cellToEdges, + arrayView2d const &edgeToNodes, + arrayView2d const &nodesCoord) +{ + array1d parentEdgeIndex(elemNodesLocations.size()); + for (localIndex i = 0; i < elemNodesLocations.size(); ++i) + { + + auto cell3dIndex = elemTo3dElem.toCellIndex[i][0]; + double elemMinDistance = std::numeric_limits::max(); + localIndex targetEdgeIndex; + + 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 const &nodesCoord = nodeManager.referencePosition(); + arrayView2d const edgeToNodes = edgeManager.nodeList(); + arrayView1d const edgeGhostRank = edgeManager.ghostRank(); + arrayView1d const &edgeLocalToGlobal = edgeManager.localToGlobalMap(); + + localIndex const numElems = embeddedSurfaceBlock.numEmbeddedSurfElem(); + + resize(embeddedSurfaceBlock.numEmbeddedSurfElem()); + + ToCellRelation> elemTo3dElem = embeddedSurfaceBlock.getEmbeddedSurfElemTo3dElem(); + ArrayOfArrays elemToNodes = embeddedSurfaceBlock.getEmbeddedSurfElemToNodes(); + ArrayOfArrays elemNodesLocations = embeddedSurfaceBlock.getEmbeddedSurfElemNodesCoords(); + ArrayOfArrays elemNormalVectors = embeddedSurfaceBlock.getEmbeddedSurfElemNormalVectors(); + ArrayOfArrays elemTangentialVectors1 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialWidthVectors(); + ArrayOfArrays elemTangentialVectors2 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialLengthVectors(); + + localIndex numEdfmNodes = elemNodesLocations.size(); + array1d allPointsParentIndex = getEdfmNodeParentEdgeIndex(elemNodesLocations.toView(),elemTo3dElem, cellToEdges, edgeToNodes, nodesCoord); + + // EDFM nodes get their gost rank and global parent index + // from their parent 3d element local edge index. + array1d allPointsGhostRank(numEdfmNodes); + array1d 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 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; + +} bool EmbeddedSurfaceSubRegion::copyFromCellBlock(EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock){ @@ -281,6 +403,8 @@ bool EmbeddedSurfaceSubRegion::copyFromCellBlock(EmbeddedSurfaceBlockABC const & EmbeddedSurfaceNodeManager & embSurfNodeManager = meshLevel.getEmbSurfNodeManager(); EdgeManager & edgeManager = meshLevel.getEdgeManager(); + ElementRegionManager & elemManager = meshLevel.getElemManager(); + localIndex const numElems= embeddedSurfaceBlock.numEmbeddedSurfElem(); resize( embeddedSurfaceBlock.numEmbeddedSurfElem() ); @@ -296,6 +420,7 @@ bool EmbeddedSurfaceSubRegion::copyFromCellBlock(EmbeddedSurfaceBlockABC const & for(auto i = 0; i(m_tangentVector2[i], elemTangentialVectors2[i]); this->calculateElementGeometricQuantities(elemNodeCoords.toViewConst(), i); } + + // add nodes to embNodeManager + for(auto i =0; i < elemNodesLocations.size(); ++i) + { + + array2d elemNodeCoords(1,3); + auto ghostRank = -2; + + // // Add the point to the node Manager + // globalIndex parentEdgeID = edgeLocalToGlobal[ pointParentIndex[ originalIndices[ j ] ] ]; + // nodeIndex = embSurfNodeManager.size(); + embSurfNodeManager.appendNode(elemNodesLocations[i].toSliceConst(), ghostRank); + + // arrayView1d< localIndex > const & parentIndex = + // embSurfNodeManager.getField< fields::parentEdgeIndex >(); + + // parentIndex[nodeIndex] = pointParentIndex[ originalIndices[ j ] ]; + + // array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); + // parentEdgeGlobalIndex[nodeIndex] = parentEdgeID; + } + + + elemManager.forElementSubRegionsComplete< CellElementSubRegion >( + [&] ( localIndex, localIndex, ElementRegionBase const & region, CellElementSubRegion & subRegion ) + { + + for(auto i = 0; i > > const & cellGhostRank ) diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp index 60abe76f7f5..e0c01511dea 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp @@ -30,6 +30,10 @@ //Do we really need this include Rectangle? #include "simpleGeometricObjects/Rectangle.hpp" +//#include "mainInterface/GeosxState.hpp" +//#include "mainInterface/ProblemManager.hpp" +//#include "mesh/DomainPartition.hpp" + namespace geos { @@ -147,6 +151,20 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion FixedOneToManyRelation const & cellToEdges, PlanarGeometricObject const * fracture ); +array1d getEdfmNodeParentEdgeIndex(ArrayOfArraysView const &elemNodesLocations, + ToCellRelation> const &elemTo3dElem, + FixedOneToManyRelation const &cellToEdges, + arrayView2d const &edgeToNodes, + arrayView2d const &nodesCoord); + +bool addAllEmbeddedSurfaces( localIndex const regionIndex, + localIndex const subRegionIndex, + NodeManager const & nodeManager, + EmbeddedSurfaceNodeManager & embSurfNodeManager, + EdgeManager const & edgeManager, + FixedOneToManyRelation const & cellToEdges, + EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock); + /** * @brief Fill @p EmbeddedSurfaceSubRegion by copying the data from the source embedded surface block * diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index 9c523adc133..48b5c930e8e 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -56,25 +56,27 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) { - bool experiment = true; - if (experiment) - { - - EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); - if (faceBlocks.hasGroup("EmbeddedSurface")) - { - EmbeddedSurfaceBlockABC const &source = faceBlocks.getGroup("EmbeddedSurface"); - subRegion.copyFromCellBlock(source); - } - else - { - GEOS_LOG_RANK_0("No face block \"" << m_faceBlockName << "\" was found in the mesh. Empty surface region was created."); - } - } - else - { - elementSubRegions.registerGroup(m_faceBlockName); - } + // bool experiment =false; + // if (experiment) + // { + + // EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); + // if (faceBlocks.hasGroup("EmbeddedSurface")) + // { + // EmbeddedSurfaceBlockABC const &source = faceBlocks.getGroup("EmbeddedSurface"); + // subRegion.copyFromCellBlock(source); + // } + // else + // { + // GEOS_LOG_RANK_0("No face block \"" << m_faceBlockName << "\" was found in the mesh. Empty surface region was created."); + // } + // } + // else + // { + // elementSubRegions.registerGroup(m_faceBlockName); + // } + + elementSubRegions.registerGroup("EmbeddedSurface"); } else if( m_subRegionType == SurfaceSubRegionType::faceElement ) { diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index 7b13c4618d5..6cf5d925b10 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -113,6 +113,7 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() GeometricObjectManager & geometricObjManager = GeometricObjectManager::getInstance(); // Get meshLevel + MeshBody & meshBody = domain.getMeshBody( 0 ); MeshLevel & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); // Get managers @@ -129,6 +130,53 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() localIndex localNumberOfSurfaceElems = 0; NewObjectLists newObjects; + + bool bulk_edfm_fill = true; + if (bulk_edfm_fill){ + + CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager(); + CellBlockManagerABC & cellBlockManagerNoConst = const_cast(domain.getMeshBody(0).getCellBlockManager()); + CellBlockManager & cellBlockManagerConcrete = dynamic_cast(cellBlockManagerNoConst); + Group const & embSurfBlocks = cellBlockManagerConcrete.getEmbeddedSurfaceBlocks(); + if (embSurfBlocks.hasGroup("EmbeddedSurface")) + { + EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup("EmbeddedSurface"); + + + elemManager.forElementSubRegionsComplete< CellElementSubRegion >( + [&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion ) + { + arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList(); + FixedOneToManyRelation const & cellToEdges = subRegion.edgeList(); + + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + { + bool added = embeddedSurfaceSubRegion.addAllEmbeddedSurfaces( er, + esr, + nodeManager, + embSurfNodeManager, + edgeManager, + cellToEdges, + embSurf); + + 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, @@ -205,7 +253,7 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() } );// 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 >( [&]( localIndex const, localIndex const, ElementRegionBase &, CellElementSubRegion & subRegion ) From dd8d49e74afd9329eecbbb952ddc407c41d6c5b7 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Mon, 29 Apr 2024 15:30:28 -0500 Subject: [PATCH 15/32] All is running --- .../mainInterface/ProblemManager.cpp | 2 +- .../mesh/ElementRegionManager.cpp | 2 +- .../mesh/EmbeddedSurfaceSubRegion.cpp | 20 +++++++++++++++++-- .../mesh/EmbeddedSurfaceSubRegion.hpp | 1 + 4 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp index 0a134a6bcea..a4bf4f0dc37 100644 --- a/src/coreComponents/mainInterface/ProblemManager.cpp +++ b/src/coreComponents/mainInterface/ProblemManager.cpp @@ -691,7 +691,7 @@ void ProblemManager::generateMesh() else if( meshBody.hasGroup( keys::cellManager ) ) { // meshBody.deregisterGroup( keys::cellManager ); - meshBody.deregisterCellBlockManager(); + //meshBody.deregisterCellBlockManager(); } meshBody.forMeshLevels( [&]( MeshLevel & meshLevel ) diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index e6d75e41e3f..a1b0ac5de42 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -141,7 +141,7 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa // testing only // - bool experiment = false; + bool experiment = true; if (experiment) { diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index 0c4a4e1df67..e8e212b34fe 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -272,16 +272,32 @@ bool EmbeddedSurfaceSubRegion::addNewEmbeddedSurface( localIndex const cellIndex return addEmbeddedElem; } array1d EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex(ArrayOfArraysView const &elemNodesLocations, + ArrayOfArraysView const &elemToNodes, ToCellRelation> const &elemTo3dElem, FixedOneToManyRelation const &cellToEdges, arrayView2d const &edgeToNodes, arrayView2d const &nodesCoord) { + // first we build a node index to parent cell index map (we oly need one) + array1d nodeToElem(elemNodesLocations.size()); + for(int i = 0; i parentEdgeIndex(elemNodesLocations.size()); for (localIndex i = 0; i < elemNodesLocations.size(); ++i) { - auto cell3dIndex = elemTo3dElem.toCellIndex[i][0]; + auto elementIndex =nodeToElem[i]; + auto cell3dIndex = elemTo3dElem.toCellIndex[elementIndex][0]; double elemMinDistance = std::numeric_limits::max(); localIndex targetEdgeIndex; @@ -342,7 +358,7 @@ array1d EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex(ArrayOf ArrayOfArrays elemTangentialVectors2 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialLengthVectors(); localIndex numEdfmNodes = elemNodesLocations.size(); - array1d allPointsParentIndex = getEdfmNodeParentEdgeIndex(elemNodesLocations.toView(),elemTo3dElem, cellToEdges, edgeToNodes, nodesCoord); + array1d 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. diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp index e0c01511dea..aebbbaeb784 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp @@ -152,6 +152,7 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion PlanarGeometricObject const * fracture ); array1d getEdfmNodeParentEdgeIndex(ArrayOfArraysView const &elemNodesLocations, + ArrayOfArraysView const &elemToNodes, ToCellRelation> const &elemTo3dElem, FixedOneToManyRelation const &cellToEdges, arrayView2d const &edgeToNodes, From 2e43180198193b67ab368f3aa80a202f892bc64f Mon Sep 17 00:00:00 2001 From: Ouassim Date: Mon, 29 Apr 2024 15:36:26 -0500 Subject: [PATCH 16/32] uncrustify pass --- .../mesh/ElementRegionManager.cpp | 23 +- .../mesh/EmbeddedSurfaceSubRegion.cpp | 296 +++++++++--------- .../mesh/EmbeddedSurfaceSubRegion.hpp | 30 +- .../mesh/SurfaceElementRegion.cpp | 4 +- .../mesh/SurfaceElementRegion.hpp | 2 +- .../mesh/generators/CellBlockManager.cpp | 10 +- .../mesh/generators/EmbeddedSurfaceBlock.cpp | 76 +++-- .../mesh/generators/EmbeddedSurfaceBlock.hpp | 128 ++++---- .../generators/EmbeddedSurfaceBlockABC.hpp | 168 +++++----- .../EmbeddedSurfaceGenerator.cpp | 143 +++++++-- 10 files changed, 491 insertions(+), 389 deletions(-) diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index a1b0ac5de42..e32090eec47 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -140,22 +140,23 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa { // testing only // - + bool experiment = true; - if (experiment) + if( experiment ) { - - CellBlockManagerABC & cellBlockManagerNoConst = const_cast(cellBlockManager); - CellBlockManager & cellBlockManagerConcrete = dynamic_cast(cellBlockManagerNoConst); - string name = "EmbeddedSurface"; - createDummyEmbeddedSurfaceBlockOrig(name, cellBlockManagerConcrete); - elemRegion.generateMesh( cellBlockManager.getEmbeddedSurfaceBlocks() ); + CellBlockManagerABC & cellBlockManagerNoConst = const_cast< CellBlockManagerABC & >(cellBlockManager); + CellBlockManager & cellBlockManagerConcrete = dynamic_cast< CellBlockManager & >(cellBlockManagerNoConst); + + string name = "EmbeddedSurface"; + createDummyEmbeddedSurfaceBlockOrig( name, cellBlockManagerConcrete ); + elemRegion.generateMesh( cellBlockManager.getEmbeddedSurfaceBlocks() ); } - else{ - elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); + else + { + elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); } - + } ); // 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 e8e212b34fe..c2b906b7463 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -271,237 +271,239 @@ bool EmbeddedSurfaceSubRegion::addNewEmbeddedSurface( localIndex const cellIndex } return addEmbeddedElem; } -array1d EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex(ArrayOfArraysView const &elemNodesLocations, - ArrayOfArraysView const &elemToNodes, - ToCellRelation> const &elemTo3dElem, - FixedOneToManyRelation const &cellToEdges, - arrayView2d const &edgeToNodes, - arrayView2d const &nodesCoord) +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 oly need one) - array1d nodeToElem(elemNodesLocations.size()); - for(int i = 0; i nodeToElem( elemNodesLocations.size()); + for( int i = 0; i parentEdgeIndex(elemNodesLocations.size()); - for (localIndex i = 0; i < elemNodesLocations.size(); ++i) + array1d< localIndex > parentEdgeIndex( elemNodesLocations.size()); + for( localIndex i = 0; i < elemNodesLocations.size(); ++i ) { - auto elementIndex =nodeToElem[i]; + auto elementIndex =nodeToElem[i]; auto cell3dIndex = elemTo3dElem.toCellIndex[elementIndex][0]; - double elemMinDistance = std::numeric_limits::max(); + double elemMinDistance = std::numeric_limits< double >::max(); localIndex targetEdgeIndex; R1Tensor nodexyz; - LvArray::tensorOps::copy<3>(nodexyz, elemNodesLocations[i]); - for (localIndex ke = 0; ke < cellToEdges.size(1); ke++) + 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 >( 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::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); + 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) + if( dist < elemMinDistance ) { elemMinDistance = dist; targetEdgeIndex = edgeIndex; } } - parentEdgeIndex[i] = targetEdgeIndex; + parentEdgeIndex[i] = targetEdgeIndex; } - return parentEdgeIndex; + 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) - { +bool EmbeddedSurfaceSubRegion::addAllEmbeddedSurfaces( + localIndex const regionIndex, + localIndex const subRegionIndex, + NodeManager const & nodeManager, + EmbeddedSurfaceNodeManager & embSurfNodeManager, + EdgeManager const & edgeManager, + FixedOneToManyRelation const & cellToEdges, + EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock ) +{ - arrayView2d const &nodesCoord = nodeManager.referencePosition(); - arrayView2d const edgeToNodes = edgeManager.nodeList(); - arrayView1d const edgeGhostRank = edgeManager.ghostRank(); - arrayView1d const &edgeLocalToGlobal = edgeManager.localToGlobalMap(); + 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(); + localIndex const numElems = embeddedSurfaceBlock.numEmbeddedSurfElem(); - resize(embeddedSurfaceBlock.numEmbeddedSurfElem()); + resize( embeddedSurfaceBlock.numEmbeddedSurfElem()); - ToCellRelation> elemTo3dElem = embeddedSurfaceBlock.getEmbeddedSurfElemTo3dElem(); - ArrayOfArrays elemToNodes = embeddedSurfaceBlock.getEmbeddedSurfElemToNodes(); - ArrayOfArrays elemNodesLocations = embeddedSurfaceBlock.getEmbeddedSurfElemNodesCoords(); - ArrayOfArrays elemNormalVectors = embeddedSurfaceBlock.getEmbeddedSurfElemNormalVectors(); - ArrayOfArrays elemTangentialVectors1 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialWidthVectors(); - ArrayOfArrays elemTangentialVectors2 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialLengthVectors(); + 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]]; - localIndex numEdfmNodes = elemNodesLocations.size(); - array1d allPointsParentIndex = getEdfmNodeParentEdgeIndex(elemNodesLocations.toView(),elemToNodes.toView(),elemTo3dElem, cellToEdges, edgeToNodes, nodesCoord); + embSurfNodeManager.appendNode( elemNodesLocations[i].toSliceConst(), pointGhostRank ); - // EDFM nodes get their gost rank and global parent index - // from their parent 3d element local edge index. - array1d allPointsGhostRank(numEdfmNodes); - array1d allPointsGlobalParenIndex(numEdfmNodes); - for (auto i = 0; i < numEdfmNodes; ++i) - { - integer pointGhostRank = edgeGhostRank[allPointsParentIndex[i]]; - localIndex pointLocalParentEdgeIndex = allPointsParentIndex[i]; - globalIndex pointGlobalParentEdgeIndex = edgeLocalToGlobal[allPointsParentIndex[i]]; + arrayView1d< localIndex > const & parentIndex = embSurfNodeManager.getField< fields::parentEdgeIndex >(); - embSurfNodeManager.appendNode(elemNodesLocations[i].toSliceConst(), pointGhostRank); + parentIndex[i] = pointLocalParentEdgeIndex; - arrayView1d< localIndex > const & parentIndex = embSurfNodeManager.getField< fields::parentEdgeIndex >(); + array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); + parentEdgeGlobalIndex[i] = pointGlobalParentEdgeIndex; + } - parentIndex[i] = pointLocalParentEdgeIndex; + 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 ); - array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); - parentEdgeGlobalIndex[i] = pointGlobalParentEdgeIndex; + // 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] ); } - 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 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 ); + } - 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; - + } -bool EmbeddedSurfaceSubRegion::copyFromCellBlock(EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock){ - - +bool EmbeddedSurfaceSubRegion::copyFromCellBlock( EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock ) +{ + + DomainPartition & domain = getGlobalState().getProblemManager().getDomainPartition(); MeshLevel & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); EmbeddedSurfaceNodeManager & embSurfNodeManager = meshLevel.getEmbSurfNodeManager(); EdgeManager & edgeManager = meshLevel.getEdgeManager(); - + ElementRegionManager & elemManager = meshLevel.getElemManager(); - - localIndex const numElems= embeddedSurfaceBlock.numEmbeddedSurfElem(); + + localIndex const numElems= embeddedSurfaceBlock.numEmbeddedSurfElem(); resize( embeddedSurfaceBlock.numEmbeddedSurfElem() ); - - ToCellRelation> elemTo3dElem = embeddedSurfaceBlock.getEmbeddedSurfElemTo3dElem(); + + ToCellRelation< ArrayOfArrays< localIndex > > elemTo3dElem = embeddedSurfaceBlock.getEmbeddedSurfElemTo3dElem(); ArrayOfArrays< localIndex > elemToNodes = embeddedSurfaceBlock.getEmbeddedSurfElemToNodes(); - ArrayOfArrays elemNodesLocations = embeddedSurfaceBlock.getEmbeddedSurfElemNodesCoords(); - ArrayOfArrays elemNormalVectors = embeddedSurfaceBlock.getEmbeddedSurfElemNormalVectors(); - ArrayOfArrays elemTangentialVectors1 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialWidthVectors(); - ArrayOfArrays elemTangentialVectors2 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialLengthVectors(); - - - - - for(auto i = 0; i elemNodesLocations = embeddedSurfaceBlock.getEmbeddedSurfElemNodesCoords(); + ArrayOfArrays< real64 > elemNormalVectors = embeddedSurfaceBlock.getEmbeddedSurfElemNormalVectors(); + ArrayOfArrays< real64 > elemTangentialVectors1 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialWidthVectors(); + ArrayOfArrays< real64 > elemTangentialVectors2 = embeddedSurfaceBlock.getEmbeddedSurfElemTangentialLengthVectors(); + + + for( auto i = 0; i elemNodeCoords(4,3); + 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]); + 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); + + 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 ); } - + // add nodes to embNodeManager - for(auto i =0; i < elemNodesLocations.size(); ++i) + for( auto i =0; i < elemNodesLocations.size(); ++i ) { - array2d elemNodeCoords(1,3); + array2d< real64 > elemNodeCoords( 1, 3 ); auto ghostRank = -2; - // // Add the point to the node Manager - // globalIndex parentEdgeID = edgeLocalToGlobal[ pointParentIndex[ originalIndices[ j ] ] ]; - // nodeIndex = embSurfNodeManager.size(); - embSurfNodeManager.appendNode(elemNodesLocations[i].toSliceConst(), ghostRank); + // // Add the point to the node Manager + // globalIndex parentEdgeID = edgeLocalToGlobal[ pointParentIndex[ originalIndices[ j ] ] ]; + // nodeIndex = embSurfNodeManager.size(); + embSurfNodeManager.appendNode( elemNodesLocations[i].toSliceConst(), ghostRank ); - // arrayView1d< localIndex > const & parentIndex = - // embSurfNodeManager.getField< fields::parentEdgeIndex >(); + // arrayView1d< localIndex > const & parentIndex = + // embSurfNodeManager.getField< fields::parentEdgeIndex >(); + + // parentIndex[nodeIndex] = pointParentIndex[ originalIndices[ j ] ]; + + // array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); + // parentEdgeGlobalIndex[nodeIndex] = parentEdgeID; + } - // parentIndex[nodeIndex] = pointParentIndex[ originalIndices[ j ] ]; - // array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); - // parentEdgeGlobalIndex[nodeIndex] = parentEdgeID; - } - - elemManager.forElementSubRegionsComplete< CellElementSubRegion >( [&] ( localIndex, localIndex, ElementRegionBase const & region, CellElementSubRegion & subRegion ) { - for(auto i = 0; i > > const & cellGhostRank ) diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp index aebbbaeb784..7599210d76d 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp @@ -151,26 +151,26 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion FixedOneToManyRelation const & cellToEdges, PlanarGeometricObject const * fracture ); -array1d getEdfmNodeParentEdgeIndex(ArrayOfArraysView const &elemNodesLocations, - ArrayOfArraysView const &elemToNodes, - ToCellRelation> const &elemTo3dElem, - FixedOneToManyRelation const &cellToEdges, - arrayView2d const &edgeToNodes, - arrayView2d const &nodesCoord); - -bool addAllEmbeddedSurfaces( localIndex const regionIndex, - localIndex const subRegionIndex, - NodeManager const & nodeManager, - EmbeddedSurfaceNodeManager & embSurfNodeManager, - EdgeManager const & edgeManager, - FixedOneToManyRelation const & cellToEdges, - EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock); + 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 ); + + bool addAllEmbeddedSurfaces( localIndex const regionIndex, + localIndex const subRegionIndex, + NodeManager const & nodeManager, + EmbeddedSurfaceNodeManager & embSurfNodeManager, + EdgeManager const & edgeManager, + FixedOneToManyRelation const & cellToEdges, + EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock ); /** * @brief Fill @p EmbeddedSurfaceSubRegion by copying the data from the source embedded surface block * * **/ - bool copyFromCellBlock(EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock); + bool copyFromCellBlock( EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock ); /** * @brief inherit ghost rank from cell elements. diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index 48b5c930e8e..5f00d2bca1c 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -75,8 +75,8 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) // { // elementSubRegions.registerGroup(m_faceBlockName); // } - - elementSubRegions.registerGroup("EmbeddedSurface"); + + elementSubRegions.registerGroup< EmbeddedSurfaceSubRegion >( "EmbeddedSurface" ); } else if( m_subRegionType == SurfaceSubRegionType::faceElement ) { diff --git a/src/coreComponents/mesh/SurfaceElementRegion.hpp b/src/coreComponents/mesh/SurfaceElementRegion.hpp index feb4d4633eb..e2b27838db9 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.hpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.hpp @@ -224,7 +224,7 @@ ENUM_STRINGS( SurfaceElementRegion::SurfaceSubRegionType, "faceElement", "embeddedElement" ); -EmbeddedSurfaceBlockABC & createDummyEmbeddedSurfaceBlock(string embeddedSurfFace); +EmbeddedSurfaceBlockABC & createDummyEmbeddedSurfaceBlock( string embeddedSurfFace ); } /* namespace geos */ #endif /* CORECOMPONENTS_MESH_SurfaceElementRegion_HPP_ */ diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index b24188abc09..ec193806f7c 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.cpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp @@ -669,7 +669,7 @@ void CellBlockManager::buildMaps() buildNodeToEdges(); fillElementToEdgesOfCellBlocks( m_faceToEdges.toViewConst(), this->getCellBlocks() ); - + //buildEmbeddedSurfaceMaps(); } @@ -690,7 +690,8 @@ ToCellRelation< array2d< localIndex > > CellBlockManager::getFaceToElements() co return m_faceToCells; } -ToCellRelation< localIndex > CellBlockManager::getEmbeddedSurfaceToElements() const { +ToCellRelation< localIndex > CellBlockManager::getEmbeddedSurfaceToElements() const +{ return m_embeddedSurfToCells; } @@ -809,8 +810,9 @@ LineBlock & CellBlockManager::registerLineBlock( string const & name ) return this->getLineBlocks().registerGroup< LineBlock >( name ); } -EmbeddedSurfaceBlock & CellBlockManager::registerEmbeddedSurfaceBlock( string const & name ){ - +EmbeddedSurfaceBlock & CellBlockManager::registerEmbeddedSurfaceBlock( string const & name ) +{ + return this->getEmbeddedSurfaceBlocks().registerGroup< EmbeddedSurfaceBlock >( name ); } diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp index 0b86848bd34..98d7e3c9a50 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp @@ -17,59 +17,73 @@ #include "EmbeddedSurfaceBlock.hpp" namespace geos -{ +{ -localIndex EmbeddedSurfaceBlock::numEmbeddedSurfElem() const{ - return m_numEmbeddedSurfaces; +localIndex EmbeddedSurfaceBlock::numEmbeddedSurfElem() const +{ + return m_numEmbeddedSurfaces; } -ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemToNodes() const { - return m_embeddedSurfElemToNodes; +ArrayOfArrays< localIndex > EmbeddedSurfaceBlock::getEmbeddedSurfElemToNodes() const +{ + return m_embeddedSurfElemToNodes; } -ToCellRelation> EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const { - return m_embeddedSurfElemTo3dElem; +ToCellRelation< ArrayOfArrays< localIndex > > EmbeddedSurfaceBlock::getEmbeddedSurfElemTo3dElem() const +{ + return m_embeddedSurfElemTo3dElem; } -ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNodesCoords() const { - return m_embeddedSurfElemNodesCoords; +ArrayOfArrays< real64 > EmbeddedSurfaceBlock::getEmbeddedSurfElemNodesCoords() const +{ + return m_embeddedSurfElemNodesCoords; } -ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemNormalVectors() const { - return m_embeddedSurfElemNormals; +ArrayOfArrays< real64 > EmbeddedSurfaceBlock::getEmbeddedSurfElemNormalVectors() const +{ + return m_embeddedSurfElemNormals; } -ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemTangentialLengthVectors() const{ - return m_embeddedSurfElemLengthVectors; +ArrayOfArrays< real64 > EmbeddedSurfaceBlock::getEmbeddedSurfElemTangentialLengthVectors() const +{ + return m_embeddedSurfElemLengthVectors; } -ArrayOfArrays EmbeddedSurfaceBlock::getEmbeddedSurfElemTangentialWidthVectors() const{ - return m_embeddedSurfElemWidthVectors; +ArrayOfArrays< real64 > EmbeddedSurfaceBlock::getEmbeddedSurfElemTangentialWidthVectors() const +{ + return m_embeddedSurfElemWidthVectors; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemNormalVectors(ArrayOfArrays&& _normals){ - m_embeddedSurfElemNormals = _normals; +void EmbeddedSurfaceBlock::setEmbeddedSurfElemNormalVectors( ArrayOfArrays< real64 > && _normals ) +{ + m_embeddedSurfElemNormals = _normals; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialLengthVectors(ArrayOfArrays && _lengthVectors){ - m_embeddedSurfElemLengthVectors = _lengthVectors; +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialLengthVectors( ArrayOfArrays< real64 > && _lengthVectors ) +{ + m_embeddedSurfElemLengthVectors = _lengthVectors; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialWidthVectors(ArrayOfArrays && _widthVectors){ - m_embeddedSurfElemWidthVectors= _widthVectors; +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialWidthVectors( ArrayOfArrays< real64 > && _widthVectors ) +{ + m_embeddedSurfElemWidthVectors= _widthVectors; } -void EmbeddedSurfaceBlock::setNumEmbeddedSurfElem(localIndex _numEmbeddedSurfaces){ - - m_numEmbeddedSurfaces = _numEmbeddedSurfaces; +void EmbeddedSurfaceBlock::setNumEmbeddedSurfElem( localIndex _numEmbeddedSurfaces ) +{ + + m_numEmbeddedSurfaces = _numEmbeddedSurfaces; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemToNodes(ArrayOfArrays && _embeddedSurfElemToNodes){ - m_embeddedSurfElemToNodes = _embeddedSurfElemToNodes; +void EmbeddedSurfaceBlock::setEmbeddedSurfElemToNodes( ArrayOfArrays< localIndex > && _embeddedSurfElemToNodes ) +{ + m_embeddedSurfElemToNodes = _embeddedSurfElemToNodes; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem(ToCellRelation< ArrayOfArrays> && _embeddedSurfElemTo3dElem){ - m_embeddedSurfElemTo3dElem = _embeddedSurfElemTo3dElem; +void EmbeddedSurfaceBlock::setEmbeddedSurfElemTo3dElem( ToCellRelation< ArrayOfArrays< localIndex > > && _embeddedSurfElemTo3dElem ) +{ + m_embeddedSurfElemTo3dElem = _embeddedSurfElemTo3dElem; } -void EmbeddedSurfaceBlock::setEmbeddedSurfElemNodes(ArrayOfArrays && _embeddedSurfElemNodes){ - m_embeddedSurfElemNodesCoords = _embeddedSurfElemNodes; +void EmbeddedSurfaceBlock::setEmbeddedSurfElemNodes( ArrayOfArrays< real64 > && _embeddedSurfElemNodes ) +{ + m_embeddedSurfElemNodesCoords = _embeddedSurfElemNodes; } -} \ No newline at end of file +} diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp index 44f6b2bf632..31de2d1abc8 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -19,13 +19,13 @@ namespace geos { - + /** * @brief Simple implementation of the @p EmbeddedSurfaceBlockABC contract - * -*/ + * + */ -class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC +class EmbeddedSurfaceBlock : public EmbeddedSurfaceBlockABC { public: /** @@ -34,73 +34,73 @@ class EmbeddedSurfaceBlock: public EmbeddedSurfaceBlockABC * @param[in] parent Parent group. */ EmbeddedSurfaceBlock( string const & name, - Group * const parent ) + Group * const parent ) : EmbeddedSurfaceBlockABC( name, parent ) { } - localIndex numEmbeddedSurfElem() const override; - ArrayOfArrays getEmbeddedSurfElemToNodes() const override; - ToCellRelation> getEmbeddedSurfElemTo3dElem() const override; - ArrayOfArrays getEmbeddedSurfElemNodesCoords() const override; - ArrayOfArrays getEmbeddedSurfElemNormalVectors() const override; - ArrayOfArrays getEmbeddedSurfElemTangentialLengthVectors() const override; - ArrayOfArrays 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 && _embeddedSurfElemToNodes); - - /** - * @brief Sets the embedded elements to 1 3d element mapping - * @param _embeddedSurfElemTo3dElem the input mapping. - */ - void setEmbeddedSurfElemTo3dElem(ToCellRelation> && _embeddedSurfElemTo3dElem); - - /** - * @brief Sets the embedded elements nodes coordinates - * @param _embeddedSurfElemNodes the coordinates of embedded elements nodes. - */ - void setEmbeddedSurfElemNodes(ArrayOfArrays && _embeddedSurfElemNodes); - - /** - * @brief Sets the embedded elements normals - * @param _normals the coordinates of embedded elements normals. - */ - void setEmbeddedSurfElemNormalVectors(ArrayOfArrays && _normals); - - /** - * @brief Sets the embedded elements tangential length vectors (horizontal plan) - * @param _lengthVectors the coordinates of embedded elements tangential length vectors. - */ - void setEmbeddedSurfElemTangentialLengthVectors(ArrayOfArrays && _lengthVectors); - - /** - * @brief Sets the embedded elements tangential width vectors (vertical plan) - * @param _widthVectors the coordinates of embedded elements tangential width vectors. - */ - void setEmbeddedSurfElemTangentialWidthVectors(ArrayOfArrays && _widthVectors); + 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 ); private: - - localIndex m_numEmbeddedSurfaces; - ArrayOfArrays< localIndex > m_embeddedSurfElemToNodes; - ToCellRelation> m_embeddedSurfElemTo3dElem; - ArrayOfArrays< real64 > m_embeddedSurfElemNodesCoords; - ArrayOfArrays< real64 > m_embeddedSurfElemNormals; - ArrayOfArrays< real64 > m_embeddedSurfElemLengthVectors; - ArrayOfArrays< real64 > m_embeddedSurfElemWidthVectors; + + 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; }; } -#endif //inlcude guard \ No newline at end of file +#endif //inlcude guard diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp index 57da75c2564..7ced1c7c089 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlockABC.hpp @@ -26,99 +26,99 @@ 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 + * + * @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. + * 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 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> 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 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 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 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 getEmbeddedSurfElemTangentialWidthVectors() const = 0; + /** + * @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 \ No newline at end of file +#endif // include guard diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index 6cf5d925b10..bd65bf38eca 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -130,40 +130,41 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() localIndex localNumberOfSurfaceElems = 0; NewObjectLists newObjects; - + bool bulk_edfm_fill = true; - if (bulk_edfm_fill){ - + if( bulk_edfm_fill ) + { + CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager(); - CellBlockManagerABC & cellBlockManagerNoConst = const_cast(domain.getMeshBody(0).getCellBlockManager()); - CellBlockManager & cellBlockManagerConcrete = dynamic_cast(cellBlockManagerNoConst); + CellBlockManagerABC & cellBlockManagerNoConst = const_cast< CellBlockManagerABC & >(domain.getMeshBody( 0 ).getCellBlockManager()); + CellBlockManager & cellBlockManagerConcrete = dynamic_cast< CellBlockManager & >(cellBlockManagerNoConst); Group const & embSurfBlocks = cellBlockManagerConcrete.getEmbeddedSurfaceBlocks(); - if (embSurfBlocks.hasGroup("EmbeddedSurface")) - { - EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup("EmbeddedSurface"); - - - elemManager.forElementSubRegionsComplete< CellElementSubRegion >( - [&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion ) + if( embSurfBlocks.hasGroup( "EmbeddedSurface" )) { - arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList(); - FixedOneToManyRelation const & cellToEdges = subRegion.edgeList(); + EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( "EmbeddedSurface" ); - arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); - { - bool added = embeddedSurfaceSubRegion.addAllEmbeddedSurfaces( er, - esr, - nodeManager, - embSurfNodeManager, - edgeManager, - cellToEdges, - embSurf); - if( added ) + elemManager.forElementSubRegionsComplete< CellElementSubRegion >( + [&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion ) + { + arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList(); + FixedOneToManyRelation const & cellToEdges = subRegion.edgeList(); + + arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + { + bool added = embeddedSurfaceSubRegion.addAllEmbeddedSurfaces( er, + esr, + nodeManager, + embSurfNodeManager, + edgeManager, + cellToEdges, + embSurf ); + + if( added ) + { + // Add all the fracture information to the CellElementSubRegion + for( localIndex edfmIndex=0; edfmIndex < embSurf.numEmbeddedSurfElem(); ++edfmIndex ) { - // 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 ); @@ -250,9 +251,91 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() } } } - } );// end loop over cells - } );// end loop over subregions - } );// end loop over planes + } );// 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 ) + { + 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 ) + { + 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 ) + { + isPositive = 1; + } + else if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) < 0 ) + { + 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_RANK_0( 2, "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 } // Launch kernel to compute connectivity index of each fractured element. elemManager.forElementSubRegionsComplete< CellElementSubRegion >( From cd555e0af92708352d89bf5f12f425505c652022 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Tue, 4 Jun 2024 15:59:36 -0500 Subject: [PATCH 17/32] start VTK edfm support --- .../mesh/ElementRegionManager.cpp | 2 +- .../VTKEmbeddedSurfaceBlockUtilities.cpp | 45 +++++++++++++++++++ .../VTKEmbeddedSurfaceBlockUtilities.hpp | 41 +++++++++++++++++ 3 files changed, 87 insertions(+), 1 deletion(-) create mode 100644 src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp create mode 100644 src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index e32090eec47..5afc690844f 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -141,7 +141,7 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa // testing only // - bool experiment = true; + bool experiment = false; if( experiment ) { diff --git a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp new file mode 100644 index 00000000000..87b78cec081 --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp @@ -0,0 +1,45 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "CellBlockManager.hpp" + +namespace geos::vtk +{ + +void importFractureNetwork( string const & embeddedSurfaceBlockName, + vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, + vtkSmartPointer< vtkDataSet > mesh, + CellBlockManager & cellBlockManager ){ + +// prepare the data +vtkIdType const numEdfmSurf = embeddedSurfaceMesh->GetNumberOfCells(); + + +EmbeddedSurfaceBlock & embeddedSurfBlock = cellBlockManager.registerEmbeddedSurfaceBlock(embeddedSurfaceBlockName); + + +embeddedSurfBlock.setNumEmbeddedSurfElem(2828); +embeddedSurfBlock.setEmbeddedSurfElemNodes(2828); +embeddedSurfBlock.setEmbeddedSurfElemNormalVectors(2828); +embeddedSurfBlock.setEmbeddedSurfElemTangentialLengthVectors(2828); +embeddedSurfBlock.setEmbeddedSurfElemTangentialWidthVectors(2828); +embeddedSurfBlock.setEmbeddedSurfElemTo3dElem(2828); + + + + } +} + diff --git a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp new file mode 100644 index 00000000000..9dcc0a39010 --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp @@ -0,0 +1,41 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 mesh[in] The vtk volumic mesh. + * @param cellBlockManager[inout] The cell block manager that will receive the embedded surface block information. + */ +void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName, + vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, + vtkSmartPointer< vtkDataSet > mesh, + CellBlockManager & cellBlockManager ); +} + +#endif // include guard From 041d2905ffec0ab3f24a9389e8ecdf2807f5d18a Mon Sep 17 00:00:00 2001 From: Ouassim Date: Mon, 24 Jun 2024 15:51:52 -0500 Subject: [PATCH 18/32] Finish EDFM loading from VTK (no parallel) --- src/coreComponents/mesh/CMakeLists.txt | 2 + .../mesh/ElementRegionManager.cpp | 23 ++-- .../mesh/SurfaceElementRegion.cpp | 31 ++--- .../mesh/generators/EmbeddedSurfaceBlock.cpp | 10 ++ .../mesh/generators/EmbeddedSurfaceBlock.hpp | 14 +++ .../VTKEmbeddedSurfaceBlockUtilities.cpp | 116 +++++++++++++++--- .../mesh/generators/VTKMeshGenerator.cpp | 19 ++- .../mesh/generators/VTKMeshGenerator.hpp | 8 ++ .../mesh/generators/VTKUtilities.cpp | 35 +++++- .../mesh/generators/VTKUtilities.hpp | 31 ++++- .../mesh/generators/VTKWellGenerator.cpp | 2 +- src/coreComponents/schema/docs/VTKMesh.rst | 27 ++++ src/coreComponents/schema/schema.xsd | 2 + 13 files changed, 259 insertions(+), 61 deletions(-) create mode 100644 src/coreComponents/schema/docs/VTKMesh.rst diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt index 6375ef6fcd5..0f6f3822e0f 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -160,6 +160,7 @@ if( ENABLE_VTK ) set( mesh_headers ${mesh_headers} generators/CollocatedNodes.hpp generators/VTKFaceBlockUtilities.hpp + generators/VTKEmbeddedSurfaceBlockUtilities.hpp generators/VTKMeshGenerator.hpp generators/VTKMeshGeneratorTools.hpp generators/VTKWellGenerator.hpp @@ -168,6 +169,7 @@ if( ENABLE_VTK ) set( mesh_sources ${mesh_sources} generators/CollocatedNodes.cpp generators/VTKFaceBlockUtilities.cpp + generators/VTKEmbeddedSurfaceBlockUtilities.cpp generators/VTKMeshGenerator.cpp generators/VTKMeshGeneratorTools.cpp generators/VTKWellGenerator.cpp diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index 5afc690844f..1376f1672b2 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -138,24 +138,17 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa this->forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion & elemRegion ) { - // testing only - // - - bool experiment = false; - if( experiment ) + + if( elemRegion.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::faceElement ) + { + Group const * fracBlocks = &cellBlockManager.getFaceBlocks(); + elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); + } + else if( elemRegion.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::embeddedElement ) { - - CellBlockManagerABC & cellBlockManagerNoConst = const_cast< CellBlockManagerABC & >(cellBlockManager); - CellBlockManager & cellBlockManagerConcrete = dynamic_cast< CellBlockManager & >(cellBlockManagerNoConst); - - string name = "EmbeddedSurface"; - createDummyEmbeddedSurfaceBlockOrig( name, cellBlockManagerConcrete ); + Group const * edfmBlocks = &cellBlockManager.getEmbeddedSurfaceBlocks(); elemRegion.generateMesh( cellBlockManager.getEmbeddedSurfaceBlocks() ); } - else - { - elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); - } } ); diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index 5f00d2bca1c..a8362d24dde 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -56,27 +56,16 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) { - // bool experiment =false; - // if (experiment) - // { - - // EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); - // if (faceBlocks.hasGroup("EmbeddedSurface")) - // { - // EmbeddedSurfaceBlockABC const &source = faceBlocks.getGroup("EmbeddedSurface"); - // subRegion.copyFromCellBlock(source); - // } - // else - // { - // GEOS_LOG_RANK_0("No face block \"" << m_faceBlockName << "\" was found in the mesh. Empty surface region was created."); - // } - // } - // else - // { - // elementSubRegions.registerGroup(m_faceBlockName); - // } - - elementSubRegions.registerGroup< EmbeddedSurfaceSubRegion >( "EmbeddedSurface" ); + EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); + if (faceBlocks.hasGroup(m_faceBlockName)) + { + EmbeddedSurfaceBlockABC const &source = faceBlocks.getGroup(m_faceBlockName); + subRegion.copyFromCellBlock(source); + } + else + { + GEOS_LOG_RANK_0("No face block \"" << m_faceBlockName << "\" was found in the mesh. Empty embedded surface region was created."); + } } else if( m_subRegionType == SurfaceSubRegionType::faceElement ) { diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp index 98d7e3c9a50..900615f78ac 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.cpp @@ -65,6 +65,16 @@ void EmbeddedSurfaceBlock::setEmbeddedSurfElemTangentialWidthVectors( ArrayOfArr 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 ) { diff --git a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp index 31de2d1abc8..4b4abff3c6c 100644 --- a/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp +++ b/src/coreComponents/mesh/generators/EmbeddedSurfaceBlock.hpp @@ -90,6 +90,18 @@ class EmbeddedSurfaceBlock : public EmbeddedSurfaceBlockABC */ 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 ); + private: localIndex m_numEmbeddedSurfaces; @@ -99,6 +111,8 @@ class EmbeddedSurfaceBlock : public EmbeddedSurfaceBlockABC ArrayOfArrays< real64 > m_embeddedSurfElemNormals; ArrayOfArrays< real64 > m_embeddedSurfElemLengthVectors; ArrayOfArrays< real64 > m_embeddedSurfElemWidthVectors; + array1d< real64 > m_embeddedSurfElemApertures; + array1d< real64 > m_embeddedSurfElemPermeability; }; } diff --git a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp index 87b78cec081..7317d9181c4 100644 --- a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp @@ -14,32 +14,120 @@ #include "VTKEmbeddedSurfaceBlockUtilities.hpp" +#include +#include +#include #include "CellBlockManager.hpp" namespace geos::vtk { -void importFractureNetwork( string const & embeddedSurfaceBlockName, +void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName, vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, vtkSmartPointer< vtkDataSet > mesh, CellBlockManager & cellBlockManager ){ -// prepare the data -vtkIdType const numEdfmSurf = embeddedSurfaceMesh->GetNumberOfCells(); + // Ouassim change + vtkIdType numEdfmFracs = embeddedSurfaceMesh->GetNumberOfCells(); + //1.Get EDFM vertices + vtkUnstructuredGrid * grid = vtkUnstructuredGrid::SafeDownCast(embeddedSurfaceMesh); + vtkPoints* const nodes = grid->GetPoints(); + vtkIdType numNodes = nodes->GetNumberOfPoints(); + ArrayOfArrays< localIndex> elem2dToNodes( numEdfmFracs ); + for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ + auto 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 * normals =grid->GetCellData()->GetArray(normal_str); + for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ + auto 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 * tangential_len =grid->GetCellData()->GetArray(tl_str); + for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ + auto 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 * tangential_width =grid->GetCellData()->GetArray(tw_str); + for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ + auto 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 * apertures =grid->GetCellData()->GetArray(aperture_str); + for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ + auto val =apertures->GetTuple1(i); + embeddedSurfaceAperture[i] = val; + } -EmbeddedSurfaceBlock & embeddedSurfBlock = cellBlockManager.registerEmbeddedSurfaceBlock(embeddedSurfaceBlockName); + //5. Get EDFM permeability + array1d< real64 > embeddedSurfacePermeability( numEdfmFracs ); + const char * const perm_str = "permeability"; + vtkDataArray * 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 * fid_mid =grid->GetCellData()->GetArray(fid_to_mid_mapping); + + ArrayOfArrays< localIndex > toBlockIndex( numEdfmFracs ); + ArrayOfArrays< localIndex > toCellIndex( numEdfmFracs ); + for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ + auto fidmid =fid_mid->GetTuple2(i); + toBlockIndex.emplaceBack(fidmid[0], 0 );// cell block is set to 0 for now + toCellIndex.emplaceBack( fidmid[0], fidmid[1]); + } + ToCellRelation< ArrayOfArrays< localIndex > > elem2dTo3d( std::move(toBlockIndex), std::move(toCellIndex) ); +EmbeddedSurfaceBlock & embeddedSurfBlock = cellBlockManager.registerEmbeddedSurfaceBlock(embeddedSurfaceBlockName); -embeddedSurfBlock.setNumEmbeddedSurfElem(2828); -embeddedSurfBlock.setEmbeddedSurfElemNodes(2828); -embeddedSurfBlock.setEmbeddedSurfElemNormalVectors(2828); -embeddedSurfBlock.setEmbeddedSurfElemTangentialLengthVectors(2828); -embeddedSurfBlock.setEmbeddedSurfElemTangentialWidthVectors(2828); -embeddedSurfBlock.setEmbeddedSurfElemTo3dElem(2828); - - - - } +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/VTKMeshGenerator.cpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp index d1402d21821..18514013dd5 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp @@ -20,6 +20,7 @@ #include "VTKMeshGenerator.hpp" #include "mesh/generators/VTKFaceBlockUtilities.hpp" +#include "mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp" #include "mesh/generators/VTKMeshGeneratorTools.hpp" #include "mesh/generators/CellBlockManager.hpp" #include "common/DataTypes.hpp" @@ -57,6 +58,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 ). @@ -89,12 +95,13 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager GEOS_LOG_RANK_0( GEOS_FMT( "{} '{}': reading mesh from {}", catalogName(), getName(), m_filePath ) ); { GEOS_LOG_LEVEL_RANK_0( 2, " reading the dataset..." ); - vtk::AllMeshes allMeshes = vtk::loadAllMeshes( m_filePath, m_mainBlockName, m_faceBlockNames ); + vtk::AllMeshes allMeshes = vtk::loadAllMeshes( m_filePath, m_mainBlockName, m_faceBlockNames, m_embeddedSurfaceBlockNames); GEOS_LOG_LEVEL_RANK_0( 2, " redistributing mesh..." ); 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, " finding neighbor ranks..." ); std::vector< vtkBoundingBox > boxes = vtk::exchangeBoundingBoxes( *m_vtkMesh, comm ); std::vector< int > const neighbors = vtk::findNeighborRanks( std::move( boxes ) ); @@ -121,7 +128,15 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager for( auto const & [name, mesh]: m_faceBlockMeshes ) { + //TODO: Ouassim change for experimentation vtk::importFractureNetwork( name, mesh, m_vtkMesh, cellBlockManager ); + + } + for( auto const & [name, mesh]: m_embeddedSurfaceBlockMeshes ) + { + //TODO: Ouassim change for experimentation + vtk::importEmbeddedFractureNetwork( name, mesh, m_vtkMesh, cellBlockManager ); + } GEOS_LOG_LEVEL_RANK_0( 2, " done!" ); diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp index 6d2eabe911f..51bea1ab1e7 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp @@ -106,6 +106,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"; } @@ -138,8 +139,15 @@ 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 72045afa4e5..fcc84ed887b 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 ); } @@ -687,8 +693,10 @@ AllMeshes redistributeByCellGraph( AllMeshes & input, vtkSmartPointer< vtkUnstructuredGrid > const finalFracMesh = vtk::redistribute( *splitFracMesh, MPI_COMM_GEOS ); 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 ); } /** @@ -900,6 +908,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,6 +922,12 @@ 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 ) ); @@ -922,6 +937,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 +964,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 +983,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, ", " ) ); diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp index 6ebf4e6e89d..db08081a39d 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.hpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.hpp @@ -78,9 +78,11 @@ class AllMeshes * @param faceBlocks The fractures meshes. */ AllMeshes( vtkSmartPointer< vtkDataSet > const & main, - std::map< string, vtkSmartPointer< vtkDataSet > > const & faceBlocks ) + std::map< string, vtkSmartPointer< vtkDataSet > > const & faceBlocks, + std::map> const & edfmSurfBlocks ) : m_main( main ), - m_faceBlocks( faceBlocks ) + m_faceBlocks( faceBlocks ), + m_embeddedSurfaceBlocks(edfmSurfBlocks) { } /** @@ -99,6 +101,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 +127,23 @@ class AllMeshes m_faceBlocks = faceBlocks; } + /** + * @brief Defines the face blocks/fractures. + * @param edfmSurBlocks 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 +151,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 const & edfmSurfBlockNames); /** * @brief Compute the rank neighbor candidate list. @@ -149,6 +172,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 +183,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 c6294c35683..03f5d6f657d 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/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 6c137c9f37f..a8d7ddfd229 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -1697,6 +1697,8 @@ stress - traction is applied to the faces as specified by the inner product of i + + From 3bda862b8cc3839cdba39ab2301bf2096e404633 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Thu, 31 Oct 2024 16:26:27 -0500 Subject: [PATCH 19/32] builds after master rebase issues --- .../mesh/EmbeddedSurfaceSubRegion.cpp | 56 ++++++------ .../mesh/generators/CellBlockManagerABC.hpp | 2 +- .../EmbeddedSurfaceGenerator.cpp | 91 ++----------------- 3 files changed, 35 insertions(+), 114 deletions(-) diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index c2b906b7463..560855bd8e1 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -418,12 +418,12 @@ bool EmbeddedSurfaceSubRegion::copyFromCellBlock( EmbeddedSurfaceBlockABC const { - DomainPartition & domain = getGlobalState().getProblemManager().getDomainPartition(); - MeshLevel & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); - EmbeddedSurfaceNodeManager & embSurfNodeManager = meshLevel.getEmbSurfNodeManager(); - EdgeManager & edgeManager = meshLevel.getEdgeManager(); + // DomainPartition & domain = getGlobalState().getProblemManager().getDomainPartition(); + // MeshLevel & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); + // EmbeddedSurfaceNodeManager & embSurfNodeManager = meshLevel.getEmbSurfNodeManager(); + // EdgeManager & edgeManager = meshLevel.getEdgeManager(); - ElementRegionManager & elemManager = meshLevel.getElemManager(); + // ElementRegionManager & elemManager = meshLevel.getElemManager(); localIndex const numElems= embeddedSurfaceBlock.numEmbeddedSurfElem(); @@ -465,40 +465,40 @@ bool EmbeddedSurfaceSubRegion::copyFromCellBlock( EmbeddedSurfaceBlockABC const } // add nodes to embNodeManager - for( auto i =0; i < elemNodesLocations.size(); ++i ) - { + /// for( auto i =0; i < elemNodesLocations.size(); ++i ) + /// { - array2d< real64 > elemNodeCoords( 1, 3 ); - auto ghostRank = -2; + /// array2d< real64 > elemNodeCoords( 1, 3 ); + /// auto ghostRank = -2; - // // Add the point to the node Manager - // globalIndex parentEdgeID = edgeLocalToGlobal[ pointParentIndex[ originalIndices[ j ] ] ]; - // nodeIndex = embSurfNodeManager.size(); - embSurfNodeManager.appendNode( elemNodesLocations[i].toSliceConst(), ghostRank ); + /// // // Add the point to the node Manager + /// // globalIndex parentEdgeID = edgeLocalToGlobal[ pointParentIndex[ originalIndices[ j ] ] ]; + /// // nodeIndex = embSurfNodeManager.size(); + /// embSurfNodeManager.appendNode( elemNodesLocations[i].toSliceConst(), ghostRank ); - // arrayView1d< localIndex > const & parentIndex = - // embSurfNodeManager.getField< fields::parentEdgeIndex >(); + /// // arrayView1d< localIndex > const & parentIndex = + /// // embSurfNodeManager.getField< fields::parentEdgeIndex >(); - // parentIndex[nodeIndex] = pointParentIndex[ originalIndices[ j ] ]; + /// // parentIndex[nodeIndex] = pointParentIndex[ originalIndices[ j ] ]; - // array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); - // parentEdgeGlobalIndex[nodeIndex] = parentEdgeID; - } + /// // array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); + /// // parentEdgeGlobalIndex[nodeIndex] = parentEdgeID; + /// } - elemManager.forElementSubRegionsComplete< CellElementSubRegion >( - [&] ( localIndex, localIndex, ElementRegionBase const & region, CellElementSubRegion & subRegion ) - { + // elemManager.forElementSubRegionsComplete< CellElementSubRegion >( + // [&] ( localIndex, localIndex, ElementRegionBase const & region, CellElementSubRegion & subRegion ) + // { - for( auto i = 0; i(domain.getMeshBody( 0 ).getCellBlockManager()); - CellBlockManager & cellBlockManagerConcrete = dynamic_cast< CellBlockManager & >(cellBlockManagerNoConst); - Group const & embSurfBlocks = cellBlockManagerConcrete.getEmbeddedSurfaceBlocks(); + // CellBlockManagerABC & cellBlockManagerNoConst = const_cast< CellBlockManagerABC & >(domain.getMeshBody( 0 ).getCellBlockManager()); + // CellBlockManager & cellBlockManagerConcrete = dynamic_cast< CellBlockManager & >(cellBlockManagerNoConst); + //Group const & embSurfBlocks = cellBlockManagerConcrete.getEmbeddedSurfaceBlocks(); + Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks(); if( embSurfBlocks.hasGroup( "EmbeddedSurface" )) { EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( "EmbeddedSurface" ); @@ -253,89 +254,9 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() } } );// 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 ) - { - 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 ) - { - 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 ) - { - isPositive = 1; - } - else if( LvArray::tensorOps::AiBi< 3 >( distVec, normalVector ) < 0 ) - { - 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_RANK_0( 2, "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 + ); } // Launch kernel to compute connectivity index of each fractured element. elemManager.forElementSubRegionsComplete< CellElementSubRegion >( From 256efee6fac389a9e77ac3012972b2a72cc1f9f5 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Tue, 5 Nov 2024 10:16:24 -0600 Subject: [PATCH 20/32] working serial run --- .../surfaceGeneration/EmbeddedSurfaceGenerator.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index 03a906f39c5..d49892e4931 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -140,9 +140,10 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() // CellBlockManager & cellBlockManagerConcrete = dynamic_cast< CellBlockManager & >(cellBlockManagerNoConst); //Group const & embSurfBlocks = cellBlockManagerConcrete.getEmbeddedSurfaceBlocks(); Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks(); - if( embSurfBlocks.hasGroup( "EmbeddedSurface" )) + //if( embSurfBlocks.hasGroup( "EmbeddedSurface" )) + if( embSurfBlocks.hasGroup( "fracture" )) { - EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( "EmbeddedSurface" ); + EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( "fracture" ); elemManager.forElementSubRegionsComplete< CellElementSubRegion >( From 2d6583d7a81bb2ea5e0bd15d83b67dd238d42a0f Mon Sep 17 00:00:00 2001 From: Ouassim Date: Thu, 7 Nov 2024 16:25:37 -0600 Subject: [PATCH 21/32] clearning --- .../mainInterface/ProblemManager.cpp | 2 + src/coreComponents/mesh/CMakeLists.txt | 4 - .../mesh/ElementRegionManager.cpp | 2 - .../mesh/EmbeddedSurfaceSubRegion.cpp | 99 +------------------ .../mesh/EmbeddedSurfaceSubRegion.hpp | 31 ++++-- .../mesh/SurfaceElementRegion.cpp | 20 ++-- .../mesh/generators/CellBlockManager.cpp | 11 --- .../mesh/generators/CellBlockManager.hpp | 7 -- .../mesh/generators/CellBlockManagerABC.hpp | 12 --- 9 files changed, 31 insertions(+), 157 deletions(-) diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp index a4bf4f0dc37..bbf09842795 100644 --- a/src/coreComponents/mainInterface/ProblemManager.cpp +++ b/src/coreComponents/mainInterface/ProblemManager.cpp @@ -691,6 +691,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 0f6f3822e0f..b758dbf2b9c 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -210,7 +210,3 @@ endif() if( GEOS_ENABLE_TESTS ) add_subdirectory( unitTests ) endif( ) - - -#add_dependencies(mesh mainInterface) -#add_dependencies(mesh linearAlgebra) \ No newline at end of file diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index 1376f1672b2..e7d7b2ac117 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -141,12 +141,10 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa if( elemRegion.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::faceElement ) { - Group const * fracBlocks = &cellBlockManager.getFaceBlocks(); elemRegion.generateMesh( cellBlockManager.getFaceBlocks() ); } else if( elemRegion.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::embeddedElement ) { - Group const * edfmBlocks = &cellBlockManager.getEmbeddedSurfaceBlocks(); elemRegion.generateMesh( cellBlockManager.getEmbeddedSurfaceBlocks() ); } diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index 560855bd8e1..25b37432e5c 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -29,13 +29,8 @@ #include "mainInterface/GeosxState.hpp" -//#include "common/DataTypes.hpp" -//#include "linearAlgebra/DofManager.hpp" -//#include "mainInterface/initialization.hpp" #include "mainInterface/ProblemManager.hpp" #include "mesh/DomainPartition.hpp" -//#include "mesh/MeshManager.hpp" -//#include "mesh/mpiCommunications/CommunicationTools.hpp" namespace geos { @@ -278,7 +273,7 @@ array1d< localIndex > EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex( Arra 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 oly need one) + // 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 > 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(); - - - - for( auto i = 0; i 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 ); - } - - // add nodes to embNodeManager - /// for( auto i =0; i < elemNodesLocations.size(); ++i ) - /// { - - /// array2d< real64 > elemNodeCoords( 1, 3 ); - /// auto ghostRank = -2; - - /// // // Add the point to the node Manager - /// // globalIndex parentEdgeID = edgeLocalToGlobal[ pointParentIndex[ originalIndices[ j ] ] ]; - /// // nodeIndex = embSurfNodeManager.size(); - /// embSurfNodeManager.appendNode( elemNodesLocations[i].toSliceConst(), ghostRank ); - - /// // arrayView1d< localIndex > const & parentIndex = - /// // embSurfNodeManager.getField< fields::parentEdgeIndex >(); - - /// // parentIndex[nodeIndex] = pointParentIndex[ originalIndices[ j ] ]; - - /// // array1d< globalIndex > & parentEdgeGlobalIndex = embSurfNodeManager.getParentEdgeGlobalIndex(); - /// // parentEdgeGlobalIndex[nodeIndex] = parentEdgeID; - /// } - - - // elemManager.forElementSubRegionsComplete< CellElementSubRegion >( - // [&] ( localIndex, localIndex, ElementRegionBase const & region, CellElementSubRegion & subRegion ) - // { - - // for( auto i = 0; i > > const & cellGhostRank ) { arrayView1d< integer > const & ghostRank = this->ghostRank(); diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp index 7599210d76d..e72ef6a91a6 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.hpp @@ -30,10 +30,6 @@ //Do we really need this include Rectangle? #include "simpleGeometricObjects/Rectangle.hpp" -//#include "mainInterface/GeosxState.hpp" -//#include "mainInterface/ProblemManager.hpp" -//#include "mesh/DomainPartition.hpp" - namespace geos { @@ -151,6 +147,16 @@ 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, @@ -158,6 +164,18 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion 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, @@ -166,11 +184,6 @@ class EmbeddedSurfaceSubRegion : public SurfaceElementSubRegion FixedOneToManyRelation const & cellToEdges, EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock ); - /** - * @brief Fill @p EmbeddedSurfaceSubRegion by copying the data from the source embedded surface block - * - * **/ - bool copyFromCellBlock( EmbeddedSurfaceBlockABC const & embeddedSurfaceBlock ); /** * @brief inherit ghost rank from cell elements. diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index a8362d24dde..25a3a409896 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -54,20 +54,12 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) Group & elementSubRegions = this->getGroup( viewKeyStruct::elementSubRegions() ); - if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) - { - EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); - if (faceBlocks.hasGroup(m_faceBlockName)) - { - EmbeddedSurfaceBlockABC const &source = faceBlocks.getGroup(m_faceBlockName); - subRegion.copyFromCellBlock(source); - } - else - { - GEOS_LOG_RANK_0("No face block \"" << m_faceBlockName << "\" was found in the mesh. Empty embedded surface region was created."); - } - } - else if( m_subRegionType == SurfaceSubRegionType::faceElement ) + if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) + { + // We just register the subregion copying of data is done at the EmbeddedSurfaceGenerator + EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); + } + if( m_subRegionType == SurfaceSubRegionType::faceElement ) { FaceElementSubRegion & subRegion = elementSubRegions.registerGroup< FaceElementSubRegion >( m_faceBlockName ); if( faceBlocks.hasGroup( m_faceBlockName ) ) diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index ec193806f7c..11b896b6d6d 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.cpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp @@ -670,7 +670,6 @@ void CellBlockManager::buildMaps() fillElementToEdgesOfCellBlocks( m_faceToEdges.toViewConst(), this->getCellBlocks() ); - //buildEmbeddedSurfaceMaps(); } ArrayOfArrays< localIndex > CellBlockManager::getFaceToNodes() const @@ -690,12 +689,6 @@ ToCellRelation< array2d< localIndex > > CellBlockManager::getFaceToElements() co return m_faceToCells; } -ToCellRelation< localIndex > CellBlockManager::getEmbeddedSurfaceToElements() const -{ - - return m_embeddedSurfToCells; -} - const Group & CellBlockManager::getCellBlocks() const { return this->getGroup( viewKeyStruct::cellBlocks() ); @@ -755,10 +748,6 @@ localIndex CellBlockManager::numFaces() const return m_numFaces; } -localIndex CellBlockManager::numEmbeddedSurfaces() const -{ - return m_numEmbeddedSurfElem; -} ArrayOfArrays< localIndex > CellBlockManager::getEdgeToFaces() const { diff --git a/src/coreComponents/mesh/generators/CellBlockManager.hpp b/src/coreComponents/mesh/generators/CellBlockManager.hpp index 64a80422d26..c97cdee825e 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.hpp @@ -82,8 +82,6 @@ class CellBlockManager : public CellBlockManagerABC ToCellRelation< array2d< localIndex > > getFaceToElements() const override; - ToCellRelation< localIndex > getEmbeddedSurfaceToElements() const override; - array1d< globalIndex > getNodeLocalToGlobal() const override; /** @@ -128,8 +126,6 @@ class CellBlockManager : public CellBlockManagerABC localIndex numFaces() const override; - localIndex numEmbeddedSurfaces() const override; - using Group::resize; /** @@ -297,8 +293,6 @@ class CellBlockManager : public CellBlockManagerABC ArrayOfArrays< localIndex > m_faceToEdges; ToCellRelation< array2d< localIndex > > m_faceToCells; - ToCellRelation< localIndex > m_embeddedSurfToCells; - array1d< globalIndex > m_nodeLocalToGlobal; std::map< string, SortedArray< localIndex > > m_nodeSets; @@ -310,7 +304,6 @@ class CellBlockManager : public CellBlockManagerABC localIndex m_numNodes; localIndex m_numFaces; localIndex m_numEdges; - localIndex m_numEmbeddedSurfElem; }; } diff --git a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp index 1069e680cc2..c48cfdde452 100644 --- a/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp +++ b/src/coreComponents/mesh/generators/CellBlockManagerABC.hpp @@ -160,12 +160,6 @@ class CellBlockManagerABC : public dataRepository::Group */ virtual localIndex numFaces() const = 0; - /** - * @brief Total number of embedded surfaces across all the cell blocks. - * @return The total number of embedded surfaces. - */ - virtual localIndex numEmbeddedSurfaces() const = 0; - /** * @brief Returns the node coordinates in a (numNodes, 3) 2d array. * @return A const view to the array. @@ -221,12 +215,6 @@ class CellBlockManagerABC : public dataRepository::Group */ virtual ToCellRelation< array2d< localIndex > > getFaceToElements() const = 0; - /** - * @brief Returns the embedded surface to elements mapping. - * @return A 1 to 1 relationship. The result is mapping from 1 fracture element to its parent matrix element. - * - */ - virtual ToCellRelation< localIndex > getEmbeddedSurfaceToElements() const = 0; /** * @brief The node to global mapping for nodes. * @return The mapping as an array of size numNodes. From 8d52abcd44940b9d7d96a3b4b767c0a3c4434c5c Mon Sep 17 00:00:00 2001 From: Ouassim Date: Fri, 8 Nov 2024 14:01:22 -0600 Subject: [PATCH 22/32] cleraning --- src/coreComponents/mesh/SurfaceElementRegion.cpp | 2 +- src/coreComponents/mesh/SurfaceElementRegion.hpp | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index 25a3a409896..dced747acc5 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -59,7 +59,7 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) // We just register the subregion copying of data is done at the EmbeddedSurfaceGenerator EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); } - if( m_subRegionType == SurfaceSubRegionType::faceElement ) + else if( m_subRegionType == SurfaceSubRegionType::faceElement ) { FaceElementSubRegion & subRegion = elementSubRegions.registerGroup< FaceElementSubRegion >( m_faceBlockName ); if( faceBlocks.hasGroup( m_faceBlockName ) ) diff --git a/src/coreComponents/mesh/SurfaceElementRegion.hpp b/src/coreComponents/mesh/SurfaceElementRegion.hpp index e2b27838db9..182d4f730b0 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.hpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.hpp @@ -224,7 +224,6 @@ ENUM_STRINGS( SurfaceElementRegion::SurfaceSubRegionType, "faceElement", "embeddedElement" ); -EmbeddedSurfaceBlockABC & createDummyEmbeddedSurfaceBlock( string embeddedSurfFace ); } /* namespace geos */ #endif /* CORECOMPONENTS_MESH_SurfaceElementRegion_HPP_ */ From c718d96496237e8a1b6003f8e38ccc9608944937 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Fri, 8 Nov 2024 14:11:03 -0600 Subject: [PATCH 23/32] bug --- src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp index 25b37432e5c..a9c0832bcfd 100644 --- a/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp +++ b/src/coreComponents/mesh/EmbeddedSurfaceSubRegion.cpp @@ -298,7 +298,7 @@ array1d< localIndex > EmbeddedSurfaceSubRegion::getEdfmNodeParentEdgeIndex( Arra auto elementIndex =nodeToElem[i]; auto cell3dIndex = elemTo3dElem.toCellIndex[elementIndex][0]; double elemMinDistance = std::numeric_limits< double >::max(); - localIndex targetEdgeIndex; + localIndex targetEdgeIndex = -1; R1Tensor nodexyz; LvArray::tensorOps::copy< 3 >( nodexyz, elemNodesLocations[i] ); From ba87177c3d9f12ff41d57c9674e32f3f5530df4e Mon Sep 17 00:00:00 2001 From: Ouassim Date: Fri, 8 Nov 2024 14:49:36 -0600 Subject: [PATCH 24/32] uncrustify pass --- .../mesh/ElementRegionManager.cpp | 10 +- .../mesh/SurfaceElementRegion.cpp | 10 +- .../VTKEmbeddedSurfaceBlockUtilities.cpp | 148 +++++++++--------- .../VTKEmbeddedSurfaceBlockUtilities.hpp | 6 +- .../mesh/generators/VTKMeshGenerator.cpp | 5 +- .../mesh/generators/VTKMeshGenerator.hpp | 4 +- .../mesh/generators/VTKUtilities.cpp | 8 +- .../mesh/generators/VTKUtilities.hpp | 6 +- .../EmbeddedSurfaceGenerator.cpp | 140 +++++++++-------- 9 files changed, 174 insertions(+), 163 deletions(-) diff --git a/src/coreComponents/mesh/ElementRegionManager.cpp b/src/coreComponents/mesh/ElementRegionManager.cpp index e7d7b2ac117..03079b41faf 100644 --- a/src/coreComponents/mesh/ElementRegionManager.cpp +++ b/src/coreComponents/mesh/ElementRegionManager.cpp @@ -138,11 +138,11 @@ void ElementRegionManager::generateMesh( CellBlockManagerABC const & cellBlockMa this->forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion & elemRegion ) { - - if( elemRegion.subRegionType() == SurfaceElementRegion::SurfaceSubRegionType::faceElement ) - { - 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() ); diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index dced747acc5..f5918bb78f9 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -54,11 +54,11 @@ 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 - EmbeddedSurfaceSubRegion &subRegion = elementSubRegions.registerGroup(m_faceBlockName); - } + if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) + { + // We just register the subregion copying of data is done at the EmbeddedSurfaceGenerator + EmbeddedSurfaceSubRegion & subRegion = elementSubRegions.registerGroup< EmbeddedSurfaceSubRegion >( m_faceBlockName ); + } else if( m_subRegionType == SurfaceSubRegionType::faceElement ) { FaceElementSubRegion & subRegion = elementSubRegions.registerGroup< FaceElementSubRegion >( m_faceBlockName ); diff --git a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp index 7317d9181c4..a44c6c636a4 100644 --- a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp @@ -23,111 +23,119 @@ namespace geos::vtk { void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName, - vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, - vtkSmartPointer< vtkDataSet > mesh, - CellBlockManager & cellBlockManager ){ - + vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, + vtkSmartPointer< vtkDataSet > mesh, + CellBlockManager & cellBlockManager ) +{ + // Ouassim change vtkIdType numEdfmFracs = embeddedSurfaceMesh->GetNumberOfCells(); //1.Get EDFM vertices - vtkUnstructuredGrid * grid = vtkUnstructuredGrid::SafeDownCast(embeddedSurfaceMesh); - vtkPoints* const nodes = grid->GetPoints(); + vtkUnstructuredGrid * grid = vtkUnstructuredGrid::SafeDownCast( embeddedSurfaceMesh ); + vtkPoints * const nodes = grid->GetPoints(); vtkIdType numNodes = nodes->GetNumberOfPoints(); - ArrayOfArrays< localIndex> elem2dToNodes( numEdfmFracs ); - for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ - auto 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< localIndex > elem2dToNodes( numEdfmFracs ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto 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 ){ + 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]); + 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 ); + ArrayOfArrays< real64 > embeddedSurfaceNormalVectors( numEdfmFracs ); const char * const normal_str = "normal_vectors"; - vtkDataArray * normals =grid->GetCellData()->GetArray(normal_str); - for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ - auto val = normals->GetTuple3(i); - embeddedSurfaceNormalVectors.emplaceBack(i, val[0]); - embeddedSurfaceNormalVectors.emplaceBack(i, val[1]); - embeddedSurfaceNormalVectors.emplaceBack(i, val[2]); + vtkDataArray * normals =grid->GetCellData()->GetArray( normal_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto 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 * tangential_len =grid->GetCellData()->GetArray(tl_str); - for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ - auto val = tangential_len->GetTuple3(i); - embeddedSurfaceTangentialLengthVectors.emplaceBack(i, val[0]); - embeddedSurfaceTangentialLengthVectors.emplaceBack(i, val[1]); - embeddedSurfaceTangentialLengthVectors.emplaceBack(i, val[2]); + vtkDataArray * tangential_len =grid->GetCellData()->GetArray( tl_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto 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 * tangential_width =grid->GetCellData()->GetArray(tw_str); - for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ - auto val = tangential_width->GetTuple3(i); - embeddedSurfaceTangentialWidthVectors.emplaceBack(i, val[0]); - embeddedSurfaceTangentialWidthVectors.emplaceBack(i, val[1]); - embeddedSurfaceTangentialWidthVectors.emplaceBack(i, val[2]); + vtkDataArray * tangential_width =grid->GetCellData()->GetArray( tw_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto 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); + array1d< real64 > embeddedSurfaceAperture( numEdfmFracs ); const char * const aperture_str = "aperture"; - vtkDataArray * apertures =grid->GetCellData()->GetArray(aperture_str); - for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ - auto val =apertures->GetTuple1(i); + vtkDataArray * apertures =grid->GetCellData()->GetArray( aperture_str ); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto val =apertures->GetTuple1( i ); embeddedSurfaceAperture[i] = val; } //5. Get EDFM permeability array1d< real64 > embeddedSurfacePermeability( numEdfmFracs ); const char * const perm_str = "permeability"; - vtkDataArray * perms =grid->GetCellData()->GetArray(perm_str); - for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ - auto val =perms->GetTuple1(i); + vtkDataArray * 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 * fid_mid =grid->GetCellData()->GetArray(fid_to_mid_mapping); + vtkDataArray * fid_mid =grid->GetCellData()->GetArray( fid_to_mid_mapping ); ArrayOfArrays< localIndex > toBlockIndex( numEdfmFracs ); ArrayOfArrays< localIndex > toCellIndex( numEdfmFracs ); - for (vtkIdType i = 0; i < numEdfmFracs; ++i ){ - auto fidmid =fid_mid->GetTuple2(i); - toBlockIndex.emplaceBack(fidmid[0], 0 );// cell block is set to 0 for now - toCellIndex.emplaceBack( fidmid[0], fidmid[1]); + for( vtkIdType i = 0; i < numEdfmFracs; ++i ) + { + auto fidmid =fid_mid->GetTuple2( i ); + toBlockIndex.emplaceBack( fidmid[0], 0 );// cell block is set to 0 for now + toCellIndex.emplaceBack( fidmid[0], fidmid[1] ); } - ToCellRelation< ArrayOfArrays< localIndex > > elem2dTo3d( std::move(toBlockIndex), std::move(toCellIndex) ); - -EmbeddedSurfaceBlock & embeddedSurfBlock = cellBlockManager.registerEmbeddedSurfaceBlock(embeddedSurfaceBlockName); - -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)); + ToCellRelation< ArrayOfArrays< localIndex > > elem2dTo3d( std::move( toBlockIndex ), std::move( toCellIndex ) ); + + EmbeddedSurfaceBlock & embeddedSurfBlock = cellBlockManager.registerEmbeddedSurfaceBlock( embeddedSurfaceBlockName ); + + 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 index 9dcc0a39010..212835ade65 100644 --- a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.hpp @@ -33,9 +33,9 @@ namespace geos::vtk * @param cellBlockManager[inout] The cell block manager that will receive the embedded surface block information. */ void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName, - vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, - vtkSmartPointer< vtkDataSet > mesh, - CellBlockManager & cellBlockManager ); + vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, + vtkSmartPointer< vtkDataSet > mesh, + CellBlockManager & cellBlockManager ); } #endif // include guard diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp index 18514013dd5..4f8fd4a729d 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp @@ -95,10 +95,11 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager GEOS_LOG_RANK_0( GEOS_FMT( "{} '{}': reading mesh from {}", catalogName(), getName(), m_filePath ) ); { GEOS_LOG_LEVEL_RANK_0( 2, " reading the dataset..." ); - vtk::AllMeshes allMeshes = vtk::loadAllMeshes( m_filePath, m_mainBlockName, m_faceBlockNames, m_embeddedSurfaceBlockNames); + vtk::AllMeshes allMeshes = vtk::loadAllMeshes( m_filePath, m_mainBlockName, m_faceBlockNames, m_embeddedSurfaceBlockNames ); GEOS_LOG_LEVEL_RANK_0( 2, " redistributing mesh..." ); vtk::AllMeshes redistributedMeshes = - vtk::redistributeMeshes( getLogLevel(), allMeshes.getMainMesh(), allMeshes.getFaceBlocks(), allMeshes.getEmbeddedSurfaceBlocks(), 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(); diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp index 51bea1ab1e7..86712d33f9f 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp @@ -144,9 +144,9 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase /// 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 diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index fcc84ed887b..8d7d87da33b 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -566,7 +566,7 @@ loadMesh( Path const & filePath, AllMeshes loadAllMeshes( Path const & filePath, string const & mainBlockName, array1d< string > const & faceBlockNames, - array1d< string > const & edfmSurfBlockNames) + array1d< string > const & edfmSurfBlockNames ) { int const lastRank = MpiWrapper::commSize() - 1; vtkSmartPointer< vtkDataSet > main = loadMesh( filePath, mainBlockName ); @@ -693,9 +693,9 @@ AllMeshes redistributeByCellGraph( AllMeshes & input, vtkSmartPointer< vtkUnstructuredGrid > const finalFracMesh = vtk::redistribute( *splitFracMesh, MPI_COMM_GEOS ); finalFractures[fractureName] = finalFracMesh; } - + // Ouassim: just add the edfm mesh at the moment and see. - auto edfmMesh = input.getEmbeddedSurfaceBlocks(); + auto edfmMesh = input.getEmbeddedSurfaceBlocks(); return AllMeshes( finalMesh, finalFractures, edfmMesh ); } @@ -971,7 +971,7 @@ redistributeMeshes( integer const logLevel, { result.setMainMesh( mesh ); result.setFaceBlocks( namesToFractures ); - result.setEmbeddedSurfaceBlocks( namesToEdfmFractures); + result.setEmbeddedSurfaceBlocks( namesToEdfmFractures ); } // Logging some information about the redistribution. diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp index db08081a39d..e2b2fe3ed15 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.hpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.hpp @@ -79,10 +79,10 @@ class AllMeshes */ AllMeshes( vtkSmartPointer< vtkDataSet > const & main, std::map< string, vtkSmartPointer< vtkDataSet > > const & faceBlocks, - std::map> const & edfmSurfBlocks ) + std::map< string, vtkSmartPointer< vtkDataSet > > const & edfmSurfBlocks ) : m_main( main ), m_faceBlocks( faceBlocks ), - m_embeddedSurfaceBlocks(edfmSurfBlocks) + m_embeddedSurfaceBlocks( edfmSurfBlocks ) { } /** @@ -157,7 +157,7 @@ class AllMeshes AllMeshes loadAllMeshes( Path const & filePath, string const & mainBlockName, array1d< string > const & faceBlockNames, - array1d const & edfmSurfBlockNames); + array1d< string > const & edfmSurfBlockNames ); /** * @brief Compute the rank neighbor candidate list. diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index d49892e4931..5a5a2ab985d 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -170,94 +170,96 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() localIndex cellIndex = embSurf.getEmbeddedSurfElemTo3dElem().toCellIndex[edfmIndex][0]; subRegion.addFracturedElement( cellIndex, edfmIndex ); newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( edfmIndex ); - } } } - } );// end loop over subregions - - } - + } + } );// end loop over subregions - }else{ + } - // Loop over all the fracture planes - geometricObjManager.forGeometricObject< PlanarGeometricObject >( m_targetObjectsName, [&]( localIndex const, - PlanarGeometricObject & fracture ) + } + else { - /* 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 ) - { - 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 ) + // 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" ); + bool added = embeddedSurfaceSubRegion.addNewEmbeddedSurface( cellIndex, + er, + esr, + nodeManager, + embSurfNodeManager, + edgeManager, + cellToEdges, + &fracture ); - // Add the information to the CellElementSubRegion - subRegion.addFracturedElement( cellIndex, localNumberOfSurfaceElems ); + if( added ) + { + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::SurfaceGenerator, "Element " << cellIndex << " is fractured" ); - newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( localNumberOfSurfaceElems ); + // Add the information to the CellElementSubRegion + subRegion.addFracturedElement( cellIndex, localNumberOfSurfaceElems ); - localNumberOfSurfaceElems++; + newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( localNumberOfSurfaceElems ); + + localNumberOfSurfaceElems++; + } } } - } - } );// end loop over subregions + } );// end loop over subregions - }); - } - ); + } ); + } + ); } // Launch kernel to compute connectivity index of each fractured element. elemManager.forElementSubRegionsComplete< CellElementSubRegion >( From 6d70f054024a50c289d12f56908eb9d0d2c5a1a9 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Mon, 11 Nov 2024 10:11:01 -0600 Subject: [PATCH 25/32] remove comments --- src/coreComponents/mesh/generators/VTKMeshGenerator.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp index 4f8fd4a729d..7e05c026aa1 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp @@ -129,13 +129,11 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager for( auto const & [name, mesh]: m_faceBlockMeshes ) { - //TODO: Ouassim change for experimentation vtk::importFractureNetwork( name, mesh, m_vtkMesh, cellBlockManager ); } for( auto const & [name, mesh]: m_embeddedSurfaceBlockMeshes ) { - //TODO: Ouassim change for experimentation vtk::importEmbeddedFractureNetwork( name, mesh, m_vtkMesh, cellBlockManager ); } From c4a913b03b5f4e3ffc0d396f70e7781008d0e3ab Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 11 Nov 2024 10:17:57 -0600 Subject: [PATCH 26/32] use faceBlockName + some fixes for build (#3436) Co-authored-by: Pavel Tomin --- src/coreComponents/mesh/SurfaceElementRegion.cpp | 2 +- src/coreComponents/mesh/SurfaceElementRegion.hpp | 2 ++ .../mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp | 2 +- .../surfaceGeneration/EmbeddedSurfaceGenerator.cpp | 9 +++++---- 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/coreComponents/mesh/SurfaceElementRegion.cpp b/src/coreComponents/mesh/SurfaceElementRegion.cpp index f5918bb78f9..492a5f25603 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.cpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.cpp @@ -57,7 +57,7 @@ void SurfaceElementRegion::generateMesh( Group const & faceBlocks ) if( m_subRegionType == SurfaceSubRegionType::embeddedElement ) { // We just register the subregion copying of data is done at the EmbeddedSurfaceGenerator - EmbeddedSurfaceSubRegion & subRegion = elementSubRegions.registerGroup< EmbeddedSurfaceSubRegion >( m_faceBlockName ); + 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 182d4f730b0..46de8eb6bf8 100644 --- a/src/coreComponents/mesh/SurfaceElementRegion.hpp +++ b/src/coreComponents/mesh/SurfaceElementRegion.hpp @@ -134,6 +134,8 @@ class SurfaceElementRegion : public ElementRegionBase */ SurfaceSubRegionType subRegionType() const { return m_subRegionType; } + 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/VTKEmbeddedSurfaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp index a44c6c636a4..90103124d98 100644 --- a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp @@ -24,7 +24,7 @@ namespace geos::vtk void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName, vtkSmartPointer< vtkDataSet > embeddedSurfaceMesh, - vtkSmartPointer< vtkDataSet > mesh, + vtkSmartPointer< vtkDataSet > GEOS_UNUSED_PARAM( mesh ), CellBlockManager & cellBlockManager ) { diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index 5a5a2ab985d..efaf5f1ff9e 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -131,8 +131,9 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() NewObjectLists newObjects; - bool bulk_edfm_fill = true; - if( bulk_edfm_fill ) + string const & faceBlockName = embeddedSurfaceRegion.getFaceBlockName(); + + if( !faceBlockName.empty() ) { CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager(); @@ -141,9 +142,9 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() //Group const & embSurfBlocks = cellBlockManagerConcrete.getEmbeddedSurfaceBlocks(); Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks(); //if( embSurfBlocks.hasGroup( "EmbeddedSurface" )) - if( embSurfBlocks.hasGroup( "fracture" )) + if( embSurfBlocks.hasGroup( faceBlockName )) { - EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( "fracture" ); + EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( faceBlockName ); elemManager.forElementSubRegionsComplete< CellElementSubRegion >( From 584a601646416c59c7dce4b2e25aa274e9ec38a1 Mon Sep 17 00:00:00 2001 From: Ouassim Date: Wed, 13 Nov 2024 14:27:10 -0600 Subject: [PATCH 27/32] add edfm vtk example --- .../edfm_vtk_fractures.vtu | 75 +++++++++++ .../edfm_vtk_matrix.vtu | 122 ++++++++++++++++++ ...ed.xml => fractureMatrixFlow_edfm_vtk.xml} | 17 +-- ... => fractureMatrixFlow_edfm_vtk_smoke.xml} | 36 ++---- .../main_edfm_vtk.vtm | 6 + .../EmbeddedSurfaceGenerator.cpp | 4 - 6 files changed, 224 insertions(+), 36 deletions(-) create mode 100644 inputFiles/singlePhaseFlowFractures/edfm_vtk_fractures.vtu create mode 100644 inputFiles/singlePhaseFlowFractures/edfm_vtk_matrix.vtu rename inputFiles/singlePhaseFlowFractures/{fractureMatrixFlow_edfm_base_modified.xml => fractureMatrixFlow_edfm_vtk.xml} (90%) rename inputFiles/singlePhaseFlowFractures/{fractureMatrixFlow_edfm_horizontalFrac_smoke_modified.xml => fractureMatrixFlow_edfm_vtk_smoke.xml} (52%) create mode 100644 inputFiles/singlePhaseFlowFractures/main_edfm_vtk.vtm diff --git a/inputFiles/singlePhaseFlowFractures/edfm_vtk_fractures.vtu b/inputFiles/singlePhaseFlowFractures/edfm_vtk_fractures.vtu new file mode 100644 index 00000000000..129af955834 --- /dev/null +++ b/inputFiles/singlePhaseFlowFractures/edfm_vtk_fractures.vtu @@ -0,0 +1,75 @@ + + + + + + + 0 + 1 + 2 + 3 + 4 + 5 + + + + + 0 4 + 1 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_base_modified.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk.xml similarity index 90% rename from inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base_modified.xml rename to inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk.xml index d928a856304..48519910732 100644 --- a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_base_modified.xml +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk.xml @@ -15,14 +15,14 @@ directParallel="0"/> - + targetObjects="{ FracturePlane }" + mpiCommOrder="1"/> @@ -40,11 +40,12 @@ @@ -100,7 +101,7 @@ name="Porosity" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/RockMatrix/cb1" + objectPath="ElementRegions/RockMatrix/hexahedra" fieldName="rockPorosity_referencePorosity" scale="0.1"/> @@ -116,20 +117,20 @@ name="initialPressure" initialCondition="1" setNames="{ all }" - objectPath="ElementRegions/RockMatrix/cb1" + objectPath="ElementRegions/RockMatrix/hexahedra" fieldName="pressure" scale="0.0"/> diff --git a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke_modified.xml b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk_smoke.xml similarity index 52% rename from inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke_modified.xml rename to inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk_smoke.xml index 45a837b1436..bc1a60fa14e 100644 --- a/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_horizontalFrac_smoke_modified.xml +++ b/inputFiles/singlePhaseFlowFractures/fractureMatrixFlow_edfm_vtk_smoke.xml @@ -2,57 +2,45 @@ - + - + + xMax="{ 0.34, 0.34, 0.51 }"/> - - + target="/Solvers/SurfaceGenerator"/> - + target="/Outputs/vtkOutput"/> + + + + + diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index efaf5f1ff9e..51d90859904 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -137,11 +137,7 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() { CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager(); - // CellBlockManagerABC & cellBlockManagerNoConst = const_cast< CellBlockManagerABC & >(domain.getMeshBody( 0 ).getCellBlockManager()); - // CellBlockManager & cellBlockManagerConcrete = dynamic_cast< CellBlockManager & >(cellBlockManagerNoConst); - //Group const & embSurfBlocks = cellBlockManagerConcrete.getEmbeddedSurfaceBlocks(); Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks(); - //if( embSurfBlocks.hasGroup( "EmbeddedSurface" )) if( embSurfBlocks.hasGroup( faceBlockName )) { EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( faceBlockName ); From c84e27a76b655d8114bbb9f2c938fc2894630edb Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 9 Dec 2024 17:52:26 -0600 Subject: [PATCH 28/32] build fix --- .../mesh/generators/CellBlockManager.cpp | 2 - .../EmbeddedSurfaceGenerator.cpp | 49 +++++++------------ src/docs/doxygen/GeosxConfig.hpp | 2 +- 3 files changed, 20 insertions(+), 33 deletions(-) diff --git a/src/coreComponents/mesh/generators/CellBlockManager.cpp b/src/coreComponents/mesh/generators/CellBlockManager.cpp index 11b896b6d6d..c7a4dd57958 100644 --- a/src/coreComponents/mesh/generators/CellBlockManager.cpp +++ b/src/coreComponents/mesh/generators/CellBlockManager.cpp @@ -669,7 +669,6 @@ void CellBlockManager::buildMaps() buildNodeToEdges(); fillElementToEdgesOfCellBlocks( m_faceToEdges.toViewConst(), this->getCellBlocks() ); - } ArrayOfArrays< localIndex > CellBlockManager::getFaceToNodes() const @@ -748,7 +747,6 @@ localIndex CellBlockManager::numFaces() const return m_numFaces; } - ArrayOfArrays< localIndex > CellBlockManager::getEdgeToFaces() const { return m_edgeToFaces; diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index 1b94505daa3..f3609564748 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -135,51 +135,40 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() if( !faceBlockName.empty() ) { - CellBlockManagerABC const & cellBlockManager = meshBody.getCellBlockManager(); Group const & embSurfBlocks = cellBlockManager.getEmbeddedSurfaceBlocks(); if( embSurfBlocks.hasGroup( faceBlockName )) { EmbeddedSurfaceBlockABC const & embSurf = embSurfBlocks.getGroup< EmbeddedSurfaceBlockABC >( faceBlockName ); - elemManager.forElementSubRegionsComplete< CellElementSubRegion >( [&]( localIndex const er, localIndex const esr, ElementRegionBase &, CellElementSubRegion & subRegion ) { - arrayView2d< localIndex const, cells::NODE_MAP_USD > const cellToNodes = subRegion.nodeList(); FixedOneToManyRelation const & cellToEdges = subRegion.edgeList(); - arrayView1d< integer const > const ghostRank = subRegion.ghostRank(); + bool added = embeddedSurfaceSubRegion.addAllEmbeddedSurfaces( er, + esr, + nodeManager, + embSurfNodeManager, + edgeManager, + cellToEdges, + embSurf ); + + if( added ) { - bool added = embeddedSurfaceSubRegion.addAllEmbeddedSurfaces( er, - esr, - nodeManager, - embSurfNodeManager, - edgeManager, - cellToEdges, - embSurf ); - - if( added ) + // Add all the fracture information to the CellElementSubRegion + for( localIndex edfmIndex=0; edfmIndex < embSurf.numEmbeddedSurfElem(); ++edfmIndex ) { - // 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 ); - } + localIndex cellIndex = embSurf.getEmbeddedSurfElemTo3dElem().toCellIndex[edfmIndex][0]; + subRegion.addFracturedElement( cellIndex, edfmIndex ); + newObjects.newElements[ {embeddedSurfaceRegion.getIndexInParent(), embeddedSurfaceSubRegion.getIndexInParent()} ].insert( edfmIndex ); } } - } );// end loop over subregions - + } ); // end loop over subregions } - - } else { - - // Loop over all the fracture planes geometricObjManager.forGeometricObject< PlanarGeometricObject >( m_targetObjectsName, [&]( localIndex const, PlanarGeometricObject & fracture ) @@ -228,6 +217,7 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() isNegative = 1; } } // end loop over nodes + if( isPositive * isNegative == 1 ) { bool added = embeddedSurfaceSubRegion.addNewEmbeddedSurface( cellIndex, @@ -252,12 +242,11 @@ void EmbeddedSurfaceGenerator::initializePostSubGroups() } } } - } );// end loop over subregions - + } ); // end loop over subregions } ); - } - ); + } ); } + // Launch kernel to compute connectivity index of each fractured element. elemManager.forElementSubRegionsComplete< CellElementSubRegion >( [&]( localIndex const, localIndex const, ElementRegionBase &, CellElementSubRegion & subRegion ) diff --git a/src/docs/doxygen/GeosxConfig.hpp b/src/docs/doxygen/GeosxConfig.hpp index f264cce1d54..44a7b9be831 100644 --- a/src/docs/doxygen/GeosxConfig.hpp +++ b/src/docs/doxygen/GeosxConfig.hpp @@ -171,7 +171,7 @@ #define fmt_VERSION 10.0.0 /// Version information for python -#define Python3_VERSION 3.10.8 +#define Python3_VERSION 3.10.6 /// Version information for CUDAToolkit /* #undef CUDAToolkit_VERSION */ From b4a5df4a2af75b43ea90e46ce12c2cff6a9534b6 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Tue, 10 Dec 2024 09:24:13 -0600 Subject: [PATCH 29/32] remove --- .../surfaceGeneration/EmbeddedSurfaceGenerator.cpp | 3 --- .../surfaceGeneration/EmbeddedSurfaceGenerator.hpp | 2 -- 2 files changed, 5 deletions(-) diff --git a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp index f3609564748..1812c5d4dbd 100644 --- a/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp +++ b/src/coreComponents/physicsSolvers/surfaceGeneration/EmbeddedSurfaceGenerator.cpp @@ -308,9 +308,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 043530d94ad..38496bcf178 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" ); From 3e0357e2e4c03bf55d893d66a16a7fa8e505c6e5 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 16 Dec 2024 16:17:45 -0600 Subject: [PATCH 30/32] doxygen fix --- src/coreComponents/mesh/generators/VTKUtilities.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp index 864d0903388..d064595f6dd 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.hpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.hpp @@ -129,7 +129,7 @@ class AllMeshes /** * @brief Defines the face blocks/fractures. - * @param edfmSurBlocks A map which connects each name of the face block to its mesh. + * @param edfmSurfBlocks A map which connects each name of the face block to its mesh. */ void setEmbeddedSurfaceBlocks( std::map< string, vtkSmartPointer< vtkDataSet > > const & edfmSurfBlocks ) { From 13889707847a94603a489391474bc5f7131c97e4 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 16 Dec 2024 17:44:08 -0600 Subject: [PATCH 31/32] fix one crash and introduce another --- .../mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp index 90103124d98..514306d37d8 100644 --- a/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKEmbeddedSurfaceBlockUtilities.cpp @@ -33,7 +33,7 @@ void importEmbeddedFractureNetwork( string const & embeddedSurfaceBlockName, //1.Get EDFM vertices vtkUnstructuredGrid * grid = vtkUnstructuredGrid::SafeDownCast( embeddedSurfaceMesh ); vtkPoints * const nodes = grid->GetPoints(); - vtkIdType numNodes = nodes->GetNumberOfPoints(); + vtkIdType numNodes = nodes ? nodes->GetNumberOfPoints() : 0; ArrayOfArrays< localIndex > elem2dToNodes( numEdfmFracs ); for( vtkIdType i = 0; i < numEdfmFracs; ++i ) From ef54e9c62652bff419e6489445534917e7096916 Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Mon, 16 Dec 2024 17:52:40 -0600 Subject: [PATCH 32/32] just to see if test fail because of this --- src/coreComponents/mainInterface/ProblemManager.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coreComponents/mainInterface/ProblemManager.cpp b/src/coreComponents/mainInterface/ProblemManager.cpp index 18c55131919..7cc0ecfc9d0 100644 --- a/src/coreComponents/mainInterface/ProblemManager.cpp +++ b/src/coreComponents/mainInterface/ProblemManager.cpp @@ -696,7 +696,7 @@ void ProblemManager::generateMesh() // meshBody.deregisterGroup( keys::cellManager ); // Cell block manager is needed for EDFM loading from VTK // Can't be deregistred here - //meshBody.deregisterCellBlockManager(); + meshBody.deregisterCellBlockManager(); } meshBody.forMeshLevels( [&]( MeshLevel & meshLevel )