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 
45 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 5 )
46 #include <dune/geometry/type.hh>
47 #else
48 #include <dune/geometry/genericgeometry/geometrytraits.hh>
49 #include <dune/geometry/genericgeometry/matrixhelper.hh>
50 #endif
51 
52 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
53 
54 #include <opm/grid/utility/ErrorMacros.hpp>
55 
56 namespace Dune
57 {
58  namespace cpgrid
59  {
60 
69  template <int mydim, int cdim>
70  class Geometry
71  {
72  };
73 
74 
75 
76 
78  template <int cdim>
79  class Geometry<3, cdim>
80  {
81  static_assert(cdim == 3, "");
82  public:
84  enum { dimension = 3 };
86  enum { mydimension = 3 };
88  enum { coorddimension = cdim };
90  enum { dimensionworld = 3 };
91 
93  typedef double ctype;
94 
96  typedef FieldVector<ctype, mydimension> LocalCoordinate;
98  typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
99 
101  typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
103  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
105  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
106 
107 #if DUNE_VERSION_NEWER(DUNE_GRID,2,5)
108  typedef Dune::Impl::FieldMatrixHelper< double > MatrixHelperType;
109 #else
110  typedef Dune::GenericGeometry::MatrixHelper< Dune::GenericGeometry::DuneCoordTraits<double> > MatrixHelperType;
111 #endif
112 
122  ctype vol,
123  const GlobalCoordinate* allcorners,
124  const int* corner_indices)
125  : pos_(pos), vol_(vol), allcorners_(allcorners), cor_idx_(corner_indices)
126  {
127  assert(allcorners && corner_indices);
128  }
129 
140  ctype vol)
141  : pos_(pos), vol_(vol)
142  {
143  }
144 
147  : pos_(0.0), vol_(0.0), allcorners_(0), cor_idx_(0)
148  {
149  }
150 
155  GlobalCoordinate global(const LocalCoordinate& local_coord) const
156  {
157  static_assert(mydimension == 3, "");
158  static_assert(coorddimension == 3, "");
159  // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
160  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
161  uvw[0] -= local_coord;
162  // Access pattern for uvw matching ordering of corners.
163  const int pat[8][3] = { { 0, 0, 0 },
164  { 1, 0, 0 },
165  { 0, 1, 0 },
166  { 1, 1, 0 },
167  { 0, 0, 1 },
168  { 1, 0, 1 },
169  { 0, 1, 1 },
170  { 1, 1, 1 } };
171  GlobalCoordinate xyz(0.0);
172  for (int i = 0; i < 8; ++i) {
173  GlobalCoordinate corner_contrib = corner(i);
174  double factor = 1.0;
175  for (int j = 0; j < 3; ++j) {
176  factor *= uvw[pat[i][j]][j];
177  }
178  corner_contrib *= factor;
179  xyz += corner_contrib;
180  }
181  return xyz;
182  }
183 
187  {
188  static_assert(mydimension == 3, "");
189  static_assert(coorddimension == 3, "");
190  // This code is modified from dune/grid/genericgeometry/mapping.hh
191  // \todo: Implement direct computation.
192  const ctype epsilon = 1e-12;
193  const ReferenceElement< ctype , 3 > & refElement =
194  ReferenceElements< ctype, 3 >::general(type());
195  LocalCoordinate x = refElement.position(0,0);
196  LocalCoordinate dx;
197  do {
198  // DF^n dx^n = F^n, x^{n+1} -= dx^n
199  JacobianTransposed JT = jacobianTransposed(x);
200  GlobalCoordinate z = global(x);
201  z -= y;
202  MatrixHelperType::template xTRightInvA<3, 3>(JT, z, dx );
203  x -= dx;
204  } while (dx.two_norm2() > epsilon*epsilon);
205  return x;
206  }
207 
212  double integrationElement(const LocalCoordinate& local_coord) const
213  {
214  JacobianTransposed Jt = jacobianTransposed(local_coord);
215  return MatrixHelperType::template sqrtDetAAT<3, 3>(Jt);
216  }
217 
220  GeometryType type() const
221  {
222  GeometryType t;
223  t.makeCube(mydimension);
224  return t;
225  }
226 
229  int corners() const
230  {
231  return 8;
232  }
233 
235  GlobalCoordinate corner(int cor) const
236  {
237  assert(allcorners_ && cor_idx_);
238  return allcorners_[cor_idx_[cor]];
239  }
240 
242  ctype volume() const
243  {
244  return vol_;
245  }
246 
248  const GlobalCoordinate& center() const
249  {
250  return pos_;
251  }
252 
257  const JacobianTransposed
258  jacobianTransposed(const LocalCoordinate& local_coord) const
259  {
260  static_assert(mydimension == 3, "");
261  static_assert(coorddimension == 3, "");
262 
263  // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
264  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local_coord };
265  uvw[0] -= local_coord;
266  // Access pattern for uvw matching ordering of corners.
267  const int pat[8][3] = { { 0, 0, 0 },
268  { 1, 0, 0 },
269  { 0, 1, 0 },
270  { 1, 1, 0 },
271  { 0, 0, 1 },
272  { 1, 0, 1 },
273  { 0, 1, 1 },
274  { 1, 1, 1 } };
275  JacobianTransposed Jt(0.0);
276  for (int i = 0; i < 8; ++i) {
277  for (int deriv = 0; deriv < 3; ++deriv) {
278  // This part contributing to dg/du_{deriv}
279  double factor = 1.0;
280  for (int j = 0; j < 3; ++j) {
281  factor *= (j != deriv) ? uvw[pat[i][j]][j]
282  : (pat[i][j] == 0 ? -1.0 : 1.0);
283  }
284  GlobalCoordinate corner_contrib = corner(i);
285  corner_contrib *= factor;
286  Jt[deriv] += corner_contrib; // using FieldMatrix row access.
287  }
288  }
289  return Jt;
290  }
291 
293  const JacobianInverseTransposed
294  jacobianInverseTransposed(const LocalCoordinate& local_coord) const
295  {
296  JacobianInverseTransposed Jti = jacobianTransposed(local_coord);
297  Jti.invert();
298  return Jti;
299  }
300 
302  bool affine() const
303  {
304  return false;
305  }
306 
307  private:
308  GlobalCoordinate pos_;
309  double vol_;
310  const GlobalCoordinate* allcorners_; // For dimension 3 only
311  const int* cor_idx_; // For dimension 3 only
312  };
313 
314 
315 
316 
317 
320  template <int cdim> // GridImp arg never used
321  class Geometry<2, cdim>
322  {
323  static_assert(cdim == 3, "");
324  public:
326  enum { dimension = 3 };
328  enum { mydimension = 2 };
330  enum { coorddimension = cdim };
332  enum { dimensionworld = 3 };
333 
335  typedef double ctype;
336 
338  typedef FieldVector<ctype, mydimension> LocalCoordinate;
340  typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
341 
343  typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
345  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
347  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
348 
353  ctype vol)
354  : pos_(pos), vol_(vol)
355  {
356  }
357 
360  : pos_(0.0), vol_(0.0)
361  {
362  }
363 
366  {
367  OPM_THROW(std::runtime_error, "Geometry::global() meaningless on singular geometry.");
368  }
369 
372  {
373  OPM_THROW(std::runtime_error, "Geometry::local() meaningless on singular geometry.");
374  }
375 
378  double integrationElement(const LocalCoordinate&) const
379  {
380  return vol_;
381  }
382 
384  GeometryType type() const
385  {
386  GeometryType t;
387  t.makeNone(mydimension);
388  return t;
389  }
390 
393  int corners() const
394  {
395  return 0;
396  }
397 
399  GlobalCoordinate corner(int /* cor */) const
400  {
401  // Meaningless call to cpgrid::Geometry::corner(int):
402  //"singular geometry has no corners.
403  // But the DUNE tests assume at least one corner.
404  return GlobalCoordinate( 0.0 );
405  }
406 
408  ctype volume() const
409  {
410  return vol_;
411  }
412 
414  const GlobalCoordinate& center() const
415  {
416  return pos_;
417  }
418 
420  const FieldMatrix<ctype, mydimension, coorddimension>&
421  jacobianTransposed(const LocalCoordinate& /* local */) const
422  {
423  OPM_THROW(std::runtime_error, "Meaningless to call jacobianTransposed() on singular geometries.");
424  }
425 
427  const FieldMatrix<ctype, coorddimension, mydimension>&
429  {
430  OPM_THROW(std::runtime_error, "Meaningless to call jacobianInverseTransposed() on singular geometries.");
431  }
432 
434  bool affine() const
435  {
436  return true;
437  }
438 
439  private:
440  GlobalCoordinate pos_;
441  ctype vol_;
442  };
443 
444 
445 
446 
448  template <int cdim> // GridImp arg never used
449  class Geometry<0, cdim>
450  {
451  static_assert(cdim == 3, "");
452  public:
454  enum { dimension = 3 };
456  enum { mydimension = 0};
458  enum { coorddimension = cdim };
460  enum { dimensionworld = 3 };
461 
463  typedef double ctype;
464 
466  typedef FieldVector<ctype, mydimension> LocalCoordinate;
468  typedef FieldVector<ctype, coorddimension> GlobalCoordinate;
469 
471  typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian;
473  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
475  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
476 
477 
481  : pos_(pos)
482  {
483  }
484 
487  : pos_(0.0)
488  {
489  }
490 
493  {
494  return pos_;
495  }
496 
499  {
500  // return 0 to make the geometry check happy.
501  return LocalCoordinate(0.0);
502  }
503 
505  double integrationElement(const LocalCoordinate&) const
506  {
507  return volume();
508  }
509 
511  GeometryType type() const
512  {
513  GeometryType t;
514  t.makeCube(mydimension);
515  return t;
516  }
517 
519  int corners() const
520  {
521  return 1;
522  }
523 
525  GlobalCoordinate corner(int cor) const
526  {
527  static_cast<void>(cor);
528  assert(cor == 0);
529  return pos_;
530  }
531 
533  ctype volume() const
534  {
535  return 1.0;
536  }
537 
539  const GlobalCoordinate& center() const
540  {
541  return pos_;
542  }
543 
545  FieldMatrix<ctype, mydimension, coorddimension>
546  jacobianTransposed(const LocalCoordinate& /* local */) const
547  {
548 
549  // Meaningless to call jacobianTransposed() on singular geometries. But we need to make DUNE happy.
550  return FieldMatrix<ctype, mydimension, coorddimension>();
551  }
552 
554  FieldMatrix<ctype, coorddimension, mydimension>
556  {
557  // Meaningless to call jacobianInverseTransposed() on singular geometries. But we need to make DUNE happy.
558  return FieldMatrix<ctype, coorddimension, mydimension>();
559  }
560 
562  bool affine() const
563  {
564  return true;
565  }
566 
567  private:
568  GlobalCoordinate pos_;
569  };
570 
571 
572 
573 
574 
575  } // namespace cpgrid
576 } // namespace Dune
577 
578 #endif // OPM_GEOMETRY_HEADER
const GlobalCoordinate & global(const LocalCoordinate &) const
Returns the position of the vertex.
Definition: Geometry.hpp:492
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:343
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:220
ctype volume() const
Volume of vertex is arbitrarily set to 1.
Definition: Geometry.hpp:533
Geometry(const GlobalCoordinate &pos, ctype vol)
Construct from centroid and volume (1- and 0-moments).
Definition: Geometry.hpp:139
int corners() const
The number of corners of this convex polytope.
Definition: Geometry.hpp:229
Geometry(const GlobalCoordinate &pos, ctype vol, const GlobalCoordinate *allcorners, const int *corner_indices)
Construct from centroid, volume (1- and 0-moments) and corners.
Definition: Geometry.hpp:121
double ctype
Coordinate element type.
Definition: Geometry.hpp:463
LocalCoordinate local(const GlobalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:371
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:347
Geometry(const GlobalCoordinate &pos, ctype vol)
Construct from centroid and volume (1- and 0-moments).
Definition: Geometry.hpp:352
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:146
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:471
const JacobianTransposed jacobianTransposed(const LocalCoordinate &local_coord) const
Jacobian transposed.
Definition: Geometry.hpp:258
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:475
Holds the implementation of the CpGrid as a pimple.
Definition: OpmParserIncludes.hpp:42
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:486
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:98
double ctype
Coordinate element type.
Definition: Geometry.hpp:93
This class encapsulates geometry for both vertices, intersections and cells.
Definition: CpGridData.hpp:92
double integrationElement(const LocalCoordinate &) const
Returns 1 for the vertex geometry.
Definition: Geometry.hpp:505
bool affine() const
The mapping implemented by this geometry is constant, therefore affine.
Definition: Geometry.hpp:562
GlobalCoordinate global(const LocalCoordinate &local_coord) const
Provide a trilinear mapping.
Definition: Geometry.hpp:155
double integrationElement(const LocalCoordinate &local_coord) const
Equal to {{J^T J}} where J is the Jacobian.
Definition: Geometry.hpp:212
ctype volume() const
Volume (area, actually) of intersection.
Definition: Geometry.hpp:408
int corners() const
A vertex is defined by a single corner.
Definition: Geometry.hpp:519
const FieldMatrix< ctype, mydimension, coorddimension > & jacobianTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:421
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:466
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:248
Geometry()
Default constructor, giving a non-valid geometry.
Definition: Geometry.hpp:359
const JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local_coord) const
Inverse of Jacobian transposed.
Definition: Geometry.hpp:294
int corners() const
The number of corners of this convex polytope.
Definition: Geometry.hpp:393
const FieldMatrix< ctype, coorddimension, mydimension > & jacobianInverseTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:428
double integrationElement(const LocalCoordinate &) const
For the singular geometry, we return a constant integration element equal to the volume.
Definition: Geometry.hpp:378
FieldMatrix< ctype, mydimension, coorddimension > jacobianTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:546
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:539
LocalCoordinate local(const GlobalCoordinate &y) const
Mapping from the cell to the reference domain.
Definition: Geometry.hpp:186
bool affine() const
The mapping implemented by this geometry is not generally affine.
Definition: Geometry.hpp:302
GeometryType type() const
Using the cube type for vertices.
Definition: Geometry.hpp:511
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:468
GlobalCoordinate corner(int cor) const
Returns the single corner: the vertex itself.
Definition: Geometry.hpp:525
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:338
LocalCoordinate local(const GlobalCoordinate &) const
Meaningless for the vertex geometry.
Definition: Geometry.hpp:498
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:345
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Type of the inverse of the transposed Jacobian matrix.
Definition: Geometry.hpp:105
const GlobalCoordinate & center() const
Returns the centroid of the geometry.
Definition: Geometry.hpp:414
FieldMatrix< ctype, coorddimension, mydimension > jacobianInverseTransposed(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:555
ctype volume() const
Cell volume.
Definition: Geometry.hpp:242
Geometry(const GlobalCoordinate &pos)
Construct from vertex position.
Definition: Geometry.hpp:480
T volume(const Point< T, 3 > *c)
Computes the volume of a 3D simplex (embedded i 3D space).
Definition: Volumes.hpp:137
const GlobalCoordinate & global(const LocalCoordinate &) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:365
bool affine() const
Since integrationElement() is constant, returns true.
Definition: Geometry.hpp:434
GlobalCoordinate corner(int) const
This method is meaningless for singular geometries.
Definition: Geometry.hpp:399
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:473
FieldVector< ctype, mydimension > LocalCoordinate
Domain type of.
Definition: Geometry.hpp:96
GeometryType type() const
We use the singular type (None) for intersections.
Definition: Geometry.hpp:384
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition: Geometry.hpp:103
FieldVector< ctype, coorddimension > GlobalCoordinate
Range type of.
Definition: Geometry.hpp:340
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition: Geometry.hpp:101
double ctype
Coordinate element type.
Definition: Geometry.hpp:335
GlobalCoordinate corner(int cor) const
The 8 corners of the hexahedral base cell.
Definition: Geometry.hpp:235