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 
113 
114  typedef Dune::Impl::FieldMatrixHelper< ctype > MatrixHelperType;
115 
116  explicit PolyhedralGridBasicGeometry ( ExtraData data )
117  : storage_( data )
118  {}
119 
120  PolyhedralGridBasicGeometry ( ExtraData data, const EntitySeed& seed )
121  : storage_( data, seed )
122  {
123  GeometryType myType = type();
124  if( ! myType.isNone() && storage_.isValid() )
125  {
126  geometryImpl_.reset( new MultiLinearGeometryType(myType, storage_) );
127  }
128  //std::cout << myType << " " << storage_.corners() << std::endl;
129  }
130 
131  GeometryType type () const { return data()->geometryType( storage_.seed() ); }
132  bool affine () const { return (geometryImpl_) ? geometryImpl_->affine() : false; }
133 
134  int corners () const { return storage_.corners(); }
135  GlobalCoordinate corner ( const int i ) const { return storage_.corner( i ); }
136  GlobalCoordinate center () const
137  {
138  if( type().isNone() )
139  {
140  return storage_.center();
141  }
142  else
143  {
144 #if DUNE_VERSION_NEWER(DUNE_GRID,2,7)
145  const auto refElem = Dune::referenceElement< ctype, mydim > ( type() );
146 #else
147  const auto& refElem = Dune::ReferenceElements< ctype, mydim >::general( type() );
148 #endif
149  return global( refElem.position(0,0) );
150  }
151  }
152 
153  GlobalCoordinate global(const LocalCoordinate& local) const
154  {
155  if( geometryImpl_ )
156  {
157  return geometryImpl_->global( local );
158  }
159 
160  return center();
161  }
162 
165  LocalCoordinate local(const GlobalCoordinate& global) const
166  {
167  if( geometryImpl_ )
168  {
169  return geometryImpl_->local( global );
170  }
171 
172  // if no geometry type return a vector filled with 1
173  return LocalCoordinate( 1 );
174  }
175 
176  ctype integrationElement ( const LocalCoordinate &local ) const
177  {
178  if( geometryImpl_ )
179  {
180  return geometryImpl_->integrationElement( local );
181  }
182  return volume();
183  }
184 
185  ctype volume () const
186  {
187  if( geometryImpl_ )
188  {
189  return geometryImpl_->volume();
190  }
191  return storage_.volume();
192  }
193 
194  JacobianTransposed jacobianTransposed ( const LocalCoordinate & local ) const
195  {
196  if( geometryImpl_ )
197  {
198  return geometryImpl_->jacobianTransposed( local );
199  }
200 
201  DUNE_THROW(NotImplemented,"jacobianTransposed not implemented");
202  return JacobianTransposed( 0 );
203  }
204 
205  JacobianInverseTransposed jacobianInverseTransposed ( const LocalCoordinate & local ) const
206  {
207  if( geometryImpl_ )
208  {
209  return geometryImpl_->jacobianInverseTransposed( local );
210  }
211 
212  DUNE_THROW(NotImplemented,"jacobianInverseTransposed not implemented");
213  return JacobianInverseTransposed( 0 );
214  }
215 
216  ExtraData data() const { return storage_.data(); }
217 
218  protected:
219  CornerStorageType storage_;
220  std::shared_ptr< MultiLinearGeometryType > geometryImpl_;
221  };
222 
223 
224  // PolyhedralGridGeometry
225  // --------------
226 
227  template< int mydim, int cdim, class Grid >
229  : public PolyhedralGridBasicGeometry< mydim, cdim, Grid >
230  {
232 
233  public:
234  typedef typename Base::ExtraData ExtraData;
235  typedef typename Base::EntitySeed EntitySeed;
236 
237  explicit PolyhedralGridGeometry ( ExtraData data )
238  : Base( data )
239  {}
240 
241  PolyhedralGridGeometry ( ExtraData data, const EntitySeed& seed )
242  : Base( data, seed )
243  {}
244  };
245 
246  template< int mydim, int cdim, class Grid >
248  : public PolyhedralGridBasicGeometry< mydim, cdim, Grid >
249  {
251 
252  public:
253  typedef typename Base::ExtraData ExtraData;
254 
255  explicit PolyhedralGridLocalGeometry ( ExtraData data )
256  : Base( data )
257  {}
258  };
259 
260 
261 } // namespace Dune
262 
263 #endif // #ifndef DUNE_POLYHEDRALGRID_GEOMETRY_HH
Definition: geometry.hh:230
Definition: geometry.hh:249
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
Definition: geometry.hh:31
FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed
type of jacobian inverse transposed
Definition: geometry.hh:108
FieldMatrix< ctype, mydim, cdim > JacobianTransposed
type of jacobian transposed
Definition: geometry.hh:111
LocalCoordinate local(const GlobalCoordinate &global) const
Mapping from the cell to the reference domain.
Definition: geometry.hh:165