Dune Core Modules (2.4.1)

geometrycache.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ALBERTA_GEOMETRYCACHE_HH
4#define DUNE_ALBERTA_GEOMETRYCACHE_HH
5
6#include <dune/grid/albertagrid/misc.hh>
7#include <dune/grid/albertagrid/algebra.hh>
8
9#if HAVE_ALBERTA
10
11namespace Dune
12{
13
14 namespace Alberta
15 {
16
17 // GeometryCache
18 // -------------
19
20 template< int dim >
21 class GeometryCache
22 {
23 static const unsigned int flagIntegrationElement = (1 << 0);
24 static const unsigned int flagJacobianTransposed = (1 << 1);
25 static const unsigned int flagJacobianInverseTransposed = (1 << 2);
26
27 public:
28 typedef FieldMatrix< Real, dimWorld, dim > JacobianInverseTransposed;
29 typedef FieldMatrix< Real, dim, dimWorld > JacobianTransposed;
30
31 GeometryCache ()
32 : flags_( 0 )
33 {}
34
35 const Real &integrationElement ( const ALBERTA EL_INFO &elInfo )
36 {
37 if( (flags_ & flagIntegrationElement) == 0 )
38 {
39 integrationElement_ = std::abs( determinant( jacobianTransposed( elInfo ) ) );
40 assert( integrationElement_ > 1e-14 );
41 flags_ |= flagIntegrationElement;
42 }
43 return integrationElement_;
44 }
45
46 const JacobianTransposed &jacobianTransposed ( const ALBERTA EL_INFO &elInfo )
47 {
48 if( (flags_ & flagJacobianTransposed) == 0 )
49 {
50 assert( (elInfo.fill_flag & FillFlags< dim >::coords) != 0 );
51 const GlobalVector &x = elInfo.coord[ 0 ];
52 for( int i = 0; i < dim; ++i )
53 {
54 const GlobalVector &y = elInfo.coord[ i+1 ];
55 for( int j = 0; j < dimWorld; ++j )
56 jacobianTransposed_[ i ][ j ] = y[ j ] - x[ j ];
57 }
58 flags_ |= flagJacobianTransposed;
59 }
60 return jacobianTransposed_;
61 }
62
63 const JacobianInverseTransposed &
64 jacobianInverseTransposed ( const ALBERTA EL_INFO &elInfo )
65 {
66 if( (flags_ & flagJacobianInverseTransposed) == 0 )
67 {
68 integrationElement_ = std::abs( invert( jacobianTransposed( elInfo ), jacobianInverseTransposed_ ) );
69 assert( integrationElement_ > 1e-14 );
70 flags_ |= flagIntegrationElement | flagJacobianInverseTransposed;
71 }
72 return jacobianInverseTransposed_;
73 }
74
75 private:
76 unsigned int flags_;
77 Real integrationElement_;
78 FieldMatrix< Real, dim, dimWorld > jacobianTransposed_;
79 FieldMatrix< Real, dimWorld, dim > jacobianInverseTransposed_;
80 };
81
82
83
84 // GeometryCacheProxy
85 // ------------------
86
87 template< int dim >
88 struct GeometryCacheProxy
89 {
90 typedef FieldMatrix< Real, dimWorld, dim > JacobianInverseTransposed;
91 typedef FieldMatrix< Real, dim, dimWorld > JacobianTransposed;
92
93 GeometryCacheProxy ( GeometryCache< dim > &geometryCache, const ALBERTA EL_INFO &elInfo )
94 : geometryCache_( geometryCache ),
95 elInfo_( elInfo )
96 {}
97
98 const Real &integrationElement ()
99 {
100 return geometryCache_.integrationElement( elInfo_ );
101 }
102
103 const JacobianTransposed &jacobianTransposed ()
104 {
105 return geometryCache_.jacobianTransposed( elInfo_ );
106 }
107
108 const JacobianInverseTransposed &jacobianInverseTransposed ()
109 {
110 return geometryCache_.jacobianInverseTransposed( elInfo_ );
111 }
112
113 private:
114 GeometryCache< dim > &geometryCache_;
115 const ALBERTA EL_INFO &elInfo_;
116 };
117
118 } // namespace Alberta
119
120} // namespace Dune
121
122#endif // #if HAVE_ALBERTA
123
124#endif // #ifndef DUNE_ALBERTA_GEOMETRYCACHE_HH
Dune namespace.
Definition: alignment.hh:10
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)