v0.8.23
Templates.hpp
Go to the documentation of this file.
1 /** \file Templates.hpp
2  * \brief Templates declarations
3  *
4  * MoFEM is free software: you can redistribute it and/or modify it under
5  * the terms of the GNU Lesser General Public License as published by the
6  * Free Software Foundation, either version 3 of the License, or (at your
7  * option) any later version.
8  *
9  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12  * License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
16  */
17 
18 #ifndef __TEMPLATES_HPP__
19 #define __TEMPLATES_HPP__
20 
21 namespace MoFEM {
22 
23 /**
24  * @brief Get Vector adaptor
25  *
26  * \code
27  *
28  * double *a;
29  * CHKERR VecGetArray(v,&a);
30  *
31  * for(int n = 0; n != nodes; ++n) {
32  *
33  * auto a = getVectorAdaptor(&a[3*n], 3);
34  * double dot = inner_prod(a, a);
35  *
36  * }
37  *
38  * CHKERR VecRetsoreArray(v,&a);
39  * \endcode
40  *
41  */
42 template <typename T1>
43 inline auto getVectorAdaptor(T1 ptr, const size_t n) {
44  typedef typename std::remove_pointer<T1>::type T;
46  ublas::shallow_array_adaptor<T>(n, ptr));
47 };
48 
49 /**
50  * @brief Get Matrix adaptor
51  *
52  * \code
53  *
54  * double *a;
55  * CHKERR VecGetArray(v,&a);
56  *
57  * for(int n = 0; n != nodes; ++n) {
58  *
59  * auto F = getMatrixAdaptor(&a[3*3*n], 3, 3);
60  * MatrixDouble C = prod(F, trans(F));
61  *
62  * }
63  *
64  * CHKERR VecRetsoreArray(v,&a);
65  * \endcode
66  *
67  */
68 template <typename T1>
69 inline auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m) {
70  typedef typename std::remove_pointer<T1>::type T;
72  n, m, ublas::shallow_array_adaptor<T>(n * m, ptr));
73 };
74 
75 /**
76  * This small utility that cascades two key extractors will be
77  * used throughout the boost example
78  * <a
79  * href=http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp>
80  * http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp
81  * </a>
82  */
83 template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
84 public:
85  typedef typename KeyExtractor1::result_type result_type;
86 
87  KeyFromKey(const KeyExtractor1 &key1_ = KeyExtractor1(),
88  const KeyExtractor2 &key2_ = KeyExtractor2())
89  : key1(key1_), key2(key2_) {}
90 
91  template <typename Arg> result_type operator()(Arg &arg) const {
92  return key1(key2(arg));
93  }
94 
95 private:
96  KeyExtractor1 key1;
97  KeyExtractor2 key2;
98 };
99 
100 template <typename id_type> struct LtBit {
101  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
102  return valueA.to_ulong() < valueB.to_ulong();
103  }
104 };
105 
106 template <typename id_type> struct EqBit {
107  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
108  return valueA.to_ulong() == valueB.to_ulong();
109  }
110 };
111 
112 template <typename id_type> struct HashBit {
113  inline unsigned int operator()(const id_type &value) const {
114  return value.to_ulong();
115  }
116 };
117 
118 template <class X> inline std::string toString(X x) {
119  std::ostringstream buffer;
120  buffer << x;
121  return buffer.str();
122 }
123 
124 /**
125 * \brief Get tensor rank 0 (scalar) form data vector
126 * \ingroup mofem_forces_and_sources_user_data_operators
127 
128 Example how to use it.
129 \code
130 VectorDouble vec;
131 vec.resize(nb_gauss_pts,false);
132 vec.clear();
133 auto t0 = getFTensor0FromData(data);
134 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
135 
136  ++t0;
137 }
138 \endcode
139 
140 */
141 template <class T, class A>
143 getFTensor0FromVec(ublas::vector<T, A> &data) {
144  static_assert(!std::is_same<T, T>::value, "not implemented");
145 }
146 
147 template <>
150  ublas::vector<double, DoubleAllocator> &data) {
151  return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&*data.data().begin());
152 }
153 
154 /**
155  * \brief Get tensor rank 1 (vector) form data matrix
156  * \ingroup mofem_forces_and_sources_user_data_operators
157  */
158 template <int Tensor_Dim, class T, class L, class A>
159 static inline FTensor::Tensor1<FTensor::PackPtr<T *, 1>, Tensor_Dim>
160 getFTensor1FromMat(ublas::matrix<T, L, A> &data) {
161  static_assert(!std::is_same<T, T>::value, "not implemented");
162 }
163 
164 /**
165  * \brief Get tensor rank 1 (vector) form data matrix (specialization)
166  * \ingroup mofem_forces_and_sources_user_data_operators
167  */
168 template <int Tensor_Dim>
169 static inline FTensor::Tensor1<FTensor::PackPtr<double *, 1>, Tensor_Dim>
171  return getFTensor1FromMat<Tensor_Dim, double, ublas::row_major,
172  DoubleAllocator>(data);
173 }
174 
175 template <>
177 getFTensor1FromMat<3, double, ublas::row_major, DoubleAllocator>(
178  MatrixDouble &data) {
179  if (data.size1() != 3) {
180  THROW_MESSAGE("Wrong size of data matrix");
181  }
183  &data(0, 0), &data(1, 0), &data(2, 0));
184 }
185 
186 template <>
188 getFTensor1FromMat<2, double, ublas::row_major, DoubleAllocator>(
189  MatrixDouble &data) {
190  if (data.size1() != 2) {
191  THROW_MESSAGE("Wrong size of data matrix");
192  }
193  return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2>(&data(0, 0),
194  &data(1, 0));
195 }
196 
197 // /**
198 // * @deprecated Name change to getFTensor1FromMat
199 // */
200 // template <int Tensor_Dim>
201 // DEPRECATED static inline FTensor::Tensor1<FTensor::PackPtr<double *, 1>,
202 // Tensor_Dim>
203 // getTensor1FormData(MatrixDouble &data) {
204 // return getFTensor1FromMat<Tensor_Dim>(data);
205 // }
206 
207 /**
208  * \brief Get tensor rank 2 (matrix) form data matrix
209  * \ingroup mofem_forces_and_sources_user_data_operators
210  */
211 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
212 static inline FTensor::Tensor2<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
213  Tensor_Dim1>
214 getFTensor2FromMat(ublas::matrix<T, L, A> &data) {
215  static_assert(!std::is_same<T, T>::value, "not implemented");
216 }
217 
218 /**
219  * Template specialization for getFTensor2FromMat
220  * \ingroup mofem_forces_and_sources_user_data_operators
221  *
222  */
223 template <>
226  if (data.size1() != 9) {
227  THROW_MESSAGE("Wrong size of data matrix; numer of rows is " +
228  boost::lexical_cast<std::string>(data.size1()));
229  }
231  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
232  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
233 }
234 
235 /**
236  * Template specialization for getFTensor2FromMat
237  */
238 template <>
241  if (data.size1() != 6) {
242  THROW_MESSAGE("Wrong size of data matrix");
243  }
245  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
246  &data(5, 0));
247 }
248 
249 /**
250  * \brief Get tensor rank 2 (matrix) form data matrix (specialization)
251  * \ingroup mofem_forces_and_sources_user_data_operators
252  */
253 template <int Tensor_Dim0, int Tensor_Dim1>
254 static inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim0,
255  Tensor_Dim1>
257  return getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
258  DoubleAllocator>(data);
259 }
260 
261 // /**
262 // * @deprecated Name change to getFTensor1FromMat
263 // */
264 // template <int Tensor_Dim0, int Tensor_Dim1>
265 // static inline DEPRECATED
266 // FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim0, Tensor_Dim1>
267 // getTensor2FormData(MatrixDouble &data) {
268 // return getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(data);
269 // }
270 
271 /**
272  * \brief Get symmetric tensor rank 2 (matrix) form data matrix
273  * \ingroup mofem_forces_and_sources_user_data_operators
274  */
275 template <int Tensor_Dim, class T, class L, class A>
277 getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
278  static_assert(!std::is_same<T, T>::value, "not implemented");
279 }
280 
281 /**
282  * @brief Get symmetric tensor rank 2 form matrix of for dimension 3
283  *
284  * Specialisation for symmetric tensor 2
285  *
286  * @tparam
287  * @param data
288  * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, 3>
289  */
290 template <>
293  if (data.size1() != 6) {
294  THROW_MESSAGE("Wrong size of data matrix");
295  }
297  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
298  &data(5, 0));
299 }
300 
301 /**
302  * @brief Get symmetric tensor rank 2 form matrix
303  *
304  * Specialisation for symmetric tensor 2
305  *
306  * @tparam Tensor_Dim
307  * @param data
308  * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, Tensor_Dim>
309  */
310 template <int Tensor_Dim>
312  Tensor_Dim>
315  DoubleAllocator>(data);
316 }
317 
318 /**
319  * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
320  *
321  * @tparam T
322  * @param t
323  * @return double
324  */
325 template <class T> static inline double dEterminant(T &t) {
326  return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
327  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
328  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
329 }
330 
331 /**
332  * \brief Calculate inverse of tensor rank 2 at integration points
333 
334  * \ingroup mofem_forces_and_sources
335  */
336 template <int Tensor_Dim, class T, class L, class A>
337 inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
338  ublas::vector<T, A> &det_data,
339  ublas::matrix<T, L, A> &inv_jac_data) {
341  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
342  "Specialization for this template not yet implemented");
344 }
345 
346 template <>
347 inline MoFEMErrorCode
348 invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
349  MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
350 
351 /**
352  * \brief Calculate determinant
353 
354  * \ingroup mofem_forces_and_sources
355  */
356 template <class T1, class T2>
357 inline MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det) {
359  det = +t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
360  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
361  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
363 }
364 
365 /**
366  * \brief Calculate matrix inverse
367 
368  * \ingroup mofem_forces_and_sources
369  */
370 template <class T1, class T2, class T3>
371 inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
373  const auto inv_det = 1. / det;
374  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
375  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
376  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
377  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
378  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
379  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
380  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
381  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
382  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
384 }
385 
386 #ifdef WITH_ADOL_C
387 
388 /**
389  * \brief Calculate matrix inverse, specialization for adouble tensor
390 
391  * \ingroup mofem_forces_and_sources
392  */
393 template <>
394 inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
399  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
400  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
401  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
402  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
403  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
404  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
405  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
406  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
407  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
409 }
410 
411 #endif
412 
413 /**
414  * \brief Calculate matrix inverse, specialization for symmetric tensor
415 
416  * \ingroup mofem_forces_and_sources
417  */
418 template <>
419 inline MoFEMErrorCode
420 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
425  const auto inv_det = 1. / det;
426  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
427  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
428  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
429  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
430  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
431  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
433 }
434 
435 #ifdef WITH_ADOL_C
436 
437 /**
438  * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
439 
440  * \ingroup mofem_forces_and_sources
441  */
442 template <>
443 inline MoFEMErrorCode
444 invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
449  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
450  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
451  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
452  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
453  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
454  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
456 }
457 
458 #endif
459 
460 /**
461  * \brief Calculate matrix inverse, specialization for symmetric (pointer)
462  tensor
463 
464  * \ingroup mofem_forces_and_sources
465  */
466 template <>
467 inline MoFEMErrorCode
468 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
473  const auto inv_det = 1. / det;
474  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
475  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
476  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
477  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
478  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
479  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
481 }
482 
483 } // namespace MoFEM
484 
485 #endif //__TEMPLATES_HPP__
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:337
static FTensor::Tensor2< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:214
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:500
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition: Types.hpp:110
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant.
Definition: Templates.hpp:357
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec< double, DoubleAllocator >(ublas::vector< double, DoubleAllocator > &data)
Definition: Templates.hpp:149
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition: Templates.hpp:87
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:619
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vectorExample how to use it.
Definition: Templates.hpp:143
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
static FTensor::Tensor2_symmetric< FTensor::PackPtr< T *, 1 >, Tensor_Dim > getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:277
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
KeyExtractor1 key1
Definition: Templates.hpp:96
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:107
KeyExtractor2 key2
Definition: Templates.hpp:97
unsigned int operator()(const id_type &value) const
Definition: Templates.hpp:113
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:325
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
std::string toString(X x)
Definition: Templates.hpp:118
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:101
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition: Types.hpp:103
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:69
result_type operator()(Arg &arg) const
Definition: Templates.hpp:91
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:76
KeyExtractor1::result_type result_type
Definition: Templates.hpp:85
static FTensor::Tensor1< FTensor::PackPtr< T *, 1 >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
Definition: Templates.hpp:160
std::vector< double, std::allocator< double > > DoubleAllocator
Definition: Types.hpp:70