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