DUNE-FEM (unstable)

piolatransformation.hh
1#ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_PIOLATRANSFORMATION_HH
2#define DUNE_FEM_SPACE_BASISFUNCTIONSET_PIOLATRANSFORMATION_HH
3
4#include <algorithm>
5
6#include <dune/common/diagonalmatrix.hh>
7#include <dune/common/fmatrix.hh>
8#include <dune/common/fvector.hh>
9#include <dune/common/diagonalmatrix.hh>
10
11#include <dune/fem/common/fmatrixcol.hh>
12
13namespace Dune
14{
15
16 namespace Fem
17 {
18
19
20 // usefull implementations
21 // -----------------------
22
23 template< class Mat >
24 double determinante ( const Mat &m )
25 { return m.det(); }
26
27 template< class F, int d, int l >
28 double determinante ( const FieldMatrix< F, d, l > &m )
29 { return m.determinant(); }
30
31 template< class F, int d >
32 double determinante ( const DiagonalMatrix< F, d > &m )
33 { return m.determinant(); }
34
35
36 // internal forward declaration
37 // ----------------------------
38
39 template< class Geometry, int dimRange >
40 class InversePiolaTransformation;
41
42
43 // PiolaTransformation
44 // -------------------
45
46 template< class Geometry, int dimRange >
47 class PiolaTransformation
48 {
49 typedef PiolaTransformation< Geometry, dimRange > ThisType;
50
51 static const int dimDomain = Geometry::GlobalCoordinate::dimension;
52 typedef typename Geometry::JacobianTransposed JacobianTransposed;
53
54 static_assert( dimRange % dimDomain == 0, "PiolaTransformation can only be applied if dimRange is a multiple of dimDomain." );
55
56 typedef typename FieldTraits< JacobianTransposed >::real_type real_type;
57 static const int blocks = dimRange / dimDomain;
58
59 public:
60 typedef InversePiolaTransformation< Geometry, dimRange > InverseTransformationType;
61
62 template< class Point >
63 PiolaTransformation ( const Geometry &geo, const Point &p )
64 : gjt_( geo.jacobianTransposed( p ) ),
65 detInv_( 1.0 / determinante( gjt_ ) )
66 {
67 }
68
69 template< class F >
70 FieldVector< F, dimRange > apply ( const FieldVector< F, dimRange > &d ) const
71 {
72 FieldVector< F, dimRange > ret;
73 if constexpr ( dimDomain == dimRange )
74 {
75 gjt_.mtv( d, ret );
76 }
77 else
78 {
79 FieldVector< F, dimDomain > arg, dest;
80 for( int r = 0; r < blocks; ++r )
81 {
82 std::copy_n( d.begin() + r*dimDomain, dimDomain, arg.begin() );
83 gjt_.mtv( arg, dest );
84 std::copy_n( dest.begin(), dimDomain, ret.begin() + r*dimDomain );
85 }
86 }
87 ret *= detInv_;
88 return ret;
89 }
90
91 template< class F >
92 FieldVector< F, dimRange > apply_t ( const FieldVector< F, dimRange > &d ) const
93 {
94 FieldVector< F, dimRange > ret;
95 if constexpr ( dimDomain == dimRange )
96 {
97 gjt_.mv( d, ret );
98 }
99 else
100 {
101 FieldVector< F, dimDomain > arg, dest;
102 for( int r = 0; r < blocks; ++r )
103 {
104 std::copy_n( ret.begin() + r*dimDomain, dimDomain, arg.begin() );
105 gjt_.mv( arg, dest );
106 std::copy_n( dest.begin(), dimDomain, ret.begin() + r*dimDomain );
107 }
108 }
109 ret *= detInv_;
110 return ret;
111 }
112
113 template< class F >
114 FieldMatrix< F, dimRange, dimDomain > apply ( const FieldMatrix< F, dimRange, dimDomain > &d ) const
115 {
116 FieldMatrix< F, dimRange, dimDomain > ret( d );
117 FieldVector< F, dimDomain > arg, dest;
118
119 for( std::size_t r = 0; r < dimDomain; ++r )
120 {
121 FieldMatrixColumn< FieldMatrix< F, dimRange, dimDomain > > col( ret, r );
122 for( std::size_t b = 0; b < blocks; ++b )
123 {
124 std::copy_n( col.begin() + b*dimDomain, dimDomain, arg.begin() );
125 gjt_.mtv( arg, dest );
126 std::copy_n( dest.begin(), dimDomain, col.begin() + b*dimDomain );
127 }
128 }
129
130 ret *= detInv_;
131 return ret;
132 }
133
134 template< class F >
135 FieldMatrix< F, dimRange, dimDomain > apply_t ( const FieldMatrix< F, dimRange, dimDomain > &d ) const
136 {
137 FieldMatrix< F, dimRange, dimDomain > ret( d );
138 FieldVector< F, dimDomain > arg, dest;
139 for( std::size_t r = 0; r < dimDomain; ++r )
140 {
141 FieldMatrixColumn< FieldMatrix< F, dimRange, dimDomain > > col( ret, r );
142 for( std::size_t b = 0; b < blocks; ++b )
143 {
144 std::copy_n( col.begin() + b*dimDomain, dimDomain, arg.begin() );
145 gjt_.mv( arg, dest );
146 std::copy_n( dest.begin(), dimDomain, col.begin() + b*dimDomain );
147 }
148 }
149 ret *= detInv_;
150 return ret;
151 }
152
153 protected:
154 JacobianTransposed gjt_;
155 real_type detInv_;
156 };
157
158
159
160 // InverseTransformationType
161 // -------------------------
162
163 template< class Geometry, int dimRange >
164 class InversePiolaTransformation
165 {
166 typedef InversePiolaTransformation< Geometry, dimRange > ThisType;
167
168 static const int dimDomain = Geometry::GlobalCoordinate::dimension;
169 typedef typename Geometry::JacobianInverseTransposed JacobianInverseTransposed;
170
171 static_assert( dimRange % dimDomain == 0, "PiolaTransformation can only be applied if dimRange is a multiple of dimDomain." );
172
173 typedef typename FieldTraits< JacobianInverseTransposed >::real_type real_type;
174 static const int blocks = dimRange / dimDomain;
175
176 public:
177 typedef PiolaTransformation< Geometry, dimRange > InverseTransformationType;
178
179 template< class Point >
180 InversePiolaTransformation ( const Geometry &geo, const Point &p )
181 : gjit_( geo.jacobianInverseTransposed( p ) ),
182 detInv_( 1.0 / determinante( gjit_ ) )
183 {}
184
185 template< class F >
186 FieldVector< F, dimRange > apply ( const FieldVector< F, dimRange > &d ) const
187 {
188 FieldVector< F, dimRange > ret( d );
189 FieldVector< F, dimDomain > arg, dest;
190 for( std::size_t r = 0; r < blocks; ++r )
191 {
192 std::copy_n( d.begin() + r*dimDomain, dimDomain, arg.begin() );
193 gjit_.mtv( arg, dest );
194 std::copy_n( dest.begin(), dimDomain, ret.begin() + r*dimDomain );
195 }
196 ret *= detInv_;
197 return ret;
198 }
199
200 template< class F >
201 FieldVector< F, dimRange > apply_t ( const FieldVector< F, dimRange > &d ) const
202 {
203 FieldVector< F, dimRange > ret( d );
204 FieldVector< F, dimDomain > arg, dest;
205 for( std::size_t r = 0; r < blocks; ++r )
206 {
207 std::copy_n( ret.begin() + r*dimDomain, dimDomain, arg.begin() );
208 gjit_.mv( arg, dest );
209 std::copy_n( dest.begin(), dimDomain, ret.begin() + r*dimDomain );
210 }
211 ret *= detInv_;
212 return ret;
213 }
214
215 template< class F >
216 FieldMatrix< F, dimRange, dimDomain > apply ( const FieldMatrix< F, dimRange, dimDomain > &d ) const
217 {
218 FieldMatrix< F, dimRange, dimDomain > ret( d );
219 FieldVector< F, dimDomain > arg, dest;
220
221 for( std::size_t r = 0; r < dimDomain; ++r )
222 {
223 FieldMatrixColumn< FieldMatrix< F, dimRange, dimDomain > > col( ret, r );
224 for( std::size_t b = 0; b < blocks; ++b )
225 {
226 std::copy_n( col.begin() + b*dimDomain, dimDomain, arg.begin() );
227 gjit_.mtv( arg, dest );
228 std::copy_n( dest.begin(), dimDomain, col.begin() + b*dimDomain );
229 }
230 }
231
232 ret *= detInv_;
233 return ret;
234 }
235
236 template< class F >
237 FieldMatrix< F, dimRange, dimDomain > apply_t ( const FieldMatrix< F, dimRange, dimDomain > &d ) const
238 {
239 FieldMatrix< F, dimRange, dimDomain > ret( d );
240 FieldVector< F, dimDomain > arg, dest;
241 for( std::size_t r = 0; r < dimDomain; ++r )
242 {
243 FieldMatrixColumn< FieldMatrix< F, dimRange, dimDomain > > col( ret, r );
244 for( std::size_t b = 0; b < blocks; ++b )
245 {
246 std::copy_n( col.begin() + b*dimDomain, dimDomain, arg.begin() );
247 gjit_.mv( arg, dest );
248 std::copy_n( dest.begin(), dimDomain, col.begin() + b*dimDomain );
249 }
250 }
251 ret *= detInv_;
252 return ret;
253 }
254
255 protected:
256 JacobianInverseTransposed gjit_;
257 real_type detInv_;
258 };
259
260 } // namespace Fem
261
262} // namespace Dune
263
264#endif // #ifndef DUNE_FEM_SPACE_BASISFUNCTIONSET_PIOLATRANSFORMATION_HH
Dune::FieldVector< ctype, cdim >::dimension
static constexpr int dimension
The size of this vector.
Definition: fvector.hh:106
Dune::Geometry::JacobianTransposed
Implementation::JacobianTransposed JacobianTransposed
type of jacobian transposed
Definition: geometry.hh:131
Dune::Geometry::JacobianInverseTransposed
Implementation::JacobianInverseTransposed JacobianInverseTransposed
type of jacobian inverse transposed
Definition: geometry.hh:120
diagonalmatrix.hh
This file implements a quadratic diagonal matrix of fixed size.
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 27, 22:40, 2025)