36#ifndef OPM_GEOMETRY_HEADER
37#define OPM_GEOMETRY_HEADER
42#include <opm/grid/utility/platform_dependent/disable_warnings.h>
44#include <dune/common/version.hh>
45#include <dune/geometry/referenceelements.hh>
46#include <dune/grid/common/geometry.hh>
48#include <dune/geometry/type.hh>
50#include <opm/grid/cpgrid/EntityRep.hpp>
51#include <opm/grid/common/Volumes.hpp>
52#include <opm/grid/utility/platform_dependent/reenable_warnings.h>
54#include <opm/grid/utility/ErrorMacros.hpp>
69 template <
int mydim,
int cdim>
81 static_assert(cdim == 3,
"");
84 enum { dimension = 3 };
86 enum { mydimension = 0};
88 enum { coorddimension = cdim };
90 enum { dimensionworld = 3 };
101 typedef FieldMatrix< ctype, coorddimension, mydimension >
Jacobian;
145 return Dune::GeometryTypes::cube(mydimension);
157 static_cast<void>(cor);
184 JacobianInverseTransposed
212 GlobalCoordinate pos_;
222 static_assert(cdim == 3,
"");
225 enum { dimension = 3 };
227 enum { mydimension = 3 };
229 enum { coorddimension = cdim };
231 enum { dimensionworld = 3 };
242 typedef FieldMatrix< ctype, coorddimension, mydimension >
Jacobian;
250 typedef Dune::Impl::FieldMatrixHelper< double > MatrixHelperType;
263 const int* corner_indices)
264 : pos_(pos), vol_(vol), allcorners_(allcorners.data()), cor_idx_(corner_indices)
266 assert(allcorners_ && corner_indices);
280 : pos_(pos), vol_(vol)
286 : pos_(0.0), vol_(0.0), allcorners_(0), cor_idx_(0)
296 static_assert(mydimension == 3,
"");
297 static_assert(coorddimension == 3,
"");
300 uvw[0] -= local_coord;
302 const int pat[8][3] = { { 0, 0, 0 },
311 for (
int i = 0; i < 8; ++i) {
314 for (
int j = 0; j < 3; ++j) {
315 factor *= uvw[pat[i][j]][j];
317 corner_contrib *= factor;
318 xyz += corner_contrib;
327 static_assert(mydimension == 3,
"");
328 static_assert(coorddimension == 3,
"");
331 const ctype epsilon = 1e-12;
332 auto refElement = Dune::ReferenceElements<ctype, 3>::cube();
340 MatrixHelperType::template xTRightInvA<3, 3>(JT, z, dx );
342 }
while (dx.two_norm2() > epsilon*epsilon);
353 return MatrixHelperType::template sqrtDetAAT<3, 3>(Jt);
360 return Dune::GeometryTypes::cube(mydimension);
373 assert(allcorners_ && cor_idx_);
374 return allcorners_[cor_idx_[cor]].center();
383 void set_volume(ctype
volume) {
397 const JacobianTransposed
400 static_assert(mydimension == 3,
"");
401 static_assert(coorddimension == 3,
"");
405 uvw[0] -= local_coord;
407 const int pat[8][3] = { { 0, 0, 0 },
416 for (
int i = 0; i < 8; ++i) {
417 for (
int deriv = 0; deriv < 3; ++deriv) {
420 for (
int j = 0; j < 3; ++j) {
421 factor *= (j != deriv) ? uvw[pat[i][j]][j]
422 : (pat[i][j] == 0 ? -1.0 : 1.0);
425 corner_contrib *= factor;
426 Jt[deriv] += corner_contrib;
433 const JacobianInverseTransposed
445 return jacobianTransposed(local_coord).transposed();
452 return jacobianInverseTransposed(local_coord).transposed();
477 std::vector<Geometry<3, cdim>>
refine(
const std::array<int, 3>& cells_per_dim,
479 std::vector<std::array<int, 8>>& indices_storage)
497 const int face_corner_indices[6][4] = {
512 const int tetra_edge_indices[6][4][2] = {
513 {{0, 1}, {0, 2}, {1, 3}, {2, 3}},
514 {{0, 1}, {0, 4}, {1, 5}, {4, 5}},
515 {{0, 2}, {0, 4}, {2, 6}, {4, 6}},
516 {{1, 3}, {1, 5}, {3, 7}, {5, 7}},
517 {{2, 3}, {2, 6}, {3, 7}, {6, 7}},
518 {{4, 5}, {4, 6}, {5, 7}, {6, 7}},
522 std::vector<Geometry<3, cdim>> result;
523 result.reserve(cells_per_dim[0] * cells_per_dim[1] * cells_per_dim[2]);
525 auto pcs = corner_storage.begin();
526 auto pis = indices_storage.begin();
529 for (
int k = 0; k < cells_per_dim[2]; k++) {
533 refined_center[2] = (parent_center[2] + k) / cells_per_dim[2];
534 for (
int h = 0; h < 8; h++) {
535 refined_corners[h][2] = (parent_corners[h][2] + k) / cells_per_dim[2];
537 for (
int j = 0; j < cells_per_dim[1]; j++) {
538 refined_center[1] = (parent_center[1] + j) / cells_per_dim[1];
539 for (
int h = 0; h < 8; h++) {
540 refined_corners[h][1] = (parent_corners[h][1] + j) / cells_per_dim[1];
542 for (
int i = 0; i < cells_per_dim[0]; i++) {
543 refined_center[0] = (parent_center[0] + i) / cells_per_dim[0];
544 for (
int h = 0; h < 8; h++) {
545 refined_corners[h][0] = (parent_corners[h][0] + i) / cells_per_dim[0];
548 auto& global_refined_corners = *pcs++;
549 global_refined_corners.reserve(8);
550 for (
const auto& corner : refined_corners) {
551 global_refined_corners.push_back(
Geometry<0, 3>(this->global(corner)));
557 auto& indices = *pis++;
558 indices = {0, 1, 2, 3, 4, 5, 6, 7};
562 this->global(refined_center));
565 const auto& hex_corners = global_refined_corners.data();
567 for (
int f = 0; f < 6; f++) {
568 face_centers[f] = hex_corners[face_corner_indices[f][0]].center();
569 face_centers[f] += hex_corners[face_corner_indices[f][1]].center();
570 face_centers[f] += hex_corners[face_corner_indices[f][2]].center();
571 face_centers[f] += hex_corners[face_corner_indices[f][3]].center();
572 face_centers[f] /= 4;
577 for (
int f = 0; f < 6; f++) {
578 for (
int e = 0; e < 4; e++) {
580 = {hex_corners[tetra_edge_indices[f][e][0]].center(),
581 hex_corners[tetra_edge_indices[f][e][1]].center(),
583 global_refined_center};
590 global_refined_center,
volume, global_refined_corners, indices.data()));
596 if (std::fabs(total_volume - this->
volume())
599 for (
auto& r : result) {
600 r.set_volume(r.volume() * correction);
608 GlobalCoordinate pos_;
623 static_assert(cdim == 3,
"");
626 enum { dimension = 3 };
628 enum { mydimension = 2 };
630 enum { coorddimension = cdim };
632 enum { dimensionworld = 3 };
643 typedef FieldMatrix< ctype, coorddimension, mydimension >
Jacobian;
656 : pos_(pos), vol_(vol)
662 : pos_(0.0), vol_(0.0)
669 OPM_THROW(std::runtime_error,
"Geometry::global() meaningless on singular geometry.");
675 OPM_THROW(std::runtime_error,
"Geometry::local() meaningless on singular geometry.");
688 return Dune::GeometryTypes::none(mydimension);
720 const FieldMatrix<ctype, mydimension, coorddimension>&
723 OPM_THROW(std::runtime_error,
"Meaningless to call jacobianTransposed() on singular geometries.");
727 const FieldMatrix<ctype, coorddimension, mydimension>&
730 OPM_THROW(std::runtime_error,
"Meaningless to call jacobianInverseTransposed() on singular geometries.");
737 return jacobianTransposed({}).transposed();
744 return jacobianInverseTransposed({}).transposed();
754 GlobalCoordinate pos_;
764 template<
int mydim,
int cdim >
765 auto referenceElement(
const cpgrid::Geometry<mydim,cdim>& geo) ->
decltype(referenceElement<double,mydim>(geo.type()))
767 return referenceElement<double,mydim>(geo.type());
A class design to hold a variable with a value for each entity of the given codimension,...
Definition: EntityRep.hpp:264
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:98
bool affine() const
The mapping implemented by this geometry is constant, therefore affine.
Definition: Geometry.hpp:206
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:101
JacobianInverse jacobianInverse(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:200
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:96
Jacobian jacobian(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:193
double integrationElement(const LocalCoordinate &) const
Returns 1 for the vertex geometry.
Definition: Geometry.hpp:137
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:107
Geometry(const GlobalCoordinate &pos)
Construct from vertex position.
Definition: Geometry.hpp:112
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse
Type of inverse of Jacobian matrix.
Definition: Geometry.hpp:103
JacobianTransposed jacobianTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:176
GeometryType type() const
Using the cube type for vertices.
Definition: Geometry.hpp:143
ctype volume() const
Volume of vertex is arbitrarily set to 1.
Definition: Geometry.hpp:163
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:105
GlobalCoordinate corner(int cor) const
Returns the single corner: the vertex itself.
Definition: Geometry.hpp:155
int corners() const
A vertex is defined by a single corner.
Definition: Geometry.hpp:149
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:169
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:118
const GlobalCoordinate & global(const LocalCoordinate &) const
Returns the position of the vertex.
Definition: Geometry.hpp:124
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:185
LocalCoordinate local(const GlobalCoordinate &) const
Meaningless for the vertex geometry.
Definition: Geometry.hpp:130
double ctype
Coordinate element type.
Definition: Geometry.hpp:93
JacobianInverse jacobianInverse(const LocalCoordinate &) const
The inverse of the jacobian.
Definition: Geometry.hpp:742
const FieldMatrix< ctype, mydimension, coorddimension > & jacobianTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:721
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:714
int corners() const
The number of corners of this convex polytope.
Definition: Geometry.hpp:693
LocalCoordinate local(const GlobalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:673
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:643
bool affine() const
Since integrationElement() is constant, returns true.
Definition: Geometry.hpp:748
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:638
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:649
ctype volume() const
Volume (area, actually) of intersection.
Definition: Geometry.hpp:708
const FieldMatrix< ctype, coorddimension, mydimension > & jacobianInverseTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:728
GeometryType type() const
We use the singular type (None) for intersections.
Definition: Geometry.hpp:686
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:647
double integrationElement(const LocalCoordinate &) const
For the singular geometry, we return a constant integration element equal to the volume.
Definition: Geometry.hpp:680
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:661
const GlobalCoordinate & global(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:667
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse
Type of Jacobian matrix.
Definition: Geometry.hpp:645
Geometry(const GlobalCoordinate &pos, ctype vol)
Construct from centroid and volume (1- and 0-moments).
Definition: Geometry.hpp:654
double ctype
Coordinate element type.
Definition: Geometry.hpp:635
Jacobian jacobian(const LocalCoordinate &) const
The jacobian.
Definition: Geometry.hpp:735
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:640
GlobalCoordinate corner(int) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:699
Specialization for 3-dimensional geometries, i.e. cells.
Definition: Geometry.hpp:221
const JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local_coord) const
Inverse of Jacobian transposed.
Definition: Geometry.hpp:434
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:248
bool affine() const
The mapping implemented by this geometry is not generally affine.
Definition: Geometry.hpp:456
double ctype
Coordinate element type.
Definition: Geometry.hpp:234
Geometry(const GlobalCoordinate &pos, ctype vol, const EntityVariable< cpgrid::Geometry< 0, 3 >, 3 > &allcorners, const int *corner_indices)
Construct from centroid, volume (1- and 0-moments) and corners.
Definition: Geometry.hpp:260
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:239
GlobalCoordinate corner(int cor) const
The 8 corners of the hexahedral base cell.
Definition: Geometry.hpp:371
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:246
GeometryType type() const
Using the cube type for all entities now (cells and vertices), but we use the singular type for inter...
Definition: Geometry.hpp:358
double integrationElement(const LocalCoordinate &local_coord) const
Equal to \sqrt{\det{J^T J}} where J is the Jacobian.
Definition: Geometry.hpp:350
const JacobianTransposed jacobianTransposed(const LocalCoordinate &local_coord) const
Jacobian transposed.
Definition: Geometry.hpp:398
GlobalCoordinate global(const LocalCoordinate &local_coord) const
Provide a trilinear mapping.
Definition: Geometry.hpp:294
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:285
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse
Type of inverse of Jacobian matrix.
Definition: Geometry.hpp:244
LocalCoordinate local(const GlobalCoordinate &y) const
Mapping from the cell to the reference domain.
Definition: Geometry.hpp:325
int corners() const
The number of corners of this convex polytope.
Definition: Geometry.hpp:365
ctype volume() const
Cell volume.
Definition: Geometry.hpp:378
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:388
Jacobian jacobian(const LocalCoordinate &local_coord) const
The jacobian.
Definition: Geometry.hpp:443
Geometry(const GlobalCoordinate &pos, ctype vol)
Construct from centroid and volume (1- and 0-moments).
Definition: Geometry.hpp:278
std::vector< Geometry< 3, cdim > > refine(const std::array< int, 3 > &cells_per_dim, std::vector< EntityVariable< Geometry< 0, 3 >, 3 > > &corner_storage, std::vector< std::array< int, 8 > > &indices_storage)
Refine a single cell with regular intervals.
Definition: Geometry.hpp:477
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:237
JacobianInverse jacobianInverse(const LocalCoordinate &local_coord) const
The inverse of the jacobian.
Definition: Geometry.hpp:450
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:242
This class encapsulates geometry for both vertices, intersections and cells.
Definition: Geometry.hpp:71
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
T volume(const Point< T, 3 > *c)
Computes the volume of a 3D simplex (embedded i 3D space).
Definition: Volumes.hpp:137
T simplex_volume(const Point< T, Dim > *a)
Computes the volume of a simplex consisting of (Dim+1) vertices embedded in Euclidean space of dimens...
Definition: Volumes.hpp:104