Skip to content

Return geometry by value (alternative 2) #884

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: master
Choose a base branch
from
Draft

Return geometry by value (alternative 2) #884

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
b59a352
Return geometry objects by value
SoilRos Jun 13, 2025
f2e1a87
Store a variant in the geometry
SoilRos Jun 13, 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
8 changes: 4 additions & 4 deletions opm/grid/CpGrid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1743,28 +1743,28 @@ namespace Dune
/// \brief Get the Position of a vertex.
/// \param cell The index identifying the cell.
/// \return The coordinates of the vertex.
const Vector& vertexPosition(int vertex) const;
Vector vertexPosition(int vertex) const;

/// \brief Get the area of a face.
/// \param cell The index identifying the face.
double faceArea(int face) const;

/// \brief Get the coordinates of the center of a face.
/// \param cell The index identifying the face.
const Vector& faceCentroid(int face) const;
Vector faceCentroid(int face) const;

/// \brief Get the unit normal of a face.
/// \param cell The index identifying the face.
/// \see faceCell
const Vector& faceNormal(int face) const;
Vector faceNormal(int face) const;

/// \brief Get the volume of the cell.
/// \param cell The index identifying the cell.
double cellVolume(int cell) const;

/// \brief Get the coordinates of the center of a cell.
/// \param cell The index identifying the face.
const Vector& cellCentroid(int cell) const;
Vector cellCentroid(int cell) const;

/// \brief An iterator over the centroids of the geometry of the entities.
/// \tparam codim The co-dimension of the entities.
Expand Down
40 changes: 19 additions & 21 deletions opm/grid/cpgrid/CpGrid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1755,7 +1755,7 @@ const Dune::FieldVector<double,3> CpGrid::faceAreaNormalEcl(int face) const
}
}

const Dune::FieldVector<double,3>& CpGrid::vertexPosition(int vertex) const
Dune::FieldVector<double,3> CpGrid::vertexPosition(int vertex) const
{
return current_view_data_->geomVector<3>()[cpgrid::EntityRep<3>(vertex, true)].center();
}
Expand All @@ -1765,12 +1765,12 @@ double CpGrid::faceArea(int face) const
return current_view_data_->geomVector<1>()[cpgrid::EntityRep<1>(face, true)].volume();
}

const Dune::FieldVector<double,3>& CpGrid::faceCentroid(int face) const
Dune::FieldVector<double,3> CpGrid::faceCentroid(int face) const
{
return current_view_data_->geomVector<1>()[cpgrid::EntityRep<1>(face, true)].center();
}

const Dune::FieldVector<double,3>& CpGrid::faceNormal(int face) const
Dune::FieldVector<double,3> CpGrid::faceNormal(int face) const
{
return current_view_data_->face_normals_.get(face);
}
Expand All @@ -1780,7 +1780,7 @@ double CpGrid::cellVolume(int cell) const
return current_view_data_->geomVector<0>()[cpgrid::EntityRep<0>(cell, true)].volume();
}

