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/fmatrix.hh>
9#include <dune/grid/common/geometry.hh>
10
11#include <dune/geometry/referenceelements.hh>
12#include <dune/geometry/type.hh>
13#include <dune/geometry/multilineargeometry.hh>
14
15
16namespace Dune
17{
18
19 // Internal Forward Declarations
20 // -----------------------------
21
22 template< int, int, class > class PolyhedralGridGeometry;
23 template< int, int, class > class PolyhedralGridLocalGeometry;
24
25 // PolyhedralGridBasicGeometry
26 // -------------------
27
28 template< int mydim, int cdim, class Grid >
30 {
31 static const int dimension = Grid::dimension;
32 static const int mydimension = mydim;
33 static const int codimension = dimension - mydimension;
34
35 static const int dimensionworld = Grid::dimensionworld;
36 static const int coorddimension = dimensionworld;
37
38 typedef typename Grid::ctype ctype;
39 typedef Dune::FieldVector< ctype, coorddimension > GlobalCoordinate;
40 typedef Dune::FieldVector< ctype, mydimension > LocalCoordinate;
41
42 typedef typename Grid::Traits::ExtraData ExtraData;
43 typedef typename Grid::Traits::template Codim<codimension>::EntitySeed EntitySeed;
44
45 template< class ct >
47 : public Dune::MultiLinearGeometryTraits< ct >
48 {
49 struct Storage
50 {
51 struct Iterator
52 : public Dune::ForwardIteratorFacade< Iterator, GlobalCoordinate, GlobalCoordinate >
53 {
54 const Storage* data_;
55 int count_;
56 explicit Iterator( const Storage* ptr, int count ) : data_( ptr ), count_( count ) {}
57
58 GlobalCoordinate dereference() const { return data_->corner( count_ ); }
59 void increment() { ++count_; }
60
61 bool equals( const Iterator& other ) const { return count_ == other.count_; }
62 };
63
64 ExtraData data_;
65 // host geometry object
66 EntitySeed seed_;
67
68 Storage( ExtraData data, EntitySeed seed )
69 : data_( data ), seed_( seed )
70 {}
71
72 Storage( ExtraData data )
73 : data_( data ), seed_()
74 {}
75
76 ExtraData data() const { return data_; }
77 bool isValid () const { return seed_.isValid(); }
78
79 GlobalCoordinate operator [] (const int i) const { return corner( i ); }
80
81 Iterator begin() const { return Iterator(this, 0); }
82 Iterator end () const { return Iterator(this, corners()); }
83
84 int corners () const { return data()->corners( seed_ ); }
85 GlobalCoordinate corner ( const int i ) const { return data()->corner( seed_, i ); }
86 GlobalCoordinate center () const { return data()->centroids( seed_ ); }
87
88 ctype volume() const { return data()->volumes( seed_ ); }
89
90 const EntitySeed& seed () const { return seed_; }
91 };
92
93 template <int mdim, int cordim>
95 {
96 typedef Storage Type;
97 };
98 };
99
100 typedef Dune::MultiLinearGeometry< ctype, mydimension, coorddimension, PolyhedralMultiLinearGeometryTraits<ctype> >
101 MultiLinearGeometryType;
102
103 typedef typename PolyhedralMultiLinearGeometryTraits< ctype > ::template
104 CornerStorage<mydimension, coorddimension >::Type CornerStorageType;
105
107 typedef FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed;
108
110 typedef FieldMatrix< ctype, mydim, cdim > JacobianTransposed;
111
113 using JacobianInverse = FieldMatrix< ctype, mydim, cdim >;
114
116 using Jacobian = FieldMatrix< ctype, cdim, mydim >;
117
118 typedef Dune::Impl::FieldMatrixHelper< ctype > MatrixHelperType;
119
120 explicit PolyhedralGridBasicGeometry ( ExtraData data )
121 : storage_( data )
122 {}
123
124 PolyhedralGridBasicGeometry ( ExtraData data, const EntitySeed& seed )
125 : storage_( data, seed )
126 {
127 GeometryType myType = type();
128 if( ! myType.isNone() && storage_.isValid() )
129 {
130 geometryImpl_.reset( new MultiLinearGeometryType(myType, storage_) );
131 }
132 //std::cout << myType << " " << storage_.corners() << std::endl;
133 }
134
135 GeometryType type () const { return data()->geometryType( storage_.seed() ); }
136 bool affine () const { return (geometryImpl_) ? geometryImpl_->affine() : false; }
137
138 int corners () const { return storage_.corners(); }
139 GlobalCoordinate corner ( const int i ) const { return storage_.corner( i ); }
140 GlobalCoordinate center () const
141 {
142 if( type().isNone() )
143 {
144 return storage_.center();
145 }
146 else
147 {
148 const auto refElem = Dune::referenceElement< ctype, mydim > ( type() );
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
219 jacobian(const LocalCoordinate& local ) const
220 {
221 return jacobianTransposed(local).transposed();
222 }
223
226 jacobianInverse(const LocalCoordinate& local) const
227 {
228 return jacobianInverseTransposed(local).transposed();
229 }
230
231 ExtraData data() const { return storage_.data(); }
232
233 protected:
234 CornerStorageType storage_;
235 std::shared_ptr< MultiLinearGeometryType > geometryImpl_;
236 };
237
238
239 // PolyhedralGridGeometry
240 // --------------
241
242 template< int mydim, int cdim, class Grid >
244 : public PolyhedralGridBasicGeometry< mydim, cdim, Grid >
245 {
247
248 public:
249 typedef typename Base::ExtraData ExtraData;
250 typedef typename Base::EntitySeed EntitySeed;
251
252 explicit PolyhedralGridGeometry ( ExtraData data )
253 : Base( data )
254 {}
255
256 PolyhedralGridGeometry ( ExtraData data, const EntitySeed& seed )
257 : Base( data, seed )
258 {}
259 };
260
261 template< int mydim, int cdim, class Grid >
263 : public PolyhedralGridBasicGeometry< mydim, cdim, Grid >
264 {
266
267 public:
268 typedef typename Base::ExtraData ExtraData;
269
270 explicit PolyhedralGridLocalGeometry ( ExtraData data )
271 : Base( data )
272 {}
273 };
274
275
276} // namespace Dune
277
278#endif // #ifndef DUNE_POLYHEDRALGRID_GEOMETRY_HH
Definition: geometry.hh:245
Definition: geometry.hh:264
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
Definition: geometry.hh:30
FieldMatrix< ctype, mydim, cdim > JacobianInverse
type of jacobian inverse transposed
Definition: geometry.hh:113
Jacobian jacobian(const LocalCoordinate &local) const
The jacobian.
Definition: geometry.hh:219
FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed
type of jacobian inverse transposed
Definition: geometry.hh:107
FieldMatrix< ctype, cdim, mydim > Jacobian
type of jacobian transposed
Definition: geometry.hh:116
FieldMatrix< ctype, mydim, cdim > JacobianTransposed
type of jacobian transposed
Definition: geometry.hh:110
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
The inverse of the jacobian.
Definition: geometry.hh:226
LocalCoordinate local(const GlobalCoordinate &global) const
Mapping from the cell to the reference domain.
Definition: geometry.hh:165