Skip to content

Return geometry by value (alternative 4) #1055

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

Open
SoilRos wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Return geometry by value (alternative 4) #1055

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
81162fb
Use `Dune::AxisAlignedCubeGeometry` for entity local geometries
SoilRos Jul 30, 2026
649ab52
Return geometry objects by value
SoilRos Jun 13, 2025
fe28d09
Refactor Geometry to use non-owning pointers for corner data
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 @@ -1242,28 +1242,28 @@ namespace Dune
/// \brief Get the Position of a vertex.
/// \param vertex The index identifying the vertex.
/// \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 face The index identifying the face.
double faceArea(int face) const;

/// \brief Get the coordinates of the center of a face.
/// \param face 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 face 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
8 changes: 4 additions & 4 deletions opm/grid/cpgrid/CpGrid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1563,7 +1563,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_data_->back()->geomVector<3>()[cpgrid::EntityRep<3>(vertex, true)].center();
}
Expand All @@ -1573,12 +1573,12 @@ double CpGrid::faceArea(int face) const
return current_data_->back()->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_data_->back()->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_data_->back()->face_normals_.get(face);
}
Expand All @@ -1588,7 +1588,7 @@ double CpGrid::cellVolume(int cell) const
return current_data_->back()->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_data_->back()->geomVector<0>()[cpgrid::EntityRep<0>(cell, true)].center();
}
Expand Down
2 changes: 1 addition & 1 deletion opm/grid/cpgrid/CpGridData.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -536,7 +536,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
42 changes: 10 additions & 32 deletions opm/grid/cpgrid/Entity.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@
#include <dune/common/version.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/axisalignedcubegeometry.hh>
#include <dune/grid/common/gridenums.hh>

#include "PartitionTypeIndicator.hpp"
Expand Down Expand Up @@ -130,8 +131,8 @@ class Entity : public EntityRep<codim>
template <int cd>
using Codim = typename Impl::CodimTraits<cd>;

typedef cpgrid::Geometry<3-codim,3> Geometry;
typedef Geometry LocalGeometry;
using Geometry = cpgrid::Geometry<3-codim,3>;
using LocalGeometry = Dune::AxisAlignedCubeGeometry<double,3,3>;

typedef cpgrid::IntersectionIterator LeafIntersectionIterator;
typedef cpgrid::IntersectionIterator LevelIntersectionIterator;
Expand Down Expand Up @@ -193,7 +194,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 @@ -284,7 +285,7 @@ class Entity : public EntityRep<codim>
/// Currently, LGR is built via refinement of a block-shaped patch from the coarse grid. So the LocalGeometry
/// of an entity coming from the LGR is one of the refined cells of the unit cube, with suitable amount of cells
/// in each direction.
Dune::cpgrid::Geometry<3,3> geometryInFather() const;
Dune::AxisAlignedCubeGeometry<double,3,3> geometryInFather() const;

/// Returns true if any of my intersections are on the boundary.
/// Implementation note:
Expand Down Expand Up @@ -422,7 +423,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 @@ -559,41 +560,18 @@ int Dune::cpgrid::Entity<codim>::getIdxInParentCell() const


template<int codim>
Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() const
Dune::AxisAlignedCubeGeometry<double,3,3> Dune::cpgrid::Entity<codim>::geometryInFather() const
{
if (!(this->hasFather())){
OPM_THROW(std::logic_error, "Entity has no father.");
}

// Indices of corners in entity's geometry in father reference element.
static constexpr std::array<int,8> in_father_reference_elem_corner_indices = {0,1,2,3,4,5,6,7};
// 'static': The returned object Geometry<3,3> stores a pointer to in_father_reference_elem_corner_indices. Therefore,
// this variable is declared static to prolongate its lifetime beyond this function (static storage duration).

auto idx_in_parent_cell = pgrid_ -> cell_to_idxInParentCell_[this->index()];
if (idx_in_parent_cell !=-1) {
const auto& cells_per_dim = (*(pgrid_ -> level_data_ptr_))[this->level()] -> cells_per_dim_;
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;
// 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,
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.};
for (int corn = 0; corn < 8; ++corn) {
for (int c = 0; c < 3; ++c)
{
center_in_father_reference_elem[c] += corners_in_father_reference_elem_temp[corn][c]/8.;
}
}
// Compute the volume of the 'local-entity'.
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());
const auto& corners = pgrid_ -> getReferenceRefinedCorners(idx_in_parent_cell, cells_per_dim);
// Create new geometry from lower left and upper right corners of the refined cell in the reference element of the father cell.
return {corners[0], corners[7]};
}
else {
OPM_THROW(std::logic_error, "Entity has no father.");
Expand Down
34 changes: 19 additions & 15 deletions opm/grid/cpgrid/Geometry.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,7 @@ namespace Dune
}

