dune-fem  2.8-git
piolatransformation.hh
Go to the documentation of this file.
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 
12 
13 namespace 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 >
48  {
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:
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  {
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  {
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 >
150  {
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:
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  {
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  {
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
Definition: bindguard.hh:11
double determinante(const Mat &m)
Definition: piolatransformation.hh:24
Definition: fmatrixcol.hh:16
Definition: piolatransformation.hh:150
InversePiolaTransformation(const Geometry &geo, const Point &p)
Definition: piolatransformation.hh:165
FieldVector< F, dimRange > apply_t(const FieldVector< F, dimRange > &d) const
Definition: piolatransformation.hh:186
real_type detInv_
Definition: piolatransformation.hh:242
PiolaTransformation< Geometry, dimRange > InverseTransformationType
Definition: piolatransformation.hh:162
FieldMatrix< F, dimRange, dimDomain > apply(const FieldMatrix< F, dimRange, dimDomain > &d) const
Definition: piolatransformation.hh:201
FieldVector< F, dimRange > apply(const FieldVector< F, dimRange > &d) const
Definition: piolatransformation.hh:171
FieldMatrix< F, dimRange, dimDomain > apply_t(const FieldMatrix< F, dimRange, dimDomain > &d) const
Definition: piolatransformation.hh:222
JacobianInverseTransposed gjit_
Definition: piolatransformation.hh:241
Definition: piolatransformation.hh:48
FieldMatrix< F, dimRange, dimDomain > apply_t(const FieldMatrix< F, dimRange, dimDomain > &d) const
Definition: piolatransformation.hh:120
FieldVector< F, dimRange > apply_t(const FieldVector< F, dimRange > &d) const
Definition: piolatransformation.hh:84
PiolaTransformation(const Geometry &geo, const Point &p)
Definition: piolatransformation.hh:63
JacobianTransposed gjt_
Definition: piolatransformation.hh:139
InversePiolaTransformation< Geometry, dimRange > InverseTransformationType
Definition: piolatransformation.hh:60
FieldMatrix< F, dimRange, dimDomain > apply(const FieldMatrix< F, dimRange, dimDomain > &d) const
Definition: piolatransformation.hh:99
FieldVector< F, dimRange > apply(const FieldVector< F, dimRange > &d) const
Definition: piolatransformation.hh:69
real_type detInv_
Definition: piolatransformation.hh:140