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 template< class F >
69 FieldVector< F, dimRange > apply ( const FieldVector< F, dimRange > &d ) const
70 {
71 FieldVector< F, dimRange > ret( d );
72 FieldVector< F, dimDomain > arg, dest;
73 for( std::size_t r = 0; r < blocks; ++r )
74 {
75 std::copy_n( d.begin() + r*dimDomain, dimDomain, arg.begin() );
76 gjt_.mtv( arg, dest );
77 std::copy_n( dest.begin(), dimDomain, ret.begin() + r*dimDomain );
78 }
79 ret *= detInv_;
80 return ret;
81 }
82
83 template< class F >
84 FieldVector< F, dimRange > apply_t ( const FieldVector< F, dimRange > &d ) const
85 {
86 FieldVector< F, dimRange > ret( d );
87 FieldVector< F, dimDomain > arg, dest;
88 for( std::size_t r = 0; r < blocks; ++r )
89 {
90 std::copy_n( ret.begin() + r*dimDomain, dimDomain, arg.begin() );
91 gjt_.mv( arg, dest );
92 std::copy_n( dest.begin(), dimDomain, ret.begin() + r*dimDomain );
93 }
94 ret *= detInv_;
95 return ret;
96 }
97
98 template< class F >
99 FieldMatrix< F, dimRange, dimDomain > apply ( const FieldMatrix< F, dimRange, dimDomain > &d ) const
100 {
101 FieldMatrix< F, dimRange, dimDomain > ret( d );
102 FieldVector< F, dimDomain > arg, dest;
103
104 for( std::size_t r = 0; r < dimDomain; ++r )
105 {
106 FieldMatrixColumn< FieldMatrix< F, dimRange, dimDomain > > col( ret, r );
107 for( std::size_t b = 0; b < blocks; ++b )
108 {
109 std::copy_n( col.begin() + b*dimDomain, dimDomain, arg.begin() );
110 gjt_.mv( arg, dest );
111 std::copy_n( dest.begin(), dimDomain, col.begin() + b*dimDomain );
112 }
113 }
114
115 ret *= detInv_;
116 return ret;
117 }
118
119 template< class F >
120 FieldMatrix< F, dimRange, dimDomain > apply_t ( const FieldMatrix< F, dimRange, dimDomain > &d ) const
121 {
122 FieldMatrix< F, dimRange, dimDomain > ret( d );
123 FieldVector< F, dimDomain > arg, dest;
124 for( std::size_t r = 0; r < dimDomain; ++r )
125 {
126 FieldMatrixColumn< FieldMatrix< F, dimRange, dimDomain > > col( ret, r );
127 for( std::size_t b = 0; b < blocks; ++b )
128 {
129 std::copy_n( col.begin() + b*dimDomain, dimDomain, arg.begin() );
130 gjt_.mtv( arg, dest );
131 std::copy_n( dest.begin(), dimDomain, col.begin() + b*dimDomain );
132 }
133 }
134 ret *= detInv_;
135 return ret;
136 }
137
138 protected:
139 JacobianTransposed gjt_;
140 real_type detInv_;
141 };
142
143
144
145 // InverseTransformationType
146 // -------------------------
147
148 template< class Geometry, int dimRange >
149 class InversePiolaTransformation
150 {
151 typedef InversePiolaTransformation< Geometry, dimRange > ThisType;
152
153 static const int dimDomain = Geometry::GlobalCoordinate::dimension;
154 typedef typename Geometry::JacobianInverseTransposed JacobianInverseTransposed;
155
156 static_assert( dimRange % dimDomain == 0, "PiolaTransformation can only be applied if dimRange is a multiple of dimDomain." );
157
158 typedef typename FieldTraits< JacobianInverseTransposed >::real_type real_type;
159 static const int blocks = dimRange / dimDomain;
160
161 public:
162 typedef PiolaTransformation< Geometry, dimRange > InverseTransformationType;
163
164 template< class Point >
165 InversePiolaTransformation ( const Geometry &geo, const Point &p )
166 : gjit_( geo.jacobianInverseTransposed( p ) ),
167 detInv_( 1.0 / determinante( gjit_ ) )
168 {}
169
170 template< class F >
171 FieldVector< F, dimRange > apply ( const FieldVector< F, dimRange > &d ) const
172 {
173 FieldVector< F, dimRange > ret( d );
174 FieldVector< F, dimDomain > arg, dest;
175 for( std::size_t r = 0; r < blocks; ++r )
176 {
177 std::copy_n( d.begin() + r*dimDomain, dimDomain, arg.begin() );
178 gjit_.mtv( arg, dest );
179 std::copy_n( dest.begin(), dimDomain, ret.begin() + r*dimDomain );
180 }
181 ret *= detInv_;
182 return ret;
183 }
184
185 template< class F >
186 FieldVector< F, dimRange > apply_t ( 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( ret.begin() + r*dimDomain, dimDomain, arg.begin() );
193 gjit_.mv( 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 FieldMatrix< F, dimRange, dimDomain > apply ( const FieldMatrix< F, dimRange, dimDomain > &d ) const
202 {
203 FieldMatrix< F, dimRange, dimDomain > ret( d );
204 FieldVector< F, dimDomain > arg, dest;
205
206 for( std::size_t r = 0; r < dimDomain; ++r )
207 {
208 FieldMatrixColumn< FieldMatrix< F, dimRange, dimDomain > > col( ret, r );
209 for( std::size_t b = 0; b < blocks; ++b )
210 {
211 std::copy_n( col.begin() + b*dimDomain, dimDomain, arg.begin() );
212 gjit_.mv( arg, dest );
213 std::copy_n( dest.begin(), dimDomain, col.begin() + b*dimDomain );
214 }
215 }
216
217 ret *= detInv_;
218 return ret;
219 }
220
221 template< class F >
222 FieldMatrix< F, dimRange, dimDomain > apply_t ( const FieldMatrix< F, dimRange, dimDomain > &d ) const
223 {
224 FieldMatrix< F, dimRange, dimDomain > ret( d );
225 FieldVector< F, dimDomain > arg, dest;
226 for( std::size_t r = 0; r < dimDomain; ++r )
227 {
228 FieldMatrixColumn< FieldMatrix< F, dimRange, dimDomain > > col( ret, r );
229 for( std::size_t b = 0; b < blocks; ++b )
230 {
231 std::copy_n( col.begin() + b*dimDomain, dimDomain, arg.begin() );
232 gjit_.mtv( arg, dest );
233 std::copy_n( dest.begin(), dimDomain, col.begin() + b*dimDomain );
234 }
235 }
236 ret *= detInv_;
237 return ret;
238 }
239
240 protected:
241 JacobianInverseTransposed gjit_;
242 real_type detInv_;
243 };
244
245 } // namespace Fem
246
247} // namespace Dune
248
249#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:96
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  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)