My Project
Loading...
Searching...
No Matches
Geometry.hpp
1//===========================================================================
2//
3// File: Geometry.hpp
4//
5// Created: Fri May 29 23:29:24 2009
6//
7// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8// B�rd Skaflestad <bard.skaflestad@sintef.no>
9// Antonella Ritorto <antonella.ritorto@opm-op.com>
10//
11// $Date$
12//
13// $Revision$
14//
15//===========================================================================
16
17/*
18 Copyright 2009, 2010, 2011 SINTEF ICT, Applied Mathematics.
19 Copyright 2009, 2010, 2011, 2022 Equinor ASA.
20
21 This file is part of the Open Porous Media project (OPM).
22
23 OPM is free software: you can redistribute it and/or modify
24 it under the terms of the GNU General Public License as published by
25 the Free Software Foundation, either version 3 of the License, or
26 (at your option) any later version.
27
28 OPM is distributed in the hope that it will be useful,
29 but WITHOUT ANY WARRANTY; without even the implied warranty of
30 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 GNU General Public License for more details.
32
33 You should have received a copy of the GNU General Public License
34 along with OPM. If not, see <http://www.gnu.org/licenses/>.
35*/
36
37#ifndef OPM_GEOMETRY_HEADER
38#define OPM_GEOMETRY_HEADER
39
40#include <cmath>
41
42// Warning suppression for Dune includes.
43#include <opm/grid/utility/platform_dependent/disable_warnings.h>
44
45#include <dune/common/version.hh>
46#include <dune/geometry/referenceelements.hh>
47#include <dune/grid/common/geometry.hh>
48
49#include <dune/geometry/type.hh>
50
51#include <opm/grid/cpgrid/EntityRep.hpp>
52#include <opm/grid/cpgrid/DefaultGeometryPolicy.hpp>
53#include <opm/grid/cpgrid/OrientedEntityTable.hpp>
54#include <opm/grid/common/Volumes.hpp>
55#include <opm/grid/utility/platform_dependent/reenable_warnings.h>
56#include <opm/grid/utility/SparseTable.hpp>
57
58#include <opm/common/ErrorMacros.hpp>
59
60namespace Dune
61{
62 namespace cpgrid
63 {
64
73 template <int mydim, int cdim>
75 {
76 };
77
78
79
80
82 template <int cdim> // GridImp arg never used
83 class Geometry<0, cdim>
84 {
85 static_assert(cdim == 3, "");
86 public:
88 enum { dimension = 3 };
90 enum { mydimension = 0};
92 enum { coorddimension = cdim };
94 enum { dimensionworld = 3 };
95
97 typedef double ctype;
98
100 typedef FieldVector<ctype, mydimension> LocalCoordinate;
102 typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
103
105 typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
107 typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse;
109 typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
111 typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
112
113
117 : pos_(pos)
118 {
119 }
120
123 : pos_(0.0)
124 {
125 }
126
129 {
130 return pos_;
131 }
132
135 {
136 // return 0 to make the geometry check happy.
137 return LocalCoordinate(0.0);
138 }
139
142 {
143 return volume();
144 }
145
147 GeometryType type() const
148 {
149 return Dune::GeometryTypes::cube(mydimension);
150 }
151
153 int corners() const
154 {
155 return 1;
156 }
157
160 {
161 static_cast<void>(cor);
162 assert(cor == 0);
163 return pos_;
164 }
165
167 ctype volume() const
168 {
169 return 1.0;
170 }
171
174 {
175 return pos_;
176 }
177
179 JacobianTransposed
180 jacobianTransposed(const LocalCoordinate& /* local */) const
181 {
182
183 // Meaningless to call jacobianTransposed() on singular geometries. But we need to make DUNE happy.
184 return {};
185 }
186
188 JacobianInverseTransposed
190 {
191 // Meaningless to call jacobianInverseTransposed() on singular geometries. But we need to make DUNE happy.
192 return {};
193 }
194
196 Jacobian
197 jacobian(const LocalCoordinate& /*local*/) const
198 {
199 return {};
200 }
201
203 JacobianInverse
204 jacobianInverse(const LocalCoordinate& /*local*/) const
205 {
206 return {};
207 }
208
210 bool affine() const
211 {
212 return true;
213 }
214
215 private:
216 GlobalCoordinate pos_;
217 }; // class Geometry<0,cdim>
218
219
220
221
224 template <int cdim> // GridImp arg never used
225 class Geometry<2, cdim>
226 {
227 static_assert(cdim == 3, "");
228 public:
230 enum { dimension = 3 };
232 enum { mydimension = 2 };
234 enum { coorddimension = cdim };
236 enum { dimensionworld = 3 };
237
239 typedef double ctype;
240
242 typedef FieldVector<ctype, mydimension> LocalCoordinate;
244 typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
245
247 typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
249 typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse;
251 typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
253 typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
254
259 ctype vol)
260 : pos_(pos), vol_(vol)
261 {
262 }
263
266 : pos_(0.0), vol_(0.0)
267 {
268 }
269
272 {
273 OPM_THROW(std::runtime_error, "Geometry::global() meaningless on singular geometry.");
274 }
275
278 {
279 OPM_THROW(std::runtime_error, "Geometry::local() meaningless on singular geometry.");
280 }
281
285 {
286 return vol_;
287 }
288
290 GeometryType type() const
291 {
292 return Dune::GeometryTypes::none(mydimension);
293 }
294
297 int corners() const
298 {
299 return 0;
300 }
301
303 GlobalCoordinate corner(int /* cor */) const
304 {
305 // Meaningless call to cpgrid::Geometry::corner(int):
306 //"singular geometry has no corners.
307 // But the DUNE tests assume at least one corner.
308 return GlobalCoordinate( 0.0 );
309 }
310
312 ctype volume() const
313 {
314 return vol_;
315 }
316
319 {
320 return pos_;
321 }
322
324 const FieldMatrix<ctype, mydimension, coorddimension>&
325 jacobianTransposed(const LocalCoordinate& /* local */) const
326 {
327 OPM_THROW(std::runtime_error, "Meaningless to call jacobianTransposed() on singular geometries.");
328 }
329
331 const FieldMatrix<ctype, coorddimension, mydimension>&
333 {
334 OPM_THROW(std::runtime_error, "Meaningless to call jacobianInverseTransposed() on singular geometries.");
335 }
336
338 Jacobian
339 jacobian(const LocalCoordinate& /*local*/) const
340 {
341 return jacobianTransposed({}).transposed();
342 }
343
345 JacobianInverse
346 jacobianInverse(const LocalCoordinate& /*local*/) const
347 {
348 return jacobianInverseTransposed({}).transposed();
349 }
350
352 bool affine() const
353 {
354 return true;
355 }
356
357 private:
358 GlobalCoordinate pos_;
359 ctype vol_;
360 };
361
362
363
364
365
367 template <int cdim>
368 class Geometry<3, cdim>
369 {
370 static_assert(cdim == 3, "");
371 public:
373 enum { dimension = 3 };
375 enum { mydimension = 3 };
377 enum { coorddimension = cdim };
379 enum { dimensionworld = 3 };
380
382 typedef double ctype;
383
385 typedef FieldVector<ctype, mydimension> LocalCoordinate;
387 typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
388
390 typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
392 typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse;
394 typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
396 typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
397
398 typedef Dune::Impl::FieldMatrixHelper< double > MatrixHelperType;
399
409 ctype vol,
410 std::shared_ptr<const EntityVariable<cpgrid::Geometry<0, 3>, 3>> allcorners_ptr,
411 const int* corner_indices)
412 : pos_(pos), vol_(vol),
413 allcorners_(allcorners_ptr), cor_idx_(corner_indices)
414 {
415 assert(allcorners_ && corner_indices);
416 }
417
428 ctype vol)
429 : pos_(pos), vol_(vol)
430 {
431 }
432
435 : pos_(0.0), vol_(0.0), allcorners_(0), cor_idx_(0)
436 {
437 }
438
444 GlobalCoordinate global(const LocalCoordinate& local_coord) const
445 {
446 static_assert(mydimension == 3, "");
447 static_assert(coorddimension == 3, "");
448 // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
449 LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
450 uvw[0] -= local_coord;
451 // Access pattern for uvw matching ordering of corners.
452 const int pat[8][3] = { { 0, 0, 0 },
453 { 1, 0, 0 },
454 { 0, 1, 0 },
455 { 1, 1, 0 },
456 { 0, 0, 1 },
457 { 1, 0, 1 },
458 { 0, 1, 1 },
459 { 1, 1, 1 } };
460 GlobalCoordinate xyz(0.0);
461 for (int i = 0; i < 8; ++i) {
462 GlobalCoordinate corner_contrib = corner(i);
463 double factor = 1.0;
464 for (int j = 0; j < 3; ++j) {
465 factor *= uvw[pat[i][j]][j];
466 }
467 corner_contrib *= factor;
468 xyz += corner_contrib;
469 }
470 return xyz;
471 }
472
476 {
477 static_assert(mydimension == 3, "");
478 static_assert(coorddimension == 3, "");
479 // This code is modified from dune/grid/genericgeometry/mapping.hh
480 // \todo: Implement direct computation.
481 const ctype epsilon = 1e-12;
482 auto refElement = Dune::ReferenceElements<ctype, 3>::cube();
483 LocalCoordinate x = refElement.position(0,0);
485 do {
486 // DF^n dx^n = F^n, x^{n+1} -= dx^n
487 JacobianTransposed JT = jacobianTransposed(x);
488 GlobalCoordinate z = global(x);
489 z -= y;
490 MatrixHelperType::template xTRightInvA<3, 3>(JT, z, dx );
491 x -= dx;
492 } while (dx.two_norm2() > epsilon*epsilon);
493 return x;
494 }
495
500 double integrationElement(const LocalCoordinate& local_coord) const
501 {
502 JacobianTransposed Jt = jacobianTransposed(local_coord);
503 return MatrixHelperType::template sqrtDetAAT<3, 3>(Jt);
504 }
505
508 GeometryType type() const
509 {
510 return Dune::GeometryTypes::cube(mydimension);
511 }
512
515 int corners() const
516 {
517 return 8;
518 }
519
522 {
523 assert(allcorners_ && cor_idx_);
524 return (allcorners_->data())[cor_idx_[cor]].center();
525 }
526
528 ctype volume() const
529 {
530 return vol_;
531 }
532
533 void set_volume(ctype volume) {
534 vol_ = volume;
535 }
536
539 {
540 return pos_;
541 }
542
549 const JacobianTransposed
550 jacobianTransposed(const LocalCoordinate& local_coord) const
551 {
552 static_assert(mydimension == 3, "");
553 static_assert(coorddimension == 3, "");
554
555 // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
556 LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
557 uvw[0] -= local_coord;
558 // Access pattern for uvw matching ordering of corners.
559 const int pat[8][3] = { { 0, 0, 0 },
560 { 1, 0, 0 },
561 { 0, 1, 0 },
562 { 1, 1, 0 },
563 { 0, 0, 1 },
564 { 1, 0, 1 },
565 { 0, 1, 1 },
566 { 1, 1, 1 } };
567 JacobianTransposed Jt(0.0);
568 for (int i = 0; i < 8; ++i) {
569 for (int deriv = 0; deriv < 3; ++deriv) {
570 // This part contributing to dg/du_{deriv}
571 double factor = 1.0;
572 for (int j = 0; j < 3; ++j) {
573 factor *= (j != deriv) ? uvw[pat[i][j]][j]
574 : (pat[i][j] == 0 ? -1.0 : 1.0);
575 }
576 GlobalCoordinate corner_contrib = corner(i);
577 corner_contrib *= factor;
578 Jt[deriv] += corner_contrib; // using FieldMatrix row access.
579 }
580 }
581 return Jt;
582 }
583
585 const JacobianInverseTransposed
587 {
588 JacobianInverseTransposed Jti = jacobianTransposed(local_coord);
589 Jti.invert();
590 return Jti;
591 }
592
594 Jacobian
595 jacobian(const LocalCoordinate& local_coord) const
596 {
597 return jacobianTransposed(local_coord).transposed();
598 }
599
601 JacobianInverse
602 jacobianInverse(const LocalCoordinate& local_coord) const
603 {
604 return jacobianInverseTransposed(local_coord).transposed();
605 }
606
608 bool affine() const
609 {
610 return false;
611 }
612
629 typedef Dune::FieldVector<double,3> PointType;
630 void refine(const std::array<int,3>& cells_per_dim,
631 DefaultGeometryPolicy& all_geom,
632 std::vector<std::array<int,8>>& refined_cell_to_point,
633 cpgrid::OrientedEntityTable<0,1>& refined_cell_to_face,
634 Opm::SparseTable<int>& refined_face_to_point,
635 cpgrid::OrientedEntityTable<1,0>& refined_face_to_cell,
637 cpgrid::SignedEntityVariable<PointType, 1>& refined_face_normals) const
638 {
640 *(all_geom.geomVector(std::integral_constant<int,3>()));
642 *(all_geom.geomVector(std::integral_constant<int,1>()));
644 *(all_geom.geomVector(std::integral_constant<int,0>()));
645 EntityVariableBase<enum face_tag>& mutable_face_tags = refined_face_tags;
646 EntityVariableBase<PointType>& mutable_face_normals = refined_face_normals;
648 // The strategy is to compute the local refined corners
649 // of the unit/reference cube, and then apply the map global().
650 // Determine the size of the vector containing all the corners
651 // of all the global refined cells (children cells).
652 refined_corners.resize((cells_per_dim[0] + 1) *(cells_per_dim[1] + 1) * (cells_per_dim[2] + 1));
653 // The nummbering starts at the botton, so k=0 (z-axis), and j=0 (y-axis), i=0 (x-axis).
654 // Then, increasing k ('going up'), followed by increasing i ('going right->'),
655 // and finally, increasing j ('going back'). This order criteria for corners
656 // 'Up [increasing k]- Right [incresing i]- Back [increasing j]'
657 // is consistant with cpgrid numbering.
658 for (int j = 0; j < cells_per_dim[1] + 1; ++j) {
659 for (int i = 0; i < cells_per_dim[0] + 1; ++i) {
660 for (int k = 0; k < cells_per_dim[2] + 1; ++k) {
661 // Compute the index of each global refined corner associated with 'jik'.
662 int refined_corner_idx =
663 (j*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (i*(cells_per_dim[2]+1)) +k;
664 // Compute the local refined corner of the unit/reference cube associated with 'jik'.
665 const LocalCoordinate& local_refined_corner = {
666 double(i)/cells_per_dim[0], double(j)/cells_per_dim[1], double(k)/cells_per_dim[2] };
667 // Compute the global refined corner 'jik' and add it in its corresponfing entry in "refined_corners".
668 refined_corners[refined_corner_idx] = Geometry<0, 3>(this->global(local_refined_corner));
669 } // end k-for-loop
670 } // end i-for-loop
671 } // end j-for-loop
673 //
675 // We want to populate "refined_faces". The size of "refined_faces" is
676 const int refined_faces_size =
677 (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2]+1)) // 'bottom/top faces'
678 + ((cells_per_dim[0]+1)*cells_per_dim[1]*cells_per_dim[2]) // 'left/right faces'
679 + (cells_per_dim[0]*(cells_per_dim[1]+1)*cells_per_dim[2]); // 'front/back faces'
680 refined_faces.resize(refined_faces_size);
681 refined_face_tags.resize(refined_faces_size);
682 refined_face_normals.resize(refined_faces_size);
683 //
684 // To create a face as a Geometry<2,3> type object we need its CENTROID and its VOLUME(area).
685 // We store the centroids/areas in the following order:
686 // - Bottom-top faces -> 3rd coordinate constant in each face.
687 // - Left-right faces -> 1st coordinate constant in each face.
688 // - Front-back faces -> 2nd coordinate constant in each face.
689 //
690 // REFINED FACE AREAS
691 // To compute the area of each face, we divide it in 4 triangles,
692 // compute the area of those with "area()", where the arguments
693 // are the 3 corners of each triangle. Then, sum them up to get the area
694 // of the global refined face.
695 // For each face, we construct 4 triangles with
696 // (1) centroid of the face,
697 // (2) one of the edges of the face.
698 //
699 // A triangle has 3 edges. Once we choose a face to base a triangle on,
700 // we choose an edge of that face as one of the edges of the triangle.
701 // The other 2 edges are fixed, since the centroid of the face the triangle
702 // is based on is one of its corners. That's why to identify
703 // a triangle we only need two things:
704 // (1) the face it's based on and
705 // (2) one of the four edges of that face.
706 //
707 // For each face, we need
708 // 1. index of the refined face
709 // [needed to access indices of the 4 edges of the face in "refined_face_to_edges"]
710 // 2. centroid of the face (common corner of the 4 triangles based on that face).
711 // [available via "['face'].center()"
712 // 3. container of 4 entries (the 4 edges of the face).
713 // Each entry consists in the 2 indices defining each edge of the face.
714 // [available in "refined_face_to_edges"].
715 //
716 // Populate "mutable_face_tags/normals", "refined_face_to_point/cell",
717 // "refined_faces".
718 //
719 for (int constant_direction = 0; constant_direction < 3; ++constant_direction){
720 // adding %3 and constant_direction, we go through the 3 type of faces.
721 // 0 -> 3rd coordinate constant: l('k') < cells_per_dim[2]+1, m('j') < cells_per_dim[1], n('i') < cells_per_dim[0]
722 // 1 -> 1rt coordinate constant: l('i') < cells_per_dim[0]+1, m('k') < cells_per_dim[2], n('j') < cells_per_dim[1]
723 // 2 -> 2nd coordinate constant: l('j') < cells_per_dim[1]+1, m('i') < cells_per_dim[0], n('k') < cells_per_dim[2]
724 std::array<int,3> cells_per_dim_mixed = {
725 cells_per_dim[(2+constant_direction)%3],
726 cells_per_dim[(1+constant_direction)%3],
727 cells_per_dim[constant_direction % 3] };
728 for (int l = 0; l < cells_per_dim_mixed[0] + 1; ++l) {
729 for (int m = 0; m < cells_per_dim_mixed[1]; ++m) {
730 for (int n = 0; n < cells_per_dim_mixed[2]; ++n) {
731 // Compute the face data.
732 auto [face_type, idx, face4corners,
733 neighboring_cells_of_one_face, local_refined_face_centroid] =
734 getIndicesFace(l, m, n, constant_direction, cells_per_dim);
735 // Add the tag to "refined_face_tags".
736 mutable_face_tags[idx]= face_type;
737 // Add the 4 corners of the face to "refined_face_to_point".
738 refined_face_to_point.appendRow(face4corners.begin(), face4corners.end());
739 // Add the neighboring cells of the face to "refined_face_to_cell".
740 refined_face_to_cell.appendRow(neighboring_cells_of_one_face.begin(),
741 neighboring_cells_of_one_face.end());
742 // Construct global face normal(s) (only one 'needed') and add it to "mutable_face_normals"
743 // Construct two vectors in the face, e.g. difference of two conners with the centroid,
744 // then obtain an orthogonal vector to both of them. Finally, normalize.
745 // Auxuliary vectors on the face:
746 GlobalCoordinate face_vector0 =
747 refined_corners[face4corners[0]].center() - global(local_refined_face_centroid);
748 GlobalCoordinate face_vector1 =
749 refined_corners[face4corners[1]].center() - global(local_refined_face_centroid);
750 mutable_face_normals[idx] = {
751 (face_vector0[1]*face_vector1[2]) - (face_vector0[2]*face_vector1[1]),
752 (face_vector0[2]*face_vector1[0]) - (face_vector0[0]*face_vector1[2]),
753 (face_vector0[0]*face_vector1[1]) - (face_vector0[1]*face_vector1[0])};
754 mutable_face_normals[idx] /= mutable_face_normals[idx].two_norm();
755 if (face_type == J_FACE) {
756 mutable_face_normals[idx] *= -1;
757 }
758 // Construct "refined_face_to_edges"
759 // with the {edge_indix[0], edge_index[1]} for each edge of the refined face.
760 std::vector<std::array<int,2>> refined_face_to_edges = {
761 { face4corners[0], face4corners[1]},
762 { face4corners[0], face4corners[2]},
763 { face4corners[1], face4corners[3]},
764 { face4corners[2], face4corners[3]}};
765 // Calculate the AREA of each face of a global refined face,
766 // by adding the 4 areas of the triangles partitioning each face.
767 double refined_face_area = 0.0;
768 for (int edge = 0; edge < 4; ++edge) {
769 // Construction of each triangle on the current face with one
770 // of its edges equal to "edge".
771 Geometry<0,3>::GlobalCoordinate trian_corners[3] = {
772 refined_corners[refined_face_to_edges[edge][0]].center(),
773 refined_corners[refined_face_to_edges[edge][1]].center(),
774 global(local_refined_face_centroid)};
775 refined_face_area += std::fabs(area(trian_corners));
776 } // end edge-for-loop
777 //
778 //
779 // Construct the Geometry<2,3> of the global refined face.
780 refined_faces[idx] = Geometry<2,cdim>(this->global(local_refined_face_centroid),
781 refined_face_area);
782 } // end n-for-loop
783 } // end m-for-loop
784 } // end l-for-loop
785 } // end r-for-loop
787 //
789 // We need to populate "refined_cells"
790 // "refined_cells"'s size is cells_per_dim[0] * cells_per_dim[1] * cells_per_dim[2].
791 // To build each global refined cell, we need
792 // 1. its global refined CENTER
793 // 2. its VOLUME
794 // 3. all global refined corners [available in "refined_corners"]
795 // 4. indices of its 8 corners.
796 //
797 refined_cells.resize(cells_per_dim[0] * cells_per_dim[1] * cells_per_dim[2]);
798 // Vector to store, in each entry, the 8 indices of the 8 corners
799 // of each global refined cell. Determine its size.
800 refined_cell_to_point.resize(cells_per_dim[0] * cells_per_dim[1] * cells_per_dim[2]);
801 // The numembering starts with index 0 for the refined cell with corners
802 // {0,0,0}, ...,{1/cells_per_dim[0], 1/cells_per_dim[1], 1/cells_per_dim[2]},
803 // then the indices grow first picking the cells in the x-axis (Right, i), then y-axis (Back, j), and
804 // finally, z-axis (Up, k).
805 //
806 // CENTERS
807 // REFINED CELL CENTERS
808 // The strategy is to compute the centers of the refined local
809 // unit/reference cube, and then apply the map global().
810 //
811 // VOLUMES OF THE REFINED CELLS
812 // REMARK: Each global refined 'cell' is a hexahedron since it may not be cube-shaped
813 // since its a 'deformation' of unit/reference cube. We may use 'hexahedron' to refer
814 // to the global refined cell in the computation of its volume.
815 //
816 // The strategy is to construct 24 tetrahedra in each hexahedron.
817 // Each tetrahedron is built with
818 // (1) the center of the hexahedron,
819 // (2) the middle point of the face the tetrahedron is based on, and
820 // (3) one of the edges of the face mentioned in 2.
821 // Each face 'supports' 4 tetrahedra, and we have 6 faces per hexahedron, which
822 // gives us the 24 tetrahedra per 'cell' (hexahedron).
823 //
824 // To compute the volume of each tetrahedron, we use "simplex_volume()" with
825 // the 6 corners of the tetrahedron as arguments. Summing up the 24 volumes,
826 // we get the volumne of the hexahedorn (global refined 'cell').
827 //
828 // Sum of all the volumes of all the (children) global refined cells.
829 double sum_all_refined_cell_volumes = 0.0;
830 //
831 // For each (global refined 'cell') hexahedron, to create 24 tetrahedra and their volumes,
832 // we introduce
833 // Vol1. "hexa_to_face" (needed to access face centroids).
834 // Vol2. "hexa_face_centroids" (one of the 6 corners of all 4 tetrahedra based on that face).
835 // Vol3. the center of the global refined 'cell' (hexahedron)
836 // (common corner of the 24 tetrahedra).
837 // Vol4. "tetra_edge_indices" indices of the 4x6 tetrahedra per 'cell',
838 // grouped by the face they are based on.
839 // Then we construct and compute the volume of the 24 tetrahedra with mainly
840 // "hexa_face_centroids" (Vol2.), global refined cell center (Vol3.), and "tetra_edge_indices" (Vol4.).
841 //
842 for (int k = 0; k < cells_per_dim[2]; ++k) {
843 for (int j = 0; j < cells_per_dim[1]; ++j) {
844 for (int i = 0; i < cells_per_dim[0]; ++i) {
845 // INDEX of the global refined cell associated with 'kji'.
846 int refined_cell_idx = (k*cells_per_dim[0]*cells_per_dim[1]) + (j*cells_per_dim[0]) +i;
847 // 1. CENTER of the global refined cell associated with 'kji' (Vol3.)
848 // Compute the center of the local refined unit/reference cube associated with 'kji'.
849 const LocalCoordinate& local_refined_cell_center = {
850 (.5 + i)/cells_per_dim[0], (.5 + j)/cells_per_dim[1], (.5 + k)/cells_per_dim[2]};
851 // Obtain the global refined center with 'this->global(local_refined_cell_center)'.
852 // 2. VOLUME of the global refined 'kji' cell
853 double refined_cell_volume = 0.0; // (computed below!)
854 // 3. All Global refined corners ("refined_corners")
855 // 4. Indices of the 8 corners of the global refined cell associated with 'kji'.
856 std::array<int,8> cell8corners_indices = { //
857 (j*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (i*(cells_per_dim[2]+1)) +k, // fake '0' {0,0,0}
858 (j*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((i+1)*(cells_per_dim[2]+1)) +k, // fake '1' {1,0,0}
859 ((j+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (i*(cells_per_dim[2]+1)) +k, // fake '2' {0,1,0}
860 ((j+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((i+1)*(cells_per_dim[2]+1)) +k, // fake '3' {1,1,0}
861 (j*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (i*(cells_per_dim[2]+1)) +k+1, // fake '4' {0,0,1}
862 (j*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((i+1)*(cells_per_dim[2]+1)) +k+1, // fake '5' {1,0,1}
863 ((j+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (i*(cells_per_dim[2]+1)) +k+1, // fake '6' {0,1,1}
864 ((j+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((i+1)*(cells_per_dim[2]+1)) +k+1 // fake '7' {1,1,1}
865 };
866 // Add this 8 corners to the corresponding entry of "refined_cell_to_point".
867 refined_cell_to_point[refined_cell_idx] = cell8corners_indices;
868 //
869 // VOLUME HEXAHEDRON (GLOBAL REFINED 'CELL')
870 // Vol1. INDICES ('from 0 to 5') of the faces of the hexahedron (needed to access face centroids).
871 std::vector<int> hexa_to_face = { //hexa_face_0to5_indices = {
872 // index face '0' bottom
873 (k*cells_per_dim[0]*cells_per_dim[1]) + (j*cells_per_dim[0]) + i,
874 // index face '1' front
875 (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2]+1))
876 + ((cells_per_dim[0]+1)*cells_per_dim[1]*cells_per_dim[2])
877 + (j*cells_per_dim[0]*cells_per_dim[2]) + (i*cells_per_dim[2]) + k,
878 // index face '2' left
879 (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2]+1))
880 + (i*cells_per_dim[1]*cells_per_dim[2]) + (k*cells_per_dim[1]) + j,
881 // index face '3' right
882 (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2]+1))
883 + ((i+1)*cells_per_dim[1]*cells_per_dim[2]) + (k*cells_per_dim[1]) + j,
884 // index face '4' back
885 (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2]+1)) +
886 ((cells_per_dim[0]+1)*cells_per_dim[1]*cells_per_dim[2])
887 + ((j+1)*cells_per_dim[0]*cells_per_dim[2]) + (i*cells_per_dim[2]) + k,
888 // index face '5' top
889 ((k+1)*cells_per_dim[0]*cells_per_dim[1]) + (j*cells_per_dim[0]) + i};
890 //
891 // We add the 6 faces of the cell into "refined_cell_to_face".
892 using cpgrid::EntityRep;
893 // First value is index, Second is orientation.
894 // Still have to find out what the orientation should be.
895 // right face ('3') outer normal points 'from left to right' -> orientation true
896 // back face ('4') outer normal points 'from front to back' -> orientation true
897 // top face ('5') outer normal points 'from bottom to top' -> orientation true
898 // (the other cases are false)
899 std::vector<cpgrid::EntityRep<1>> faces_of_one_cell = {
900 { hexa_to_face[0], false}, {hexa_to_face[1], false},
901 { hexa_to_face[2], false}, {hexa_to_face[3], true},
902 { hexa_to_face[4], true}, {hexa_to_face[5], true} };
903 refined_cell_to_face.appendRow(faces_of_one_cell.begin(), faces_of_one_cell.end());
904 //
905 // Vol2. CENTROIDS of the faces of the hexahedron.
906 // (one of the 6 corners of all 4 tetrahedra based on that face).
907 std::vector<Geometry<0,3>::GlobalCoordinate> hexa_face_centroids;
908 for (auto& idx : hexa_to_face) {
909 hexa_face_centroids.push_back(refined_faces[idx].center());
910 }
911 // Indices of the 4 edges of each face of the hexahedron.
912 // A tetrahedron has six edges. Once we choose a face to base a
913 // tetrahedron on, we choose an edge of that face as one of the
914 // edges of the tetrahedron. The other five edges are fixed, since
915 // the center of the hexahedron and the center of the face are
916 // the reminder 2 corners of the tetrahedron. That's why to identify
917 // a tetrahedron we only need two things:
918 // (1) the face it's based on and
919 // (2) one of the four edges of that face.
920 //
921 // Container with 6 entries, one per face. Each entry has the
922 // 4 indices of the 4 corners of each face.
923 std::vector<std::array<int,4>> cell_face4corners;
924 cell_face4corners.reserve(6);
925 for (int face = 0; face < 6; ++face) {
926 cell_face4corners.push_back({ refined_face_to_point[hexa_to_face[face]][0],
927 refined_face_to_point[hexa_to_face[face]][1],
928 refined_face_to_point[hexa_to_face[face]][2],
929 refined_face_to_point[hexa_to_face[face]][3]});
930 }
931 // Vol4. Container with indices of the edges of the 4 tetrahedra per face
932 // [according to description above]
933 std::vector<std::vector<std::array<int,2>>> tetra_edge_indices;
934 tetra_edge_indices.reserve(6);
935 for (auto& face_indices : cell_face4corners)
936 {
937 std::vector<std::array<int,2>> face4edges_indices = {
938 { face_indices[0], face_indices[1]}, // fake '{0,1}'/'{4,5}'
939 { face_indices[0], face_indices[2]}, // fake '{0,2}'/'{4,6}'
940 { face_indices[1], face_indices[3]}, // fake '{1,3}'/'{5,7}'
941 { face_indices[2], face_indices[3]} }; // fake '{2,3}'/'{6,7}'
942 tetra_edge_indices.push_back(face4edges_indices);
943 }
944 // Sum of the 24 volumes to get the volume of the hexahedron,
945 // stored in "refined_cell_volume".
946 // Calculate the volume of each hexahedron, by adding
947 // the 4 tetrahedra at each face (4x6 = 24 tetrahedra).
948 for (int face = 0; face < 6; ++face) {
949 for (int edge = 0; edge < 4; ++edge) {
950 // Construction of each tetrahedron based on "face" with one
951 // of its edges equal to "edge".
952 const Geometry<0, 3>::GlobalCoordinate tetra_corners[4] = {
953 refined_corners[tetra_edge_indices[face][edge][0]].center(), // (see Vol4.)
954 refined_corners[tetra_edge_indices[face][edge][1]].center(), // (see Vol4.)
955 hexa_face_centroids[face], // (see Vol2.)
956 // global_refined_cell_center
957 this->global(local_refined_cell_center)}; // (see Vol3.)
958 refined_cell_volume += std::fabs(simplex_volume(tetra_corners));
959 } // end edge-for-loop
960 } // end face-for-loop
961 // Add the volume of the hexahedron (global refined 'cell')
962 // to the container with of all volumes of all the refined cells.
963 sum_all_refined_cell_volumes += refined_cell_volume;
964 // Create a pointer to the first element of "refined_cell_to_point"
965 // (required as the fourth argement to construct a Geometry<3,3> type object).
966 int* indices_storage_ptr = refined_cell_to_point[refined_cell_idx].data();
967 // Construct the Geometry of the refined cell associated with 'kji'.
968 refined_cells[refined_cell_idx] =
969 Geometry<3,cdim>(this->global(local_refined_cell_center),
970 refined_cell_volume,
971 all_geom.geomVector(std::integral_constant<int,3>()),
972 indices_storage_ptr);
973 } // end i-for-loop
974 } // end j-for-loop
975 } // end k-for-loop
976 // Rescale all volumes if the sum of volume of all the global refined 'cells' does not match the
977 // volume of the 'parent cell'.
978 // Compare the sum of all the volumes of all refined cells with 'parent cell' volume.
979 if (std::fabs(sum_all_refined_cell_volumes - this->volume())
980 > std::numeric_limits<Geometry<3, cdim>::ctype>::epsilon()) {
981 Geometry<3, cdim>::ctype correction = this->volume() / sum_all_refined_cell_volumes;
982 for(auto& cell: refined_cells){
983 cell.vol_ *= correction;
984 }
985 } // end if-statement
987 }
988
989 private:
990 GlobalCoordinate pos_;
991 double vol_;
992 std::shared_ptr<const EntityVariable<Geometry<0, 3>,3>> allcorners_; // For dimension 3 only
993 const int* cor_idx_; // For dimension 3 only
994
1008 const std::tuple< enum face_tag, int,
1009 std::array<int, 4>, std::vector<cpgrid::EntityRep<0>>,
1010 LocalCoordinate>
1011 getIndicesFace(int l, int m, int n, int constant_direction, const std::array<int, 3>& cells_per_dim) const
1012 {
1013 using cpgrid::EntityRep;
1014 std::vector<cpgrid::EntityRep<0>> neighboring_cells_of_one_face; // {index, orientation}
1015 switch(constant_direction) {
1016 case 0: // {l,m,n} = {k,j,i}, constant in z-direction
1017 // Orientation true when outer normal points 'from bottom to top'
1018 // Orientation false when outer normal points 'from top to bottom'
1019 if (l != 0) {
1020 neighboring_cells_of_one_face.push_back({((l-1)*cells_per_dim[0]*cells_per_dim[1])
1021 + (m*cells_per_dim[0]) + n, true});
1022 }
1023 if (l != cells_per_dim[2]) {
1024 neighboring_cells_of_one_face.push_back({ (l*cells_per_dim[0]*cells_per_dim[1])
1025 + (m*cells_per_dim[0]) + n, false});
1026 }
1027 return { face_tag::K_FACE, (l*cells_per_dim[0]*cells_per_dim[1]) + (m*cells_per_dim[0]) + n,
1028 {(m*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (n*(cells_per_dim[2]+1)) +l,
1029 (m*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((n+1)*(cells_per_dim[2]+1)) +l,
1030 ((m+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (n*(cells_per_dim[2]+1)) +l,
1031 ((m+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((n+1)*(cells_per_dim[2]+1)) +l},
1032 neighboring_cells_of_one_face,
1033 {(.5 + n)/cells_per_dim[0], (.5 + m)/cells_per_dim[1], double(l)/cells_per_dim[2]}};
1034 case 1: // {l,m,n} = {i,k,j}, constant in the x-direction
1035 // Orientation true when outer normal points 'from left to right'
1036 // Orientation false when outer normal points 'from right to left'
1037 if (l != 0) {
1038 neighboring_cells_of_one_face.push_back({(m*cells_per_dim[0]*cells_per_dim[1])
1039 + (n*cells_per_dim[0]) +l-1, true});
1040 }
1041 if (l != cells_per_dim[0]) {
1042 neighboring_cells_of_one_face.push_back({ (m*cells_per_dim[0]*cells_per_dim[1])
1043 + (n*cells_per_dim[0]) + l, false});
1044 }
1045 return { face_tag::I_FACE, (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2]+1))
1046 + (l*cells_per_dim[1]*cells_per_dim[2]) + (m*cells_per_dim[1]) + n,
1047 {(n*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (l*(cells_per_dim[2]+1)) +m,
1048 ((n+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (l*(cells_per_dim[2]+1)) +m,
1049 (n*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (l*(cells_per_dim[2]+1)) +m+1,
1050 ((n+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (l*(cells_per_dim[2]+1)) +m+1},
1051 neighboring_cells_of_one_face,
1052 { double(l)/cells_per_dim[0], (.5 + n)/cells_per_dim[1], (.5 + m)/cells_per_dim[2]}};
1053 case 2: // {l,m,n} = {j,i,k}, constant in the y-direction
1054 // Orientation true when outer normal points 'from front to back'
1055 // Orientation false when outer normal points 'from back to front'
1056 if (l != 0) {
1057 neighboring_cells_of_one_face.push_back({(n*cells_per_dim[0]*cells_per_dim[1])
1058 + ((l-1)*cells_per_dim[0]) +m, true});
1059 }
1060 if (l != cells_per_dim[1]) {
1061 neighboring_cells_of_one_face.push_back({(n*cells_per_dim[0]*cells_per_dim[1])
1062 + (l*cells_per_dim[0]) + m, false});
1063 }
1064 return { face_tag::J_FACE, (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2] +1))
1065 + ((cells_per_dim[0]+1)*cells_per_dim[1]*cells_per_dim[2])
1066 + (l*cells_per_dim[0]*cells_per_dim[2]) + (m*cells_per_dim[2]) + n,
1067 {(l*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (m*(cells_per_dim[2]+1)) +n,
1068 (l*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((m+1)*(cells_per_dim[2]+1)) +n,
1069 (l*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (m*(cells_per_dim[2]+1)) +n+1,
1070 (l*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((m+1)*(cells_per_dim[2]+1)) +n+1},
1071 neighboring_cells_of_one_face,
1072 {(.5 + m)/cells_per_dim[0], double(l)/cells_per_dim[1], (.5 + n)/cells_per_dim[2]}};
1073 default:
1074 // Should never be reached, but prevents compiler warning
1075 OPM_THROW(std::logic_error, "Unhandled dimension. This should never happen!");
1076 }
1077 }
1078 };
1079 } // namespace cpgrid
1080
1081 template< int mydim, int cdim >
1082 auto referenceElement(const cpgrid::Geometry<mydim,cdim>& geo) -> decltype(referenceElement<double,mydim>(geo.type()))
1083 {
1084 return referenceElement<double,mydim>(geo.type());
1085 }
1086
1087} // namespace Dune
1088
1089#endif // OPM_GEOMETRY_HEADER
Definition DefaultGeometryPolicy.hpp:53
const EntityVariable< cpgrid::Geometry< 3 - codim, 3 >, codim > & geomVector() const
Definition DefaultGeometryPolicy.hpp:86
Represents an entity of a given codim, with positive or negative orientation.
Definition EntityRep.hpp:99
Base class for EntityVariable and SignedEntityVariable.
Definition EntityRep.hpp:219
A class design to hold a variable with a value for each entity of the given codimension,...
Definition EntityRep.hpp:267
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition Geometry.hpp:102
bool affine() const
The mapping implemented by this geometry is constant, therefore affine.
Definition Geometry.hpp:210
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition Geometry.hpp:105
JacobianInverse jacobianInverse(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:204
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition Geometry.hpp:100
Jacobian jacobian(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:197
double integrationElement(const LocalCoordinate &) const
Returns 1 for the vertex geometry.
Definition Geometry.hpp:141
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition Geometry.hpp:111
Geometry(const GlobalCoordinate &pos)
Construct from vertex position.
Definition Geometry.hpp:116
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse
Type of inverse of Jacobian matrix.
Definition Geometry.hpp:107
JacobianTransposed jacobianTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:180
GeometryType type() const
Using the cube type for vertices.
Definition Geometry.hpp:147
ctype volume() const
Volume of vertex is arbitrarily set to 1.
Definition Geometry.hpp:167
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition Geometry.hpp:109
GlobalCoordinate corner(int cor) const
Returns the single corner: the vertex itself.
Definition Geometry.hpp:159
int corners() const
A vertex is defined by a single corner.
Definition Geometry.hpp:153
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition Geometry.hpp:173
Geometry()
Default constructor, giving a non-valid geometry.
Definition Geometry.hpp:122
const GlobalCoordinate & global(const LocalCoordinate &) const
Returns the position of the vertex.
Definition Geometry.hpp:128
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:189
LocalCoordinate local(const GlobalCoordinate &) const
Meaningless for the vertex geometry.
Definition Geometry.hpp:134
double ctype
Coordinate element type.
Definition Geometry.hpp:97
JacobianInverse jacobianInverse(const LocalCoordinate &) const
The inverse of the jacobian.
Definition Geometry.hpp:346
const FieldMatrix< ctype, mydimension, coorddimension > & jacobianTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:325
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition Geometry.hpp:318
int corners() const
The number of corners of this convex polytope.
Definition Geometry.hpp:297
LocalCoordinate local(const GlobalCoordinate &) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:277
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition Geometry.hpp:247
bool affine() const
Since integrationElement() is constant, returns true.
Definition Geometry.hpp:352
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition Geometry.hpp:242
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition Geometry.hpp:253
ctype volume() const
Volume (area, actually) of intersection.
Definition Geometry.hpp:312
const FieldMatrix< ctype, coorddimension, mydimension > & jacobianInverseTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:332
GeometryType type() const
We use the singular type (None) for intersections.
Definition Geometry.hpp:290
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition Geometry.hpp:251
double integrationElement(const LocalCoordinate &) const
For the singular geometry, we return a constant integration element equal to the volume.
Definition Geometry.hpp:284
Geometry()
Default constructor, giving a non-valid geometry.
Definition Geometry.hpp:265
const GlobalCoordinate & global(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:271
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse
Type of inverse of Jacobian matrix.
Definition Geometry.hpp:249
Geometry(const GlobalCoordinate &pos, ctype vol)
Construct from centroid and volume (1- and 0-moments).
Definition Geometry.hpp:258
double ctype
Coordinate element type.
Definition Geometry.hpp:239
Jacobian jacobian(const LocalCoordinate &) const
The jacobian.
Definition Geometry.hpp:339
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition Geometry.hpp:244
GlobalCoordinate corner(int) const
This method is meaningless for singular geometries.
Definition Geometry.hpp:303
const JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local_coord) const
Inverse of Jacobian transposed.
Definition Geometry.hpp:586
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition Geometry.hpp:396
bool affine() const
The mapping implemented by this geometry is not generally affine.
Definition Geometry.hpp:608
double ctype
Coordinate element type.
Definition Geometry.hpp:382
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition Geometry.hpp:387
GlobalCoordinate corner(int cor) const
Get the cor-th of 8 corners of the hexahedral base cell.
Definition Geometry.hpp:521
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition Geometry.hpp:394
Geometry(const GlobalCoordinate &pos, ctype vol, std::shared_ptr< const EntityVariable< cpgrid::Geometry< 0, 3 >, 3 > > allcorners_ptr, const int *corner_indices)
Construct from center, volume (1- and 0-moments) and corners.
Definition Geometry.hpp:408
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:508
double integrationElement(const LocalCoordinate &local_coord) const
Equal to \sqrt{\det{J^T J}} where J is the Jacobian.
Definition Geometry.hpp:500
const JacobianTransposed jacobianTransposed(const LocalCoordinate &local_coord) const
Jacobian transposed.
Definition Geometry.hpp:550
GlobalCoordinate global(const LocalCoordinate &local_coord) const
Provide a trilinear mapping.
Definition Geometry.hpp:444
Geometry()
Default constructor, giving a non-valid geometry.
Definition Geometry.hpp:434
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverse
Type of inverse of Jacobian matrix.
Definition Geometry.hpp:392
LocalCoordinate local(const GlobalCoordinate &y) const
Mapping from the cell to the reference domain.
Definition Geometry.hpp:475
int corners() const
The number of corners of this convex polytope.
Definition Geometry.hpp:515
ctype volume() const
Cell volume.
Definition Geometry.hpp:528
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition Geometry.hpp:538
Jacobian jacobian(const LocalCoordinate &local_coord) const
The jacobian.
Definition Geometry.hpp:595
Geometry(const GlobalCoordinate &pos, ctype vol)
Construct from centroid and volume (1- and 0-moments).
Definition Geometry.hpp:427
void refine(const std::array< int, 3 > &cells_per_dim, DefaultGeometryPolicy &all_geom, std::vector< std::array< int, 8 > > &refined_cell_to_point, cpgrid::OrientedEntityTable< 0, 1 > &refined_cell_to_face, Opm::SparseTable< int > &refined_face_to_point, cpgrid::OrientedEntityTable< 1, 0 > &refined_face_to_cell, cpgrid::EntityVariable< enum face_tag, 1 > &refined_face_tags, cpgrid::SignedEntityVariable< PointType, 1 > &refined_face_normals) const
Definition Geometry.hpp:630
Dune::FieldVector< double, 3 > PointType
Refine a single cell with regular intervals.
Definition Geometry.hpp:629
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition Geometry.hpp:385
JacobianInverse jacobianInverse(const LocalCoordinate &local_coord) const
The inverse of the jacobian.
Definition Geometry.hpp:602
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition Geometry.hpp:390
This class encapsulates geometry for vertices, intersections, and cells.
Definition Geometry.hpp:75
Represents the topological relationships between sets of entities, for example cells and faces.
Definition OrientedEntityTable.hpp:139
void appendRow(DataIter row_beg, DataIter row_end)
Appends a row to the table.
Definition SparseTable.hpp:108
A class design to hold a variable with a value for each entity of the given codimension,...
Definition EntityRep.hpp:302
A SparseTable stores a table with rows of varying size as efficiently as possible.
Definition SparseTable.hpp:55
void appendRow(DataIter row_beg, DataIter row_end)
Appends a row to the table.
Definition SparseTable.hpp:108
Copyright 2019 Equinor AS.
Definition CartesianIndexMapper.hpp:10
T area(const Point< T, 2 > *c)
Computes the area of a 2-dimensional triangle.
Definition Volumes.hpp:118
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
face_tag
Connection taxonomy.
Definition preprocess.h:66
@ K_FACE
Connection topologically normal to I-J plane.
Definition preprocess.h:69
@ J_FACE
Connection topologically normal to I-K plane.
Definition preprocess.h:68
@ I_FACE
Connection topologically normal to J-K plane.
Definition preprocess.h:67