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
420 : pos_(0.0), vol_(0.0), allcorners_(0), cor_idx_(0)
421 {
422 }
423
429 GlobalCoordinate global(const LocalCoordinate& local_coord) const
430 {
431 static_assert(mydimension == 3, "");
432 static_assert(coorddimension == 3, "");
433 // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
434 LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
435 uvw[0] -= local_coord;
436 // Access pattern for uvw matching ordering of corners.
437 const int pat[8][3] = { { 0, 0, 0 },
438 { 1, 0, 0 },
439 { 0, 1, 0 },
440 { 1, 1, 0 },
441 { 0, 0, 1 },
442 { 1, 0, 1 },
443 { 0, 1, 1 },
444 { 1, 1, 1 } };
445 GlobalCoordinate xyz(0.0);
446 for (int i = 0; i < 8; ++i) {
447 GlobalCoordinate corner_contrib = corner(i);
448 double factor = 1.0;
449 for (int j = 0; j < 3; ++j) {
450 factor *= uvw[pat[i][j]][j];
451 }
452 corner_contrib *= factor;
453 xyz += corner_contrib;
454 }
455 return xyz;
456 }
457
461 {
462 static_assert(mydimension == 3, "");
463 static_assert(coorddimension == 3, "");
464 // This code is modified from dune/grid/genericgeometry/mapping.hh
465 // \todo: Implement direct computation.
466 const ctype epsilon = 1e-12;
467 auto refElement = Dune::ReferenceElements<ctype, 3>::cube();
468 LocalCoordinate x = refElement.position(0,0);
470 do {
471 // DF^n dx^n = F^n, x^{n+1} -= dx^n
472 JacobianTransposed JT = jacobianTransposed(x);
473 GlobalCoordinate z = global(x);
474 z -= y;
475 MatrixHelperType::template xTRightInvA<3, 3>(JT, z, dx );
476 x -= dx;
477 } while (dx.two_norm2() > epsilon*epsilon);
478 return x;
479 }
480
485 double integrationElement(const LocalCoordinate& local_coord) const
486 {
487 JacobianTransposed Jt = jacobianTransposed(local_coord);
488 return MatrixHelperType::template sqrtDetAAT<3, 3>(Jt);
489 }
490
493 GeometryType type() const
494 {
495 return Dune::GeometryTypes::cube(mydimension);
496 }
497
500 int corners() const
501 {
502 return 8;
503 }
504
507 {
508 assert(allcorners_ && cor_idx_);
509 return (allcorners_->data())[cor_idx_[cor]].center();
510 }
511
513 ctype volume() const
514 {
515 return vol_;
516 }
517
518 void set_volume(ctype volume) {
519 vol_ = volume;
520 }
521
524 {
525 return pos_;
526 }
527
534 const JacobianTransposed
535 jacobianTransposed(const LocalCoordinate& local_coord) const
536 {
537 static_assert(mydimension == 3, "");
538 static_assert(coorddimension == 3, "");
539
540 // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
541 LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
542 uvw[0] -= local_coord;
543 // Access pattern for uvw matching ordering of corners.
544 const int pat[8][3] = { { 0, 0, 0 },
545 { 1, 0, 0 },
546 { 0, 1, 0 },
547 { 1, 1, 0 },
548 { 0, 0, 1 },
549 { 1, 0, 1 },
550 { 0, 1, 1 },
551 { 1, 1, 1 } };
552 JacobianTransposed Jt(0.0);
553 for (int i = 0; i < 8; ++i) {
554 for (int deriv = 0; deriv < 3; ++deriv) {
555 // This part contributing to dg/du_{deriv}
556 double factor = 1.0;
557 for (int j = 0; j < 3; ++j) {
558 factor *= (j != deriv) ? uvw[pat[i][j]][j]
559 : (pat[i][j] == 0 ? -1.0 : 1.0);
560 }
561 GlobalCoordinate corner_contrib = corner(i);
562 corner_contrib *= factor;
563 Jt[deriv] += corner_contrib; // using FieldMatrix row access.
564 }
565 }
566 return Jt;
567 }
568
570 const JacobianInverseTransposed
572 {
573 JacobianInverseTransposed Jti = jacobianTransposed(local_coord);
574 Jti.invert();
575 return Jti;
576 }
577
579 Jacobian
580 jacobian(const LocalCoordinate& local_coord) const
581 {
582 return jacobianTransposed(local_coord).transposed();
583 }
584
586 JacobianInverse
587 jacobianInverse(const LocalCoordinate& local_coord) const
588 {
589 return jacobianInverseTransposed(local_coord).transposed();
590 }
591
593 bool affine() const
594 {
595 return false;
596 }
597
616 typedef Dune::FieldVector<double,3> PointType;
617 void refineCellifiedPatch(const std::array<int,3>& cells_per_dim,
618 DefaultGeometryPolicy& all_geom,
619 std::vector<std::array<int,8>>& refined_cell_to_point,
620 cpgrid::OrientedEntityTable<0,1>& refined_cell_to_face,
621 Opm::SparseTable<int>& refined_face_to_point,
622 cpgrid::OrientedEntityTable<1,0>& refined_face_to_cell,
624 cpgrid::SignedEntityVariable<PointType, 1>& refined_face_normals,
625 const std::array<int,3>& patch_dim,
626 const std::vector<double>& widthsX,
627 const std::vector<double>& lengthsY,
628 const std::vector<double>& heightsZ) const
629 {
631 *(all_geom.geomVector(std::integral_constant<int,3>()));
633 *(all_geom.geomVector(std::integral_constant<int,1>()));
635 *(all_geom.geomVector(std::integral_constant<int,0>()));
636 EntityVariableBase<enum face_tag>& mutable_face_tags = refined_face_tags;
637 EntityVariableBase<PointType>& mutable_face_normals = refined_face_normals;
638
640 // The strategy is to compute the local refined corners
641 // of the unit/reference cube, and then apply the map global().
642 // Determine the size of the vector containing all the corners
643 // of all the global refined cells (children cells).
644 // For easier notation:
645 const std::array<int,3>& refined_dim = { cells_per_dim[0]*patch_dim[0],
646 cells_per_dim[1]*patch_dim[1],
647 cells_per_dim[2]*patch_dim[2]};
648 refined_corners.resize((refined_dim[0] + 1)*(refined_dim[1] + 1)*(refined_dim[2] + 1));
649 // The nummbering starts at the botton, so k=0 (z-axis), and j=0 (y-axis), i=0 (x-axis).
650 // Then, increasing k ('going up'), followed by increasing i ('going right->'),
651 // and finally, increasing j ('going back'). This order criteria for corners
652 // 'Up [increasing k]- Right [incresing i]- Back [increasing j]'
653 // is consistant with cpgrid numbering.
654 //
655 assert(static_cast<int>(widthsX.size()) == patch_dim[0]);
656 assert(static_cast<int>(lengthsY.size()) == patch_dim[1]);
657 assert(static_cast<int>(heightsZ.size()) == patch_dim[2]);
658 const auto localCoordNumerator = []( const std::vector<double>& vec, int sumLimit, double multiplier) {
659 double lcn = 0;
660 assert(!vec.empty());
661 assert(sumLimit < static_cast<int>(vec.size()));
662 lcn += multiplier*vec[sumLimit];
663 for (int m = 0; m < sumLimit; ++m) {
664 lcn += vec[m];
665 }
666 return lcn;
667 };
668 // E.g. localCoordNumerator( dx, 3, 0.25) = x0 + x1 + x2 + 0.25.x3
669 //
670 const double sumWidths = std::accumulate(widthsX.begin(), widthsX.end(), double(0));
671 // x0 + x1 + ... + xL, if dx = {x0, x1, ..., xL}
672 const double sumLengths = std::accumulate(lengthsY.begin(), lengthsY.end(), double(0));
673 // y0 + y1 + ... + yM, if dy = {y0, y1, ..., yM}
674 const double sumHeights = std::accumulate(heightsZ.begin(), heightsZ.end(), double(0));
675 // z0 + z1 + ... + zN, if dz = {z0, z1, ..., zN}
676
677 for (int j = 0; j < refined_dim[1] +1; ++j) {
678 double local_y = 0;
679 for (int i = 0; i < refined_dim[0] +1; ++i) {
680 double local_x = 0.;
681 for (int k = 0; k < refined_dim[2] +1; ++k) {
682 double local_z = 0.;
683
684 // Compute the index of each global refined corner associated with 'jik'.
685 int refined_corner_idx =
686 (j*(refined_dim[0]+1)*(refined_dim[2]+1)) + (i*(refined_dim[2]+1)) + k;
687
688 // Compute the local refined corner of the unit/reference cube associated with 'jik'.
689 if ( i == refined_dim[0]) { // last corner in the x-direction
690 local_x = sumWidths;
691 } else {
692 local_x = localCoordNumerator(widthsX, i/cells_per_dim[0], double((i % cells_per_dim[0])) / cells_per_dim[0]);
693 }
694 if ( j == refined_dim[1]) { // last corner in the y-direction
695 local_y = sumLengths;
696 } else {
697 local_y = localCoordNumerator(lengthsY, j/cells_per_dim[1], double((j % cells_per_dim[1])) / cells_per_dim[1]);
698 }
699 if ( k == refined_dim[2]) { // last corner in the z-direction
700 local_z = sumHeights;
701 } else {
702 local_z = localCoordNumerator(heightsZ, k/cells_per_dim[2], double((k % cells_per_dim[2])) /cells_per_dim[2]);
703 }
704
705 const LocalCoordinate& local_refined_corner = { local_x/sumWidths, local_y/sumLengths, local_z/sumHeights };
706 assert(local_x/sumWidths <= 1.);
707 assert(local_y/sumLengths <= 1.);
708 assert(local_z/sumHeights <= 1.);
709
710 // Compute the global refined corner 'jik' and add it in its corresponfing entry in "refined_corners".
711 refined_corners[refined_corner_idx] = Geometry<0, 3>(this->global(local_refined_corner));
712 } // end k-for-loop
713 } // end i-for-loop
714 } // end j-for-loop
716 //
718 // We want to populate "refined_faces". The size of "refined_faces" is
719 const int refined_faces_size =
720 (refined_dim[0]*refined_dim[1]*(refined_dim[2]+1)) // 'bottom/top faces'
721 + ((refined_dim[0]+1)*refined_dim[1]*refined_dim[2]) // 'left/right faces'
722 + (refined_dim[0]*(refined_dim[1]+1)*refined_dim[2]); // 'front/back faces'
723 refined_faces.resize(refined_faces_size);
724 refined_face_tags.resize(refined_faces_size);
725 refined_face_normals.resize(refined_faces_size);
726 //
727 // To create a face as a Geometry<2,3> type object we need its CENTROID and its VOLUME(area).
728 // We store the centroids/areas in the following order:
729 // - Bottom-top faces -> 3rd coordinate constant in each face.
730 // - Left-right faces -> 1st coordinate constant in each face.
731 // - Front-back faces -> 2nd coordinate constant in each face.
732 //
733 // REFINED FACE AREAS
734 // To compute the area of each face, we divide it in 4 triangles,
735 // compute the area of those with "area()", where the arguments
736 // are the 3 corners of each triangle. Then, sum them up to get the area
737 // of the global refined face.
738 // For each face, we construct 4 triangles with
739 // (1) centroid of the face,
740 // (2) one of the edges of the face.
741 //
742 // A triangle has 3 edges. Once we choose a face to base a triangle on,
743 // we choose an edge of that face as one of the edges of the triangle.
744 // The other 2 edges are fixed, since the centroid of the face the triangle
745 // is based on is one of its corners. That's why to identify
746 // a triangle we only need two things:
747 // (1) the face it's based on and
748 // (2) one of the four edges of that face.
749 //
750 // For each face, we need
751 // 1. index of the refined face
752 // [needed to access indices of the 4 edges of the face in "refined_face_to_edges"]
753 // 2. centroid of the face (common corner of the 4 triangles based on that face).
754 // [available via "['face'].center()"
755 // 3. container of 4 entries (the 4 edges of the face).
756 // Each entry consists in the 2 indices defining each edge of the face.
757 // [available in "refined_face_to_edges"].
758 //
759 // Populate "mutable_face_tags/normals", "refined_face_to_point/cell",
760 // "refined_faces".
761 //
762 for (int constant_direction = 0; constant_direction < 3; ++constant_direction){
763 // adding %3 and constant_direction, we go through the 3 type of faces.
764 // 0 -> 3rd coordinate constant: l('k') < cells_per_dim[2]+1, m('j') < cells_per_dim[1], n('i') < cells_per_dim[0]
765 // 1 -> 1rt coordinate constant: l('i') < cells_per_dim[0]+1, m('k') < cells_per_dim[2], n('j') < cells_per_dim[1]
766 // 2 -> 2nd coordinate constant: l('j') < cells_per_dim[1]+1, m('i') < cells_per_dim[0], n('k') < cells_per_dim[2]
767 std::array<int,3> refined_dim_mixed = {
768 refined_dim[(2+constant_direction)%3],
769 refined_dim[(1+constant_direction)%3],
770 refined_dim[constant_direction % 3] };
771 for (int l = 0; l < refined_dim_mixed[0] + 1; ++l) {
772 for (int m = 0; m < refined_dim_mixed[1]; ++m) {
773 for (int n = 0; n < refined_dim_mixed[2]; ++n) {
774 // Compute the face data.
775 auto [face_tag, idx, face_to_point, face_to_cell, wrong_local_centroid] =
776 getIndicesFace(l, m, n, constant_direction, refined_dim);
777 // Add the tag to "refined_face_tags".
778 mutable_face_tags[idx]= face_tag;
779 // Add the 4 corners of the face to "refined_face_to_point".
780 refined_face_to_point.appendRow(face_to_point.begin(), face_to_point.end());
781 // Add the neighboring cells of the face to "refined_face_to_cell".
782 refined_face_to_cell.appendRow(face_to_cell.begin(), face_to_cell.end());
783 // Compute the centroid as the average of the 4 corners of the face
784 GlobalCoordinate face_center = { 0., 0., 0.};
785 for (int corn = 0; corn < 4; ++corn){
786 face_center += refined_corners[face_to_point[corn]].center();
787 }
788 face_center /= 4.;
789 // Construct global face normal(s) (only one 'needed') and add it to "mutable_face_normals"
790 // Construct two vectors in the face, e.g. difference of two conners with the centroid,
791 // then obtain an orthogonal vector to both of them. Finally, normalize.
792 // Auxuliary vectors on the face:
793 GlobalCoordinate face_vector0 = refined_corners[face_to_point[0]].center() - face_center;
794 GlobalCoordinate face_vector1 = refined_corners[face_to_point[1]].center() - face_center;
795 mutable_face_normals[idx] = {
796 (face_vector0[1]*face_vector1[2]) - (face_vector0[2]*face_vector1[1]),
797 (face_vector0[2]*face_vector1[0]) - (face_vector0[0]*face_vector1[2]),
798 (face_vector0[0]*face_vector1[1]) - (face_vector0[1]*face_vector1[0])};
799 mutable_face_normals[idx] /= mutable_face_normals[idx].two_norm();
800 if (face_tag == J_FACE) {
801 mutable_face_normals[idx] *= -1;
802 }
803 // Construct "refined_face_to_edges"
804 // with the {edge_indix[0], edge_index[1]} for each edge of the refined face.
805 std::vector<std::array<int,2>> refined_face_to_edges = {
806 { face_to_point[0], face_to_point[1] },
807 { face_to_point[0], face_to_point[2] },
808 { face_to_point[1], face_to_point[3] },
809 { face_to_point[2], face_to_point[3] }
810 };
811 // Calculate the AREA of each face of a global refined face,
812 // by adding the 4 areas of the triangles partitioning each face.
813 double refined_face_area = 0.0;
814 for (int edge = 0; edge < 4; ++edge) {
815 // Construction of each triangle on the current face with one
816 // of its edges equal to "edge".
817 Geometry<0,3>::GlobalCoordinate trian_corners[3] = {
818 refined_corners[refined_face_to_edges[edge][0]].center(),
819 refined_corners[refined_face_to_edges[edge][1]].center(),
820 face_center };
821 refined_face_area += std::fabs(area(trian_corners));
822 } // end edge-for-loop
823 //
824 //
825 // Construct the Geometry<2,3> of the global refined face.
826 refined_faces[idx] = Geometry<2,cdim>(face_center, refined_face_area);
827 } // end n-for-loop
828 } // end m-for-loop
829 } // end l-for-loop
830 } // end r-for-loop
832 //
834 // We need to populate "refined_cells"
835 // "refined_cells"'s size is cells_per_dim[0] * cells_per_dim[1] * cells_per_dim[2].
836 // To build each global refined cell, we need
837 // 1. its global refined CENTER
838 // 2. its VOLUME
839 // 3. all global refined corners [available in "refined_corners"]
840 // 4. indices of its 8 corners.
841 //
842 refined_cells.resize(refined_dim[0] * refined_dim[1] * refined_dim[2]);
843 // Vector to store, in each entry, the 8 indices of the 8 corners
844 // of each global refined cell. Determine its size.
845 refined_cell_to_point.resize(refined_dim[0] * refined_dim[1] * refined_dim[2]);
846 // The numembering starts with index 0 for the refined cell with corners
847 // {0,0,0}, ...,{1/cells_per_dim[0], 1/cells_per_dim[1], 1/cells_per_dim[2]},
848 // then the indices grow first picking the cells in the x-axis (Right, i), then y-axis (Back, j), and
849 // finally, z-axis (Up, k).
850 //
851 // CENTERS
852 // REFINED CELL CENTERS
853 // The strategy is to compute the centers of the refined local
854 // unit/reference cube, and then apply the map global().
855 //
856 // VOLUMES OF THE REFINED CELLS
857 // REMARK: Each global refined 'cell' is a hexahedron since it may not be cube-shaped
858 // since its a 'deformation' of unit/reference cube. We may use 'hexahedron' to refer
859 // to the global refined cell in the computation of its volume.
860 //
861 // The strategy is to construct 24 tetrahedra in each hexahedron.
862 // Each tetrahedron is built with
863 // (1) the center of the hexahedron,
864 // (2) the middle point of the face the tetrahedron is based on, and
865 // (3) one of the edges of the face mentioned in 2.
866 // Each face 'supports' 4 tetrahedra, and we have 6 faces per hexahedron, which
867 // gives us the 24 tetrahedra per 'cell' (hexahedron).
868 //
869 // To compute the volume of each tetrahedron, we use "simplex_volume()" with
870 // the 6 corners of the tetrahedron as arguments. Summing up the 24 volumes,
871 // we get the volumne of the hexahedorn (global refined 'cell').
872 //
873 // Sum of all the volumes of all the (children) global refined cells.
874 double sum_all_refined_cell_volumes = 0.0;
875 //
876 // For each (global refined 'cell') hexahedron, to create 24 tetrahedra and their volumes,
877 // we introduce
878 // Vol1. "hexa_to_face" (needed to access face centroids).
879 // Vol2. "hexa_face_centroids" (one of the 6 corners of all 4 tetrahedra based on that face).
880 // Vol3. the center of the global refined 'cell' (hexahedron)
881 // (common corner of the 24 tetrahedra).
882 // Vol4. "tetra_edge_indices" indices of the 4x6 tetrahedra per 'cell',
883 // grouped by the face they are based on.
884 // Then we construct and compute the volume of the 24 tetrahedra with mainly
885 // "hexa_face_centroids" (Vol2.), global refined cell center (Vol3.), and "tetra_edge_indices" (Vol4.).
886 //
887 for (int k = 0; k < refined_dim[2]; ++k) {
888 for (int j = 0; j < refined_dim[1]; ++j) {
889 for (int i = 0; i < refined_dim[0]; ++i) {
890 // INDEX of the global refined cell associated with 'kji'.
891 int refined_cell_idx = (k*refined_dim[0]*refined_dim[1]) + (j*refined_dim[0]) +i;
892 // Obtain the global refined center with 'this->global(local_refined_cell_center)'.
893 // 2. VOLUME of the global refined 'kji' cell
894 double refined_cell_volume = 0.0; // (computed below!)
895 // 3. All Global refined corners ("refined_corners")
896 // 4. Indices of the 8 corners of the global refined cell associated with 'kji'.
897 std::array<int,8> cell_to_point = { //
898 (j*(refined_dim[0]+1)*(refined_dim[2]+1)) + (i*(refined_dim[2]+1)) +k, // fake '0' {0,0,0}
899 (j*(refined_dim[0]+1)*(refined_dim[2]+1)) + ((i+1)*(refined_dim[2]+1)) +k, // fake '1' {1,0,0}
900 ((j+1)*(refined_dim[0]+1)*(refined_dim[2]+1)) + (i*(refined_dim[2]+1)) +k, // fake '2' {0,1,0}
901 ((j+1)*(refined_dim[0]+1)*(refined_dim[2]+1)) + ((i+1)*(refined_dim[2]+1)) +k, // fake '3' {1,1,0}
902 (j*(refined_dim[0]+1)*(refined_dim[2]+1)) + (i*(refined_dim[2]+1)) +k+1, // fake '4' {0,0,1}
903 (j*(refined_dim[0]+1)*(refined_dim[2]+1)) + ((i+1)*(refined_dim[2]+1)) +k+1, // fake '5' {1,0,1}
904 ((j+1)*(refined_dim[0]+1)*(refined_dim[2]+1)) + (i*(refined_dim[2]+1)) +k+1, // fake '6' {0,1,1}
905 ((j+1)*(refined_dim[0]+1)*(refined_dim[2]+1)) + ((i+1)*(refined_dim[2]+1)) +k+1 // fake '7' {1,1,1}
906 };
907 // Add this 8 corners to the corresponding entry of "refined_cell_to_point".
908 refined_cell_to_point[refined_cell_idx] = cell_to_point;
909 // 1. CENTER of the global refined cell associated with 'kji' (Vol3.)
910 // Compute the center of the local refined unit/reference cube associated with 'kji'.
911 GlobalCoordinate refined_cell_center = {0., 0., 0.};
912 for (int corn = 0; corn < 8; ++corn) {
913 refined_cell_center += refined_corners[cell_to_point[corn]].center();
914 }
915 refined_cell_center /= 8.;
916 //
917 // VOLUME HEXAHEDRON (GLOBAL REFINED 'CELL')
918 // Vol1. INDICES ('from 0 to 5') of the faces of the hexahedron (needed to access face centroids).
919 std::vector<int> hexa_to_face = { //hexa_face_0to5_indices = {
920 // index face '0' bottom
921 (k*refined_dim[0]*refined_dim[1]) + (j*refined_dim[0]) + i,
922 // index face '1' front
923 (refined_dim[0]*refined_dim[1]*(refined_dim[2]+1))
924 + ((refined_dim[0]+1)*refined_dim[1]*refined_dim[2])
925 + (j*refined_dim[0]*refined_dim[2]) + (i*refined_dim[2]) + k,
926 // index face '2' left
927 (refined_dim[0]*refined_dim[1]*(refined_dim[2]+1))
928 + (i*refined_dim[1]*refined_dim[2]) + (k*refined_dim[1]) + j,
929 // index face '3' right
930 (refined_dim[0]*refined_dim[1]*(refined_dim[2]+1))
931 + ((i+1)*refined_dim[1]*refined_dim[2]) + (k*refined_dim[1]) + j,
932 // index face '4' back
933 (refined_dim[0]*refined_dim[1]*(refined_dim[2]+1)) +
934 ((refined_dim[0]+1)*refined_dim[1]*refined_dim[2])
935 + ((j+1)*refined_dim[0]*refined_dim[2]) + (i*refined_dim[2]) + k,
936 // index face '5' top
937 ((k+1)*refined_dim[0]*refined_dim[1]) + (j*refined_dim[0]) + i};
938 //
939 // We add the 6 faces of the cell into "refined_cell_to_face".
940 using cpgrid::EntityRep;
941 // First value is index, Second is orientation.
942 // Still have to find out what the orientation should be.
943 // right face ('3') outer normal points 'from left to right' -> orientation true
944 // back face ('4') outer normal points 'from front to back' -> orientation true
945 // top face ('5') outer normal points 'from bottom to top' -> orientation true
946 // (the other cases are false)
947 std::vector<cpgrid::EntityRep<1>> cell_to_face = {
948 { hexa_to_face[0], false}, {hexa_to_face[1], false},
949 { hexa_to_face[2], false}, {hexa_to_face[3], true},
950 { hexa_to_face[4], true}, {hexa_to_face[5], true} };
951 refined_cell_to_face.appendRow(cell_to_face.begin(), cell_to_face.end());
952 //
953 // Vol2. CENTROIDS of the faces of the hexahedron.
954 // (one of the 6 corners of all 4 tetrahedra based on that face).
955 std::vector<Geometry<0,3>::GlobalCoordinate> hexa_face_centroids;
956 for (auto& idx : hexa_to_face) {
957 hexa_face_centroids.push_back(refined_faces[idx].center());
958 }
959 // Indices of the 4 edges of each face of the hexahedron.
960 // A tetrahedron has six edges. Once we choose a face to base a
961 // tetrahedron on, we choose an edge of that face as one of the
962 // edges of the tetrahedron. The other five edges are fixed, since
963 // the center of the hexahedron and the center of the face are
964 // the reminder 2 corners of the tetrahedron. That's why to identify
965 // a tetrahedron we only need two things:
966 // (1) the face it's based on and
967 // (2) one of the four edges of that face.
968 //
969 // Container with 6 entries, one per face. Each entry has the
970 // 4 indices of the 4 corners of each face.
971 std::vector<std::array<int,4>> cell_face4corners;
972 cell_face4corners.reserve(6);
973 for (int face = 0; face < 6; ++face) {
974 cell_face4corners.push_back({
975 refined_face_to_point[hexa_to_face[face]][0],
976 refined_face_to_point[hexa_to_face[face]][1],
977 refined_face_to_point[hexa_to_face[face]][2],
978 refined_face_to_point[hexa_to_face[face]][3] });
979 }
980 // Vol4. Container with indices of the edges of the 4 tetrahedra per face
981 // [according to description above]
982 std::vector<std::vector<std::array<int,2>>> tetra_edge_indices;
983 tetra_edge_indices.reserve(6);
984 for (auto& face_indices : cell_face4corners)
985 {
986 std::vector<std::array<int,2>> face4edges_indices = {
987 { face_indices[0], face_indices[1]}, // fake '{0,1}'/'{4,5}'
988 { face_indices[0], face_indices[2]}, // fake '{0,2}'/'{4,6}'
989 { face_indices[1], face_indices[3]}, // fake '{1,3}'/'{5,7}'
990 { face_indices[2], face_indices[3]} }; // fake '{2,3}'/'{6,7}'
991 tetra_edge_indices.push_back(face4edges_indices);
992 }
993 // Sum of the 24 volumes to get the volume of the hexahedron,
994 // stored in "refined_cell_volume".
995 // Calculate the volume of each hexahedron, by adding
996 // the 4 tetrahedra at each face (4x6 = 24 tetrahedra).
997 for (int face = 0; face < 6; ++face) {
998 for (int edge = 0; edge < 4; ++edge) {
999 // Construction of each tetrahedron based on "face" with one
1000 // of its edges equal to "edge".
1001 const Geometry<0, 3>::GlobalCoordinate tetra_corners[4] = {
1002 refined_corners[tetra_edge_indices[face][edge][0]].center(), // (see Vol4.)
1003 refined_corners[tetra_edge_indices[face][edge][1]].center(), // (see Vol4.)
1004 hexa_face_centroids[face], // (see Vol2.)
1005 refined_cell_center }; // (see Vol3.)
1006 refined_cell_volume += std::fabs(simplex_volume(tetra_corners));
1007 } // end edge-for-loop
1008 } // end face-for-loop
1009 // Add the volume of the hexahedron (global refined 'cell')
1010 // to the container with of all volumes of all the refined cells.
1011 sum_all_refined_cell_volumes += refined_cell_volume;
1012 // Create a pointer to the first element of "refined_cell_to_point"
1013 // (required as the fourth argement to construct a Geometry<3,3> type object).
1014 int* indices_storage_ptr = refined_cell_to_point[refined_cell_idx].data();
1015 // Construct the Geometry of the refined cell associated with 'kji'.
1016 refined_cells[refined_cell_idx] =
1017 Geometry<3,cdim>(refined_cell_center,
1018 refined_cell_volume,
1019 all_geom.geomVector(std::integral_constant<int,3>()),
1020 indices_storage_ptr);
1021 } // end i-for-loop
1022 } // end j-for-loop
1023 } // end k-for-loop
1024 // Rescale all volumes if the sum of volume of all the global refined 'cells' does not match the
1025 // volume of the 'parent cell'.
1026 // Compare the sum of all the volumes of all refined cells with 'parent cell' volume.
1027 if (std::fabs(sum_all_refined_cell_volumes - this->volume())
1028 > std::numeric_limits<Geometry<3, cdim>::ctype>::epsilon()) {
1029 Geometry<3, cdim>::ctype correction = this->volume() / sum_all_refined_cell_volumes;
1030 for(auto& cell: refined_cells){
1031 cell.vol_ *= correction;
1032 }
1033 } // end if-statement
1035 }
1036
1037 private:
1038 GlobalCoordinate pos_;
1039 double vol_;
1040 std::shared_ptr<const EntityVariable<Geometry<0, 3>,3>> allcorners_; // For dimension 3 only
1041 const int* cor_idx_; // For dimension 3 only
1042
1056 const std::tuple< enum face_tag, int,
1057 std::array<int, 4>, std::vector<cpgrid::EntityRep<0>>,
1058 LocalCoordinate>
1059 getIndicesFace(int l, int m, int n, int constant_direction, const std::array<int, 3>& cells_per_dim) const
1060 {
1061 using cpgrid::EntityRep;
1062 std::vector<cpgrid::EntityRep<0>> neighboring_cells_of_one_face; // {index, orientation}
1063 switch(constant_direction) {
1064 case 0: // {l,m,n} = {k,j,i}, constant in z-direction
1065 // Orientation true when outer normal points 'from bottom to top'
1066 // Orientation false when outer normal points 'from top to bottom'
1067 if (l != 0) {
1068 neighboring_cells_of_one_face.push_back({((l-1)*cells_per_dim[0]*cells_per_dim[1])
1069 + (m*cells_per_dim[0]) + n, true});
1070 }
1071 if (l != cells_per_dim[2]) {
1072 neighboring_cells_of_one_face.push_back({ (l*cells_per_dim[0]*cells_per_dim[1])
1073 + (m*cells_per_dim[0]) + n, false});
1074 }
1075 return { face_tag::K_FACE, (l*cells_per_dim[0]*cells_per_dim[1]) + (m*cells_per_dim[0]) + n,
1076 {(m*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (n*(cells_per_dim[2]+1)) +l,
1077 (m*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((n+1)*(cells_per_dim[2]+1)) +l,
1078 ((m+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (n*(cells_per_dim[2]+1)) +l,
1079 ((m+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((n+1)*(cells_per_dim[2]+1)) +l},
1080 neighboring_cells_of_one_face,
1081 {(.5 + n)/cells_per_dim[0], (.5 + m)/cells_per_dim[1], double(l)/cells_per_dim[2]}};
1082 case 1: // {l,m,n} = {i,k,j}, constant in the x-direction
1083 // Orientation true when outer normal points 'from left to right'
1084 // Orientation false when outer normal points 'from right to left'
1085 if (l != 0) {
1086 neighboring_cells_of_one_face.push_back({(m*cells_per_dim[0]*cells_per_dim[1])
1087 + (n*cells_per_dim[0]) +l-1, true});
1088 }
1089 if (l != cells_per_dim[0]) {
1090 neighboring_cells_of_one_face.push_back({ (m*cells_per_dim[0]*cells_per_dim[1])
1091 + (n*cells_per_dim[0]) + l, false});
1092 }
1093 return { face_tag::I_FACE, (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2]+1))
1094 + (l*cells_per_dim[1]*cells_per_dim[2]) + (m*cells_per_dim[1]) + n,
1095 {(n*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (l*(cells_per_dim[2]+1)) +m,
1096 ((n+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (l*(cells_per_dim[2]+1)) +m,
1097 (n*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (l*(cells_per_dim[2]+1)) +m+1,
1098 ((n+1)*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (l*(cells_per_dim[2]+1)) +m+1},
1099 neighboring_cells_of_one_face,
1100 { double(l)/cells_per_dim[0], (.5 + n)/cells_per_dim[1], (.5 + m)/cells_per_dim[2]}};
1101 case 2: // {l,m,n} = {j,i,k}, constant in the y-direction
1102 // Orientation true when outer normal points 'from front to back'
1103 // Orientation false when outer normal points 'from back to front'
1104 if (l != 0) {
1105 neighboring_cells_of_one_face.push_back({(n*cells_per_dim[0]*cells_per_dim[1])
1106 + ((l-1)*cells_per_dim[0]) +m, true});
1107 }
1108 if (l != cells_per_dim[1]) {
1109 neighboring_cells_of_one_face.push_back({(n*cells_per_dim[0]*cells_per_dim[1])
1110 + (l*cells_per_dim[0]) + m, false});
1111 }
1112 return { face_tag::J_FACE, (cells_per_dim[0]*cells_per_dim[1]*(cells_per_dim[2] +1))
1113 + ((cells_per_dim[0]+1)*cells_per_dim[1]*cells_per_dim[2])
1114 + (l*cells_per_dim[0]*cells_per_dim[2]) + (m*cells_per_dim[2]) + n,
1115 {(l*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (m*(cells_per_dim[2]+1)) +n,
1116 (l*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((m+1)*(cells_per_dim[2]+1)) +n,
1117 (l*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + (m*(cells_per_dim[2]+1)) +n+1,
1118 (l*(cells_per_dim[0]+1)*(cells_per_dim[2]+1)) + ((m+1)*(cells_per_dim[2]+1)) +n+1},
1119 neighboring_cells_of_one_face,
1120 {(.5 + m)/cells_per_dim[0], double(l)/cells_per_dim[1], (.5 + n)/cells_per_dim[2]}};
1121 default:
1122 // Should never be reached, but prevents compiler warning
1123 OPM_THROW(std::logic_error, "Unhandled dimension. This should never happen!");
1124 }
1125 }
1126 };
1127 } // namespace cpgrid
1128
1129 template< int mydim, int cdim >
1130 auto referenceElement(const cpgrid::Geometry<mydim,cdim>& geo) -> decltype(referenceElement<double,mydim>(geo.type()))
1131 {
1132 return referenceElement<double,mydim>(geo.type());
1133 }
1134
1135} // namespace Dune
1136
1137#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:272
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:571
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:593
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:506
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:493
void refineCellifiedPatch(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 std::array< int, 3 > &patch_dim, const std::vector< double > &widthsX, const std::vector< double > &lengthsY, const std::vector< double > &heightsZ) const
Definition Geometry.hpp:617
double integrationElement(const LocalCoordinate &local_coord) const
Equal to \sqrt{\det{J^T J}} where J is the Jacobian.
Definition Geometry.hpp:485
const JacobianTransposed jacobianTransposed(const LocalCoordinate &local_coord) const
Jacobian transposed.
Definition Geometry.hpp:535
GlobalCoordinate global(const LocalCoordinate &local_coord) const
Provide a trilinear mapping.
Definition Geometry.hpp:429
Geometry()
Default constructor, giving a non-valid geometry.
Definition Geometry.hpp:419
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:460
int corners() const
The number of corners of this convex polytope.
Definition Geometry.hpp:500
ctype volume() const
Cell volume.
Definition Geometry.hpp:513
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition Geometry.hpp:523
Jacobian jacobian(const LocalCoordinate &local_coord) const
The jacobian.
Definition Geometry.hpp:580
Dune::FieldVector< double, 3 > PointType
Refine a single cell considering different widths, lengths, and heights.
Definition Geometry.hpp:616
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:587
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:307
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