diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt index 3b7e7a825cb..c3f5c23d540 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -29,6 +29,7 @@ set( mesh_headers CellElementRegion.hpp CellElementRegionSelector.hpp CellElementSubRegion.hpp + CellElementSubRegionVariable.hpp DomainPartition.hpp EdgeManager.hpp ElementRegionBase.hpp @@ -126,6 +127,7 @@ set( mesh_sources CellElementRegion.cpp CellElementRegionSelector.cpp CellElementSubRegion.cpp + CellElementSubRegionVariable.cpp DomainPartition.cpp EdgeManager.cpp ElementRegionBase.cpp diff --git a/src/coreComponents/mesh/CellElementSubRegionVariable.cpp b/src/coreComponents/mesh/CellElementSubRegionVariable.cpp new file mode 100644 index 00000000000..bed60438bfb --- /dev/null +++ b/src/coreComponents/mesh/CellElementSubRegionVariable.cpp @@ -0,0 +1,298 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +#include "CellElementSubRegionVariable.hpp" + +#include "common/TypeDispatch.hpp" +#include "mesh/MeshLevel.hpp" +#include "mesh/generators/CellBlockUtilities.hpp" + +namespace geos +{ +using namespace dataRepository; + +CellElementSubRegionVariable::CellElementSubRegionVariable( string const & name, Group * const parent ): + ElementSubRegionBase( name, parent ) +{ + registerWrapper( viewKeyStruct::nodeListString(), &m_toNodesRelation ); + registerWrapper( viewKeyStruct::edgeListString(), &m_toEdgesRelation ); + registerWrapper( viewKeyStruct::faceListString(), &m_toFacesRelation ); + + excludeWrappersFromPacking( { viewKeyStruct::nodeListString(), + viewKeyStruct::edgeListString(), + viewKeyStruct::faceListString() } ); +} + + +void CellElementSubRegionVariable::copyFromCellBlock( CellBlockABC const & cellBlock ) +{ + // Defines the (unique) element type of this cell element region, + // and its associated number of nodes, edges, faces. + m_elementType = cellBlock.getElementType(); + m_numNodesPerElement = cellBlock.numNodesPerElement(); + m_numEdgesPerElement = cellBlock.numEdgesPerElement(); + m_numFacesPerElement = cellBlock.numFacesPerElement(); + + // We call the `resize` member function of the cell to (nodes, edges, faces) relations, + // before calling the `CellElementSubRegionVariable::resize` in order to keep the first dimension. + // Be careful when refactoring. + m_toNodesRelation.resize( this->size(), m_numNodesPerElement ); + m_toEdgesRelation.resize( this->size(), m_numEdgesPerElement ); + m_toFacesRelation.resize( this->size(), m_numFacesPerElement ); + this->resize( cellBlock.numElements() ); + + // this->nodeList() = cellBlock.getElemToNodes(); + // this->edgeList() = cellBlock.getElemToEdges(); + // this->faceList() = cellBlock.getElemToFaces(); + + this->m_localToGlobalMap = cellBlock.localToGlobalMap(); + + this->constructGlobalToLocalMap(); + cellBlock.forExternalProperties( [&]( WrapperBase const & wrapper ) + { + types::dispatch( types::ListofTypeList< types::StandardArrays >{}, [&]( auto tupleOfTypes ) + { + using ArrayType = camp::first< decltype( tupleOfTypes ) >; + auto const src = Wrapper< ArrayType >::cast( wrapper ).reference().toViewConst(); + this->registerWrapper( wrapper.getName(), std::make_unique< ArrayType >( &src ) ); + }, wrapper ); + } ); +} + + + +// localIndex CellElementSubRegionVariable::packUpDownMapsSize( arrayView1d< localIndex const > const & packList ) const +// { +// buffer_unit_type * junk = nullptr; +// return packUpDownMapsImpl< false >( junk, packList ); +// } + + +// localIndex CellElementSubRegionVariable::packUpDownMaps( buffer_unit_type * & buffer, +// arrayView1d< localIndex const > const & packList ) const +// { +// return packUpDownMapsImpl< true >( buffer, packList ); +// } + +// template< bool DO_PACKING > +// localIndex CellElementSubRegionVariable::packUpDownMapsImpl( buffer_unit_type * & buffer, +// arrayView1d< localIndex const > const & packList ) const +// { + +// arrayView1d< globalIndex const > const localToGlobal = this->localToGlobalMap(); +// arrayView1d< globalIndex const > nodeLocalToGlobal = nodeList().relatedObjectLocalToGlobal(); +// arrayView1d< globalIndex const > edgeLocalToGlobal = edgeList().relatedObjectLocalToGlobal(); +// arrayView1d< globalIndex const > faceLocalToGlobal = faceList().relatedObjectLocalToGlobal(); + + +// localIndex packedSize = bufferOps::Pack< DO_PACKING >( buffer, +// nodeList().base().toViewConst(), +// m_unmappedGlobalIndicesInNodelist, +// packList, +// localToGlobal, +// nodeLocalToGlobal ); + +// packedSize += bufferOps::Pack< DO_PACKING >( buffer, +// edgeList().base().toViewConst(), +// m_unmappedGlobalIndicesInEdgelist, +// packList, +// localToGlobal, +// edgeLocalToGlobal ); + + +// packedSize += bufferOps::Pack< DO_PACKING >( buffer, +// faceList().base().toViewConst(), +// m_unmappedGlobalIndicesInFacelist, +// packList, +// localToGlobal, +// faceLocalToGlobal ); + +// return packedSize; +// } + +// localIndex CellElementSubRegionVariable::unpackUpDownMaps( buffer_unit_type const * & buffer, +// localIndex_array & packList, +// bool const GEOS_UNUSED_PARAM( overwriteUpMaps ), +// bool const GEOS_UNUSED_PARAM( overwriteDownMaps ) ) +// { +// localIndex unPackedSize = 0; +// unPackedSize += bufferOps::Unpack( buffer, +// nodeList().base().toView(), +// packList, +// m_unmappedGlobalIndicesInNodelist, +// this->globalToLocalMap(), +// nodeList().relatedObjectGlobalToLocal() ); + +// unPackedSize += bufferOps::Unpack( buffer, +// edgeList().base(), +// packList, +// m_unmappedGlobalIndicesInEdgelist, +// this->globalToLocalMap(), +// edgeList().relatedObjectGlobalToLocal() ); + +// unPackedSize += bufferOps::Unpack( buffer, +// faceList().base(), +// packList, +// m_unmappedGlobalIndicesInFacelist, +// this->globalToLocalMap(), +// faceList().relatedObjectGlobalToLocal() ); + + +// GEOS_ERROR_IF_NE( m_unmappedGlobalIndicesInNodelist.size(), 0 ); +// GEOS_ERROR_IF_NE( m_unmappedGlobalIndicesInEdgelist.size(), 0 ); +// GEOS_ERROR_IF_NE( m_unmappedGlobalIndicesInFacelist.size(), 0 ); + +// return unPackedSize; +// } + + +localIndex CellElementSubRegionVariable::getFaceNodes( localIndex const elementIndex, + localIndex const localFaceIndex, + Span< localIndex > const nodeIndices ) const +{ +// return geos::getFaceNodes( m_elementType, elementIndex, localFaceIndex, m_toNodesRelation, nodeIndices ); + return 0; // To be implemented +} + +void CellElementSubRegionVariable:: + calculateElementCenterAndVolume( localIndex const k, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X ) const +{ + auto getElementCoordinatesaAndComputeElementCenter = [k, X, this]( auto & XLocal ) + { + LvArray::tensorOps::fill< 3 >( m_elementCenter[k], 0 ); + for( localIndex a = 0; a < m_numNodesPerElement; ++a ) + { + LvArray::tensorOps::copy< 3 >( XLocal[a], X[m_toNodesRelation( k, a )] ); + LvArray::tensorOps::add< 3 >( m_elementCenter[k], XLocal[a] ); + } + LvArray::tensorOps::scale< 3 >( m_elementCenter[k], 1.0 / m_numNodesPerElement ); + }; + + switch( m_elementType ) + { + case ElementType::Hexahedron: + { + real64 Xlocal[8][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::hexahedronVolume( Xlocal ); + break; + } + case ElementType::Tetrahedron: + { + real64 Xlocal[4][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::tetrahedronVolume( Xlocal ); + break; + } + case ElementType::Wedge: + { + real64 Xlocal[6][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::wedgeVolume( Xlocal ); + break; + } + case ElementType::Pyramid: + { + real64 Xlocal[5][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::pyramidVolume( Xlocal ); + break; + } + case ElementType::Prism5: + { + real64 Xlocal[10][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::prismVolume< 5 >( Xlocal ); + break; + } + case ElementType::Prism6: + { + real64 Xlocal[12][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::prismVolume< 6 >( Xlocal ); + break; + } + case ElementType::Prism7: + { + real64 Xlocal[14][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::prismVolume< 7 >( Xlocal ); + break; + } + case ElementType::Prism8: + { + real64 Xlocal[16][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::prismVolume< 8 >( Xlocal ); + break; + } + case ElementType::Prism9: + { + real64 Xlocal[18][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::prismVolume< 9 >( Xlocal ); + break; + } + case ElementType::Prism10: + { + real64 Xlocal[20][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::prismVolume< 10 >( Xlocal ); + break; + } + case ElementType::Prism11: + { + real64 Xlocal[22][3]; + getElementCoordinatesaAndComputeElementCenter( Xlocal ); + m_elementVolume[k] = computationalGeometry::prismVolume< 11 >( Xlocal ); + break; + } + default: + { + GEOS_ERROR( GEOS_FMT( "Volume calculation not supported for element type {} in subregion {}", + m_elementType, getDataContext() ) ); + } + } + + GEOS_ERROR_IF( m_elementVolume[k] <= 0.0, + GEOS_FMT( "Negative volume for element {} type {} in subregion {}", + k, m_elementType, getDataContext() ) ); +} + +void CellElementSubRegionVariable::calculateElementGeometricQuantities( NodeManager const & nodeManager, + FaceManager const & GEOS_UNUSED_PARAM( faceManager ) ) +{ + GEOS_MARK_FUNCTION; + + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition(); + + forAll< parallelHostPolicy >( this->size(), [=] ( localIndex const k ) + { + calculateElementCenterAndVolume( k, X ); + } ); +} + +void CellElementSubRegionVariable::setupRelatedObjectsInRelations( MeshLevel const & mesh ) +{ + this->m_toNodesRelation.setRelatedObject( mesh.getNodeManager() ); + this->m_toEdgesRelation.setRelatedObject( mesh.getEdgeManager() ); + this->m_toFacesRelation.setRelatedObject( mesh.getFaceManager() ); +} + +REGISTER_CATALOG_ENTRY( ObjectManagerBase, CellElementSubRegionVariable, string const &, Group * const ) + +} /* namespace geos */ diff --git a/src/coreComponents/mesh/CellElementSubRegionVariable.hpp b/src/coreComponents/mesh/CellElementSubRegionVariable.hpp new file mode 100644 index 00000000000..c8d9043f415 --- /dev/null +++ b/src/coreComponents/mesh/CellElementSubRegionVariable.hpp @@ -0,0 +1,246 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +#ifndef GEOS_MESH_CELLELEMENTSUBREGIONVARIABLE_HPP_ +#define GEOS_MESH_CELLELEMENTSUBREGIONVARIABLE_HPP_ + +#include "mesh/generators/CellBlockABC.hpp" +#include "mesh/NodeManager.hpp" +#include "mesh/FaceManager.hpp" +#include "mesh/utilities/ComputationalGeometry.hpp" +#include "ElementSubRegionBase.hpp" + + +namespace geos +{ + +class MeshLevel; + +/** + * @class CellElementSubRegionVariable + * Class specializing the element subregion for a cell element. + * This is the class used in the physics solvers to represent a collection of mesh cell elements + */ +class CellElementSubRegionVariable : public ElementSubRegionBase +{ +public: + + /// Alias for the type of the element-to-node map + using NodeMapType = InterObjectRelation< ArrayOfArrays< localIndex > >; + /// Alias for the type of the element-to-edge map + using EdgeMapType = InterObjectRelation< ArrayOfArrays< localIndex > >; + /// Alias for the type of the element-to-face map + using FaceMapType = InterObjectRelation< ArrayOfArrays< localIndex > >; + /// Type of map between cell blocks and embedded elements + using EmbSurfMapType = InterObjectRelation< ArrayOfArrays< localIndex > >; + + /** + * @brief Const getter for the catalog name. + * @return the name of this type in the catalog + */ + static string catalogName() + { return "CellElementSubRegionVariable"; } + + /** + * @copydoc catalogName() + */ + virtual string getCatalogName() const override final + { return catalogName(); } + + /** + * @name Constructor / Destructor + */ + ///@{ + + /** + * @brief Constructor for this class. + * @param[in] name the name of this object manager + * @param[in] parent the parent Group + */ + CellElementSubRegionVariable( string const & name, Group * const parent ); + + ///@} + + /** + * @name Helpers for CellElementSubRegionVariable construction + */ + ///@{ + + /** + * @brief Fill the CellElementSubRegionVariable by copying those of the source CellBlock + * @param cellBlock the CellBlock which properties (connectivity info) will be copied. + */ + void copyFromCellBlock( CellBlockABC const & cellBlock ); + + ///@} + + + /** + * @name Overriding packing / Unpacking functions + */ + ///@{ + + virtual localIndex packUpDownMapsSize( arrayView1d< localIndex const > const & packList ) const override; + + virtual localIndex packUpDownMaps( buffer_unit_type * & buffer, + arrayView1d< localIndex const > const & packList ) const override; + + virtual localIndex unpackUpDownMaps( buffer_unit_type const * & buffer, + array1d< localIndex > & packList, + bool const overwriteUpMaps, + bool const overwriteDownMaps ) override; + + + ///@} + + /** + * @brief struct to serve as a container for variable strings and keys + * @struct viewKeyStruct + */ + struct viewKeyStruct : public ElementSubRegionBase::viewKeyStruct + { + } + /// viewKey struct for the CellElementSubRegionVariable class + m_CellBlockSubRegionViewKeys; + + virtual viewKeyStruct & viewKeys() override { return m_CellBlockSubRegionViewKeys; } + virtual viewKeyStruct const & viewKeys() const override { return m_CellBlockSubRegionViewKeys; } + + /** + * @brief Get the local indices of the nodes in a face of the element. + * @param[in] elementIndex The local index of the target element. + * @param[in] localFaceIndex The local index of the target face in the element (this will be [0, numFacesInElement[) + * @param[out] nodeIndices Memory to which node indices for the face will be written, must have sufficient size. + * @return tne number of values written into @p nodeIndices + * @deprecated This method will be removed soon. + */ + localIndex getFaceNodes( localIndex const elementIndex, + localIndex const localFaceIndex, + Span< localIndex > const nodeIndices ) const; + + /** + * @brief Get the element-to-node map. + * @return a reference to the element-to-node map + */ + NodeMapType & nodeList() { return m_toNodesRelation; } + + /** + * @copydoc nodeList() + */ + NodeMapType const & nodeList() const { return m_toNodesRelation; } + + /** + * @brief Get the element-to-edge map. + * @return a reference to element-to-edge map + */ + EdgeMapType & edgeList() { return m_toEdgesRelation; } + + /** + * @copydoc edgeList() + */ + EdgeMapType const & edgeList() const { return m_toEdgesRelation; } + + /** + * @brief Get the element-to-face map. + * @return a reference to the element to face map + */ + FaceMapType & faceList() { return m_toFacesRelation; } + + /** + * @copydoc faceList() + */ + FaceMapType const & faceList() const { return m_toFacesRelation; } + + + + /** + * @brief Compute the center of each element in the subregion. + * @param[in] X an arrayView of (const) node positions + */ + void calculateElementCenters( arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X ) const + { + ElementSubRegionBase::calculateElementCenters( m_toNodesRelation, X ); + } + + void calculateElementGeometricQuantities( NodeManager const & nodeManager, + FaceManager const & faceManager ) override; + +private: + + /// Element-to-node relation + NodeMapType m_toNodesRelation; + + /// Element-to-edge relation + EdgeMapType m_toEdgesRelation; + + /// Element-to-face relation + FaceMapType m_toFacesRelation; + + /// Name of the properties registered from an external mesh + string_array m_externalPropertyNames; + + /** + * @brief Compute the volume of the k-th element in the subregion. + * @param[in] k the index of the element in the subregion + * @param[in] X an arrayView of (const) node positions + */ + void calculateCellVolumesKernel( localIndex const k, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X ) const; + + void calculateElementCenterAndVolume( localIndex const k, + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const & X ) const; + + /// The array of jacobian determinantes. + array2d< real64 > m_detJ; + + /// Map of unmapped global indices in the element-to-node map + map< localIndex, array1d< globalIndex > > m_unmappedGlobalIndicesInNodelist; + + /// Map of unmapped global indices in the element-to-face map + map< localIndex, array1d< globalIndex > > m_unmappedGlobalIndicesInEdgelist; + + /// Map of unmapped global indices in the element-to-face map + map< localIndex, array1d< globalIndex > > m_unmappedGlobalIndicesInFacelist; + + + /** + * @brief Pack element-to-node and element-to-face maps + * @tparam DO_PACKING the flag for the bufferOps::Pack function + * @param buffer the buffer used in the bufferOps::Pack function + * @param packList the packList used in the bufferOps::Pack function + * @return the pack size + */ + template< bool DO_PACKING > + localIndex packUpDownMapsImpl( buffer_unit_type * & buffer, + arrayView1d< localIndex const > const & packList ) const; + + /** + * @brief Links the managers to their mappings. + * @param[in] mesh Holds the node/edge/face managers. + * + * Defines links the element to nodes, edges and faces to their respective node/edge/face managers. + */ + void setupRelatedObjectsInRelations( MeshLevel const & mesh ) override; + + template< bool DO_PACKING > + localIndex packFracturedElementsImpl( buffer_unit_type * & buffer, + arrayView1d< localIndex const > const & packList, + arrayView1d< globalIndex const > const & embeddedSurfacesLocalToGlobal ) const; +}; + +} /* namespace geos */ + +#endif /* GEOS_MESH_CELLELEMENTSUBREGIONVARIABLE_HPP_ */