v0.10.0
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 template <typename T> using ShardVec = boost::shared_ptr<std::vector<T>>;
24 
25 /**
26  * @brief Get Vector adaptor
27  *
28  * \code
29  *
30  * double *a;
31  * CHKERR VecGetArray(v,&a);
32  *
33  * for(int n = 0; n != nodes; ++n) {
34  *
35  * auto a = getVectorAdaptor(&a[3*n], 3);
36  * double dot = inner_prod(a, a);
37  *
38  * }
39  *
40  * CHKERR VecRetsoreArray(v,&a);
41  * \endcode
42  *
43  */
44 template <typename T1> inline auto getVectorAdaptor(T1 ptr, const size_t n) {
45  typedef typename std::remove_pointer<T1>::type T;
47  ublas::shallow_array_adaptor<T>(n, ptr));
48 };
49 
50 /**
51  * @brief Get Matrix adaptor
52  *
53  * \code
54  *
55  * double *a;
56  * CHKERR VecGetArray(v,&a);
57  *
58  * for(int n = 0; n != nodes; ++n) {
59  *
60  * auto F = getMatrixAdaptor(&a[3*3*n], 3, 3);
61  * MatrixDouble C = prod(F, trans(F));
62  *
63  * }
64  *
65  * CHKERR VecRetsoreArray(v,&a);
66  * \endcode
67  *
68  */
69 template <typename T1>
70 inline auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m) {
71  typedef typename std::remove_pointer<T1>::type T;
73  n, m, ublas::shallow_array_adaptor<T>(n * m, ptr));
74 };
75 
76 /**
77  * This small utility that cascades two key extractors will be
78  * used throughout the boost example
79  * <a
80  * href=http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp>
81  * http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp
82  * </a>
83  */
84 template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
85 public:
86  typedef typename KeyExtractor1::result_type result_type;
87 
88  KeyFromKey(const KeyExtractor1 &key1_ = KeyExtractor1(),
89  const KeyExtractor2 &key2_ = KeyExtractor2())
90  : key1(key1_), key2(key2_) {}
91 
92  template <typename Arg> result_type operator()(Arg &arg) const {
93  return key1(key2(arg));
94  }
95 
96 private:
97  KeyExtractor1 key1;
98  KeyExtractor2 key2;
99 };
100 
101 template <typename id_type> struct LtBit {
102  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
103  return valueA.to_ulong() < valueB.to_ulong();
104  }
105 };
106 
107 template <typename id_type> struct EqBit {
108  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
109  return valueA.to_ulong() == valueB.to_ulong();
110  }
111 };
112 
113 template <typename id_type> struct HashBit {
114  inline unsigned int operator()(const id_type &value) const {
115  return value.to_ulong();
116  }
117 };
118 
119 template <class X> inline std::string toString(X x) {
120  std::ostringstream buffer;
121  buffer << x;
122  return buffer.str();
123 }
124 
125 /**
126 * \brief Get tensor rank 0 (scalar) form data vector
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");
146 }
147 
148 template <>
151  ublas::vector<double, DoubleAllocator> &data) {
152  return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&*data.data().begin());
153 }
154 
155 template <int Tensor_Dim, int S, class T, class L, class A>
157 
158 template <int S>
161  get(MatrixDouble &data) {
162  if (data.size1() != 3)
164  "getFTensor1FromMat<3>: wrong size of data matrix, number of "
165  "rows should be 3 but is %d" +
166  boost::lexical_cast<std::string>(data.size1()));
167 
169  &data(0, 0), &data(1, 0), &data(2, 0));
170  }
171 };
172 
173 template <int S>
176  get(MatrixDouble &data) {
177  if (data.size1() != 2)
179  "getFTensor1FromMat<2>: wrong size of data matrix, number of "
180  "rows should be 2 but is %d" +
181  boost::lexical_cast<std::string>(data.size1()));
182  return FTensor::Tensor1<FTensor::PackPtr<double *, S>, 2>(&data(0, 0),
183  &data(1, 0));
184  }
185 };
186 
187 /**
188  * \brief Get tensor rank 1 (vector) form data matrix
189  */
190 template <int Tensor_Dim, int S = 1, class T, class L, class A>
192 getFTensor1FromMat(ublas::matrix<T, L, A> &data) {
193  static_assert(!std::is_same<T, T>::value, "not implemented");
194 }
195 
196 /**
197  * \brief Get tensor rank 1 (vector) form data matrix (specialization)
198  */
199 template <int Tensor_Dim, int S = 1>
202  return getFTensor1FromMatImpl<Tensor_Dim, S, double, ublas::row_major,
203  DoubleAllocator>::get(data);
204 }
205 
206 /**
207  * \brief Get tensor rank 2 (matrix) form data matrix
208  */
209 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
210 inline FTensor::Tensor2<FTensor::PackPtr<T *, 1>, Tensor_Dim0, Tensor_Dim1>
211 getFTensor2FromMat(ublas::matrix<T, L, A> &data) {
212  static_assert(!std::is_same<T, T>::value,
213  "Such getFTensor2FromMat specialisation is not implemented");
214 }
215 
216 /**
217  * Template specialization for getFTensor2FromMat
218  *
219  */
220 template <>
223  if (data.size1() != 9)
224  THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix; numer "
225  "of rows should be 9 but is " +
226  boost::lexical_cast<std::string>(data.size1()));
227 
229  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
230  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
231 }
232 
233 /**
234  * Template specialization for getFTensor2FromMat
235  */
236 template <>
239  if (data.size1() != 6)
240  THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix, numer "
241  "of rows should be 6 but is " +
242  boost::lexical_cast<std::string>(data.size1()));
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  * Template specialization for getFTensor2FromMat
251  */
252 template <>
255  if (data.size1() != 4)
256  THROW_MESSAGE("getFTensor2FromMat<2,2>: wrong size of data matrix, numer "
257  "of rows should be 4 but is " +
258  boost::lexical_cast<std::string>(data.size1()));
259 
261  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
262 }
263 
264 /**
265  * \brief Get tensor rank 2 (matrix) form data matrix (specialization)
266  */
267 template <int Tensor_Dim0, int Tensor_Dim1>
268 static inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim0,
269  Tensor_Dim1>
271  return getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
272  DoubleAllocator>(data);
273 }
274 
275 template <int Tensor_Dim, int S, class T, class L, class A>
277 
278 template <int S, class T, class L, class A>
281  get(ublas::matrix<T, L, A> &data) {
282  if (data.size1() != 6)
284  "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
285  "of rows should be 6 but is " +
286  boost::lexical_cast<std::string>(data.size1()));
287 
289  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
290  &data(5, 0));
291  }
292 };
293 
294 template <int S, class T, class L, class A>
297  get(ublas::matrix<T, L, A> &data) {
298  if (data.size1() != 3)
300  "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
301  "of rows should be 3 but is " +
302  boost::lexical_cast<std::string>(data.size1()));
304  &data(0, 0), &data(1, 0), &data(2, 0));
305  }
306 };
307 
308 /**
309  * \brief Get symmetric tensor rank 2 (matrix) form data matrix
310  */
311 template <int Tensor_Dim, int S, class T, class L, class A>
313 getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
315 }
316 
317 template <int Tensor_Dim, int S = 1>
319  Tensor_Dim>
321  return getFTensor2SymmetricFromMat<Tensor_Dim, S, double, ublas::row_major,
322  DoubleAllocator>(data);
323 }
324 
325 template <int Tensor_Dim01, int Tensor_Dim23, int S, class T, class L, class A>
327 
328 template <int S>
330  DoubleAllocator> {
331  static inline FTensor::Ddg<FTensor::PackPtr<double *, S>, 1, 1>
332  get(MatrixDouble &data) {
333  if (data.size1() != 1)
335  "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
336  "of rows should be 36 but is " +
337  boost::lexical_cast<std::string>(data.size1()));
338 
339  return FTensor::Ddg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
340  }
341 };
342 
343 template <int S>
345  DoubleAllocator> {
346  static inline FTensor::Ddg<FTensor::PackPtr<double *, S>, 2, 2>
347  get(MatrixDouble &data) {
348  if (data.size1() != 9)
350  "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
351  "of rows should be 36 but is " +
352  boost::lexical_cast<std::string>(data.size1()));
354  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
355  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
356  }
357 };
358 
359 template <int S>
361  DoubleAllocator> {
362  static inline FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3>
363  get(MatrixDouble &data) {
364  if (data.size1() != 36)
366  "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
367  "of rows should be 36 but is " +
368  boost::lexical_cast<std::string>(data.size1()));
370  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
371  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
372  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
373  &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
374  &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
375  &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
376  &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
377  &data(35, 0)};
378  }
379 };
380 
381 /**
382  * @brief Get symmetric tensor rank 4 on first two and last indices from
383  * form data matrix
384  *
385  * @tparam Tensor_Dim01 dimension of frirst two indicies
386  * @tparam Tensor_Dim23 dimension of second two indicies
387  * @tparam T the type of object stored
388  * @tparam L the storage organization
389  * @tparam A the type of Storage array
390  * @param data data container
391  * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
392  */
393 template <int Tensor_Dim01, int Tensor_Dim23, int S = 1, class T, class L,
394  class A>
395 static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
396 getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
397  static_assert(!std::is_same<T, T>::value,
398  "Such getFTensor4DdgFromMat specialisation is not implemented");
399 }
400 
401 template <int Tensor_Dim01, int Tensor_Dim23, int S = 1>
402 static inline FTensor::Ddg<FTensor::PackPtr<double *, S>, Tensor_Dim01,
403  Tensor_Dim23>
405  return getFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, double,
407  DoubleAllocator>::get(data);
408 }
409 
410 /**
411  * @brief Make Tensor1 from pointer
412  *
413  * @tparam DIM
414  * @param ptr
415  * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
416  */
417 template <int DIM>
419 getFTensor1FromPtr(double *ptr) {
420  static_assert(DIM != 3,
421  "Such getFTensor1FromPtr specialization is not implemented");
422 };
423 
424 template <>
426 getFTensor1FromPtr<3>(double *ptr) {
428  &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
429 };
430 
431 /**
432  * @brief Make Tensor2 from pointer
433  *
434  * @tparam DIM
435  * @param ptr
436  * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
437  */
438 template <int DIM1, int DIM2>
440 getFTensor2FromPtr(double *ptr) {
441  static_assert(DIM1 != 3, "Such getFTensor2FromPtr is not implemented");
442  static_assert(DIM2 >= 2 && DIM2 <= 3,
443  "Such getFTensor2FromPtr is not implemented");
444 };
445 
446 template <>
448  2> inline getFTensor2FromPtr<3, 2>(double *ptr) {
450  &ptr[HVEC0_0], &ptr[HVEC0_1],
451 
452  &ptr[HVEC1_0], &ptr[HVEC1_1],
453 
454  &ptr[HVEC2_0], &ptr[HVEC2_1]);
455 };
456 
457 template <>
459  3> inline getFTensor2FromPtr<3, 3>(double *ptr) {
461  &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
462 
463  &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
464 
465  &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
466 };
467 
468 /**
469  * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
470  *
471  * @tparam T
472  * @param t
473  * @return double
474  */
475 template <class T> static inline double dEterminant(T &t) {
476  return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
477  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
478  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
479 }
480 
481 /**
482  * \brief Calculate inverse of tensor rank 2 at integration points
483 
484  */
485 template <int Tensor_Dim, class T, class L, class A>
486 inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
487  ublas::vector<T, A> &det_data,
488  ublas::matrix<T, L, A> &inv_jac_data) {
490  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
491  "Specialization for this template not yet implemented");
493 }
494 
495 template <>
496 inline MoFEMErrorCode
497 invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
498  MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
499 
500 /**
501  * \brief Calculate determinant 3 by 3
502 
503  */
504 template <class T1, class T2>
505 inline MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det) {
507  det = +t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
508  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
509  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
511 }
512 
513 /**
514  * \brief Calculate determinant 2 by 2
515 
516  */
517 template <class T1, class T2>
518 inline MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det) {
520  det = t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
522 }
523 
524 /**
525  * \brief Calculate matrix inverse 3 by 3
526 
527  */
528 template <class T1, class T2, class T3>
529 inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
531  const auto inv_det = 1. / det;
532  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
533  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
534  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
535  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
536  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
537  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
538  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
539  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
540  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
542 }
543 
544 /**
545  * \brief Calculate matrix inverse 2 by 2
546 
547  */
548 template <class T1, class T2, class T3>
549 inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
551  const auto inv_det = 1. / det;
552  inv_t(0, 0) = t(1, 1) * inv_det;
553  inv_t(0, 1) = -t(0, 1) * inv_det;
554  inv_t(1, 0) = -t(1, 0) * inv_det;
555  inv_t(1, 1) = t(0, 0) * inv_det;
557 }
558 
559 #ifdef WITH_ADOL_C
560 
561 /**
562  * \brief Calculate matrix inverse, specialization for adouble tensor
563 
564  */
565 template <>
566 inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
571  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
572  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
573  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
574  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
575  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
576  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
577  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
578  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
579  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
581 }
582 
583 #endif
584 
585 /**
586  * \brief Calculate matrix inverse, specialization for symmetric tensor
587 
588  */
589 template <>
590 inline MoFEMErrorCode
591 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
596  const auto inv_det = 1. / det;
597  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
598  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
599  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
600  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
601  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
602  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
604 }
605 
606 #ifdef WITH_ADOL_C
607 
608 /**
609  * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
610 
611  */
612 template <>
613 inline MoFEMErrorCode
614 invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
619  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
620  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
621  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
622  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
623  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
624  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
626 }
627 
628 #endif
629 
630 /**
631  * \brief Calculate matrix inverse, specialization for symmetric (pointer)
632  tensor
633 
634  */
635 template <>
636 inline MoFEMErrorCode
637 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
642  const auto inv_det = 1. / det;
643  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
644  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
645  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
646  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
647  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
648  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
650 }
651 
652 /**
653  * @brief Extract entity handle form multi-index container
654  *
655  */
657  template <typename Iterator>
658  static inline EntityHandle extract(const Iterator &it) {
659  return (*it)->getEnt();
660  }
661 };
662 
663 /**
664  * @brief Insert ordered mofem multi-index into range
665  *
666  * \code
667  * auto hi_rit = refEntsPtr->upper_bound(start);
668  * auto hi_rit = refEntsPtr->upper_bound(end);
669  * Range to_erase;
670  * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
671  * \endcode
672  *
673  * @tparam Iterator
674  * @param r
675  * @param begin_iter
676  * @param end_iter
677  * @return moab::Range::iterator
678  */
679 template <typename Extractor, typename Iterator>
680 moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
681  Iterator end_iter) {
682  moab::Range::iterator hint = r.begin();
683  while (begin_iter != end_iter) {
684  size_t j = 0;
685  auto bi = Extractor::extract(begin_iter);
686  Iterator pj = begin_iter;
687  while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
688  ++pj;
689  ++j;
690  }
691  hint = r.insert(hint, bi, bi + (j - 1));
692  begin_iter = pj;
693  }
694  return hint;
695 };
696 
697 /**
698  * @brief Do nothing, used to rebuild database
699  *
700  */
702  Modify_change_nothing() = default;
703  template <typename T> inline void operator()(T &e) {}
704 };
705 
706 /**
707  * @brief Template used to reconstruct multi-index
708  *
709  * @tparam MI multi-index
710  * @tparam Modifier
711  * @param mi
712  * @param mo
713  * @return MoFEMErrorCode
714  */
715 template <typename MI, typename MO = Modify_change_nothing>
717  MO &&mo = Modify_change_nothing()) {
719  for (auto it = mi.begin(); it != mi.end(); ++it) {
720  if (!const_cast<MI &>(mi).modify(it, mo))
721  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
722  "Houston we have a problem");
723  }
725 }
726 
727 } // namespace MoFEM
728 
729 #endif //__TEMPLATES_HPP__
Do nothing, used to rebuild database.
Definition: Templates.hpp:701
const double T
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:143
FTensor::Tensor2< FTensor::PackPtr< double *, DIM1 *DIM2 >, DIM1, DIM2 > getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
Definition: Templates.hpp:440
static FTensor::Ddg< FTensor::PackPtr< double *, S >, 2, 2 > get(MatrixDouble &data)
Definition: Templates.hpp:347
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:505
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
Definition: Templates.hpp:396
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
Definition: Templates.hpp:192
static FTensor::Ddg< FTensor::PackPtr< double *, S >, 3, 3 > get(MatrixDouble &data)
Definition: Templates.hpp:363
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:658
FTensor::Tensor1< FTensor::PackPtr< double *, DIM >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition: Templates.hpp:419
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition: Types.hpp:115
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
static FTensor::Tensor2_symmetric< FTensor::PackPtr< T *, S >, 2 > get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:297
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec< double, DoubleAllocator >(ublas::vector< double, DoubleAllocator > &data)
Definition: Templates.hpp:150
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
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:211
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition: Templates.hpp:88
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static FTensor::Tensor1< FTensor::PackPtr< double *, S >, 2 > get(MatrixDouble &data)
Definition: Templates.hpp:176
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
KeyExtractor1 key1
Definition: Templates.hpp:97
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:108
static Index< 'n', 3 > n
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:680
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:518
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:486
Extract entity handle form multi-index container.
Definition: Templates.hpp:656
static FTensor::Tensor2_symmetric< FTensor::PackPtr< T *, S >, 3 > get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:281
KeyExtractor2 key2
Definition: Templates.hpp:98
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
Definition: Templates.hpp:426
static Index< 'm', 3 > m
unsigned int operator()(const id_type &value) const
Definition: Templates.hpp:114
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:459
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:549
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:475
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
static Index< 'j', 3 > j
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:44
static FTensor::Ddg< FTensor::PackPtr< double *, S >, 1, 1 > get(MatrixDouble &data)
Definition: Templates.hpp:332
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:716
std::string toString(X x)
Definition: Templates.hpp:119
static FTensor::Tensor1< FTensor::PackPtr< double *, S >, 3 > get(MatrixDouble &data)
Definition: Templates.hpp:161
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:102
const double r
rate factor
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:448
Index< 'L', 3 > L
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition: Types.hpp:108
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:70
result_type operator()(Arg &arg) const
Definition: Templates.hpp:92
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
KeyExtractor1::result_type result_type
Definition: Templates.hpp:86
std::vector< double, std::allocator< double > > DoubleAllocator
Definition: Types.hpp:70
static FTensor::Tensor2_symmetric< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:313
boost::shared_ptr< std::vector< T > > ShardVec
Definition: Templates.hpp:23