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
17namespace 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
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