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 Statoil 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 // Warning suppression for Dune includes.
40 #include <opm/grid/utility/platform_dependent/disable_warnings.h>
41 
42 #include <dune/common/version.hh>
43 #include <dune/geometry/referenceelements.hh>
44 #include <dune/grid/common/geometry.hh>
45 
46 #include <dune/geometry/type.hh>
47 
48 #include <opm/grid/cpgrid/EntityRep.hpp>
49 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
50 
51 #include <opm/grid/utility/ErrorMacros.hpp>
52 
53 namespace Dune
54 {
55  namespace cpgrid
56  {
57 
66  template <int mydim, int cdim>
67  class Geometry
68  {
69  };
70 
71 
72 
73 
75  template <int cdim> // GridImp arg never used
76  class Geometry<0, cdim>
77  {
78  static_assert(cdim == 3, "");
79  public:
81  enum { dimension = 3 };
83  enum { mydimension = 0};
85  enum { coorddimension = cdim };
87  enum { dimensionworld = 3 };
88 
90  typedef double ctype;
91 
93  typedef FieldVector<ctype, mydimension> LocalCoordinate;
95  typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
96 
98  typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
100  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
102  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
103 
104 
108  : pos_(pos)
109  {
110  }
111 
114  : pos_(0.0)
115  {
116  }
117 
120  {
121  return pos_;
122  }
123 
126  {
127  // return 0 to make the geometry check happy.
128  return LocalCoordinate(0.0);
129  }
130 
132  double integrationElement(const LocalCoordinate&) const
133  {
134  return volume();
135  }
136 
138  GeometryType type() const
139  {
140  return Dune::GeometryTypes::cube(mydimension);
141  }
142 
144  int corners() const
145  {
146  return 1;
147  }
148 
150  GlobalCoordinate corner(int cor) const
151  {
152  static_cast<void>(cor);
153  assert(cor == 0);
154  return pos_;
155  }
156 
158  ctype volume() const
159  {
160  return 1.0;
161  }
162 
164  const GlobalCoordinate& center() const
165  {
166  return pos_;
167  }
168 
170  FieldMatrix<ctype, mydimension, coorddimension>
171  jacobianTransposed(const LocalCoordinate& /* local */) const
172  {
173 
174  // Meaningless to call jacobianTransposed() on singular geometries. But we need to make DUNE happy.
175  return FieldMatrix<ctype, mydimension, coorddimension>();
176  }
177 
179  FieldMatrix<ctype, coorddimension, mydimension>
181  {
182  // Meaningless to call jacobianInverseTransposed() on singular geometries. But we need to make DUNE happy.
183  return FieldMatrix<ctype, coorddimension, mydimension>();
184  }
185 
187  bool affine() const
188  {
189  return true;
190  }
191 
192  private:
193  GlobalCoordinate pos_;
194  };
195 
196 
197 
198 
200  template <int cdim>
201  class Geometry<3, cdim>
202  {
203  static_assert(cdim == 3, "");
204  public:
206  enum { dimension = 3 };
208  enum { mydimension = 3 };
210  enum { coorddimension = cdim };
212  enum { dimensionworld = 3 };
213 
215  typedef double ctype;
216 
218  typedef FieldVector<ctype, mydimension> LocalCoordinate;
220  typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
221 
223  typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
225  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
227  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
228 
229  typedef Dune::Impl::FieldMatrixHelper< double > MatrixHelperType;
230 
240  ctype vol,
241  const EntityVariable<cpgrid::Geometry<0, 3>, 3>& allcorners,
242  const int* corner_indices)
243  : pos_(pos), vol_(vol), allcorners_(allcorners.data()), cor_idx_(corner_indices)
244  {
245  assert(allcorners_ && corner_indices);
246  }
247 
258  ctype vol)
259  : pos_(pos), vol_(vol)
260  {
261  }
262 
265  : pos_(0.0), vol_(0.0), allcorners_(0), cor_idx_(0)
266  {
267  }
268 
273  GlobalCoordinate global(const LocalCoordinate& local_coord) const
274  {
275  static_assert(mydimension == 3, "");
276  static_assert(coorddimension == 3, "");
277  // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
278  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
279  uvw[0] -= local_coord;
280  // Access pattern for uvw matching ordering of corners.
281  const int pat[8][3] = { { 0, 0, 0 },
282  { 1, 0, 0 },
283  { 0, 1, 0 },
284  { 1, 1, 0 },
285  { 0, 0, 1 },
286  { 1, 0, 1 },
287  { 0, 1, 1 },
288  { 1, 1, 1 } };
289  GlobalCoordinate xyz(0.0);
290  for (int i = 0; i < 8; ++i) {
291  GlobalCoordinate corner_contrib = corner(i);
292  double factor = 1.0;
293  for (int j = 0; j < 3; ++j) {
294  factor *= uvw[pat[i][j]][j];
295  }
296  corner_contrib *= factor;
297  xyz += corner_contrib;
298  }
299  return xyz;
300  }
301 
305  {
306  static_assert(mydimension == 3, "");
307  static_assert(coorddimension == 3, "");
308  // This code is modified from dune/grid/genericgeometry/mapping.hh
309  // \todo: Implement direct computation.
310  const ctype epsilon = 1e-12;
311  auto refElement = Dune::ReferenceElements<ctype, 3>::cube();
312  LocalCoordinate x = refElement.position(0,0);
313  LocalCoordinate dx;
314  do {
315  // DF^n dx^n = F^n, x^{n+1} -= dx^n
316  JacobianTransposed JT = jacobianTransposed(x);
317  GlobalCoordinate z = global(x);
318  z -= y;
319  MatrixHelperType::template xTRightInvA<3, 3>(JT, z, dx );
320  x -= dx;
321  } while (dx.two_norm2() > epsilon*epsilon);
322  return x;
323  }
324 
329  double integrationElement(const LocalCoordinate& local_coord) const
330  {
331  JacobianTransposed Jt = jacobianTransposed(local_coord);
332  return MatrixHelperType::template sqrtDetAAT<3, 3>(Jt);
333  }
334 
337  GeometryType type() const
338  {
339  return Dune::GeometryTypes::cube(mydimension);
340  }
341 
344  int corners() const
345  {
346  return 8;
347  }
348 
350  GlobalCoordinate corner(int cor) const
351  {
352  assert(allcorners_ && cor_idx_);
353  return allcorners_[cor_idx_[cor]].center();
354  }
355 
357  ctype volume() const
358  {
359  return vol_;
360  }
361 
363  const GlobalCoordinate& center() const
364  {
365  return pos_;
366  }
367 
372  const JacobianTransposed
373  jacobianTransposed(const LocalCoordinate& local_coord) const
374  {
375  static_assert(mydimension == 3, "");
376  static_assert(coorddimension == 3, "");
377 
378  // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
379  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
380  uvw[0] -= local_coord;
381  // Access pattern for uvw matching ordering of corners.
382  const int pat[8][3] = { { 0, 0, 0 },
383  { 1, 0, 0 },
384  { 0, 1, 0 },
385  { 1, 1, 0 },
386  { 0, 0, 1 },
387  { 1, 0, 1 },
388  { 0, 1, 1 },
389  { 1, 1, 1 } };
390  JacobianTransposed Jt(0.0);
391  for (int i = 0; i < 8; ++i) {
392  for (int deriv = 0; deriv < 3; ++deriv) {
393  // This part contributing to dg/du_{deriv}
394  double factor = 1.0;
395  for (int j = 0; j < 3; ++j) {
396  factor *= (j != deriv) ? uvw[pat[i][j]][j]
397  : (pat[i][j] == 0 ? -1.0 : 1.0);
398  }
399  GlobalCoordinate corner_contrib = corner(i);
400  corner_contrib *= factor;
401  Jt[deriv] += corner_contrib; // using FieldMatrix row access.
402  }
403  }
404  return Jt;
405  }
406 
408  const JacobianInverseTransposed
409  jacobianInverseTransposed(const LocalCoordinate& local_coord) const
410  {
411  JacobianInverseTransposed Jti = jacobianTransposed(local_coord);
412  Jti.invert();
413  return Jti;
414  }
415 
417  bool affine() const
418  {
419  return false;
420  }
421 
422  private:
423  GlobalCoordinate pos_;
424  double vol_;
425  const cpgrid::Geometry<0, 3>* allcorners_; // For dimension 3 only
426  const int* cor_idx_; // For dimension 3 only
427  };
428 
429 
430 
431 
432 
435  template <int cdim> // GridImp arg never used
436  class Geometry<2, cdim>
437  {
438  static_assert(cdim == 3, "");
439  public:
441  enum { dimension = 3 };
443  enum { mydimension = 2 };
445  enum { coorddimension = cdim };
447  enum { dimensionworld = 3 };
448 
450  typedef double ctype;
451 
453  typedef FieldVector<ctype, mydimension> LocalCoordinate;
455  typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
456 
458  typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
460  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
462  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
463 
468  ctype vol)
469  : pos_(pos), vol_(vol)
470  {
471  }
472 
475  : pos_(0.0), vol_(0.0)
476  {
477  }
478 
481  {
482  OPM_THROW(std::runtime_error, "Geometry::global() meaningless on singular geometry.");
483  }
484 
487  {
488  OPM_THROW(std::runtime_error, "Geometry::local() meaningless on singular geometry.");
489  }
490 
493  double integrationElement(const LocalCoordinate&) const
494  {
495  return vol_;
496  }
497 
499  GeometryType type() const
500  {
501  return Dune::GeometryTypes::none(mydimension);
502  }
503 
506  int corners() const
507  {
508  return 0;
509  }
510 
512  GlobalCoordinate corner(int /* cor */) const
513  {
514  // Meaningless call to cpgrid::Geometry::corner(int):
515  //"singular geometry has no corners.
516  // But the DUNE tests assume at least one corner.
517  return GlobalCoordinate( 0.0 );
518  }
519 
521  ctype volume() const
522  {
523  return vol_;
524  }
525 
527  const GlobalCoordinate& center() const
528  {
529  return pos_;
530  }
531 
533  const FieldMatrix<ctype, mydimension, coorddimension>&
534  jacobianTransposed(const LocalCoordinate& /* local */) const
535  {
536  OPM_THROW(std::runtime_error, "Meaningless to call jacobianTransposed() on singular geometries.");
537  }
538 
540  const FieldMatrix<ctype, coorddimension, mydimension>&
542  {
543  OPM_THROW(std::runtime_error, "Meaningless to call jacobianInverseTransposed() on singular geometries.");
544  }
545 
547  bool affine() const
548  {
549  return true;
550  }
551 
552  private:
553  GlobalCoordinate pos_;
554  ctype vol_;
555  };
556 
557 
558 
559 
560 
561  } // namespace cpgrid
562 
563  template< int mydim, int cdim >
564  auto referenceElement(const cpgrid::Geometry<mydim,cdim>& geo) -> decltype(referenceElement<double,mydim>(geo.type()))
565  {
566  return referenceElement<double,mydim>(geo.type());
567  }
568 
569 } // namespace Dune
570 
571 #endif // OPM_GEOMETRY_HEADER
A class design to hold a variable with a value for each entity of the given codimension,...
Definition: EntityRep.hpp:264
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:95
bool affine() const
The mapping implemented by this geometry is constant, therefore affine.
Definition: Geometry.hpp:187
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:98
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:93
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:164
double integrationElement(const LocalCoordinate &) const
Returns 1 for the vertex geometry.
Definition: Geometry.hpp:132
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:102
Geometry(const GlobalCoordinate &pos)
Construct from vertex position.
Definition: Geometry.hpp:107
FieldMatrix< ctype, coorddimension, mydimension > jacobianInverseTransposed(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:138
ctype volume() const
Volume of vertex is arbitrarily set to 1.
Definition: Geometry.hpp:158
FieldMatrix< ctype, mydimension, coorddimension > jacobianTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:171
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:100
GlobalCoordinate corner(int cor) const
Returns the single corner: the vertex itself.
Definition: Geometry.hpp:150
int corners() const
A vertex is defined by a single corner.
Definition: Geometry.hpp:144
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:113
LocalCoordinate local(const GlobalCoordinate &) const
Meaningless for the vertex geometry.
Definition: Geometry.hpp:125
const GlobalCoordinate & global(const LocalCoordinate &) const
Returns the position of the vertex.
Definition: Geometry.hpp:119
double ctype
Coordinate element type.
Definition: Geometry.hpp:90
const GlobalCoordinate & global(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:480
int corners() const
The number of corners of this convex polytope.
Definition: Geometry.hpp:506
LocalCoordinate local(const GlobalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:486
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:458
bool affine() const
Since integrationElement() is constant, returns true.
Definition: Geometry.hpp:547
const FieldMatrix< ctype, coorddimension, mydimension > & jacobianInverseTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:541
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:453
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:462
ctype volume() const
Volume (area, actually) of intersection.
Definition: Geometry.hpp:521
GeometryType type() const
We use the singular type (None) for intersections.
Definition: Geometry.hpp:499
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:460
double integrationElement(const LocalCoordinate &) const
For the singular geometry, we return a constant integration element equal to the volume.
Definition: Geometry.hpp:493
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:474
const FieldMatrix< ctype, mydimension, coorddimension > & jacobianTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:534
Geometry(const GlobalCoordinate &pos, ctype vol)
Construct from centroid and volume (1- and 0-moments).
Definition: Geometry.hpp:467
double ctype
Coordinate element type.
Definition: Geometry.hpp:450
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:527
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:455
GlobalCoordinate corner(int) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:512
const JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local_coord) const
Inverse of Jacobian transposed.
Definition: Geometry.hpp:409
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:227
bool affine() const
The mapping implemented by this geometry is not generally affine.
Definition: Geometry.hpp:417
double ctype
Coordinate element type.
Definition: Geometry.hpp:215
Geometry(const GlobalCoordinate &pos, ctype vol, const EntityVariable< cpgrid::Geometry< 0, 3 >, 3 > &allcorners, const int *corner_indices)
Construct from centroid, volume (1- and 0-moments) and corners.
Definition: Geometry.hpp:239
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:220
GlobalCoordinate corner(int cor) const
The 8 corners of the hexahedral base cell.
Definition: Geometry.hpp:350
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:225
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:337
double integrationElement(const LocalCoordinate &local_coord) const
Equal to \sqrt{\det{J^T J}} where J is the Jacobian.
Definition: Geometry.hpp:329
const JacobianTransposed jacobianTransposed(const LocalCoordinate &local_coord) const
Jacobian transposed.
Definition: Geometry.hpp:373
GlobalCoordinate global(const LocalCoordinate &local_coord) const
Provide a trilinear mapping.
Definition: Geometry.hpp:273
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:264
LocalCoordinate local(const GlobalCoordinate &y) const
Mapping from the cell to the reference domain.
Definition: Geometry.hpp:304
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:363
int corners() const
The number of corners of this convex polytope.
Definition: Geometry.hpp:344
ctype volume() const
Cell volume.
Definition: Geometry.hpp:357
Geometry(const GlobalCoordinate &pos, ctype vol)
Construct from centroid and volume (1- and 0-moments).
Definition: Geometry.hpp:257
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:218
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:223
This class encapsulates geometry for both vertices, intersections and cells.
Definition: Geometry.hpp:68
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
T volume(const Point< T, 3 > *c)
Computes the volume of a 3D simplex (embedded i 3D space).
Definition: Volumes.hpp:137