/// Returns the position of the vertex.
const GlobalCoordinate& global(const LocalCoordinate&) const
GlobalCoordinate global(const LocalCoordinate&) const
{
return pos_;
}
Expand Down Expand Up @@ -173,7 +173,7 @@ namespace Dune
}

/// Returns the centroid of the geometry.
const GlobalCoordinate& center() const
GlobalCoordinate center() const
{
return pos_;
}
Expand Down Expand Up @@ -273,7 +273,7 @@ namespace Dune
}

/// This method is meaningless for singular geometries.
const GlobalCoordinate& global(const LocalCoordinate&) const
GlobalCoordinate global(const LocalCoordinate&) const
{
OPM_THROW(std::runtime_error, "Geometry::global() meaningless on singular geometry.");
}
Expand Down Expand Up @@ -320,20 +320,20 @@ namespace Dune
}

/// Returns the centroid of the geometry.
const GlobalCoordinate& center() const
GlobalCoordinate center() const
{
return pos_;
}

/// This method is meaningless for singular geometries.
const FieldMatrix<ctype, mydimension, coorddimension>&
FieldMatrix<ctype, mydimension, coorddimension>
jacobianTransposed(const LocalCoordinate& /* local */) const
{
OPM_THROW(std::runtime_error, "Meaningless to call jacobianTransposed() on singular geometries.");
}

/// This method is meaningless for singular geometries.
const FieldMatrix<ctype, coorddimension, mydimension>&
FieldMatrix<ctype, coorddimension, mydimension>
jacobianInverseTransposed(const LocalCoordinate& /*local*/) const
{
OPM_THROW(std::runtime_error, "Meaningless to call jacobianInverseTransposed() on singular geometries.");
Expand Down Expand Up @@ -406,6 +406,10 @@ namespace Dune

/// @brief Construct from center, volume (1- and 0-moments) and
/// corners.
/// @warning This constructor does not own the corners or indices,
/// thus, the pointers must remain valid for the lifetime of
/// the Geometry object.
///
/// @param pos the centroid of the entity
/// @param vol the volume(area) of the entity
/// @param allcorners_ptr pointer of all corner positions in the grid
Expand All @@ -414,17 +418,17 @@ namespace Dune
/// by (kji), i.e. i running fastest.
Geometry(const GlobalCoordinate& pos,
ctype vol,
std::shared_ptr<const EntityVariable<cpgrid::Geometry<0, 3>, 3>> allcorners_ptr,
EntityVariable<cpgrid::Geometry<0, 3>, 3> const * allcorners_ptr,
const int* corner_indices)
: pos_(pos), vol_(vol),
allcorners_(allcorners_ptr), cor_idx_(corner_indices)
{
assert(allcorners_ && corner_indices);
assert(allcorners_ptr && corner_indices);
}

/// Default constructor, giving a non-valid geometry.
Geometry()
: pos_(0.0), vol_(0.0), allcorners_(0), cor_idx_(0)
: pos_(0.0), vol_(0.0), allcorners_(nullptr), cor_idx_(nullptr)
{
}

Expand Down Expand Up @@ -527,7 +531,7 @@ namespace Dune
}

