My Project
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 <memory>
7 
8 #include <dune/common/version.hh>
9 #include <dune/common/fmatrix.hh>
10 #include <dune/grid/common/geometry.hh>
11 
12 #include <dune/geometry/referenceelements.hh>
13 #include <dune/geometry/type.hh>
14 #include <dune/geometry/multilineargeometry.hh>
15 
16 
17 namespace Dune
18 {
19 
20  // Internal Forward Declarations
21  // -----------------------------
22 
23  template< int, int, class > class PolyhedralGridGeometry;
24  template< int, int, class > class PolyhedralGridLocalGeometry;
25 
26  // PolyhedralGridBasicGeometry
27  // -------------------
28 
29  template< int mydim, int cdim, class Grid >
31  {
32  static const int dimension = Grid::dimension;
33  static const int mydimension = mydim;
34  static const int codimension = dimension - mydimension;
35 
36  static const int dimensionworld = Grid::dimensionworld;
37  static const int coorddimension = dimensionworld;
38 
39  typedef typename Grid::ctype ctype;
40  typedef Dune::FieldVector< ctype, coorddimension > GlobalCoordinate;
41  typedef Dune::FieldVector< ctype, mydimension > LocalCoordinate;
42 
43  typedef typename Grid::Traits::ExtraData ExtraData;
44  typedef typename Grid::Traits::template Codim<codimension>::EntitySeed EntitySeed;
45 
46  template< class ct >
48  : public Dune::MultiLinearGeometryTraits< ct >
49  {
50  struct Storage
51  {
52  struct Iterator
53  : public Dune::ForwardIteratorFacade< Iterator, GlobalCoordinate, GlobalCoordinate >
54  {
55  const Storage* data_;
56  int count_;
57  explicit Iterator( const Storage* ptr, int count ) : data_( ptr ), count_( count ) {}
58 
59  GlobalCoordinate dereference() const { return data_->corner( count_ ); }
60  void increment() { ++count_; }
61 
62  bool equals( const Iterator& other ) const { return count_ == other.count_; }
63  };
64 
65  ExtraData data_;
66  // host geometry object
67  EntitySeed seed_;
68 
69  Storage( ExtraData data, EntitySeed seed )
70  : data_( data ), seed_( seed )
71  {}
72 
73  Storage( ExtraData data )
74  : data_( data ), seed_()
75  {}
76 
77  ExtraData data() const { return data_; }
78  bool isValid () const { return seed_.isValid(); }
79 
80  GlobalCoordinate operator [] (const int i) const { return corner( i ); }
81 
82  Iterator begin() const { return Iterator(this, 0); }
83  Iterator end () const { return Iterator(this, corners()); }
84 
85  int corners () const { return data()->corners( seed_ ); }
86  GlobalCoordinate corner ( const int i ) const { return data()->corner( seed_, i ); }
87  GlobalCoordinate center () const { return data()->centroids( seed_ ); }
88 
89  ctype volume() const { return data()->volumes( seed_ ); }
90 
91  const EntitySeed& seed () const { return seed_; }
92  };
93 
94  template <int mdim, int cordim>
96  {
97  typedef Storage Type;
98  };
99  };
100 
101  typedef Dune::MultiLinearGeometry< ctype, mydimension, coorddimension, PolyhedralMultiLinearGeometryTraits<ctype> >
102  MultiLinearGeometryType;
103 
104  typedef typename PolyhedralMultiLinearGeometryTraits< ctype > ::template
105  CornerStorage<mydimension, coorddimension >::Type CornerStorageType;
106 
108  typedef FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed;
109 
111  typedef FieldMatrix< ctype, mydim, cdim > JacobianTransposed;
112 
114  using JacobianInverse = FieldMatrix< ctype, mydim, cdim >;
115 
117  using Jacobian = FieldMatrix< ctype, cdim, mydim >;
118 
119  typedef Dune::Impl::FieldMatrixHelper< ctype > MatrixHelperType;
120 
121  explicit PolyhedralGridBasicGeometry ( ExtraData data )
122  : storage_( data )
123  {}
124 
125  PolyhedralGridBasicGeometry ( ExtraData data, const EntitySeed& seed )
126  : storage_( data, seed )
127  {
128  GeometryType myType = type();
129  if( ! myType.isNone() && storage_.isValid() )
130  {
131  geometryImpl_.reset( new MultiLinearGeometryType(myType, storage_) );
132  }
133  //std::cout << myType << " " << storage_.corners() << std::endl;
134  }
135 
136  GeometryType type () const { return data()->geometryType( storage_.seed() ); }
137  bool affine () const { return (geometryImpl_) ? geometryImpl_->affine() : false; }
138 
139  int corners () const { return storage_.corners(); }
140  GlobalCoordinate corner ( const int i ) const { return storage_.corner( i ); }
141  GlobalCoordinate center () const
142  {
143  if( type().isNone() )
144  {
145  return storage_.center();
146  }
147  else
148  {
149 #if DUNE_VERSION_NEWER(DUNE_GRID,2,7)
150  const auto refElem = Dune::referenceElement< ctype, mydim > ( type() );
151 #else
152  const auto& refElem = Dune::ReferenceElements< ctype, mydim >::general( type() );
153 #endif
154  return global( refElem.position(0,0) );
155  }
156  }
157 
158  GlobalCoordinate global(const LocalCoordinate& local) const
159  {
160  if( geometryImpl_ )
161  {
162  return geometryImpl_->global( local );
163  }
164 
165  return center();
166  }
167 
170  LocalCoordinate local(const GlobalCoordinate& global) const
171  {
172  if( geometryImpl_ )
173  {
174  return geometryImpl_->local( global );
175  }
176 
177  // if no geometry type return a vector filled with 1
178  return LocalCoordinate( 1 );
179  }
180 
181  ctype integrationElement ( const LocalCoordinate &local ) const
182  {
183  if( geometryImpl_ )
184  {
185  return geometryImpl_->integrationElement( local );
186  }
187  return volume();
188  }
189 
190  ctype volume () const
191  {
192  if( geometryImpl_ )
193  {
194  return geometryImpl_->volume();
195  }
196  return storage_.volume();
197  }
198 
199  JacobianTransposed jacobianTransposed ( const LocalCoordinate & local ) const
200  {
201  if( geometryImpl_ )
202  {
203  return geometryImpl_->jacobianTransposed( local );
204  }
205 
206  DUNE_THROW(NotImplemented,"jacobianTransposed not implemented");
207  return JacobianTransposed( 0 );
208  }
209 
210  JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate & local ) const
211  {
212  if( geometryImpl_ )
213  {
214  return geometryImpl_->jacobianInverseTransposed( local );
215  }
216 
217  DUNE_THROW(NotImplemented,"jacobianInverseTransposed not implemented");
218  return JacobianInverseTransposed( 0 );
219  }
220 
221 
223  Jacobian
224  jacobian(const LocalCoordinate& local ) const
225  {
226  return jacobianTransposed(local).transposed();
227  }
228 
231  jacobianInverse(const LocalCoordinate& local) const
232  {
233  return jacobianInverseTransposed(local).transposed();
234  }
235 
236  ExtraData data() const { return storage_.data(); }
237 
238  protected:
239  CornerStorageType storage_;
240  std::shared_ptr< MultiLinearGeometryType > geometryImpl_;
241  };
242 
243 
244  // PolyhedralGridGeometry
245  // --------------
246 
247  template< int mydim, int cdim, class Grid >
249  : public PolyhedralGridBasicGeometry< mydim, cdim, Grid >
250  {
252 
253  public:
254  typedef typename Base::ExtraData ExtraData;
255  typedef typename Base::EntitySeed EntitySeed;
256 
257  explicit PolyhedralGridGeometry ( ExtraData data )
258  : Base( data )
259  {}
260 
261  PolyhedralGridGeometry ( ExtraData data, const EntitySeed& seed )
262  : Base( data, seed )
263  {}
264  };
265 
266  template< int mydim, int cdim, class Grid >
268  : public PolyhedralGridBasicGeometry< mydim, cdim, Grid >
269  {
271 
272  public:
273  typedef typename Base::ExtraData ExtraData;
274 
275  explicit PolyhedralGridLocalGeometry ( ExtraData data )
276  : Base( data )
277  {}
278  };
279 
280 
281 } // namespace Dune
282 
283 #endif // #ifndef DUNE_POLYHEDRALGRID_GEOMETRY_HH
Definition: geometry.hh:250
Definition: geometry.hh:269
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
Definition: geometry.hh:31
FieldMatrix< ctype, mydim, cdim > JacobianInverse
type of jacobian inverse transposed
Definition: geometry.hh:114
Jacobian jacobian(const LocalCoordinate &local) const
The jacobian.
Definition: geometry.hh:224
FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed
type of jacobian inverse transposed
Definition: geometry.hh:108
FieldMatrix< ctype, cdim, mydim > Jacobian
type of jacobian transposed
Definition: geometry.hh:117
FieldMatrix< ctype, mydim, cdim > JacobianTransposed
type of jacobian transposed
Definition: geometry.hh:111
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
The inverse of the jacobian.
Definition: geometry.hh:231
LocalCoordinate local(const GlobalCoordinate &global) const
Mapping from the cell to the reference domain.
Definition: geometry.hh:170