geometry.hh
1 // -*- mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=2 sw=2 sts=2:
3 #ifndef DUNE_POLYHEDRALGRID_GEOMETRY_HH
4 #define DUNE_POLYHEDRALGRID_GEOMETRY_HH
5 
6 #include <dune/common/version.hh>
7 #include <dune/common/fmatrix.hh>
8 #include <dune/grid/common/geometry.hh>
9 
10 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 5 )
11 #include <dune/geometry/type.hh>
12 #else
13 #include <dune/geometry/genericgeometry/geometrytraits.hh>
14 #include <dune/geometry/genericgeometry/matrixhelper.hh>
15 #endif
16 
17 
18 namespace Dune
19 {
20 
21  // Internal Forward Declarations
22  // -----------------------------
23 
24  template< int, int, class > class PolyhedralGridGeometry;
25  template< int, int, class > class PolyhedralGridLocalGeometry;
26 
27  // PolyhedralGridBasicGeometry
28  // -------------------
29 
30  template< int mydim, int cdim, class Grid >
32  {
33  static const int dimension = Grid::dimension;
34  static const int mydimension = mydim;
35  static const int codimension = dimension - mydimension;
36 
37  static const int dimensionworld = Grid::dimensionworld;
38  static const int coorddimension = dimensionworld;
39 
40  typedef typename Grid::ctype ctype;
41  typedef Dune::FieldVector< ctype, coorddimension > GlobalCoordinate;
42  typedef Dune::FieldVector< ctype, mydimension > LocalCoordinate;
43 
45  typedef FieldMatrix<ctype,cdim,mydim> JacobianInverseTransposed;
46 
48  typedef FieldMatrix< ctype, mydim, cdim > JacobianTransposed;
49 
50 
51 #if DUNE_VERSION_NEWER(DUNE_GRID,2,5)
52  typedef Dune::Impl::FieldMatrixHelper< double > MatrixHelperType;
53 #else
54  typedef Dune::GenericGeometry::MatrixHelper< Dune::GenericGeometry::DuneCoordTraits<double> > MatrixHelperType;
55 #endif
56 
57  typedef typename Grid::Traits::ExtraData ExtraData;
58  typedef typename Grid::Traits::template Codim<codimension>::EntitySeed EntitySeed;
59 
60  explicit PolyhedralGridBasicGeometry ( ExtraData data )
61  : data_( data ),
62  seed_( )
63  {}
64 
65  PolyhedralGridBasicGeometry ( ExtraData data, const EntitySeed& seed )
66  : data_( data ),
67  seed_( seed )
68  {}
69 
70  GeometryType type () const { return data()->geomTypes(codimension)[0]; }
71  bool affine () const { return false; }
72 
73  int corners () const { return data()->corners( seed_ ); }
74  GlobalCoordinate corner ( const int i ) const { return data()->corner( seed_, i ); }
75  GlobalCoordinate center () const { return data()->centroids( seed_ ); }
76 
77  GlobalCoordinate global(const LocalCoordinate& local) const
78  {
79  const GeometryType geomType = type();
80  if( geomType.isCube() )
81  {
82  assert( mydimension == 3 );
83  assert( coorddimension == 3 );
84 
85  // uvw = { (1-u, 1-v, 1-w), (u, v, w) }
86  LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
87  uvw[0] -= local;
88 
89  //const ReferenceElement< ctype , mydimension > & refElement =
90  // ReferenceElements< ctype, mydimension >::general( type() );
91 
92  // Access pattern for uvw matching ordering of corners.
93  const int pat[8][3] = { { 0, 0, 0 },
94  { 1, 0, 0 },
95  { 0, 1, 0 },
96  { 1, 1, 0 },
97  { 0, 0, 1 },
98  { 1, 0, 1 },
99  { 0, 1, 1 },
100  { 1, 1, 1 } };
101 
102  const int nCorners = corners();
103  //refElement.size( mydimension );
104 
105  GlobalCoordinate xyz(0.0);
106  for (int i = 0; i < nCorners ; ++i)
107  {
108  GlobalCoordinate cornerContrib = corner(i);
109  //LocalCoordinate refCorner = refElement.position(i,mydimension);
110  double factor = 1.0;
111  for (int j = 0; j < mydimension; ++j)
112  {
113  //factor *= uvw[ refCorner[ j ] ][ j ];
114  factor *= uvw[ pat[ i ][ j ] ][ j ];
115  }
116  cornerContrib *= factor;
117  xyz += cornerContrib;
118  }
119  return xyz;
120  }
121  else if ( geomType.isNone() )
122  {
123  // if no geometry type return the center of the element
124  return center();
125  }
126  else
127  {
128  DUNE_THROW(NotImplemented,"global for geometry type " << geomType << " is not supported!");
129  return GlobalCoordinate( 0 );
130  }
131  }
132 
135  LocalCoordinate local(const GlobalCoordinate& y) const
136  {
137  const GeometryType geomType = type();
138  if( geomType.isCube() )
139  {
140  // This code is modified from dune/grid/genericgeometry/mapping.hh
141  // \todo: Implement direct computation.
142  const ctype epsilon = 1e-12;
143 #if DUNE_VERSION_NEWER(DUNE_GRID,2,3)
144  const ReferenceElement< ctype , mydimension > & refElement =
145  ReferenceElements< ctype, mydimension >::general(type());
146 #else
147  const GenericReferenceElement< ctype , mydimension > & refElement =
148  GenericReferenceElements< ctype, mydimension >::general(type());
149 #endif
150 
151  LocalCoordinate x = refElement.position(0,0);
152  LocalCoordinate dx;
153  do {
154  // DF^n dx^n = F^n, x^{n+1} -= dx^n
155  JacobianTransposed JT = jacobianTransposed(x);
156  GlobalCoordinate z = global(x);
157  z -= y;
158  MatrixHelperType::template xTRightInvA<mydimension,coorddimension>(JT, z, dx );
159  x -= dx;
160  } while (dx.two_norm2() > epsilon*epsilon);
161  return x;
162  }
163  else if ( geomType.isNone() )
164  {
165  // if no geometry type return a vector filled with 1
166  return LocalCoordinate( 1 );
167  }
168  else
169  {
170  DUNE_THROW(NotImplemented,"local for geometry type " << geomType << " is not supported!");
171  return LocalCoordinate( 0 );
172  }
173  }
174 
175  ctype integrationElement ( const LocalCoordinate & ) const { return volume(); }
176  ctype volume () const { return data()->volumes( seed_ ); }
177 
178 #if DUNE_VERSION_NEWER(DUNE_GRID,2,4)
179  JacobianTransposed jacobianTransposed ( const LocalCoordinate & ) const
180  {
181  DUNE_THROW(NotImplemented,"jacobianTransposed not implemented");
182  return JacobianTransposed( 0 );
183  }
184 
185  JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate & ) const
186  {
187  DUNE_THROW(NotImplemented,"jacobianInverseTransposed not implemented");
188  return JacobianInverseTransposed( 0 );
189  }
190 #else
191  const JacobianTransposed& jacobianTransposed ( const LocalCoordinate &local ) const
192  {
193  DUNE_THROW(NotImplemented,"jacobianTransposed not implemented");
194  static const JacobianTransposed jac( 0 );
195  return jac;
196  }
197 
198  const JacobianInverseTransposed& jacobianInverseTransposed ( const LocalCoordinate &local ) const
199  {
200  DUNE_THROW(NotImplemented,"jacobianInverseTransposed not implemented");
201  static const JacobianInverseTransposed jac( 0 );
202  return jac;
203  }
204 #endif
205 
206  ExtraData data() const { return data_; }
207 
208  protected:
209  ExtraData data_;
210  // host geometry object
211  EntitySeed seed_;
212  };
213 
214 
215  // PolyhedralGridGeometry
216  // --------------
217 
218  template< int mydim, int cdim, class Grid >
219  class PolyhedralGridGeometry
220  : public PolyhedralGridBasicGeometry< mydim, cdim, Grid >
221  {
222  typedef PolyhedralGridBasicGeometry< mydim, cdim, Grid > Base;
223 
224  public:
225  typedef typename Base::ExtraData ExtraData;
226  typedef typename Base::EntitySeed EntitySeed;
227 
228  explicit PolyhedralGridGeometry ( ExtraData data )
229  : Base( data )
230  {}
231 
232  PolyhedralGridGeometry ( ExtraData data, const EntitySeed& seed )
233  : Base( data, seed )
234  {}
235  };
236 
237  template< int mydim, int cdim, class Grid >
238  class PolyhedralGridLocalGeometry
239  : public PolyhedralGridBasicGeometry< mydim, cdim, Grid >
240  {
241  typedef PolyhedralGridBasicGeometry< mydim, cdim, Grid > Base;
242 
243  public:
244  typedef typename Base::ExtraData ExtraData;
245 
246  explicit PolyhedralGridLocalGeometry ( ExtraData data )
247  : Base( data )
248  {}
249  };
250 
251 
252 #if ! DUNE_VERSION_NEWER(DUNE_GRID,2,4)
253  namespace FacadeOptions
254  {
255 
258  template< int mydim, int cdim, class GridImp >
259  struct StoreGeometryReference< mydim, cdim, GridImp, PolyhedralGridGeometry >
260  {
262  static const bool v = false;
263  };
264 
265  template< int mydim, int cdim, class GridImp >
266  struct StoreGeometryReference< mydim, cdim, GridImp, PolyhedralGridLocalGeometry >
267  {
269  static const bool v = false;
270  };
271 
272  }
273 #endif
274 
275 
276 } // namespace Dune
277 
278 #endif // #ifndef DUNE_POLYHEDRALGRID_GEOMETRY_HH
Definition: geometry.hh:24
LocalCoordinate local(const GlobalCoordinate &y) const
Mapping from the cell to the reference domain.
Definition: geometry.hh:135
Holds the implementation of the CpGrid as a pimple.
Definition: OpmParserIncludes.hpp:42
FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed
type of jacobian inverse transposed
Definition: geometry.hh:45
Definition: geometry.hh:31
FieldMatrix< ctype, mydim, cdim > JacobianTransposed
type of jacobian transposed
Definition: geometry.hh:48
Definition: geometry.hh:25