Skip to content

feat: add variable size subregion #3820

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 2 commits into
base: develop
Choose a base branch
from
Draft

feat: add variable size subregion #3820

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1e8b8b7
add CellElementSubRegionVariable
rrsettgast Sep 16, 2025
7590cb5
Merge branch 'develop' into feature/variableSizeSubregion
rrsettgast Sep 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions src/coreComponents/mesh/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ set( mesh_headers
CellElementRegion.hpp
CellElementRegionSelector.hpp
CellElementSubRegion.hpp
CellElementSubRegionVariable.hpp
DomainPartition.hpp
EdgeManager.hpp
ElementRegionBase.hpp
Expand Down Expand Up @@ -126,6 +127,7 @@ set( mesh_sources
CellElementRegion.cpp
CellElementRegionSelector.cpp
CellElementSubRegion.cpp
CellElementSubRegionVariable.cpp
DomainPartition.cpp
EdgeManager.cpp
ElementRegionBase.cpp
Expand Down
298 changes: 298 additions & 0 deletions src/coreComponents/mesh/CellElementSubRegionVariable.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 */
Loading
Loading