dune-alugrid 3.0.0
geometry_imp.cc
Go to the documentation of this file.
1#ifndef DUNE_ALUGRID_GEOMETRY_IMP_CC
2#define DUNE_ALUGRID_GEOMETRY_IMP_CC
3
4#include <dune/geometry/genericgeometry/topologytypes.hh>
5
6#include "grid.hh"
7#include "mappings.hh"
8#include "geometry.hh"
10
11
12namespace Dune {
13// --Geometry
14
15template< int mydim, int cdim, class GridImp>
16inline GeometryType
18{
19 return GeometryType( (elementType == tetra) ?
20 GenericGeometry :: SimplexTopology< mydim > :: type :: id :
21 GenericGeometry :: CubeTopology < mydim > :: type :: id,
22 mydim );
23}
24
25template< int mydim, int cdim, class GridImp>
26inline int
28{
29 return corners_;
30}
31
32template< int mydim, int cdim, class GridImp>
35corner (int i) const
36{
37 return geoImpl()[ i ];
38}
39
40
41template< int mydim, int cdim, class GridImp>
44global (const LocalCoordinate& local) const
45{
46 GlobalCoordinate global;
47 geoImpl().mapping().map2world(local, global);
48 return global;
49}
50
51template< int mydim, int cdim, class GridImp >
54local (const GlobalCoordinate& global) const
55{
56 LocalCoordinate local;
57 geoImpl().mapping().world2map(global, local);
58 return local;
59}
60
61template< int mydim, int cdim, class GridImp>
64integrationElement (const LocalCoordinate& local) const
65{
66 // this is the only case we need to specialize
67 if( mydim == 3 && elementType == tetra )
68 {
69 alugrid_assert ( geoImpl().valid() );
70 return 6.0 * geoImpl().volume();
71 }
72 else
73 return geoImpl().mapping().det( local );
74}
75
76template<int mydim, int cdim, class GridImp>
79volume () const
80{
81 if( mydim == 3 )
82 {
83 alugrid_assert ( geoImpl().valid() );
84 return geoImpl().volume() ;
85 }
86 else if ( mydim == 2 && elementType == tetra )
87 {
88 enum { factor = Factorial<mydim>::factorial };
89 // local vector does not affect the result
90 const LocalCoordinate dummy(0);
91 return integrationElement( dummy ) / ((ctype) factor);
92 }
93 else
94 {
95 return integrationElement(LocalCoordinate(0.5));
96 }
97}
98
99template< int mydim, int cdim, class GridImp>
100inline bool
102affine() const
103{
104 return geoImpl().mapping().affine();
105}
106
107template< int mydim, int cdim, class GridImp>
111{
112 return geoImpl().mapping().jacobianInverseTransposed( local );
113}
114
115template< int mydim, int cdim, class GridImp>
118jacobianTransposed (const LocalCoordinate& local) const
119{
120 return geoImpl().mapping().jacobianTransposed( local );
121}
122
123template <int mydim, int cdim, class GridImp>
124inline void
126print (std::ostream& ss) const
127{
128 const char* charElType = (elementType == tetra) ? "tetra" : "hexa";
129 ss << "ALU3dGridGeometry<" << mydim << ", " << cdim << ", " << charElType << "> = {\n";
130 for(int i=0; i<corners(); ++i)
131 {
132 ss << " corner " << i << " ";
133 ss << "{" << corner(i) << "}"; ss << std::endl;
134 }
135 ss << "} \n";
136}
137
138// built Geometry
139template <int mydim, int cdim, class GridImp>
140template <class GeometryType>
141inline bool
143buildGeomInFather(const GeometryType &fatherGeom , const GeometryType & myGeom)
144{
145 // update geo impl
146 geoImpl().updateInFather( fatherGeom, myGeom );
147
148 // my volume is a part of 1 for hexas, for tetra adjust with factor
149 double volume = myGeom.volume() / fatherGeom.volume() ;
150 if( elementType == tetra && mydim == 3 )
151 {
152 volume /= 6.0; // ???
153 geoImpl().setVolume( volume );
154#ifdef ALUGRIDDEBUG
155 LocalCoordinate local( 0.0 );
156 alugrid_assert ( std::abs( 6.0 * geoImpl().volume() - integrationElement( local ) ) < 1e-12 );
157#endif
158 }
159 else
160 geoImpl().setVolume( volume );
161
162 return true;
163}
164
165//--hexaBuildGeom
166template <int mydim, int cdim, class GridImp>
167inline bool
169buildGeom(const IMPLElementType& item)
170{
171 alugrid_assert( int(mydim) == 2 || int(mydim) == 3 );
172
173 if ( elementType == hexa )
174 {
175 // if this assertion is thrown, use ElementTopo::dune2aluVertex instead
176 // of number when calling myvertex
177 alugrid_assert ( ElementTopo::dune2aluVertex(0) == 0 );
178 alugrid_assert ( ElementTopo::dune2aluVertex(1) == 1 );
179 alugrid_assert ( ElementTopo::dune2aluVertex(2) == 3 );
180 alugrid_assert ( ElementTopo::dune2aluVertex(3) == 2 );
181 alugrid_assert ( ElementTopo::dune2aluVertex(4) == 4 );
182 alugrid_assert ( ElementTopo::dune2aluVertex(5) == 5 );
183 alugrid_assert ( ElementTopo::dune2aluVertex(6) == 7 );
184 alugrid_assert ( ElementTopo::dune2aluVertex(7) == 6 );
185
186 if( mydim == 3 ) // hexahedron
187 {
188 // update geo impl
189 geoImpl().update( item.myvertex(0)->Point(),
190 item.myvertex(1)->Point(),
191 item.myvertex(3)->Point(),
192 item.myvertex(2)->Point(),
193 item.myvertex(4)->Point(),
194 item.myvertex(5)->Point(),
195 item.myvertex(7)->Point(),
196 item.myvertex(6)->Point() );
197 }
198 else if ( mydim == 2 ) // quadrilateral
199 {
200 // update geo impl (drop vertex 4,5,6,7)
201 geoImpl().update( item.myvertex(0)->Point(),
202 item.myvertex(1)->Point(),
203 item.myvertex(3)->Point(),
204 item.myvertex(2)->Point() );
205 }
206 }
207 else if( elementType == tetra )
208 {
209 // if this assertion is thrown, use ElementTopo::dune2aluVertex instead
210 // of number when calling myvertex
211 alugrid_assert ( ElementTopo::dune2aluVertex(0) == 0 );
212 alugrid_assert ( ElementTopo::dune2aluVertex(1) == 1 );
213 alugrid_assert ( ElementTopo::dune2aluVertex(2) == 2 );
214 alugrid_assert ( ElementTopo::dune2aluVertex(3) == 3 );
215
216 if( mydim == 3 ) // tetrahedron
217 {
218 // update geo impl
219 geoImpl().update( item.myvertex(0)->Point(),
220 item.myvertex(1)->Point(),
221 item.myvertex(2)->Point(),
222 item.myvertex(3)->Point() );
223 }
224 else if( mydim == 2 ) // triangle
225 {
226 // update geo impl (drop vertex 0)
227 geoImpl().update( item.myvertex(1)->Point(),
228 item.myvertex(2)->Point(),
229 item.myvertex(3)->Point() );
230 }
231 }
232
233 if( mydim == 3 )
234 {
235 // get volume of element
236 geoImpl().setVolume( item.volume() );
237 }
238
239 return true;
240}
241
242// buildFaceGeom
243template <int mydim, int cdim, class GridImp>
244inline bool
246buildGeom(const HFaceType & item, int t)
247{
248 // get geo face
249 const GEOFaceType& face = static_cast<const GEOFaceType&> (item);
250
251 const int numVertices = ElementTopo::numVerticesPerFace;
252 typedef ALUTwist< numVertices, 2 > Twist;
253 const Twist twist( t );
254
255 // for all vertices of this face get rotatedIndex
256 int rotatedALUIndex[ 4 ];
257 for( int i = 0; i < numVertices; ++i )
258 rotatedALUIndex[ i ] = (elementType == tetra ? twist( i ) : twist( i ) ^ (twist( i ) >> 1));
259
260 if( elementType == hexa )
261 {
262 if( mydim == 2 ) //quadrilateral
263 {
264 // update geometry implementation
265 geoImpl().update( face.myvertex(rotatedALUIndex[0])->Point(),
266 face.myvertex(rotatedALUIndex[1])->Point(),
267 face.myvertex(rotatedALUIndex[2])->Point(),
268 face.myvertex(rotatedALUIndex[3])->Point() );
269 }
270 else if ( mydim == 1) //edge
271 {
272 //update geometry implementation
273 //we cannot use the rotatedALUIndex here, because for the codimiterator we get the wrong twist
274 geoImpl().update( face.myvertex(t < 0 ? 0 : 3)->Point(),
275 face.myvertex(t < 0 ? 3 : 0)->Point() );
276 }
277 }
278 else if ( elementType == tetra )
279 {
280 if ( mydim == 2) //triangle
281 {
282 // update geometry implementation
283 geoImpl().update( face.myvertex(rotatedALUIndex[0])->Point(),
284 face.myvertex(rotatedALUIndex[1])->Point(),
285 face.myvertex(rotatedALUIndex[2])->Point());
286 }
287 else if ( mydim == 1 ) //edge
288 {
289 //update geometry implementation (drop index 0)
290 geoImpl().update( face.myvertex( rotatedALUIndex[1])->Point(),
291 face.myvertex( rotatedALUIndex[2])->Point() );
292 }
293 }
294
295 return true;
296}
297
298// --buildFaceGeom
299template <int mydim, int cdim, class GridImp>
300template <class coord_t>
301inline bool
303buildGeom(const coord_t& p0,
304 const coord_t& p1,
305 const coord_t& p2,
306 const coord_t& p3)
307{
308 // update geometry implementation
309 geoImpl().update( p0, p1, p2, p3 );
310 return true;
311}
312
313// --buildFaceGeom
314template <int mydim, int cdim, class GridImp>
315template <class coord_t>
316inline bool
318buildGeom(const coord_t& p0,
319 const coord_t& p1,
320 const coord_t& p2)
321{
322 // update geometry implementation
323 geoImpl().update( p0, p1, p2 );
324 return true;
325}
326
327
328// --buildFaceGeom for edges
329template <int mydim, int cdim, class GridImp>
330template <class coord_t>
331inline bool
333buildGeom(const coord_t& p0,
334 const coord_t& p1)
335{
336 alugrid_assert (mydim == 1 );
337 // update geometry implementation
338 geoImpl().update( p0, p1 );
339 return true;
340}
341
342
343template <int mydim, int cdim, class GridImp> // for faces
344inline bool
346buildGeom(const FaceCoordinatesType& coords)
347{
348 if ( elementType == hexa )
349 {
350 if ( mydim == 2)
351 return buildGeom( coords[0], coords[1], coords[2], coords[3] );
352 else if ( mydim == 1 )
353 return buildGeom( coords[0], coords[1] );
354 }
355 else
356 {
357 alugrid_assert ( elementType == tetra );
358 if (mydim == 2 )
359 return buildGeom( coords[0], coords[1], coords[2] );
360 else if ( mydim == 1 )
361 return buildGeom( coords[0], coords[1] );
362 }
363}
364
365template <int mydim, int cdim, class GridImp> // for edges
366inline bool
368buildGeom(const HEdgeType & item, int twist)
369{
370 const GEOEdgeType & edge = static_cast<const GEOEdgeType &> (item);
371
372 if (mydim == 1) // edge
373 {
374 // update geometry implementation
375 geoImpl().update( edge.myvertex((twist) %2)->Point(),
376 edge.myvertex((1+twist)%2)->Point() );
377 }
378 else if ( mydim == 0) // point
379 {
380 if (elementType == hexa)
381 {
382 // update geometry implementation (drop vertex 1 as it has higher global index)
383 geoImpl().update( edge.myvertex(0)->Point() );
384 }
385 else if ( elementType == tetra)
386 {
387 // update geometry implementation (drop vertex 0)
388 geoImpl().update( edge.myvertex(1)->Point() );
389 }
390 }
391 return true;
392}
393
394template <int mydim, int cdim, class GridImp> // for Vertices ,i.e. Points
395inline bool
397buildGeom(const VertexType & item, int twist)
398{
399 // update geometry implementation
400 geoImpl().update( static_cast<const GEOVertexType &> (item).Point() );
401 return true;
402}
403
404} // end namespace Dune
405#endif // end DUNE_ALUGRID_GEOMETRY_IMP_CC
#define alugrid_assert(EX)
Definition alugrid_assert.hh:20
Definition alu3dinclude.hh:80
@ hexa
Definition topology.hh:12
@ tetra
Definition topology.hh:12
Definition geometry.hh:632
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &local) const
jacobian transposed
Definition geometry_imp.cc:118
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &local) const
Definition geometry_imp.cc:110
ALU3dImplTraitsType::GEOFaceType GEOFaceType
Definition geometry.hh:642
FieldMatrix< ctype, GridImp::dimension==3 ? EntityCount< elementType >::numVerticesPerFace :2, cdim > FaceCoordinatesType
Definition geometry.hh:678
FieldMatrix< ctype, cdim, mydim > JacobianInverseTransposed
type of jacobian inverse transposed
Definition geometry.hh:671
LocalCoordinate local(const GlobalCoordinate &global) const
Definition geometry_imp.cc:54
GlobalCoordinate corner(int i) const
access to coordinates of corners. Index is the number of the corner
Definition geometry_imp.cc:35
int corners() const
return the number of corners of this element. Corners are numbered 0..n-1
Definition geometry_imp.cc:27
FieldVector< ctype, cdim > GlobalCoordinate
type of the global coordinates
Definition geometry.hh:668
ALU3dImplTraitsType::HEdgeType HEdgeType
Definition geometry.hh:648
ALU3dImplTraitsType::VertexType VertexType
Definition geometry.hh:649
ALU3dImplTraitsType::GEOEdgeType GEOEdgeType
Definition geometry.hh:643
bool buildGeom(const IMPLElementType &item)
Methods that not belong to the Interface, but have to be public.
Definition geometry_imp.cc:169
bool buildGeomInFather(const GeometryType &fatherGeom, const GeometryType &myGeom)
build geometry of local coordinates relative to father
Definition geometry_imp.cc:143
ctype integrationElement(const LocalCoordinate &local) const
A(l) , see grid.hh.
Definition geometry_imp.cc:64
GridImp::ctype ctype
Definition geometry.hh:662
ALU3dImplTraitsType::IMPLElementType IMPLElementType
Definition geometry.hh:641
void print(std::ostream &ss) const
Definition geometry_imp.cc:126
bool affine() const
returns true if mapping is affine
Definition geometry_imp.cc:102
ctype volume() const
returns volume of geometry
Definition geometry_imp.cc:79
GlobalCoordinate global(const LocalCoordinate &local) const
Definition geometry_imp.cc:44
FieldVector< ctype, mydim > LocalCoordinate
type of local coordinates
Definition geometry.hh:665
ALU3dImplTraitsType::GEOVertexType GEOVertexType
Definition geometry.hh:644
ALU3dImplTraitsType::HFaceType HFaceType
Definition geometry.hh:647
FieldMatrix< ctype, mydim, cdim > JacobianTransposed
type of jacobian transposed
Definition geometry.hh:674
Definition twists.hh:19