DUNE-FEM (unstable)

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