v0.14.0
Templates.hpp
Go to the documentation of this file.
1 /** \file Templates.hpp
2  * \brief Templates declarations
3  */
4 
5 #ifndef __TEMPLATES_HPP__
6 #define __TEMPLATES_HPP__
7 
8 namespace MoFEM {
9 
10 template <typename T> using ShardVec = boost::shared_ptr<std::vector<T>>;
11 
12 /**
13  * @brief Get Vector adaptor
14  *
15  * \code
16  *
17  * double *a;
18  * CHKERR VecGetArray(v,&a);
19  *
20  * for(int n = 0; n != nodes; ++n) {
21  *
22  * auto a = getVectorAdaptor(&a[3*n], 3);
23  * double dot = inner_prod(a, a);
24  *
25  * }
26  *
27  * CHKERR VecRetsoreArray(v,&a);
28  * \endcode
29  *
30  */
31 template <typename T1> inline auto getVectorAdaptor(T1 ptr, const size_t n) {
32  typedef typename std::remove_pointer<T1>::type T;
34  ublas::shallow_array_adaptor<T>(n, ptr));
35 };
36 
37 /**
38  * @brief Get Matrix adaptor
39  *
40  * \code
41  *
42  * double *a;
43  * CHKERR VecGetArray(v,&a);
44  *
45  * for(int n = 0; n != nodes; ++n) {
46  *
47  * auto F = getMatrixAdaptor(&a[3*3*n], 3, 3);
48  * MatrixDouble C = prod(F, trans(F));
49  *
50  * }
51  *
52  * CHKERR VecRetsoreArray(v,&a);
53  * \endcode
54  *
55  */
56 template <typename T1>
57 inline auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m) {
58  typedef typename std::remove_pointer<T1>::type T;
60  n, m, ublas::shallow_array_adaptor<T>(n * m, ptr));
61 };
62 
63 /**
64  * This small utility that cascades two key extractors will be
65  * used throughout the boost example
66  * <a
67  * href=http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp>
68  * http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp
69  * </a>
70  */
71 template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
72 public:
73  typedef typename KeyExtractor1::result_type result_type;
74 
75  KeyFromKey(const KeyExtractor1 &key1_ = KeyExtractor1(),
76  const KeyExtractor2 &key2_ = KeyExtractor2())
77  : key1(key1_), key2(key2_) {}
78 
79  template <typename Arg> result_type operator()(Arg &arg) const {
80  return key1(key2(arg));
81  }
82 
83 private:
84  KeyExtractor1 key1;
85  KeyExtractor2 key2;
86 };
87 
88 template <typename id_type> struct LtBit {
89  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
90  return valueA.to_ulong() < valueB.to_ulong();
91  }
92 };
93 
94 template <typename id_type> struct EqBit {
95  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
96  return valueA.to_ulong() == valueB.to_ulong();
97  }
98 };
99 
100 template <typename id_type> struct HashBit {
101  inline unsigned int operator()(const id_type &value) const {
102  return value.to_ulong();
103  }
104 };
105 
106 template <class X> inline std::string toString(X x) {
107  std::ostringstream buffer;
108  buffer << x;
109  return buffer.str();
110 }
111 
112 template <int S, class T, class A> struct GetFTensor0FromVecImpl {
113  static inline auto get(ublas::vector<T, A> &data) {
114  return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(&*data.data().begin());
115  }
116 };
117 
118 /**
119 * \brief Get tensor rank 0 (scalar) form data vector
120 
121 Example how to use it.
122 \code
123 VectorDouble vec;
124 vec.resize(nb_gauss_pts,false);
125 vec.clear();
126 auto t0 = getFTensor0FromData(data);
127 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
128 
129  ++t0;
130 }
131 \endcode
132 
133 */
134 template <int S = 1, class T, class A>
135 static inline auto getFTensor0FromVec(ublas::vector<T, A> &data) {
137 }
138 
139 template <int Tensor_Dim, int S, class T, class L, class A>
141 
142 template <int S, class T, class A>
143 struct GetFTensor1FromMatImpl<3, S, T, ublas::row_major, A> {
144  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
145 #ifndef NDDEBUG
146  if (data.size1() != 3)
148  "getFTensor1FromMat<3>: wrong size of data matrix, number of "
149  "rows should be 3 but is " +
150  boost::lexical_cast<std::string>(data.size1()));
151 #endif
153  &data(0, 0), &data(1, 0), &data(2, 0));
154  }
155 };
156 
157 template <int S, class T, class A>
158 struct GetFTensor1FromMatImpl<4, S, T, ublas::row_major, A> {
159  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
160 #ifndef NDDEBUG
161  if (data.size1() != 4)
163  "getFTensor1FromMat<4>: wrong size of data matrix, number of "
164  "rows should be 4 but is " +
165  boost::lexical_cast<std::string>(data.size1()));
166 #endif
168  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
169  }
170 };
171 
172 template <int S, class T, class A>
173 struct GetFTensor1FromMatImpl<6, S, T, ublas::row_major, A> {
174  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
175 #ifndef NDDEBUG
176  if (data.size1() != 6)
178  "getFTensor1FromMat<6>: wrong size of data matrix, number of "
179  "rows should be 6 but is " +
180  boost::lexical_cast<std::string>(data.size1()));
181 #endif
183  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
184  &data(5, 0));
185  }
186 };
187 
188 template <int S, class T, class A>
189 struct GetFTensor1FromMatImpl<9, S, T, ublas::row_major, A> {
190  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
191 #ifndef NDDEBUG
192  if (data.size1() != 9)
194  "getFTensor1FromMat<9>: wrong size of data matrix, number of "
195  "rows should be 9 but is " +
196  boost::lexical_cast<std::string>(data.size1()));
197 #endif
199  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
200  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
201  }
202 };
203 
204 template <int S, class T, class A>
205 struct GetFTensor1FromMatImpl<2, S, T, ublas::row_major, A> {
206  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
207 #ifndef NDDEBUG
208  if (data.size1() != 2)
210  "getFTensor1FromMat<2>: wrong size of data matrix, number of "
211  "rows should be 2 but is " +
212  boost::lexical_cast<std::string>(data.size1()));
213 #endif
214  return FTensor::Tensor1<FTensor::PackPtr<double *, S>, 2>(&data(0, 0),
215  &data(1, 0));
216  }
217 };
218 
219 template <int S, class T, class A>
220 struct GetFTensor1FromMatImpl<1, S, T, ublas::row_major, A> {
221  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
222 #ifndef NDEBUG
223  if (data.size1() != 1)
225  "getFTensor1FromMat<1>: wrong size of data matrix, number of "
226  "rows should be 1 but is " +
227  boost::lexical_cast<std::string>(data.size1()));
228 #endif
229  return FTensor::Tensor1<FTensor::PackPtr<double *, S>, 1>(&data(0, 0));
230  }
231 };
232 
233 /**
234  * \brief Get tensor rank 1 (vector) form data matrix
235  */
236 template <int Tensor_Dim, int S = 1, class T, class L, class A>
238 getFTensor1FromMat(ublas::matrix<T, L, A> &data) {
239  static_assert(!std::is_same<T, T>::value, "not implemented");
240 }
241 
242 /**
243  * \brief Get tensor rank 1 (vector) form data matrix (specialization)
244  */
245 template <int Tensor_Dim, int S = 1>
246 inline auto getFTensor1FromMat(MatrixDouble &data) {
247  return GetFTensor1FromMatImpl<Tensor_Dim, S, double, ublas::row_major,
248  DoubleAllocator>::get(data);
249 }
250 
251 template <int Tensor_Dim1, int Tensor_Dim2, int S, class T, class L, class A>
253  static inline auto get(ublas::matrix<T, L, A> &data) {
254 #ifndef NDEBUG
255  if (data.size1() != Tensor_Dim1 * Tensor_Dim2) {
257  "getFTensor2FromMat<" +
258  boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
259  boost::lexical_cast<std::string>(Tensor_Dim2) +
260  ">: wrong size of rows in data matrix, should be " +
261  boost::lexical_cast<std::string>(Tensor_Dim1 * Tensor_Dim2) +
262  " but is " + boost::lexical_cast<std::string>(data.size1()));
263  }
264 #endif
265  std::array<double *, Tensor_Dim1 * Tensor_Dim2> ptrs;
266  for (auto i = 0; i != Tensor_Dim1 * Tensor_Dim2; ++i)
267  ptrs[i] = &data(i, 0);
269  Tensor_Dim2>(ptrs);
270  }
271 };
272 
273 /**
274  * \brief Get tensor rank 2 (matrix) form data matrix
275  */
276 template <int Tensor_Dim1, int Tensor_Dim2>
277 inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim1, Tensor_Dim2>
279  return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, 1, double,
280  ublas::row_major, DoubleAllocator>::get(data);
281 }
282 
283 template <int Tensor_Dim1, int Tensor_Dim2>
284 inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim1, Tensor_Dim2>
286 
287 /**
288  * Template specialization for getFTensor2FromMat
289  */
290 template <>
293  return FTensor::Tensor2<FTensor::PackPtr<double *, 1>, 1, 1>(&*data.begin());
294 }
295 
296 template <int Tensor_Dim, int S, class T, class L, class A>
298 
299 template <int S, class T, class L, class A>
301  static inline auto get(ublas::matrix<T, L, A> &data) {
302 #ifndef NDEBUG
303  if (data.size1() != 6)
305  "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
306  "of rows should be 6 but is " +
307  boost::lexical_cast<std::string>(data.size1()));
308 #endif
310  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
311  &data(5, 0));
312  }
313 };
314 
315 template <int S, class T, class L, class A>
317  static inline auto get(ublas::matrix<T, L, A> &data) {
318 #ifndef NDEBUG
319  if (data.size1() != 3)
321  "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
322  "of rows should be 3 but is " +
323  boost::lexical_cast<std::string>(data.size1()));
324 #endif
326  &data(0, 0), &data(1, 0), &data(2, 0));
327  }
328 };
329 
330 /**
331  * \brief Get symmetric tensor rank 2 (matrix) form data matrix
332  */
333 template <int Tensor_Dim, int S, class T, class L, class A>
334 static inline auto getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
336 }
337 
338 template <int Tensor_Dim, int S = 1>
339 static inline auto getFTensor2SymmetricFromMat(MatrixDouble &data) {
340  return getFTensor2SymmetricFromMat<Tensor_Dim, S, double, ublas::row_major,
341  DoubleAllocator>(data);
342 }
343 
344 template <int Tensor_Dim01, int Tensor_Dim23, int S, class T, class L, class A>
346 
347 template <int S, class T, class A>
348 struct GetFTensor4DdgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
349  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
350 #ifndef NDEBUG
351  if (data.size1() != 1)
353  "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
354  "of rows should be 1 but is " +
355  boost::lexical_cast<std::string>(data.size1()));
356 #endif
357  return FTensor::Ddg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
358  }
359 };
360 
361 template <int S, class T, class A>
362 struct GetFTensor4DdgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
363  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
364 #ifndef NDEBUG
365  if (data.size1() != 9) {
367  "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
368  "of rows should be 9 but is " +
369  boost::lexical_cast<std::string>(data.size1()));
370  }
371 #endif
373  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
374  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
375  }
376 };
377 
378 template <int S, class T, class A>
379 struct GetFTensor4DdgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
380  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
381 #ifndef NDEBUG
382  if (data.size1() != 36) {
383  cerr << data.size1() << endl;
385  "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
386  "of rows should be 36 but is " +
387  boost::lexical_cast<std::string>(data.size1()));
388  }
389 #endif
391  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
392  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
393  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
394  &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
395  &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
396  &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
397  &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
398  &data(35, 0)};
399  }
400 };
401 
402 /**
403  * @brief Get symmetric tensor rank 4 on first two and last indices from
404  * form data matrix
405  *
406  * @tparam Tensor_Dim01 dimension of first two indicies
407  * @tparam Tensor_Dim23 dimension of second two indicies
408  * @tparam T the type of object stored
409  * @tparam L the storage organization
410  * @tparam A the type of Storage array
411  * @param data data container
412  * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
413  */
414 template <int Tensor_Dim01, int Tensor_Dim23, int S = 1, class T, class L,
415  class A>
416 static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
417 getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
418  static_assert(!std::is_same<T, T>::value,
419  "Such getFTensor4DdgFromMat specialisation is not implemented");
420 }
421 
422 template <int Tensor_Dim01, int Tensor_Dim23, int S = 1>
423 static inline auto getFTensor4DdgFromMat(MatrixDouble &data) {
424  return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, double,
426  DoubleAllocator>::get(data);
427 }
428 
429 template <int Tensor_Dim01, int Tensor_Dim2, int S, class T, class L, class A>
431 
432 template <int S, class T, class A>
433 struct GetFTensor3DgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
434  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
435 #ifndef NDEBUG
436  if (data.size1() != 1)
438  "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
439  "of rows should be 1 but is " +
440  boost::lexical_cast<std::string>(data.size1()));
441 #endif
442  return FTensor::Dg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
443  }
444 };
445 
446 template <int S, class T, class A>
447 struct GetFTensor3DgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
448  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
449 #ifndef NDEBUG
450  if (data.size1() != 6) {
452  "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
453  "of rows should be 6 but is " +
454  boost::lexical_cast<std::string>(data.size1()));
455  }
456 #endif
458  &data(0, 0), &data(1, 0), &data(2, 0),
459  &data(3, 0), &data(4, 0), &data(5, 0)};
460  }
461 };
462 
463 template <int S, class T, class A>
464 struct GetFTensor3DgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
465  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
466 #ifndef NDEBUG
467  if (data.size1() != 18) {
468  cerr << data.size1() << endl;
470  "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
471  "of rows should be 18 but is " +
472  boost::lexical_cast<std::string>(data.size1()));
473  }
474 #endif
476  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
477  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
478  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
479  &data(15, 0), &data(16, 0), &data(17, 0)};
480  }
481 };
482 
483 /**
484  * @brief Get symmetric tensor rank 3 on the first two indices from
485  * form data matrix
486  *
487  * @tparam Tensor_Dim01 dimension of first two indicies
488  * @tparam Tensor_Dim2 dimension of last index
489  * @tparam T the type of object stored
490  * @tparam L the storage organization
491  * @tparam A the type of Storage array
492  * @param data data container
493  * @return FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
494  */
495 template <int Tensor_Dim01, int Tensor_Dim2, int S = 1, class T, class L,
496  class A>
497 static inline FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim2>
498 getFTensor3DgFromMat(ublas::matrix<T, L, A> &data) {
499  static_assert(!std::is_same<T, T>::value,
500  "Such getFTensor3DgFromMat specialisation is not implemented");
501 }
502 
503 template <int Tensor_Dim01, int Tensor_Dim2, int S = 1>
504 static inline auto getFTensor3DgFromMat(MatrixDouble &data) {
505  return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, double,
506  ublas::row_major, DoubleAllocator>::get(data);
507 }
508 
509 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
510  int S, class T, class L, class A>
512 
513 template <int S, class T, class A>
514 struct GetFTensor4FromMatImpl<1, 1, 1, 1, S, T, ublas::row_major, A> {
515  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
516 #ifndef NDEBUG
517  if (data.size1() != 1)
519  "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
520  "of rows should be 1 but is " +
521  boost::lexical_cast<std::string>(data.size1()));
522 #endif
524  &data(0, 0)};
525  }
526 };
527 
528 template <int S, class T, class A>
529 struct GetFTensor4FromMatImpl<2, 2, 2, 2, S, T, ublas::row_major, A> {
530  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
531 #ifndef NDEBUG
532  if (data.size1() != 16) {
534  "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
535  "of rows should be 16 but is " +
536  boost::lexical_cast<std::string>(data.size1()));
537  }
538 #endif
540  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
541  &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
542  &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
543  &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
544  }
545 };
546 
547 template <int S, class T, class A>
548 struct GetFTensor4FromMatImpl<3, 3, 3, 3, S, T, ublas::row_major, A> {
549  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
550 #ifndef NDEBUG
551  if (data.size1() != 81) {
552  cerr << data.size1() << endl;
554  "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
555  "of rows should be 81 but is " +
556  boost::lexical_cast<std::string>(data.size1()));
557  }
558 #endif
560  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
561  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
562  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
563  &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
564  &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
565  &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
566  &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
567  &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
568  &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
569  &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
570  &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
571  &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
572  &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
573  &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
574  &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
575  &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
576  &data(80, 0)};
577  }
578 };
579 
580 /**
581  * @brief Get tensor rank 4 (non symmetric) form data matrix
582  *
583  * @tparam Tensor_Dim0 dimension of frirst index
584  * @tparam Tensor_Dim1 dimension of second index
585  * @tparam Tensor_Dim2 dimension of third index
586  * @tparam Tensor_Dim3 dimension of fourth index
587  * @tparam T the type of object stored
588  * @tparam L the storage organization
589  * @tparam A the type of Storage array
590  * @param data data container
591  * @return FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
592  Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
593  */
594 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
595  int S = 1, class T, class L, class A>
596 static inline FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
597  Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
598 getFTensor4FromMat(ublas::matrix<T, L, A> &data) {
599  static_assert(!std::is_same<T, T>::value,
600  "Such getFTensor4FromMat specialisation is not implemented");
601 }
602 
603 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
604  int S = 1>
605 static inline auto getFTensor4FromMat(MatrixDouble &data) {
606  return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
607  Tensor_Dim3, S, double, ublas::row_major,
608  DoubleAllocator>::get(data);
609 }
610 
611 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class T,
612  class L, class A>
614 
615 template <int S, class T, class A>
616 struct GetFTensor3FromMatImpl<1, 1, 1, S, T, ublas::row_major, A> {
617  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
618 #ifndef NDEBUG
619  if (data.size1() != 1)
621  "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
622  "of rows should be 1 but is " +
623  boost::lexical_cast<std::string>(data.size1()));
624 #endif
626  &data(0, 0)};
627  }
628 };
629 
630 template <int S, class T, class A>
631 struct GetFTensor3FromMatImpl<2, 2, 2, S, T, ublas::row_major, A> {
632  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
633 #ifndef NDEBUG
634  if (data.size1() != 8)
636  "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
637  "of rows should be 8 but is " +
638  boost::lexical_cast<std::string>(data.size1()));
639 #endif
641  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
642  &data(5, 0), &data(6, 0), &data(7, 0)
643 
644  };
645  }
646 };
647 
648 template <int S, class T, class A>
649 struct GetFTensor3FromMatImpl<3, 2, 2, S, T, ublas::row_major, A> {
650  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
651 #ifndef NDEBUG
652  if (data.size1() != 12)
654  "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
655  "of rows should be 12 but is " +
656  boost::lexical_cast<std::string>(data.size1()));
657 #endif
659  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
660  &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
661  &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
662  }
663 };
664 
665 template <int S, class T, class A>
666 struct GetFTensor3FromMatImpl<2, 2, 3, S, T, ublas::row_major, A> {
667  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
668 #ifndef NDEBUG
669  if (data.size1() != 12)
671  "getFTensor3FromMat<2, 2, 3>: wrong size of data matrix, number "
672  "of rows should be 12 but is " +
673  boost::lexical_cast<std::string>(data.size1()));
674 #endif
676  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
677  &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
678  &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
679  }
680 };
681 
682 template <int S, class T, class A>
683 struct GetFTensor3FromMatImpl<3, 3, 3, S, T, ublas::row_major, A> {
684  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
685 #ifndef NDEBUG
686  if (data.size1() != 27)
688  "getFTensor3FromMat<3, 3, 3>: wrong size of data matrix, number "
689  "of rows should be 27 but is " +
690  boost::lexical_cast<std::string>(data.size1()));
691 #endif
693  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
694  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
695  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
696  &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
697  &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
698  &data(25, 0), &data(26, 0)};
699  }
700 };
701 
702 template <int S, class T, class A>
703 struct GetFTensor3FromMatImpl<6, 3, 3, S, T, ublas::row_major, A> {
704  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
705 #ifndef NDEBUG
706  if (data.size1() != 54)
708  "getFTensor3FromMat<6, 3, 3>: wrong size of data matrix, number "
709  "of rows should be 54 but is " +
710  boost::lexical_cast<std::string>(data.size1()));
711 #endif
712  std::array<double *, 54> ptrs;
713  for (auto i = 0; i != 54; ++i)
714  ptrs[i] = &data(i, 0);
715  return FTensor::Tensor3<FTensor::PackPtr<double *, S>, 6, 3, 3>(ptrs);
716  }
717 };
718 
719 template <int S, class T, class A>
720 struct GetFTensor3FromMatImpl<3, 3, 6, S, T, ublas::row_major, A> {
721  static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
722 #ifndef NDEBUG
723  if (data.size1() != 54)
725  "getFTensor3FromMat<3, 3, 6>: wrong size of data matrix, number "
726  "of rows should be 54 but is " +
727  boost::lexical_cast<std::string>(data.size1()));
728 #endif
729  std::array<double *, 54> ptrs;
730  for (auto i = 0; i != 54; ++i)
731  ptrs[i] = &data(i, 0);
732  return FTensor::Tensor3<FTensor::PackPtr<double *, S>, 3, 3, 6>(ptrs);
733  }
734 };
735 
736 /**
737  * @brief Get tensor rank 3 (non symmetries) form data matrix
738  *
739  * @tparam Tensor_Dim0 dimension of first index
740  * @tparam Tensor_Dim1 dimension of second index
741  * @tparam Tensor_Dim2 dimension of third index
742  * @tparam S shift size
743  * @tparam T the type of object stored
744  * @tparam L the storage organization
745  * @tparam A the type of Storage array
746  * @param data data container
747  * @return FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
748  Tensor_Dim1, Tensor_Dim2>
749  */
750 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1, class T,
751  class L, class A>
752 static inline FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
753  Tensor_Dim1, Tensor_Dim2>
754 getFTensor3FromMat(ublas::matrix<T, L, A> &data) {
755  static_assert(!std::is_same<T, T>::value,
756  "Such getFTensor3FromMat specialisation is not implemented");
757 }
758 
759 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1>
760 static inline auto getFTensor3FromMat(MatrixDouble &data) {
761  return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
763  DoubleAllocator>::get(data);
764 }
765 
766 template <int DIM, int S, typename T> struct GetFTensor1FromPtrImpl;
767 
768 template <int S, typename T> struct GetFTensor1FromPtrImpl<2, S, T> {
769  GetFTensor1FromPtrImpl() = delete;
770  inline static auto get(T *ptr) {
772  &ptr[HVEC1]);
773  }
774 };
775 
776 template <int S, typename T> struct GetFTensor1FromPtrImpl<3, S, T> {
777  GetFTensor1FromPtrImpl() = delete;
778  inline static auto get(T *ptr) {
780  &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
781  }
782 };
783 
784 template <int S, typename T> struct GetFTensor1FromPtrImpl<4, S, T> {
785  GetFTensor1FromPtrImpl() = delete;
786  inline static auto get(T *ptr) {
787  return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 4>(&ptr[0], &ptr[1],
788  &ptr[2], &ptr[3]);
789  }
790 };
791 
792 template <int S, typename T> struct GetFTensor1FromPtrImpl<6, S, T> {
793  GetFTensor1FromPtrImpl() = delete;
794  inline static auto get(T *ptr) {
796  &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
797  }
798 };
799 
800 /**
801  * @brief Make Tensor1 from pointer
802  *
803  * @tparam DIM
804  * @param ptr
805  * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
806  */
807 template <int DIM, int S = DIM>
809 getFTensor1FromPtr(double *ptr) {
811 };
812 
813 #ifdef WITH_ADOL_C
814 template <int DIM, int S = DIM>
818 };
819 #endif
820 
821 template <int DIM, int S = DIM>
823 getFTensor1FromPtr(std::complex<double> *ptr) {
825 };
826 
827 template <int DIM1, int DIM2, int S, typename T> struct GetFTensor2FromPtr;
828 
829 template <int S, typename T> struct GetFTensor2FromPtr<3, 2, S, T> {
830  GetFTensor2FromPtr() = delete;
831  inline static auto get(T *ptr) {
833  &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
834  }
835 };
836 
837 template <int S, typename T> struct GetFTensor2FromPtr<6, 6, S, T> {
838  GetFTensor2FromPtr() = delete;
839  inline static auto get(T *ptr) {
841  &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
842  &ptr[8], &ptr[9], &ptr[10], &ptr[11], &ptr[12], &ptr[13], &ptr[14],
843  &ptr[15], &ptr[16], &ptr[17], &ptr[18], &ptr[19], &ptr[20], &ptr[21],
844  &ptr[22], &ptr[23], &ptr[24], &ptr[25], &ptr[26], &ptr[27], &ptr[28],
845  &ptr[29], &ptr[30], &ptr[31], &ptr[32], &ptr[33], &ptr[34], &ptr[35]);
846  }
847 };
848 
849 template <int S, typename T> struct GetFTensor2FromPtr<3, 3, S, T> {
850  GetFTensor2FromPtr() = delete;
851  inline static auto get(T *ptr) {
853  &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
854  &ptr[8]);
855  }
856 };
857 
858 template <int S, typename T> struct GetFTensor2FromPtr<2, 2, S, T> {
859  GetFTensor2FromPtr() = delete;
860  inline static auto get(T *ptr) {
861  return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 2, 2>(&ptr[0], &ptr[1],
862  &ptr[2], &ptr[3]);
863  }
864 };
865 
866 template <int S, typename T> struct GetFTensor2FromPtr<1, 3, S, T> {
867  GetFTensor2FromPtr() = delete;
868  inline static auto get(T *ptr) {
869  return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 3>(&ptr[0], &ptr[1],
870  &ptr[2]);
871  }
872 };
873 
874 template <int S, typename T> struct GetFTensor2FromPtr<1, 2, S, T> {
875  GetFTensor2FromPtr() = delete;
876  inline static auto get(T *ptr) {
877  return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 2>(&ptr[0], &ptr[1]);
878  }
879 };
880 
881 template <int S, typename T> struct GetFTensor2FromPtr<1, 1, S, T> {
882  GetFTensor2FromPtr() = delete;
883  inline static auto get(T *ptr) {
884  return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 1>(&ptr[0]);
885  }
886 };
887 
888 /**
889  * @brief Make Tensor2 from pointer
890  *
891  * @tparam DIM
892  * @param ptr
893  * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
894  */
895 template <int DIM1, int DIM2, int S = DIM1 * DIM2>
896 inline auto getFTensor2FromPtr(double *ptr) {
898 };
899 
900 /**
901  * @brief Make Tensor2 from pointer
902  *
903  * @tparam DIM
904  * @param ptr
905  * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
906  */
907 template <int DIM1, int DIM2, int S = DIM1 * DIM2>
908 inline auto getFTensor2FromPtr(std::complex<double> *ptr) {
910 };
911 
912 /**
913  * @brief Make Tensor2 for HVec base from pointer
914  *
915  * @tparam DIM
916  * @param ptr
917  * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
918  */
919 template <int DIM1, int DIM2>
922  static_assert(DIM1 == DIM1 || DIM2 != DIM2,
923  "Such getFTensor2HVecFromPtr is not implemented");
924 };
925 
926 template <>
928  2> inline getFTensor2HVecFromPtr<3, 2>(double *ptr) {
930  &ptr[HVEC0_0], &ptr[HVEC0_1],
931 
932  &ptr[HVEC1_0], &ptr[HVEC1_1],
933 
934  &ptr[HVEC2_0], &ptr[HVEC2_1]);
935 };
936 
937 template <>
939  3> inline getFTensor2HVecFromPtr<3, 3>(double *ptr) {
941  &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
942 
943  &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
944 
945  &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
946 };
947 
948 /*
949  * @brief Make Tensor3 from pointer
950  *
951  * @tparam DIM
952  * @param ptr
953  * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3>, DIM1,
954  * DIM2, DIM3>
955  */
956 template <int DIM1, int DIM2, int DIM3>
958  DIM2, DIM3>
959 getFTensor3FromPtr(double *ptr) {
960  static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
961  "Such getFTensor3FromPtr is not implemented");
962 };
963 
964 template <>
968  &ptr[0], &ptr[1], &ptr[2],
969 
970  &ptr[3], &ptr[4], &ptr[5],
971 
972  &ptr[6], &ptr[7], &ptr[8],
973 
974  &ptr[9], &ptr[10], &ptr[11]
975 
976  );
977 };
978 
979 template <>
983  &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
984  &ptr[8], &ptr[9], &ptr[10], &ptr[11], &ptr[12], &ptr[13], &ptr[14],
985  &ptr[15], &ptr[16], &ptr[17], &ptr[18], &ptr[19], &ptr[20], &ptr[21],
986  &ptr[22], &ptr[23], &ptr[24], &ptr[25], &ptr[26]);
987 };
988 
989 /**
990  * @brief Make symmetric Tensor2 from pointer
991  *
992  * @tparam DIM
993  * @param ptr
994  * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
995  */
996 template <int DIM>
998  FTensor::PackPtr<double *, (DIM * (DIM + 1)) / 2>, DIM>
1000  static_assert(DIM, "Such getFTensor2SymmetricFromPtr is not implemented");
1001 }
1002 
1003 template <>
1007  &ptr[0], &ptr[1], &ptr[2],
1008 
1009  &ptr[3], &ptr[4],
1010 
1011  &ptr[5]);
1012 };
1013 
1014 template <>
1018  &ptr[0], &ptr[1], &ptr[2]);
1019 };
1020 
1021 #ifdef WITH_ADOL_C
1022 
1023 /**
1024  * @brief Make symmetric Tensor2 from pointer
1025  *
1026  * @tparam DIM
1027  * @param ptr
1028  * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
1029  */
1030 template <int DIM>
1032  FTensor::PackPtr<adouble *, (DIM * (DIM + 1)) / 2>, DIM>
1034  static_assert(DIM, "Such getFTensor2SymmetricFromPtr is not implemented");
1035 }
1036 
1037 template <>
1041  &ptr[0], &ptr[1], &ptr[2],
1042 
1043  &ptr[3], &ptr[4],
1044 
1045  &ptr[5]);
1046 };
1047 
1048 template <>
1052  &ptr[0], &ptr[1], &ptr[2]);
1053 };
1054 
1055 #endif
1056 
1057 /**
1058  * @brief Make symmetric Tensor2 from pointer, taking lower triangle of matrix
1059  *
1060  * @tparam DIM
1061  * @param ptr
1062  * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
1063  */
1064 template <int DIM>
1067  static_assert(DIM,
1068  "Such getFTensor2SymmetricLowerFromPtr is not implemented");
1069 }
1070 
1071 template <>
1075  &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
1076 
1077  &ptr[HVEC1_0], &ptr[HVEC1_1],
1078 
1079  &ptr[HVEC2_2]);
1080 };
1081 
1082 template <>
1086  &ptr[0], &ptr[1], &ptr[3]);
1087 };
1088 
1089 template <int DIM, int S> struct GetFTensor1FromArray;
1090 
1091 template <int S> struct GetFTensor1FromArray<1, S> {
1092  GetFTensor1FromArray() = delete;
1093  template <typename V> static inline auto get(V &data) {
1094  using T = typename std::remove_reference<decltype(data[0])>::type;
1095  return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 1>{&data[0]};
1096  }
1097 };
1098 
1099 template <int S> struct GetFTensor1FromArray<2, S> {
1100  GetFTensor1FromArray() = delete;
1101  template <typename V> static inline auto get(V &data) {
1102  using T = typename std::remove_reference<decltype(data[0])>::type;
1103  return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 2>{&data[0], &data[1]};
1104  }
1105 };
1106 
1107 template <int S> struct GetFTensor1FromArray<3, S> {
1108  GetFTensor1FromArray() = delete;
1109  template <typename V> static inline auto get(V &data) {
1110  using T = typename std::remove_reference<decltype(data[0])>::type;
1111  return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 3>{&data[0], &data[1],
1112  &data[2]};
1113  }
1114 };
1115 
1116 template <int S> struct GetFTensor1FromArray<4, S> {
1117  GetFTensor1FromArray() = delete;
1118  template <typename V> static inline auto get(V &data) {
1119  using T = typename std::remove_reference<decltype(data[0])>::type;
1120  return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 4>{&data[0], &data[1],
1121  &data[2], &data[3]};
1122  }
1123 };
1124 
1125 template <int S> struct GetFTensor1FromArray<6, S> {
1126  GetFTensor1FromArray() = delete;
1127  template <typename V> static inline auto get(V &data) {
1128  using T = typename std::remove_reference<decltype(data[0])>::type;
1130  &data[0], &data[1], &data[2], &data[3], &data[4], &data[5]};
1131  }
1132 };
1133 
1134 template <int S> struct GetFTensor1FromArray<9, S> {
1135  GetFTensor1FromArray() = delete;
1136  template <typename V> static inline auto get(V &data) {
1137  using T = typename std::remove_reference<decltype(data[0])>::type;
1139  &data[0], &data[1], &data[2], &data[3], &data[4],
1140  &data[5], &data[6], &data[7], &data[8]};
1141  }
1142 };
1143 
1144 /**
1145  * @brief Get FTensor1 from array
1146  *
1147  * \todo Generalise for different arrays and data types
1148  *
1149  * @tparam DIM
1150  * @param data
1151  * @return FTensor::Tensor1<FTensor::PackPtr<double *, S>, DIM>
1152  */
1153 template <int DIM, int S> inline auto getFTensor1FromArray(VectorDouble &data) {
1154  return GetFTensor1FromArray<DIM, S>::get(data);
1155 }
1156 
1157 /** @copydoc getFTensor1FromArray */
1158 template <int DIM, int S = 0>
1159 inline auto getFTensor1FromArray(VectorDouble3 &data);
1160 
1161 template <> inline auto getFTensor1FromArray<3, 0>(VectorDouble3 &data) {
1162  return GetFTensor1FromArray<3, 0>::get(data);
1163 }
1164 
1165 template <int DIM, int S>
1167 getFTensor1FromMat(MatrixDouble &data, const size_t rr);
1168 
1169 template <>
1171 getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1172  return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2>{&data(rr + 0, 0),
1173  &data(rr + 1, 0)};
1174 }
1175 
1176 template <>
1178 getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1180  &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1181 }
1182 
1183 template <>
1185 getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1187  &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0), &data(rr + 3, 0)};
1188 }
1189 
1190 template <>
1192 getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1194  &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0),
1195  &data(rr + 3, 0), &data(rr + 4, 0), &data(rr + 5, 0),
1196  &data(rr + 6, 0), &data(rr + 7, 0), &data(rr + 8, 0)};
1197 }
1198 
1199 /**
1200  * @brief Get FTensor1 from array
1201  *
1202  * \todo Generalise for diffrent arrays and data types
1203  *
1204  * @tparam DIM
1205  * @param data
1206  * @param rr
1207  * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
1208  */
1209 template <int DIM, int S>
1211 getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
1212  static_assert(DIM != DIM, "not implemented");
1214 }
1215 
1216 template <>
1218 getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
1219  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
1220  &data(rr + 1, 1)};
1221 }
1222 
1223 template <>
1225 getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
1227  &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1228 }
1229 
1230 /**
1231  * @brief Get FTensor2 from array
1232  *
1233  * \note Generalise for other data types
1234  *
1235  * @tparam DIM1
1236  * @tparam DIM2
1237  * @tparam S
1238  * @param data
1239  * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
1240  */
1241 template <int DIM1, int DIM2, int S, class T, class L, class A>
1243 
1244 template <int DIM1, int DIM2, class T, class L, class A>
1246 
1247 template <int S, class T, class L, class A>
1248 struct GetFTensor2FromArrayImpl<2, 2, S, T, L, A> {
1249  GetFTensor2FromArrayImpl() = delete;
1250  inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1251  const size_t cc) {
1253  &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1254 
1255  &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1256  }
1257 };
1258 
1259 template <int S, class T, class L, class A>
1260 struct GetFTensor2FromArrayImpl<3, 3, S, T, L, A> {
1261  GetFTensor2FromArrayImpl() = delete;
1262  inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1263  const size_t cc) {
1265  &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1266  &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1267  &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2)};
1268  }
1269 };
1270 
1271 template <class T, class L, class A>
1273  GetFTensor2FromArrayRawPtrImpl() = delete;
1274  inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1275  const size_t cc, const int ss = 0) {
1277  &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1278 
1279  &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
1280  }
1281 };
1282 
1283 template <class T, class L, class A>
1285  GetFTensor2FromArrayRawPtrImpl() = delete;
1286  inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1287  const size_t cc, const int ss = 0) {
1289  &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1290  &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1291  &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
1292  ss);
1293  }
1294 };
1295 
1296 template <int DIM1, int DIM2, int S>
1298 getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc = 0) {
1300  VecAllocator<double>>::get(data, rr, cc);
1301 }
1302 
1303 template <int DIM1, int DIM2>
1305 getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc,
1306  const int ss) {
1308  VecAllocator<double>>::get(data, rr, cc,
1309  ss);
1310 }
1311 
1312 template <int S, typename T, typename L, typename A>
1313 inline auto getFTensor2FromArray2by2(ublas::matrix<T, L, A> &data,
1314  const FTensor::Number<S> &,
1315  const size_t rr, const size_t cc = 0) {
1317 }
1318 
1319 template <int S, typename T, typename L, typename A>
1320 inline auto getFTensor2FromArray3by3(ublas::matrix<T, L, A> &data,
1321  const FTensor::Number<S> &,
1322  const size_t rr, const size_t cc = 0) {
1324 }
1325 
1326 #ifdef WITH_ADOL_C
1327 
1328 template <int DIM1, int DIM2, int S>
1329 inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
1331  VecAllocator<adouble>>::get(data, rr);
1332 }
1333 
1334 #endif
1335 
1336 // list of lapack wrappers
1337 /**
1338  * @brief compute matrix inverse with lapack dgetri
1339  *
1340  * @param mat input square matrix / output inverse matrix
1341  * @return MoFEMErrorCode
1342  */
1345 
1346  const size_t M = mat.size1();
1347  const size_t N = mat.size2();
1348 
1349  if (M != N)
1350  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1351  "The input matrix for inverse computation is not square %d != %d",
1352  M, N);
1353 
1354  int *ipv = new int[N];
1355  int lwork = N * N;
1356  double *work = new double[lwork];
1357  int info;
1358  info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1359  if (info != 0)
1360  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1361  "lapack error info = %d", info);
1362  info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1363  if (info != 0)
1364  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1365  "lapack error info = %d", info);
1366 
1367  delete[] ipv;
1368  delete[] work;
1369 
1371 }
1372 /**
1373  * @brief solve linear system with lapack dgesv
1374  *
1375  * @param mat input lhs square matrix / output L and U from the factorization
1376  * @param f input rhs vector / output solution vector
1377  * @return MoFEMErrorCode
1378  */
1381 
1382  const size_t M = mat.size1();
1383  const size_t N = mat.size2();
1384 
1385  if (M == 0 || M != N)
1386  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1387  "The input matrix for inverse computation is not square %d != %d",
1388  M, N);
1389  if (f.size() != M)
1390  f.resize(M, false);
1391 
1392  const int nrhs = 1;
1393  int info;
1394  int *ipiv = new int[M];
1395  info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1396  &*f.data().begin(), nrhs);
1397 
1398  if (info != 0) {
1399  SETERRQ1(PETSC_COMM_SELF, 1, "error lapack solve dgesv info = %d", info);
1400  }
1401 
1402  delete[] ipiv;
1404 }
1405 
1406 /**
1407  * @brief Solve linear system of equations using Lapack
1408  *
1409  * @param mat
1410  * @param f
1411  * @return MoFEMErrorCode
1412  */
1414  VectorDouble &f) {
1416  // copy matrix since on output lapack returns factorisation
1417  auto mat_copy = mat;
1418  CHKERR solveLinearSystem(mat_copy, f);
1420 }
1421 
1422 /**
1423  * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1424  *
1425  * @param mat input symmetric matrix
1426  * @param eig output eigen values sorted
1427  * @param eigen_vec output matrix of row eigen vectors
1428  * @return MoFEMErrorCode
1429  */
1431  VectorDouble &eig,
1432  MatrixDouble &eigen_vec) {
1434 
1435  const size_t M = mat.size1();
1436  const size_t N = mat.size2();
1437 
1438  if (M == 0 || M != N)
1439  SETERRQ2(
1440  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1441  "The input matrix for eigen value computation is not square %d != %d",
1442  M, N);
1443  if (eig.size() != M)
1444  eig.resize(M, false);
1445 
1446  eigen_vec = mat;
1447  const int n = M;
1448  const int lda = M;
1449  const int size = (M + 2) * M;
1450  int lwork = size;
1451  double *work = new double[size];
1452 
1453  if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
1454  &*eig.data().begin(), work, lwork) > 0)
1455  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1456  "The algorithm failed to compute eigenvalues.");
1457 
1458  delete[] work;
1460 }
1461 /**
1462  * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1463  *
1464  * @tparam DIM
1465  * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
1466  * @param eig output eigen values sorted
1467  * @return MoFEMErrorCode
1468  */
1469 template <int DIM, typename T1, typename T2>
1470 inline MoFEMErrorCode
1474 
1475  const int n = DIM;
1476  const int lda = DIM;
1477  const int lwork = (DIM + 2) * DIM;
1478  std::array<double, (DIM + 2) * DIM> work;
1479 
1480  if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
1481  lwork) > 0)
1482  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1483  "The algorithm failed to compute eigenvalues.");
1485 }
1486 
1487 /**
1488  * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1489  *
1490  * @tparam DIM
1491  * @param mat input tensor pointer of size DIM x DIM
1492  * @param eig output eigen values sorted
1493  * @param eigen_vec output matrix of row eigen vectors
1494  * @return MoFEMErrorCode
1495  */
1496 template <int DIM, typename T1, typename T2, typename T3>
1497 inline MoFEMErrorCode
1500  FTensor::Tensor2<T3, DIM, DIM> &eigen_vec) {
1502  for (int ii = 0; ii != DIM; ii++)
1503  for (int jj = 0; jj != DIM; jj++)
1504  eigen_vec(ii, jj) = mat(ii, jj);
1505 
1506  CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
1507 
1509 }
1510 
1511 /**
1512  * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
1513  *
1514  * @tparam T
1515  * @param t
1516  * @return double
1517  */
1518 template <typename T> static inline auto determinantTensor3by3(T &t) {
1519  return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
1520  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
1521  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
1522 }
1523 
1524 /**
1525  * @brief Calculate the determinant of a 2x2 matrix or a tensor of rank 2
1526  *
1527  * @tparam T
1528  * @param t
1529  * @return double
1530  */
1531 template <typename T> static inline auto determinantTensor2by2(T &t) {
1532  return t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
1533 }
1534 
1535 template <typename T, int DIM> struct DeterminantTensorImpl;
1536 
1537 template <typename T> struct DeterminantTensorImpl<T, 3> {
1538  static inline auto get(T &t) { return determinantTensor3by3(t); }
1539 };
1540 
1541 template <typename T> struct DeterminantTensorImpl<T, 2> {
1542  static auto get(T &t) { return determinantTensor2by2(t); }
1543 };
1544 
1545 /**
1546  * @brief Calculate the determinant of a tensor of rank DIM
1547  */
1548 template <typename T, int DIM>
1551 }
1552 
1553 /**
1554  * @brief Calculate the determinant of a tensor of rank DIM
1555  */
1556 template <typename T, int DIM>
1559 }
1560 
1561 /**
1562  * \brief Calculate inverse of tensor rank 2 at integration points
1563 
1564  */
1565 template <int Tensor_Dim, class T, class L, class A>
1566 inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
1567  ublas::vector<T, A> &det_data,
1568  ublas::matrix<T, L, A> &inv_jac_data) {
1570  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1571  "Specialization for this template not yet implemented");
1573 }
1574 
1575 template <>
1576 inline MoFEMErrorCode
1577 invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
1578  MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
1579 
1580 /**
1581  * \brief Calculate determinant 3 by 3
1582 
1583  */
1584 template <class T1, class T2>
1587  det = determinantTensor3by3(t);
1589 }
1590 
1591 /**
1592  * \brief Calculate determinant 2 by 2
1593 
1594  */
1595 template <class T1, class T2>
1598  det = determinantTensor2by2(t);
1600 }
1601 
1602 /**
1603  * \brief Calculate matrix inverse 3 by 3
1604 
1605  */
1606 template <class T1, class T2, class T3>
1607 inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
1609  const auto inv_det = 1. / det;
1610  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1611  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1612  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1613  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
1614  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1615  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1616  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
1617  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
1618  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1620 }
1621 
1622 /**
1623  * \brief Calculate matrix inverse 2 by 2
1624 
1625  */
1626 template <class T1, class T2, class T3>
1627 inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
1629  const auto inv_det = 1. / det;
1630  inv_t(0, 0) = t(1, 1) * inv_det;
1631  inv_t(0, 1) = -t(0, 1) * inv_det;
1632  inv_t(1, 0) = -t(1, 0) * inv_det;
1633  inv_t(1, 1) = t(0, 0) * inv_det;
1635 }
1636 
1637 #ifdef WITH_ADOL_C
1638 
1639 /**
1640  * \brief Calculate matrix inverse, specialization for adouble tensor
1641 
1642  */
1643 template <>
1644 inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
1649  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1650  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1651  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1652  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
1653  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1654  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1655  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
1656  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
1657  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1659 }
1660 
1661 #endif
1662 
1663 /**
1664  * \brief Calculate matrix inverse, specialization for symmetric tensor
1665 
1666  */
1667 template <>
1668 inline MoFEMErrorCode
1669 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1674  const auto inv_det = 1. / det;
1675  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1676  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1677  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1678  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1679  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1680  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1682 }
1683 
1684 #ifdef WITH_ADOL_C
1685 
1686 /**
1687  * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1688 
1689  */
1690 template <>
1691 inline MoFEMErrorCode
1692 invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
1697  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1698  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1699  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1700  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1701  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1702  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1704 }
1705 
1706 #endif
1707 
1708 /**
1709  * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1710  tensor
1711 
1712  */
1713 template <>
1714 inline MoFEMErrorCode
1715 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1720  const auto inv_det = 1. / det;
1721  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1722  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1723  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1724  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1725  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1726  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1728 }
1729 
1730 template <typename T1, typename T2, typename T3, int DIM>
1732 
1733 template <typename T1, typename T2, typename T3>
1734 struct InvertTensorImpl<T1, T2, T3, 3> {
1735  inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
1736  return invertTensor3by3(t, det, inv_t);
1737  }
1738 };
1739 
1740 template <typename T1, typename T2, typename T3>
1741 struct InvertTensorImpl<T1, T2, T3, 2> {
1742  inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
1743  return invertTensor2by2(t, det, inv_t);
1744  }
1745 };
1746 
1747 template <typename T1, typename T2, typename T3, int DIM>
1748 static inline MoFEMErrorCode
1752  FTensor::Tensor2<T3, DIM, DIM>, DIM>::invert(t, det,
1753  inv_t);
1754 }
1755 
1756 template <typename T1, typename T2, typename T3, int DIM>
1757 static inline MoFEMErrorCode
1762  DIM>::invert(t, det, inv_t);
1763 }
1764 
1765 /**
1766  * @brief Extract entity handle form multi-index container
1767  *
1768  */
1770  template <typename Iterator>
1771  static inline EntityHandle extract(const Iterator &it) {
1772  return (*it)->getEnt();
1773  }
1774 };
1775 
1776 /**
1777  * @brief Insert ordered mofem multi-index into range
1778  *
1779  * \note Inserted range has to be ordered.
1780  *
1781  * \code
1782  * auto hi_rit = refEntsPtr->upper_bound(start);
1783  * auto hi_rit = refEntsPtr->upper_bound(end);
1784  * Range to_erase;
1785  * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1786  * \endcode
1787  *
1788  * @tparam Iterator
1789  * @param r
1790  * @param begin_iter
1791  * @param end_iter
1792  * @return moab::Range::iterator
1793  */
1794 template <typename Extractor, typename Iterator>
1795 moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1796  Iterator end_iter) {
1797  moab::Range::iterator hint = r.begin();
1798  while (begin_iter != end_iter) {
1799  size_t j = 0;
1800  auto bi = Extractor::extract(begin_iter);
1801  Iterator pj = begin_iter;
1802  while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1803  ++pj;
1804  ++j;
1805  }
1806  hint = r.insert(hint, bi, bi + (j - 1));
1807  begin_iter = pj;
1808  }
1809  return hint;
1810 };
1811 
1812 /**
1813  * @brief Do nothing, used to rebuild database
1814  *
1815  */
1817  Modify_change_nothing() = default;
1818  template <typename T> inline void operator()(T &e) {}
1819 };
1820 
1821 /**
1822  * @brief Template used to reconstruct multi-index
1823  *
1824  * @tparam MI multi-index
1825  * @tparam Modifier
1826  * @param mi
1827  * @param mo
1828  * @return MoFEMErrorCode
1829  */
1830 template <typename MI, typename MO = Modify_change_nothing>
1832  MO &&mo = Modify_change_nothing()) {
1834  for (auto it = mi.begin(); it != mi.end(); ++it) {
1835  if (!const_cast<MI &>(mi).modify(it, mo))
1836  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1837  "Houston we have a problem");
1838  }
1840 }
1841 
1842 struct TempMeshset {
1844  rval = moab.create_meshset(MESHSET_SET, meshset);
1845  MOAB_THROW(rval);
1846  }
1847  virtual ~TempMeshset() { delete_meshset(); }
1848  operator EntityHandle() const { return meshset; }
1849  auto get_ptr() { return &meshset; }
1850 
1851 private:
1853  rval = moab.delete_entities(&meshset, 1);
1854  MOAB_THROW(rval);
1855  }
1858 };
1859 
1860 /**
1861  * @brief Create smart pointer to temporary meshset
1862  *
1863  */
1865  return boost::make_shared<TempMeshset>(moab);
1866 };
1867 
1868 inline auto id_from_handle(const EntityHandle h) {
1869  return static_cast<EntityID>(h & MB_ID_MASK);
1870 };
1871 
1872 /**
1873  * @brief get type from entity handle
1874  *
1875  */
1876 inline auto type_from_handle(const EntityHandle h) {
1877  return static_cast<EntityType>(h >> MB_ID_WIDTH);
1878 };
1879 
1880 /**
1881  * @brief get entity handle from type and id
1882  *
1883  */
1884 inline auto ent_form_type_and_id(const EntityType type, const EntityID id) {
1885  return (static_cast<EntityHandle>(type) << MB_ID_WIDTH) | id;
1886 };
1887 
1888 /**
1889  * @brief get entity dimension form handle
1890  *
1891  */
1893  return moab::CN::Dimension(type_from_handle(h));
1894 };
1895 
1896 /**
1897  * @brief get entity type name from handle
1898  *
1899  */
1901  return moab::CN::EntityTypeName(type_from_handle(h));
1902 };
1903 
1904 /**
1905  * @brief get field bit id from bit number
1906  *
1907  */
1908 inline auto field_bit_from_bit_number(const int bit_number) {
1909  return BitFieldId().set(bit_number - 1);
1910 };
1911 
1912 /**
1913  * @brief Insert ranges
1914  *
1915  * @tparam I
1916  * @param f
1917  * @param s
1918  * @param tester
1919  * @param inserter
1920  * @return auto
1921  */
1922 template <typename I>
1923 auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
1924  boost::function<MoFEMErrorCode(I f, I s)> inserter) {
1926 
1927  auto first = f;
1928  while (first != s)
1929  if (tester(first)) {
1930 
1931  auto second = first;
1932  ++second;
1933 
1934  while (second != s) {
1935  if (tester(second))
1936  ++second;
1937  else
1938  break;
1939  }
1940 
1941  CHKERR inserter(first, second);
1942 
1943  first = second;
1944  if (first != s)
1945  ++first;
1946 
1947  } else {
1948  ++first;
1949  }
1950 
1952 }
1953 
1954 /**
1955  * @brief Create Array
1956  *
1957  * See:
1958  * <a
1959  * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
1960  * stack overflow</a>
1961  *
1962  * @tparam Dest
1963  * @tparam Arg
1964  * @param arg
1965  * @return constexpr auto
1966  */
1967 template <typename Dest = void, typename... Arg>
1968 constexpr auto make_array(Arg &&...arg) {
1969  if constexpr (std::is_same<void, Dest>::value)
1970  return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
1971  {std::forward<Arg>(arg)...}};
1972  else
1973  return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
1974 }
1975 
1976 } // namespace MoFEM
1977 
1978 #endif //__TEMPLATES_HPP__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::GetFTensor2FromPtr< 3, 2, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:831
MoFEM::GetFTensor2FromMatImpl
Definition: Templates.hpp:252
MoFEM::id_from_handle
auto id_from_handle(const EntityHandle h)
Definition: Templates.hpp:1868
MoFEM::dimension_from_handle
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
Definition: Templates.hpp:1892
MoFEM::KeyFromKey::operator()
result_type operator()(Arg &arg) const
Definition: Templates.hpp:79
MoFEM::getFTensor1FromPtr
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition: Templates.hpp:809
MoFEM::TempMeshset::~TempMeshset
virtual ~TempMeshset()
Definition: Templates.hpp:1847
MoFEM::GetFTensor3FromMatImpl< 2, 2, 2, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:632
MoFEM::GetFTensor1FromArray< 4, S >::get
static auto get(V &data)
Definition: Templates.hpp:1118
MoFEM::GetFTensor2FromPtr< 1, 1, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:883
MoFEM::LtBit::operator()
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:89
MoFEM::GetFTensor1FromArray< 9, S >::get
static auto get(V &data)
Definition: Templates.hpp:1136
MoFEM::GetFTensor2FromArrayRawPtrImpl< 3, 3, T, L, A >::get
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
Definition: Templates.hpp:1286
MoFEM::GetFTensor3DgFromMatImpl< 2, 2, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:448
MoFEM::GetFTensor1FromMatImpl< 9, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:190
MoFEM::field_bit_from_bit_number
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
Definition: Templates.hpp:1908
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MoFEM::GetFTensor4DdgFromMatImpl
Definition: Templates.hpp:345
EntityHandle
HVEC0_2
@ HVEC0_2
Definition: definitions.h:198
MoFEM::getFTensor2HVecFromPtr
FTensor::Tensor2< FTensor::PackPtr< double *, DIM1 *DIM2 >, DIM1, DIM2 > getFTensor2HVecFromPtr(double *ptr)
Make Tensor2 for HVec base from pointer.
Definition: Templates.hpp:921
lapack_dgetri
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:185
MoFEM::Types::MatrixShallowArrayAdaptor
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition: Types.hpp:121
FTensor::row_major
@ row_major
Definition: Layout.hpp:13
MoFEM::GetFTensor3FromMatImpl< 3, 3, 3, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:684
MoFEM::getFTensor2SymmetricFromMat
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:334
MoFEM::GetFTensor3DgFromMatImpl
Definition: Templates.hpp:430
MoFEM::GetFTensor1FromMatImpl< 4, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:159
lapack_dgesv
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:176
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
HVEC0_1
@ HVEC0_1
Definition: definitions.h:195
MoFEM::EqBit::operator()
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:95
MoFEM::Types::BitFieldId
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:42
MoFEM::GetFTensor2FromPtr< 6, 6, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:839
MoFEM::getFTensor2FromArray
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1298
MoFEM::GetFTensor1FromArray< 3, S >::get
static auto get(V &data)
Definition: Templates.hpp:1109
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::TempMeshset::meshset
EntityHandle meshset
Definition: Templates.hpp:1856
MoFEM::GetFTensor4DdgFromMatImpl< 3, 3, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:380
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
MoFEM::solveLinearSystem
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
Definition: Templates.hpp:1379
MoFEM::Modify_change_nothing
Do nothing, used to rebuild database.
Definition: Templates.hpp:1816
MoFEM::GetFTensor2FromPtr< 2, 2, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:860
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
MoFEM::getFTensor1FromMat
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:238
MoFEM::Modify_change_nothing::Modify_change_nothing
Modify_change_nothing()=default
HVEC0_0
@ HVEC0_0
Definition: definitions.h:192
MoFEM::getFTensor2SymmetricLowerFromPtr< 3 >
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 9 >, 3 > getFTensor2SymmetricLowerFromPtr< 3 >(double *ptr)
Definition: Templates.hpp:1073
HVEC1_1
@ HVEC1_1
Definition: definitions.h:196
MoFEM::getFTensor2SymmetricFromPtr< 3 >
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricFromPtr< 3 >(double *ptr)
Definition: Templates.hpp:1005
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
MoFEM::TempMeshset::get_ptr
auto get_ptr()
Definition: Templates.hpp:1849
HVEC1
@ HVEC1
Definition: definitions.h:186
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::GetFTensor2FromPtr< 1, 2, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:876
MoFEM::reconstructMultiIndex
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1831
MoFEM::getFTensor2FromPtr
auto getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
Definition: Templates.hpp:896
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::invertTensor3by3
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:1566
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::GetFTensor2FromMatImpl::get
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:253
MoFEM::DeterminantTensorImpl< T, 3 >::get
static auto get(T &t)
Definition: Templates.hpp:1538
FTensor::Tensor1::data
T data[Tensor_Dim]
Definition: Tensor1_value.hpp:10
MoFEM::ShardVec
boost::shared_ptr< std::vector< T > > ShardVec
Definition: Templates.hpp:10
MoFEM::determinantTensor2by2
static auto determinantTensor2by2(T &t)
Calculate the determinant of a 2x2 matrix or a tensor of rank 2.
Definition: Templates.hpp:1531
MoFEM::GetFTensor2FromArrayRawPtrImpl< 2, 2, T, L, A >::get
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
Definition: Templates.hpp:1274
HVEC2_1
@ HVEC2_1
Definition: definitions.h:197
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PlasticOps::M
FTensor::Index< 'M', 3 > M
Definition: PlasticOps.hpp:117
MoFEM::GetFTensor1FromMatImpl
Definition: Templates.hpp:140
MoFEM::Types::MatrixADouble
UBlasMatrix< adouble > MatrixADouble
Definition: Types.hpp:80
MoFEM::GetFTensor2FromArrayImpl< 3, 3, S, T, L, A >::get
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
Definition: Templates.hpp:1262
MoFEM::GetFTensor3FromMatImpl< 3, 3, 6, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:721
MoFEM::GetFTensor2SymmetricFromMatImpl
Definition: Templates.hpp:297
FTensor::Number
Definition: Number.hpp:11
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
MoFEM::invertTensor2by2
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1627
HVEC2_2
@ HVEC2_2
Definition: definitions.h:200
HVEC1_2
@ HVEC1_2
Definition: definitions.h:199
MoFEM::GetFTensor1FromArray< 6, S >::get
static auto get(V &data)
Definition: Templates.hpp:1127
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::GetFTensor1FromPtrImpl< 3, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:778
MoFEM::Modify_change_nothing::operator()
void operator()(T &e)
Definition: Templates.hpp:1818
MoFEM::GetFTensor0FromVecImpl
Definition: Templates.hpp:112
MoFEM::getFTensor3FromMat
static FTensor::Tensor3< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 3 (non symmetries) form data matrix.
Definition: Templates.hpp:754
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:939
MoFEM::insertOrdered
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1795
MoFEM::getFTensor1FromArray< 3, 0 >
auto getFTensor1FromArray< 3, 0 >(VectorDouble3 &data)
Definition: Templates.hpp:1161
MoFEM::KeyFromKey
Definition: Templates.hpp:71
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1864
MoFEM::KeyFromKey::key2
KeyExtractor2 key2
Definition: Templates.hpp:85
DIM1
constexpr int DIM1
Definition: level_set.cpp:21
double
MoFEM::GetFTensor3FromMatImpl
Definition: Templates.hpp:613
MoFEM::GetFTensor3DgFromMatImpl< 3, 3, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:465
convert.type
type
Definition: convert.py:64
MoFEM::GetFTensor4FromMatImpl< 3, 3, 3, 3, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:549
MoFEM::HashBit::operator()
unsigned int operator()(const id_type &value) const
Definition: Templates.hpp:101
MoFEM::KeyFromKey::key1
KeyExtractor1 key1
Definition: Templates.hpp:84
MoFEM::TempMeshset::moab
moab::Interface & moab
Definition: Templates.hpp:1857
MoFEM::DeterminantTensorImpl< T, 2 >::get
static auto get(T &t)
Definition: Templates.hpp:1542
MoFEM::GetFTensor1FromMatImpl< 1, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:221
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::getFTensor1FromArray
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
Definition: Templates.hpp:1153
MoFEM::getFTensor2SymmetricFromPtr< 2 >
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 3 >, 2 > getFTensor2SymmetricFromPtr< 2 >(double *ptr)
Definition: Templates.hpp:1016
MoFEM::getFTensor2FromArray2by2
auto getFTensor2FromArray2by2(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1313
MoFEM::invertTensor
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
Definition: Templates.hpp:1749
h
double h
Definition: photon_diffusion.cpp:60
MoFEM::GetFTensor3FromMatImpl< 1, 1, 1, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:617
MoFEM::getFTensor2FromVec
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromVec(VectorDouble &data)
Definition: Templates.hpp:292
MoFEM::determinantTensor
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
Definition: Templates.hpp:1549
MoFEM::GetFTensor1FromPtrImpl
Definition: Templates.hpp:766
MoFEM::KeyFromKey::result_type
KeyExtractor1::result_type result_type
Definition: Templates.hpp:73
MoFEM::getMatrixAdaptor
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:57
MoFEM::getFTensor1FromArrayDiag
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr)
Get FTensor1 from array.
Definition: Templates.hpp:1211
MoFEM::GetFTensor3FromMatImpl< 3, 2, 2, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:650
MoFEM::getFTensor2SymmetricFromPtr
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(DIM *(DIM+1))/2 >, DIM > getFTensor2SymmetricFromPtr(double *ptr)
Make symmetric Tensor2 from pointer.
Definition: Templates.hpp:999
lapack_dgetrf
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:157
FTensor::PackPtr
Definition: FTensor.hpp:54
MoFEM::type_name_from_handle
auto type_name_from_handle(const EntityHandle h)
get entity type name from handle
Definition: Templates.hpp:1900
MoFEM::RefEntExtractor
Extract entity handle form multi-index container.
Definition: Templates.hpp:1769
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::getFTensor3FromPtr< 3, 2, 2 >
FTensor::Tensor3< FTensor::PackPtr< double *, 12 >, 3, 2, 2 > getFTensor3FromPtr< 3, 2, 2 >(double *ptr)
Definition: Templates.hpp:966
MoFEM::GetFTensor2FromPtr< 1, 3, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:868
MoFEM::GetFTensor4DdgFromMatImpl< 1, 1, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:349
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1876
MB_ID_WIDTH
#define MB_ID_WIDTH
Definition: definitions.h:227
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
MoFEM::GetFTensor0FromVecImpl::get
static auto get(ublas::vector< T, A > &data)
Definition: Templates.hpp:113
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::DeterminantTensorImpl
Definition: Templates.hpp:1535
MoFEM::Types::DoubleAllocator
VecAllocator< double > DoubleAllocator
Definition: Types.hpp:62
MoFEM::Types::VecAllocator
std::vector< T, std::allocator< T > > VecAllocator
Definition: Types.hpp:59
MoFEM::GetFTensor1FromArray< 2, S >::get
static auto get(V &data)
Definition: Templates.hpp:1101
MoFEM::HashBit
Definition: Templates.hpp:100
convert.n
n
Definition: convert.py:82
MoFEM::InvertTensorImpl< T1, T2, T3, 3 >::invert
static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t)
Definition: Templates.hpp:1735
MoFEM::GetFTensor2SymmetricFromMatImpl< 2, S, T, L, A >::get
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:317
MoFEM::getFTensor3FromPtr
FTensor::Tensor3< FTensor::PackPtr< double *, DIM1 *DIM2 *DIM3 >, DIM1, DIM2, DIM3 > getFTensor3FromPtr(double *ptr)
Definition: Templates.hpp:959
N
const int N
Definition: speed_test.cpp:3
MoFEM::getFTensor2FromArray3by3
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1320
MoFEM::InvertTensorImpl< T1, T2, T3, 2 >::invert
static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t)
Definition: Templates.hpp:1742
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1518
Range
MoFEM::GetFTensor3DgFromMatImpl< 1, 1, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:434
adouble
MoFEM::GetFTensor3FromMatImpl< 2, 2, 3, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:667
MoFEM::GetFTensor4FromMatImpl
Definition: Templates.hpp:511
FTensor::Dg
Definition: Dg_value.hpp:9
FTensor::Tensor0
Definition: Tensor0.hpp:16
MoFEM::GetFTensor3FromMatImpl< 6, 3, 3, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:704
MoFEM::InvertTensorImpl
Definition: Templates.hpp:1731
MoFEM::toString
std::string toString(X x)
Definition: Templates.hpp:106
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::getFTensor2SymmetricLowerFromPtr
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, DIM *DIM >, DIM > getFTensor2SymmetricLowerFromPtr(double *ptr)
Make symmetric Tensor2 from pointer, taking lower triangle of matrix.
Definition: Templates.hpp:1066
MoFEM::rangeInserter
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
Definition: Templates.hpp:1923
MoFEM::GetFTensor2FromPtr
Definition: Templates.hpp:827
MoFEM::GetFTensor1FromPtrImpl< 2, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:770
MoFEM::getFTensor2HVecFromPtr< 3, 2 >
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:928
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::GetFTensor2FromArrayImpl
Get FTensor2 from array.
Definition: Templates.hpp:1242
FTensor::Ddg
Definition: Ddg_value.hpp:7
lapack_dsyev
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:261
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::GetFTensor1FromPtrImpl< 6, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:794
MoFEM::GetFTensor4DdgFromMatImpl< 2, 2, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:363
MoFEM::computeMatrixInverse
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
Definition: Templates.hpp:1343
MoFEM::TempMeshset::delete_meshset
void delete_meshset()
Definition: Templates.hpp:1852
MoFEM::GetFTensor2FromArrayRawPtrImpl
Definition: Templates.hpp:1245
DIM2
constexpr int DIM2
Definition: level_set.cpp:22
HVEC2_0
@ HVEC2_0
Definition: definitions.h:194
MoFEM::getFTensor4DdgFromMat
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:417
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::TempMeshset
Definition: Templates.hpp:1842
MoFEM::TempMeshset::TempMeshset
TempMeshset(moab::Interface &moab)
Definition: Templates.hpp:1843
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::Types::VectorShallowArrayAdaptor
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition: Types.hpp:114
MoFEM::GetFTensor1FromMatImpl< 6, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:174
MoFEM::make_array
constexpr auto make_array(Arg &&...arg)
Create Array.
Definition: Templates.hpp:1968
MoFEM::GetFTensor1FromPtrImpl< 4, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:786
MoFEM::getFTensor2FromMat
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:278
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::GetFTensor1FromMatImpl< 2, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:206
MoFEM::getFTensor2SymmetricLowerFromPtr< 2 >
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 4 >, 2 > getFTensor2SymmetricLowerFromPtr< 2 >(double *ptr)
Definition: Templates.hpp:1084
MoFEM::LtBit
Definition: Templates.hpp:88
HVEC2
@ HVEC2
Definition: definitions.h:186
MoFEM::GetFTensor2FromPtr< 3, 3, S, T >::get
static auto get(T *ptr)
Definition: Templates.hpp:851
MoFEM::getFTensor3DgFromMat
static FTensor::Dg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim2 > getFTensor3DgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 3 on the first two indices from form data matrix.
Definition: Templates.hpp:498
MoFEM::ent_form_type_and_id
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
Definition: Templates.hpp:1884
MB_ID_MASK
#define MB_ID_MASK
Definition: definitions.h:234
MoFEM::RefEntExtractor::extract
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:1771
MoFEM::GetFTensor1FromArray
Definition: Templates.hpp:1089
MoFEM::GetFTensor1FromArray< 1, S >::get
static auto get(V &data)
Definition: Templates.hpp:1093
MoFEM::getFTensor4FromMat
static FTensor::Tensor4< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3 > getFTensor4FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 4 (non symmetric) form data matrix.
Definition: Templates.hpp:598
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::GetFTensor1FromMatImpl< 3, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:144
MoFEM::GetFTensor4FromMatImpl< 2, 2, 2, 2, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:530
MoFEM::GetFTensor4FromMatImpl< 1, 1, 1, 1, S, T, ublas::row_major, A >::get
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:515
MoFEM::EqBit
Definition: Templates.hpp:94
MoFEM::getFTensor3FromPtr< 3, 3, 3 >
FTensor::Tensor3< FTensor::PackPtr< double *, 27 >, 3, 3, 3 > getFTensor3FromPtr< 3, 3, 3 >(double *ptr)
Definition: Templates.hpp:981
MoFEM::KeyFromKey::KeyFromKey
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition: Templates.hpp:75
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1430
MoFEM::GetFTensor2FromArrayImpl< 2, 2, S, T, L, A >::get
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
Definition: Templates.hpp:1250
MoFEM::GetFTensor2SymmetricFromMatImpl< 3, S, T, L, A >::get
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:301
HVEC1_0
@ HVEC1_0
Definition: definitions.h:193
HVEC0
@ HVEC0
Definition: definitions.h:186