const Dune::FieldVector<double,3>& CpGrid::cellCentroid(int cell) const
Dune::FieldVector<double,3> CpGrid::cellCentroid(int cell) const
{
return current_view_data_->geomVector<0>()[cpgrid::EntityRep<0>(cell, true)].center();
}
Expand Down Expand Up @@ -2857,10 +2857,10 @@ void CpGrid::refineAndProvideMarkedRefinedRelations( /* Marked elements paramete
int& cell_count,
std::vector<std::vector<int>>& preAdapt_level_to_leaf_cells_vec,
/* Additional parameters */
const std::vector<std::array<int,3>>& cells_per_dim_vec) const
const std::vector<std::array<int,3>>& cells_per_dim_vec) const
{
// If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.

// Each marked element for refinement (mark equal to 1), will be refined individuality, creating its own Lgr. The element index will
// be also used to identify its lgr. Even though, in the end, all the refined entities will belong to a unique level grid.
// For this reason, we associate "-1" with those elements that are not involved in any refinement and will appear
Expand All @@ -2878,7 +2878,7 @@ void CpGrid::refineAndProvideMarkedRefinedRelations( /* Marked elements paramete
cell_count +=1;
preAdapt_level_to_leaf_cells_vec[element.level()][element.getLevelElem().index()] = cell_count;
}

// When the element is marked for refinement, we also mark its corners and faces
// since they will get replaced by refined ones.
if (getMark(element) == 1) {
Expand Down Expand Up @@ -2914,7 +2914,7 @@ void CpGrid::refineAndProvideMarkedRefinedRelations( /* Marked elements paramete
refined_cell_count_vec[shiftedLevel] +=1;

}

preAdapt_parent_to_children_cells_vec[element.level()][element.getLevelElem().index()] = std::make_pair( markedElemLevel, refinedChildrenList);
for (const auto& [markedCorner, lgrEquivCorner] : parentCorners_to_equivalentRefinedCorners) {
cornerInMarkedElemWithEquivRefinedCorner[markedCorner].push_back({elemIdx, lgrEquivCorner});
Expand All @@ -2934,7 +2934,7 @@ CpGrid::defineChildToParentAndIdxInParentCell(const std::map<std::array<int,2>,s
const int& cell_count) const
{
// If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.

// ------------------------ Refined grid parameters
// Refined child cells and their parents. Entry is {-1,-1} when cell has no father. Otherwise, {level parent cell, parent cell index}
// Each entry represents a refined level.
Expand Down Expand Up @@ -3304,7 +3304,7 @@ void CpGrid::identifyLeafGridCorners(std::map<std::array<int,2>,int>& elemLgrAnd
const std::vector<std::array<int,3>>& cells_per_dim_vec) const
{
// If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.

// Step 1. Select/store the corners from the starting grid not involved in any (new) LGR.
// Replace the corners from level zero involved in LGR by the equivalent ones, born in LGRs.
// In this case, we avoid repetition considering the last appearance of the level zero corner
Expand Down Expand Up @@ -3458,8 +3458,8 @@ void CpGrid::identifyLeafGridFaces(std::map<std::array<int,2>,int>& elemLgrAndEl
// If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.

// Max level before calling adapt.
const int& preAdaptMaxLevel = this->maxLevel();
const int& preAdaptMaxLevel = this->maxLevel();

// Step 1. Add the LGR faces, for each LGR
for (int elem = 0; elem < current_view_data_->size(0); ++elem) {
if (markedElem_to_itsLgr[elem]!=nullptr) {
Expand Down Expand Up @@ -3520,7 +3520,7 @@ void CpGrid::populateLeafGridCorners(Dune::cpgrid::EntityVariableBase<cpgrid::Ge
const std::unordered_map<int,std::array<int,2>>& adaptedCorner_to_elemLgrAndElemLgrCorner) const
{
// If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.

adapted_corners.resize(corner_count);
for (int corner = 0; corner < corner_count; ++corner) {
const auto& [elemLgr, elemLgrCorner] = adaptedCorner_to_elemLgrAndElemLgrCorner.at(corner);
Expand Down Expand Up @@ -3561,7 +3561,7 @@ void CpGrid::populateLeafGridFaces(Dune::cpgrid::EntityVariableBase<cpgrid::Geom
const int& preAdaptMaxLevel) const
{
// If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.

adapted_faces.resize(face_count);
mutable_face_tags.resize(face_count);
mutable_face_normals.resize(face_count);
Expand Down Expand Up @@ -3741,7 +3741,7 @@ void CpGrid::populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geom
const int& preAdaptMaxLevel) const
{
// If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.

// --- Adapted cells ---
// Store the adapted cells. Main difficulty: to lookup correctly the indices of the corners and faces of each cell.
adapted_cells.resize(cell_count);
Expand Down Expand Up @@ -3845,7 +3845,7 @@ void CpGrid::populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geom

// Create a pointer to the first element of "adapted_cell_to_point" (required as the fourth argement to construct a Geometry<3,3> type object).
int* indices_storage_ptr = adapted_cell_to_point[cell].data();
adapted_cells[cell] = cpgrid::Geometry<3,3>(cellGeom.center(), cellGeom.volume(), allCorners, indices_storage_ptr);
adapted_cells[cell] = cpgrid::Geometry<3,3>(cellGeom.center(), cellGeom.volume(), allCorners.get(), indices_storage_ptr);
} // adapted_cells

// Adapted/Leaf-grid-view face to cell.
Expand All @@ -3871,7 +3871,7 @@ void CpGrid::populateRefinedCells(std::vector<Dune::cpgrid::EntityVariableBase<c
const std::map<std::array<int,2>,int>& markedElemAndEquivRefinedCorn_to_corner,
const std::vector<std::vector<std::array<int,2>>>& cornerInMarkedElemWithEquivRefinedCorner,
const std::vector<std::array<int,3>>& cells_per_dim_vec) const
{
{
// --- Refined cells ---
for (std::size_t shiftedLevel = 0; shiftedLevel < refined_cell_count_vec.size(); ++shiftedLevel) {

Expand Down Expand Up @@ -3966,7 +3966,7 @@ void CpGrid::populateRefinedCells(std::vector<Dune::cpgrid::EntityVariableBase<c

// Create a pointer to the first element of "refined_cell_to_point" (required as the fourth argement to construct a Geometry<3,3> type object).
int* indices_storage_ptr = refined_cell_to_point_vec[shiftedLevel][cell].data();
refined_cells_vec[shiftedLevel][cell] = cpgrid::Geometry<3,3>(elemLgrGeom.center(), elemLgrGeom.volume(), allLevelCorners, indices_storage_ptr);
refined_cells_vec[shiftedLevel][cell] = cpgrid::Geometry<3,3>(elemLgrGeom.center(), elemLgrGeom.volume(), allLevelCorners.get(), indices_storage_ptr);
} // refined_cells
// Refined face to cell.
refined_cell_to_face_vec[shiftedLevel].makeInverseRelation(refined_face_to_cell_vec[shiftedLevel]);
Expand Down Expand Up @@ -4137,7 +4137,7 @@ void CpGrid::updateCornerHistoryLevels(const std::vector<std::vector<std::array<
currentData()[refinedLevel]->corner_history_[refinedCorner] = preAdaptGrid_corner_history.empty() ? std::array<int,2>{{0, static_cast<int>(corner)}} : preAdaptGrid_corner_history[corner];
}
}

// corner_history_ leaf grid view
for ( int leafCorner = 0; leafCorner < corner_count; ++leafCorner){
currentData().back()->corner_history_.resize(corner_count);
Expand Down Expand Up @@ -4664,5 +4664,3 @@ int CpGrid::replaceLgr1FaceIdxByLgr2FaceIdx(const std::array<int,3>& cells_per_d
}

} // namespace Dune


16 changes: 8 additions & 8 deletions opm/grid/cpgrid/CpGridData.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -535,7 +535,7 @@ struct CellGeometryHandle
buffer.read(pos[i]);

buffer.read(vol);
scatterCont_[t] = Geom(pos, vol, pointGeom_, cell2Points_[t.index()].data());
scatterCont_[t] = Geom(pos, vol, pointGeom_.get(), cell2Points_[t.index()].data());
double isAquifer;
buffer.read(isAquifer);
if (isAquifer == 1.0)
Expand Down Expand Up @@ -1577,8 +1577,8 @@ void CpGridData::distributeGlobalGrid(CpGrid& grid,

// Compute partition type for points
computePointPartitionType();
computeCommunicationInterfaces(noExistingPoints);

computeCommunicationInterfaces(noExistingPoints);
#else // #if HAVE_MPI
static_cast<void>(grid);
static_cast<void>(view_data);
Expand Down Expand Up @@ -2065,9 +2065,9 @@ int CpGridData::sharedFaceTag(const std::vector<std::array<int,3>>& startIJK_2Pa
assert(endIJK_2Patches.size() == 2);

int faceTag = -1; // 0 represents I_FACE, 1 J_FACE, and 2 K_FACE. Use -1 for no sharing face case.

if (patchesShareFace(startIJK_2Patches, endIJK_2Patches)) {

const auto& detectSharing = [](const std::vector<int>& faceIdxs, const std::vector<int>& otherFaceIdxs){
bool faceIsShared = false;
for (const auto& face : faceIdxs) {
Expand All @@ -2080,7 +2080,7 @@ int CpGridData::sharedFaceTag(const std::vector<std::array<int,3>>& startIJK_2Pa
}
return faceIsShared; // should be false here
};

const auto& [iFalse, iTrue, jFalse, jTrue, kFalse, kTrue] = this->getBoundaryPatchFaces(startIJK_2Patches[0], endIJK_2Patches[0]);
const auto& [iFalseOther, iTrueOther, jFalseOther, jTrueOther, kFalseOther, kTrueOther] =
this->getBoundaryPatchFaces(startIJK_2Patches[1], endIJK_2Patches[1]);
Expand Down Expand Up @@ -2306,7 +2306,7 @@ Geometry<3,3> CpGridData::cellifyPatch(const std::array<int,3>& startIJK, const
const int* cellifiedPatch_indices_storage_ptr = &allcorners_cellifiedPatch[0];
// Construct (and return) the Geometry<3,3> of the 'cellified patch'.
return Geometry<3,3>(cellifiedPatch_center, cellifiedPatch_volume,
cellifiedPatch_geometry.geomVector(std::integral_constant<int,3>()),
cellifiedPatch_geometry.geomVector(std::integral_constant<int,3>()).get(),
cellifiedPatch_indices_storage_ptr);
}
}
Expand Down Expand Up @@ -2699,7 +2699,7 @@ bool CpGridData::preAdapt()
if (local_empty)
mark_.resize(size(0));
}

// Detect the maximum mark across processes, and rewrite
// the local entry in mark_, i.e.,
// mark_[ element.index() ] = max{ local marks in processes where this element belongs to}.
Expand Down
11 changes: 5 additions & 6 deletions opm/grid/cpgrid/Entity.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ class Entity : public EntityRep<codim>
}

/// @brief Return the geometry of the entity (does not depend on its orientation).
const Geometry& geometry() const;
Geometry geometry() const;

/// @brief Return the level of the entity in the grid hierarchy. Level = 0 represents the coarsest grid.
int level() const;
Expand Down Expand Up @@ -387,7 +387,7 @@ unsigned int Entity<codim>::subEntities ( const unsigned int cc ) const
}

template <int codim>
const typename Entity<codim>::Geometry& Entity<codim>::geometry() const
typename Entity<codim>::Geometry Entity<codim>::geometry() const
{
return pgrid_->geomVector<codim>()[*this];
}
Expand Down Expand Up @@ -539,10 +539,9 @@ Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() cons
const auto& auxArr = pgrid_ -> getReferenceRefinedCorners(idx_in_parent_cell, cells_per_dim);
FieldVector<double, 3> corners_in_father_reference_elem_temp[8] =
{ auxArr[0], auxArr[1], auxArr[2], auxArr[3], auxArr[4], auxArr[5], auxArr[6], auxArr[7]};
auto in_father_reference_elem_corners = std::make_shared<EntityVariable<cpgrid::Geometry<0, 3>, 3>>();
EntityVariableBase<cpgrid::Geometry<0, 3>>& mutable_in_father_reference_elem_corners = *in_father_reference_elem_corners;
EntityVariable<cpgrid::Geometry<0, 3>, 3> in_father_reference_elem_corners;
// Assign the corners. Make use of the fact that pointers behave like iterators.
mutable_in_father_reference_elem_corners.assign(corners_in_father_reference_elem_temp,
in_father_reference_elem_corners.assign(corners_in_father_reference_elem_temp,
corners_in_father_reference_elem_temp + 8);
// Compute the center of the 'local-entity'.
Dune::FieldVector<double, 3> center_in_father_reference_elem = {0., 0.,0.};
Expand All @@ -556,7 +555,7 @@ Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() cons
double volume_in_father_reference_elem = double(1)/(cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2]);
// Construct (and return) the Geometry<3,3> of 'child-cell in the reference element of its father (unit cube)'.
return Dune::cpgrid::Geometry<3,3>(center_in_father_reference_elem, volume_in_father_reference_elem,
in_father_reference_elem_corners, in_father_reference_elem_corner_indices.data());
std::move(in_father_reference_elem_corners), in_father_reference_elem_corner_indices.data());
}
else {
OPM_THROW(std::logic_error, "Entity has no father.");
Expand Down
Loading