DUNE PDELab (git)

geometry.hh
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ALBERTA_GEOMETRY_HH
6#define DUNE_ALBERTA_GEOMETRY_HH
7
9#include <dune/grid/albertagrid/misc.hh>
11
12#if HAVE_ALBERTA
13
14namespace Dune
15{
16
17 // Forward Declarations
18 // --------------------
19
20 template< int dim, int dimworld >
21 class AlbertaGrid;
22
23
24
25 // AlbertaGridCoordinateReader
26 // ---------------------------
27
28 template< int codim, class GridImp >
29 struct AlbertaGridCoordinateReader
30 {
31 typedef typename std::remove_const< GridImp >::type Grid;
32
33 static constexpr int dimension = Grid::dimension;
34 static constexpr int codimension = codim;
35 static constexpr int mydimension = dimension - codimension;
36 static constexpr int coorddimension = Grid::dimensionworld;
37
38 typedef Alberta::Real ctype;
39
40 typedef Alberta::ElementInfo< dimension > ElementInfo;
41 typedef FieldVector< ctype, coorddimension > Coordinate;
42
43 AlbertaGridCoordinateReader ( const GridImp &grid,
44 const ElementInfo &elementInfo,
45 int subEntity )
46 : grid_( grid ),
47 elementInfo_( elementInfo ),
48 subEntity_( subEntity )
49 {}
50
51 const ElementInfo &elementInfo () const
52 {
53 return elementInfo_;
54 }
55
56 void coordinate ( int i, Coordinate &x ) const
57 {
58 assert( !elementInfo_ == false );
59 assert( (i >= 0) && (i <= mydimension) );
60
61 const int k = mapVertices( subEntity_, i );
62 const Alberta::GlobalVector &coord = grid_.getCoord( elementInfo_, k );
63 for( int j = 0; j < coorddimension; ++j )
64 x[ j ] = coord[ j ];
65 }
66
67 bool hasDeterminant () const
68 {
69 return false;
70 }
71
72 ctype determinant () const
73 {
74 assert( hasDeterminant() );
75 return ctype( 0 );
76 }
77
78 private:
79 static int mapVertices ( int subEntity, int i )
80 {
81 return Alberta::MapVertices< dimension, codimension >::apply( subEntity, i );
82 }
83
84 const Grid &grid_;
85 const ElementInfo &elementInfo_;
86 const int subEntity_;
87 };
88
89
90
91 // AlbertaGridGeometry
92 // -------------------
93
106 template< int mydim, int cdim, class GridImp >
108 {
110
111 // remember type of the grid
112 typedef GridImp Grid;
113
114 // dimension of barycentric coordinates
115 static constexpr int dimbary = mydim + 1;
116
117 public:
119 typedef Alberta::Real ctype;
120
121 static constexpr int dimension = Grid :: dimension;
122 static constexpr int mydimension = mydim;
123 static constexpr int codimension = dimension - mydimension;
124 static constexpr int coorddimension = cdim;
125
128
131
134
135 private:
136 static constexpr int numCorners = mydimension + 1;
137
139
140 public:
142 {
143 invalidate();
144 }
145
146 template< class CoordReader >
147 AlbertaGridGeometry ( const CoordReader &coordReader )
148 {
149 build( coordReader );
150 }
151
154 {
155 return GeometryTypes::simplex( mydimension );
156 }
157
159 bool affine () const { return true; }
160
162 int corners () const
163 {
164 return numCorners;
165 }
166
168 GlobalCoordinate corner ( const int i ) const
169 {
170 assert( (i >= 0) && (i < corners()) );
171 return coord_[ i ];
172 }
173
176 {
177 return centroid_;
178 }
179
181 GlobalCoordinate global ( const LocalCoordinate &local ) const;
182
185
192 {
193 assert( calcedDet_ );
194 return elDet_;
195 }
196
198 ctype integrationElement ( [[maybe_unused]] const LocalCoordinate &local ) const
199 {
200 return integrationElement();
201 }
202
204 ctype volume () const
205 {
206 return integrationElement() / ctype( factorial(mydimension) );
207 }
208
214 const JacobianTransposed &jacobianTransposed () const;
215
217 const JacobianTransposed &
218 jacobianTransposed ( [[maybe_unused]] const LocalCoordinate &local ) const
219 {
220 return jacobianTransposed();
221 }
222
228 const JacobianInverseTransposed &jacobianInverseTransposed () const;
229
231 const JacobianInverseTransposed &
232 jacobianInverseTransposed ( [[maybe_unused]] const LocalCoordinate &local ) const
233 {
235 }
236
238 Jacobian jacobian ( const LocalCoordinate &local ) const
239 {
240 return jacobianTransposed(local).transposed();
241 }
242
245 {
247 }
248
249 //***********************************************************************
250 // Methods that not belong to the Interface, but have to be public
251 //***********************************************************************
252
255 {
256 builtJT_ = false;
257 builtJTInv_ = false;
258 calcedDet_ = false;
259 }
260
262 template< class CoordReader >
263 void build ( const CoordReader &coordReader );
264
265 void print ( std::ostream &out ) const;
266
267 private:
268 // calculates the volume of the element
269 ctype elDeterminant () const
270 {
271 return std::abs( Alberta::determinant( jacobianTransposed() ) );
272 }
273
275 CoordMatrix coord_;
276
278 GlobalCoordinate centroid_;
279
280 // storage for the transposed of the jacobian
281 mutable JacobianTransposed jT_;
282
283 // storage for the transposed inverse of the jacboian
284 mutable JacobianInverseTransposed jTInv_;
285
286 // has jT_ been computed, yet?
287 mutable bool builtJT_;
288 // has jTInv_ been computed, yet?
289 mutable bool builtJTInv_;
290
291 mutable bool calcedDet_;
292 mutable ctype elDet_;
293 };
294
295
296
297 // AlbertaGridGlobalGeometry
298 // -------------------------
299
300 template< int mydim, int cdim, class GridImp >
301 class AlbertaGridGlobalGeometry
302 : public AlbertaGridGeometry< mydim, cdim, GridImp >
303 {
304 typedef AlbertaGridGlobalGeometry< mydim, cdim, GridImp > This;
305 typedef AlbertaGridGeometry< mydim, cdim, GridImp > Base;
306
307 public:
308 AlbertaGridGlobalGeometry ()
309 : Base()
310 {}
311
312 template< class CoordReader >
313 AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
314 : Base( coordReader )
315 {}
316 };
317
318
319#if !DUNE_ALBERTA_CACHE_COORDINATES
320 template< int dim, int cdim >
321 class AlbertaGridGlobalGeometry< dim, cdim, const AlbertaGrid< dim, cdim > >
322 {
323 typedef AlbertaGridGlobalGeometry< dim, cdim, const AlbertaGrid< dim, cdim > > This;
324
325 // remember type of the grid
326 typedef AlbertaGrid< dim, cdim > Grid;
327
328 // dimension of barycentric coordinates
329 static constexpr int dimbary = dim + 1;
330
331 typedef Alberta::ElementInfo< dim > ElementInfo;
332
333 public:
335 typedef Alberta::Real ctype;
336
337 static constexpr int dimension = Grid::dimension;
338 static constexpr int mydimension = dimension;
339 static constexpr int codimension = dimension - mydimension;
340 static constexpr int coorddimension = cdim;
341
342 typedef FieldVector< ctype, mydimension > LocalCoordinate;
343 typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
344
345 typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
346 typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
347
350
351 private:
352 static constexpr int numCorners = mydimension + 1;
353
355
356 public:
357 AlbertaGridGlobalGeometry ()
358 {
359 invalidate();
360 }
361
362 template< class CoordReader >
363 AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
364 {
365 build( coordReader );
366 }
367
369 GeometryType type () const
370 {
371 return GeometryTypes::simplex( mydimension );
372 }
373
375 int corners () const
376 {
377 return numCorners;
378 }
379
381 GlobalCoordinate corner ( const int i ) const
382 {
383 assert( (i >= 0) && (i < corners()) );
384 const Alberta::GlobalCoordinate &x = elementInfo_.coordinate( i );
385 GlobalCoordinate y;
386 for( int j = 0; j < coorddimension; ++j )
387 y[ j ] = x[ j ];
388 return y;
389 }
390
392 GlobalCoordinate center () const
393 {
394 GlobalCoordinate centroid_ = corner( 0 );
395 for( int i = 1; i < numCorners; ++i )
396 centroid_ += corner( i );
397 centroid_ *= ctype( 1 ) / ctype( numCorners );
398 return centroid_;
399 }
400
402 GlobalCoordinate global ( const LocalCoordinate &local ) const;
403
405 LocalCoordinate local ( const GlobalCoordinate &global ) const;
406
412 ctype integrationElement () const
413 {
414 return elementInfo_.geometryCache().integrationElement();
415 }
416
418 ctype integrationElement ( const LocalCoordinate &local ) const
419 {
420 return integrationElement();
421 }
422
424 ctype volume () const
425 {
426 return integrationElement() / ctype( factorial(mydimension) );
427 }
428
434 const JacobianTransposed &jacobianTransposed () const
435 {
436 return elementInfo_.geometryCache().jacobianTransposed();
437 }
438
440 const JacobianTransposed &
441 jacobianTransposed ( const LocalCoordinate &local ) const
442 {
443 return jacobianTransposed();
444 }
445
451 const JacobianInverseTransposed &jacobianInverseTransposed () const
452 {
453 return elementInfo_.geometryCache().jacobianInverseTransposed();
454 }
455
457 const JacobianInverseTransposed &
458 jacobianInverseTransposed ( const LocalCoordinate &local ) const
459 {
461 }
462
464 Jacobian jacobian ( const LocalCoordinate &local ) const
465 {
466 return jacobianTransposed(local).transposed();
467 }
468
470 JacobianInverse jacobianInverse ( const LocalCoordinate &local ) const
471 {
473 }
474
475 //***********************************************************************
476 // Methods that not belong to the Interface, but have to be public
477 //***********************************************************************
478
480 void invalidate ()
481 {
482 elementInfo_ = ElementInfo();
483 }
484
486 template< class CoordReader >
487 void build ( const CoordReader &coordReader )
488 {
489 elementInfo_ = coordReader.elementInfo();
490 }
491
492 private:
493 ElementInfo elementInfo_;
494 };
495#endif // #if !DUNE_ALBERTA_CACHE_COORDINATES
496
497
498
499 // AlbertaGridLocalGeometryProvider
500 // --------------------------------
501
502 template< class Grid >
503 class AlbertaGridLocalGeometryProvider
504 {
505 typedef AlbertaGridLocalGeometryProvider< Grid > This;
506
507 public:
508 typedef typename Grid::ctype ctype;
509
510 static constexpr int dimension = Grid::dimension;
511
512 template< int codim >
513 struct Codim
514 {
515 typedef AlbertaGridGeometry< dimension-codim, dimension, Grid > LocalGeometry;
516 };
517
518 typedef typename Codim< 0 >::LocalGeometry LocalElementGeometry;
519 typedef typename Codim< 1 >::LocalGeometry LocalFaceGeometry;
520
521 static constexpr int numChildren = 2;
522 static constexpr int numFaces = dimension + 1;
523
524 static constexpr int minFaceTwist = Alberta::Twist< dimension, dimension-1 >::minTwist;
525 static constexpr int maxFaceTwist = Alberta::Twist< dimension, dimension-1 >::maxTwist;
526 static constexpr int numFaceTwists = maxFaceTwist - minFaceTwist + 1;
527
528 private:
529 struct GeoInFatherCoordReader;
530 struct FaceCoordReader;
531
532 AlbertaGridLocalGeometryProvider ()
533 {
534 buildGeometryInFather();
535 buildFaceGeometry();
536 }
537
538 ~AlbertaGridLocalGeometryProvider ()
539 {
540 for( int child = 0; child < numChildren; ++child )
541 {
542 delete geometryInFather_[ child ][ 0 ];
543 delete geometryInFather_[ child ][ 1 ];
544 }
545
546 for( int i = 0; i < numFaces; ++i )
547 {
548 for( int j = 0; j < numFaceTwists; ++j )
549 delete faceGeometry_[ i ][ j ];
550 }
551 }
552
553 void buildGeometryInFather();
554 void buildFaceGeometry();
555
556 public:
557 const LocalElementGeometry &
558 geometryInFather ( int child, const int orientation = 1 ) const
559 {
560 assert( (child >= 0) && (child < numChildren) );
561 assert( (orientation == 1) || (orientation == -1) );
562 return *geometryInFather_[ child ][ (orientation + 1) / 2 ];
563 }
564
565 const LocalFaceGeometry &
566 faceGeometry ( int face, int twist = 0 ) const
567 {
568 assert( (face >= 0) && (face < numFaces) );
569 assert( (twist >= minFaceTwist) && (twist <= maxFaceTwist) );
570 return *faceGeometry_[ face ][ twist - minFaceTwist ];
571 }
572
573 static const This &instance ()
574 {
575 static This theInstance;
576 return theInstance;
577 }
578
579 private:
580 template< int codim >
581 static int mapVertices ( int subEntity, int i )
582 {
583 return Alberta::MapVertices< dimension, codim >::apply( subEntity, i );
584 }
585
586 const LocalElementGeometry *geometryInFather_[ numChildren ][ 2 ];
587 const LocalFaceGeometry *faceGeometry_[ numFaces ][ numFaceTwists ];
588 };
589
590} // namespace Dune
591
592#endif // #if HAVE_ALBERTA
593
594#endif // #ifndef DUNE_ALBERTA_GEOMETRY_HH
geometry implementation for AlbertaGrid
Definition: geometry.hh:108
GeometryType type() const
obtain the type of reference element
Definition: geometry.hh:153
GlobalCoordinate center() const
return center of geometry
Definition: geometry.hh:175
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &local) const
transposed of the geometry mapping's Jacobian
Definition: geometry.hh:218
const JacobianTransposed & jacobianTransposed() const
transposed of the geometry mapping's Jacobian
Definition: geometry.cc:55
int corners() const
number of corner the geometry
Definition: geometry.hh:162
const JacobianInverseTransposed & jacobianInverseTransposed() const
transposed inverse of the geometry mapping's Jacobian
Definition: geometry.cc:73
GlobalCoordinate corner(const int i) const
obtain the i-th corner of this geometry
Definition: geometry.hh:168
ctype integrationElement() const
integration element of the geometry mapping
Definition: geometry.hh:191
ctype integrationElement(const LocalCoordinate &local) const
integration element of the geometry mapping
Definition: geometry.hh:198
GlobalCoordinate global(const LocalCoordinate &local) const
map a point from the reference element to the geometry
Definition: geometry.cc:34
Alberta::Real ctype
type of coordinates
Definition: geometry.hh:119
ctype volume() const
volume of geometry
Definition: geometry.hh:204
Jacobian jacobian(const LocalCoordinate &local) const
geometry mapping's Jacobian
Definition: geometry.hh:238
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &local) const
transposed inverse of the geometry mapping's Jacobian
Definition: geometry.hh:232
void invalidate()
invalidate the geometry
Definition: geometry.hh:254
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
inverse of the geometry mapping's Jacobian
Definition: geometry.hh:244
bool affine() const
returns always true since we only have simplices
Definition: geometry.hh:159
void build(const CoordReader &coordReader)
build the geometry from a coordinate reader
Definition: geometry.cc:90
FieldMatrix< K, COLS, ROWS > transposed() const
Return transposed of the matrix as FieldMatrix.
Definition: fmatrix.hh:171
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
static constexpr int dimension
The dimension of the grid.
Definition: grid.hh:387
static constexpr int dimensionworld
The dimension of the world the grid lives in.
Definition: grid.hh:390
ct ctype
Define type used for coordinates in grid module.
Definition: grid.hh:518
Wrapper and interface classes for element geometries.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
provides a wrapper for ALBERTA's el_info structure
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:128
Dune namespace.
Definition: alignedallocator.hh:13
static constexpr T factorial(const T &n) noexcept
calculate the factorial of n as a constexpr
Definition: math.hh:111
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)