/// Returns the centroid of the geometry.
const GlobalCoordinate& center() const
GlobalCoordinate center() const
{
return pos_;
}
Expand All @@ -538,7 +542,7 @@ namespace Dune
/// and {u_i} are the reference coordinates.
/// g = g(u) = (g_1(u), g_2(u), g_3(u)), u=(u_1,u_2,u_3)
/// g = map from (local) reference domain to global cell
const JacobianTransposed
JacobianTransposed
jacobianTransposed(const LocalCoordinate& local_coord) const
{
static_assert(mydimension == 3, "");
Expand Down Expand Up @@ -574,7 +578,7 @@ namespace Dune
}

/// @brief Inverse of Jacobian transposed. \see jacobianTransposed().
const JacobianInverseTransposed
JacobianInverseTransposed
jacobianInverseTransposed(const LocalCoordinate& local_coord) const
{
JacobianInverseTransposed Jti = jacobianTransposed(local_coord);
Expand Down Expand Up @@ -1023,7 +1027,7 @@ namespace Dune
refined_cells[refined_cell_idx] =
Geometry<3,cdim>(refined_cell_center,
refined_cell_volume,
all_geom.geomVector(std::integral_constant<int,3>()),
all_geom.geomVector(std::integral_constant<int,3>()).get(),
indices_storage_ptr);
} // end i-for-loop
} // end j-for-loop
Expand All @@ -1044,7 +1048,7 @@ namespace Dune
private:
GlobalCoordinate pos_;
double vol_;
std::shared_ptr<const EntityVariable<Geometry<0, 3>,3>> allcorners_; // For dimension 3 only
const EntityVariable<Geometry<0, 3>,3>* allcorners_; // store non-owning pointer of the 8 corners
const int* cor_idx_; // For dimension 3 only

/// @brief
Expand All @@ -1058,7 +1062,7 @@ namespace Dune
/// @param [out] refined_face_tag I_FACE, J_FACE, K_FACE
/// @param [out] refined_face_index Face index of a refined cell 'lmn' generated with "refine()".
/// @param [out] refined_face_to_point Four corner indices of the corners of the refined face 'lmn'.
/// Vertex order
/// Vertex order
/// for I_FACE: jk, (j+1)k, (j+1)(k+1), j(k+1)
/// for J_FACE: (i+1)k, ik, i(k+1), (i+1)(k+1)
/// for K_FACE: ij, (i+1)j, (i+1)(j+1), (i+1)(j+1)
Expand Down
15 changes: 7 additions & 8 deletions opm/grid/cpgrid/GridHelpers.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,9 +147,9 @@ beginFaceCentroids(const Dune::CpGrid& grid)
return FaceCentroidTraits<Dune::CpGrid>::IteratorType(grid, 0);
}

const double* cellCentroid(const Dune::CpGrid& grid, int cell_index)
Vector cellCentroid(const Dune::CpGrid& grid, int cell_index)
{
return &(grid.cellCentroid(cell_index)[0]);
return grid.cellCentroid(cell_index);
}

double cellVolume(const Dune::CpGrid& grid, int cell_index)
Expand All @@ -167,8 +167,7 @@ CellVolumeIterator endCellVolumes(const Dune::CpGrid& grid)
return CellVolumeIterator(grid, numCells(grid));
}

const FaceCentroidTraits<Dune::CpGrid>::ValueType&
faceCentroid(const Dune::CpGrid& grid, int face_index)
Vector faceCentroid(const Dune::CpGrid& grid, int face_index)
{
return grid.faceCentroid(face_index);
}
Expand All @@ -190,14 +189,14 @@ face2Vertices(const Dune::CpGrid& grid)
return Dune::cpgrid::FaceVerticesContainerProxy(&grid);
}

const double* vertexCoordinates(const Dune::CpGrid& grid, int index)
Vector vertexCoordinates(const Dune::CpGrid& grid, int index)
{
return &(grid.vertexPosition(index)[0]);
return grid.vertexPosition(index);
}

const double* faceNormal(const Dune::CpGrid& grid, int face_index)
Vector faceNormal(const Dune::CpGrid& grid, int face_index)
{
return &(grid.faceNormal(face_index)[0]);
return grid.faceNormal(face_index);
}

double faceArea(const Dune::CpGrid& grid, int face_index)
Expand Down
Loading