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>
159 struct GetFTensor1FromMatImpl<3, S, double, ublas::row_major, DoubleAllocator> {
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>
174 struct GetFTensor1FromMatImpl<2, S, double, ublas::row_major, DoubleAllocator> {
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>
329 struct GetFTensor4DdgFromMatImpl<1, 1, S, double, ublas::row_major,
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>
344 struct GetFTensor4DdgFromMatImpl<2, 2, S, double, ublas::row_major,
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 9 but is " +
352  boost::lexical_cast<std::string>(data.size1()));
353  }
355  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
356  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
357  }
358 };
359 
360 template <int S>
361 struct GetFTensor4DdgFromMatImpl<3, 3, S, double, ublas::row_major,
362  DoubleAllocator> {
363  static inline FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3>
364  get(MatrixDouble &data) {
365  if (data.size1() != 36) {
366  cerr << data.size1() << endl;
368  "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
369  "of rows should be 36 but is " +
370  boost::lexical_cast<std::string>(data.size1()));
371  }
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), &data(9, 0),
375  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
376  &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
377  &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
378  &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
379  &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
380  &data(35, 0)};
381  }
382 };
383 
384 /**
385  * @brief Get symmetric tensor rank 4 on first two and last indices from
386  * form data matrix
387  *
388  * @tparam Tensor_Dim01 dimension of frirst two indicies
389  * @tparam Tensor_Dim23 dimension of second two indicies
390  * @tparam T the type of object stored
391  * @tparam L the storage organization
392  * @tparam A the type of Storage array
393  * @param data data container
394  * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
395  */
396 template <int Tensor_Dim01, int Tensor_Dim23, int S = 1, class T, class L,
397  class A>
398 static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
399 getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
400  static_assert(!std::is_same<T, T>::value,
401  "Such getFTensor4DdgFromMat specialisation is not implemented");
402 }
403 
404 template <int Tensor_Dim01, int Tensor_Dim23, int S = 1>
405 static inline FTensor::Ddg<FTensor::PackPtr<double *, S>, Tensor_Dim01,
406  Tensor_Dim23>
408  return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, double,
410  DoubleAllocator>::get(data);
411 }
412 
413 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
414  int S, class T, class L, class A>
416 
417 template <int S>
418 struct GetFTensor4FromMatImpl<1, 1, 1, 1, S, double, ublas::row_major,
419  DoubleAllocator> {
420  static inline FTensor::Tensor4<FTensor::PackPtr<double *, S>, 1, 1, 1, 1>
421  get(MatrixDouble &data) {
422  if (data.size1() != 1)
424  "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
425  "of rows should be 1 but is " +
426  boost::lexical_cast<std::string>(data.size1()));
427 
429  &data(0, 0)};
430  }
431 };
432 
433 template <int S>
434 struct GetFTensor4FromMatImpl<2, 2, 2, 2, S, double, ublas::row_major,
435  DoubleAllocator> {
436  static inline FTensor::Tensor4<FTensor::PackPtr<double *, S>, 2, 2, 2, 2>
437  get(MatrixDouble &data) {
438  if (data.size1() != 16) {
440  "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
441  "of rows should be 16 but is " +
442  boost::lexical_cast<std::string>(data.size1()));
443  }
445  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
446  &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
447  &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
448  &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
449  }
450 };
451 
452 template <int S>
453 struct GetFTensor4FromMatImpl<3, 3, 3, 3, S, double, ublas::row_major,
454  DoubleAllocator> {
455  static inline FTensor::Tensor4<FTensor::PackPtr<double *, S>, 3, 3, 3, 3>
456  get(MatrixDouble &data) {
457  if (data.size1() != 81) {
458  cerr << data.size1() << endl;
460  "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
461  "of rows should be 81 but is " +
462  boost::lexical_cast<std::string>(data.size1()));
463  }
465  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
466  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
467  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
468  &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
469  &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
470  &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
471  &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
472  &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
473  &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
474  &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
475  &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
476  &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
477  &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
478  &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
479  &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
480  &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
481  &data(80, 0)};
482  }
483 };
484 
485 /**
486  * @brief Get tensor rank 4 (non symmetric) form data matrix
487  *
488  * @tparam Tensor_Dim0 dimension of frirst index
489  * @tparam Tensor_Dim1 dimension of second index
490  * @tparam Tensor_Dim2 dimension of third index
491  * @tparam Tensor_Dim3 dimension of fourth index
492  * @tparam T the type of object stored
493  * @tparam L the storage organization
494  * @tparam A the type of Storage array
495  * @param data data container
496  * @return FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
497  Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
498  */
499 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
500  int S = 1, class T, class L, class A>
501 static inline FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
502  Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
503 getFTensor4FromMat(ublas::matrix<T, L, A> &data) {
504  static_assert(!std::is_same<T, T>::value,
505  "Such getFTensor4FromMat specialisation is not implemented");
506 }
507 
508 template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
509  int S = 1>
510 static inline FTensor::Tensor4<FTensor::PackPtr<double *, S>, Tensor_Dim0,
511  Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
513  return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
514  Tensor_Dim3, S, double, ublas::row_major,
515  DoubleAllocator>::get(data);
516 }
517 
518 /**
519  * @brief Make Tensor1 from pointer
520  *
521  * @tparam DIM
522  * @param ptr
523  * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
524  */
525 template <int DIM>
527 getFTensor1FromPtr(double *ptr) {
528  static_assert(DIM != 3 && DIM != 2,
529  "Such getFTensor1FromPtr specialization is not implemented");
530 };
531 
532 template <>
534 getFTensor1FromPtr<2>(double *ptr) {
536  &ptr[HVEC1]);
537 };
538 
539 template <>
541 getFTensor1FromPtr<3>(double *ptr) {
543  &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
544 };
545 
546 /**
547  * @brief Make Tensor2 from pointer
548  *
549  * @tparam DIM
550  * @param ptr
551  * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
552  */
553 template <int DIM1, int DIM2>
555 getFTensor2FromPtr(double *ptr) {
556  static_assert(DIM1 != 3, "Such getFTensor2FromPtr is not implemented");
557  static_assert(DIM2 >= 2 && DIM2 <= 3,
558  "Such getFTensor2FromPtr is not implemented");
559 };
560 
561 template <>
563  2> inline getFTensor2FromPtr<3, 2>(double *ptr) {
565  &ptr[HVEC0_0], &ptr[HVEC0_1],
566 
567  &ptr[HVEC1_0], &ptr[HVEC1_1],
568 
569  &ptr[HVEC2_0], &ptr[HVEC2_1]);
570 };
571 
572 template <>
574  3> inline getFTensor2FromPtr<3, 3>(double *ptr) {
576  &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
577 
578  &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
579 
580  &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
581 };
582 
583 /**
584  * @brief Get FTensor1 from array
585  *
586  * \todo Generalise for diffrent arrays and data types
587  *
588  * @tparam DIM
589  * @param data
590  * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
591  */
592 template <int DIM, int S>
595  static_assert(DIM != DIM, "not implemented");
597 }
598 
599 template <>
602  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data[0], &data[1]};
603 }
604 
605 template <>
608  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>{&data[0], &data[1],
609  &data[2]};
610 }
611 
612 /**
613  * @brief Get FTensor1 from array
614  *
615  * \todo Generalise for diffrent arrays and data types
616  *
617  * @tparam DIM
618  * @param data
619  * @param rr
620  * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
621  */
622 template <int DIM, int S>
624 getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
625  static_assert(DIM != DIM, "not implemented");
627 }
628 
629 template <>
631 getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
632  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
633  &data(rr + 1, 1)};
634 }
635 
636 template <>
638 getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
640  &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
641 }
642 
643 /**
644  * @brief Get FTensor2 from array
645  *
646  * \note Generalise for other data types
647  *
648  * @tparam DIM1
649  * @tparam DIM2
650  * @tparam S
651  * @param data
652  * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
653  */
654 template <int DIM1, int DIM2, int S>
656 getFTensor2FromArray(MatrixDouble &data, const size_t rr) {
657  static_assert(DIM1 != DIM1, "not implemented");
659 }
660 
661 template <>
663 getFTensor2FromArray(MatrixDouble &data, const size_t rr) {
665  &data(rr + 0, 0), &data(rr + 0, 1), &data(rr + 1, 0), &data(rr + 1, 1)};
666 }
667 
668 template <>
670 getFTensor2FromArray(MatrixDouble &data, const size_t rr) {
672  &data(rr + 0, 0), &data(rr + 0, 1), &data(rr + 0, 2),
673  &data(rr + 1, 0), &data(rr + 1, 1), &data(rr + 1, 2),
674  &data(rr + 2, 0), &data(rr + 2, 1), &data(rr + 2, 2)};
675 }
676 
677 // list of lapack wrappers
678 /**
679  * @brief compute matrix inverse with lapack dgetri
680  *
681  * @param mat input square matrix / output inverse matrix
682  * @return MoFEMErrorCode
683  */
686 
687  const size_t M = mat.size1();
688  const size_t N = mat.size2();
689 
690  if (M != N)
691  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
692  "The input matrix for inverse computation is not square %d != %d",
693  M, N);
694 
695  int *ipv = new int[N];
696  int lwork = N * N;
697  double *work = new double[lwork];
698  int info;
699  info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
700  if (info != 0)
701  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
702  "lapack error info = %d", info);
703  info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
704  if (info != 0)
705  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
706  "lapack error info = %d", info);
707 
708  delete[] ipv;
709  delete[] work;
710 
712 }
713 /**
714  * @brief solve linear system with lapack dgesv
715  *
716  * @param mat input lhs square matrix / output L and U from the factorization
717  * @param f input rhs vector / output solution vector
718  * @return MoFEMErrorCode
719  */
722 
723  const size_t M = mat.size1();
724  const size_t N = mat.size2();
725 
726  if (M == 0 || M != N)
727  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
728  "The input matrix for inverse computation is not square %d != %d",
729  M, N);
730  if (f.size() != M)
731  f.resize(M, false);
732 
733  const int nrhs = 1;
734  int info;
735  int *ipiv = new int[M];
736  info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
737  &*f.data().begin(), nrhs);
738 
739  if (info != 0) {
740  SETERRQ1(PETSC_COMM_SELF, 1, "error lapack solve dgesv info = %d", info);
741  }
742 
743  delete[] ipiv;
745 }
746 
748  VectorDouble &f) {
750  // copy matrix since on output lapack returns factorisation
751  auto mat_copy = mat;
752  CHKERR solveLinearSystem(mat_copy, f);
754 }
755 
756 /**
757  * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
758  *
759  * @param mat input symmetric matrix
760  * @param eig output eigen values sorted
761  * @param eigen_vec output matrix of row eigen vectors
762  * @return MoFEMErrorCode
763  */
765  VectorDouble &eig,
766  MatrixDouble &eigen_vec) {
768 
769  const size_t M = mat.size1();
770  const size_t N = mat.size2();
771 
772  if (M == 0 || M != N)
773  SETERRQ2(
774  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
775  "The input matrix for eigen value computation is not square %d != %d",
776  M, N);
777  if (eig.size() != M)
778  eig.resize(M, false);
779 
780  eigen_vec = mat;
781  const int n = M;
782  const int lda = M;
783  const int size = (M + 2) * M;
784  int lwork = size;
785  double *work = new double[size];
786 
787  if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
788  &*eig.data().begin(), work, lwork) > 0)
789  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
790  "The algorithm failed to compute eigenvalues.");
791 
792  delete[] work;
794 }
795 /**
796  * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
797  *
798  * @tparam DIM
799  * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
800  * @param eig output eigen values sorted
801  * @return MoFEMErrorCode
802  */
803 template <int DIM>
804 inline MoFEMErrorCode
808 
809  const int n = DIM;
810  const int lda = DIM;
811  const int lwork = (DIM + 2) * DIM;
812  std::array<double, (DIM + 2) * DIM> work;
813 
814  if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
815  lwork) > 0)
816  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
817  "The algorithm failed to compute eigenvalues.");
819 }
820 /**
821  * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
822  *
823  * @tparam DIM
824  * @param mat input tensor pointer of size DIM x DIM
825  * @param eig output eigen values sorted
826  * @param eigen_vec output matrix of row eigen vectors
827  * @return MoFEMErrorCode
828  */
829 template <int DIM>
835  for (int ii = 0; ii != DIM; ii++)
836  for (int jj = 0; jj != DIM; jj++)
837  eigen_vec(ii, jj) = mat(ii, jj);
838 
839  CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
840 
842 }
843 
844 /**
845  * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
846  *
847  * @tparam DIM
848  * @param mat input tensor of size DIM x DIM
849  * @param eig output eigen values sorted
850  * @param eigen_vec output matrix of row eigen vectors
851  * @return MoFEMErrorCode
852  */
853 template <int DIM>
854 inline MoFEMErrorCode
859  for (int ii = 0; ii != DIM; ii++)
860  for (int jj = 0; jj != DIM; jj++)
861  eigen_vec(ii, jj) = mat(ii, jj);
862 
863  CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
864 
866 }
867 
868 /**
869  * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
870  *
871  * @tparam T
872  * @param t
873  * @return double
874  */
875 template <class T> static inline double dEterminant(T &t) {
876  return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
877  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
878  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
879 }
880 
881 /**
882  * \brief Calculate inverse of tensor rank 2 at integration points
883 
884  */
885 template <int Tensor_Dim, class T, class L, class A>
886 inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
887  ublas::vector<T, A> &det_data,
888  ublas::matrix<T, L, A> &inv_jac_data) {
890  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
891  "Specialization for this template not yet implemented");
893 }
894 
895 template <>
896 inline MoFEMErrorCode
897 invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
898  MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
899 
900 /**
901  * \brief Calculate determinant 3 by 3
902 
903  */
904 template <class T1, class T2>
905 inline MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det) {
907  det = +t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
908  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
909  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
911 }
912 
913 /**
914  * \brief Calculate determinant 2 by 2
915 
916  */
917 template <class T1, class T2>
918 inline MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det) {
920  det = t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
922 }
923 
924 /**
925  * \brief Calculate matrix inverse 3 by 3
926 
927  */
928 template <class T1, class T2, class T3>
929 inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
931  const auto inv_det = 1. / det;
932  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
933  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
934  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
935  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
936  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
937  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
938  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
939  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
940  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
942 }
943 
944 /**
945  * \brief Calculate matrix inverse 2 by 2
946 
947  */
948 template <class T1, class T2, class T3>
949 inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
951  const auto inv_det = 1. / det;
952  inv_t(0, 0) = t(1, 1) * inv_det;
953  inv_t(0, 1) = -t(0, 1) * inv_det;
954  inv_t(1, 0) = -t(1, 0) * inv_det;
955  inv_t(1, 1) = t(0, 0) * inv_det;
957 }
958 
959 #ifdef WITH_ADOL_C
960 
961 /**
962  * \brief Calculate matrix inverse, specialization for adouble tensor
963 
964  */
965 template <>
966 inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
968  FTensor::Tensor2<adouble, 3, 3> &t, adouble &det,
971  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
972  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
973  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
974  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
975  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
976  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
977  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
978  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
979  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
981 }
982 
983 #endif
984 
985 /**
986  * \brief Calculate matrix inverse, specialization for symmetric tensor
987 
988  */
989 template <>
990 inline MoFEMErrorCode
991 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
996  const auto inv_det = 1. / det;
997  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
998  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
999  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1000  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1001  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1002  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1004 }
1005 
1006 #ifdef WITH_ADOL_C
1007 
1008 /**
1009  * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1010 
1011  */
1012 template <>
1013 inline MoFEMErrorCode
1014 invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
1016  FTensor::Tensor2_symmetric<adouble, 3> &t, adouble &det,
1019  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1020  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1021  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1022  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1023  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1024  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1026 }
1027 
1028 #endif
1029 
1030 /**
1031  * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1032  tensor
1033 
1034  */
1035 template <>
1036 inline MoFEMErrorCode
1037 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1039  FTensor::Tensor2_symmetric<double, 3> &t, double &det,
1042  const auto inv_det = 1. / det;
1043  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1044  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1045  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1046  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1047  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1048  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1050 }
1051 
1052 /**
1053  * @brief Extract entity handle form multi-index container
1054  *
1055  */
1057  template <typename Iterator>
1058  static inline EntityHandle extract(const Iterator &it) {
1059  return (*it)->getEnt();
1060  }
1061 };
1062 
1063 /**
1064  * @brief Insert ordered mofem multi-index into range
1065  *
1066  * \code
1067  * auto hi_rit = refEntsPtr->upper_bound(start);
1068  * auto hi_rit = refEntsPtr->upper_bound(end);
1069  * Range to_erase;
1070  * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1071  * \endcode
1072  *
1073  * @tparam Iterator
1074  * @param r
1075  * @param begin_iter
1076  * @param end_iter
1077  * @return moab::Range::iterator
1078  */
1079 template <typename Extractor, typename Iterator>
1080 moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1081  Iterator end_iter) {
1082  moab::Range::iterator hint = r.begin();
1083  while (begin_iter != end_iter) {
1084  size_t j = 0;
1085  auto bi = Extractor::extract(begin_iter);
1086  Iterator pj = begin_iter;
1087  while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1088  ++pj;
1089  ++j;
1090  }
1091  hint = r.insert(hint, bi, bi + (j - 1));
1092  begin_iter = pj;
1093  }
1094  return hint;
1095 };
1096 
1097 /**
1098  * @brief Do nothing, used to rebuild database
1099  *
1100  */
1103  template <typename T> inline void operator()(T &e) {}
1104 };
1105 
1106 /**
1107  * @brief Template used to reconstruct multi-index
1108  *
1109  * @tparam MI multi-index
1110  * @tparam Modifier
1111  * @param mi
1112  * @param mo
1113  * @return MoFEMErrorCode
1114  */
1115 template <typename MI, typename MO = Modify_change_nothing>
1117  MO &&mo = Modify_change_nothing()) {
1119  for (auto it = mi.begin(); it != mi.end(); ++it) {
1120  if (!const_cast<MI &>(mi).modify(it, mo))
1121  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1122  "Houston we have a problem");
1123  }
1125 }
1126 
1127 } // namespace MoFEM
1128 
1129 #endif //__TEMPLATES_HPP__
static Index< 'm', 3 > m
static Index< 'M', 3 > M
static Index< 'n', 3 > n
static Index< 'L', 3 > L
T data[Tensor_Dim]
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ HVEC0
Definition: definitions.h:255
@ HVEC1
Definition: definitions.h:255
@ HVEC2
Definition: definitions.h:255
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
@ MOFEM_INVALID_DATA
Definition: definitions.h:128
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
@ HVEC1_1
Definition: definitions.h:265
@ HVEC0_1
Definition: definitions.h:264
@ HVEC1_0
Definition: definitions.h:262
@ HVEC2_1
Definition: definitions.h:266
@ HVEC1_2
Definition: definitions.h:268
@ HVEC2_2
Definition: definitions.h:269
@ HVEC2_0
Definition: definitions.h:263
@ HVEC0_2
Definition: definitions.h:267
@ HVEC0_0
Definition: definitions.h:261
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
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:188
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:169
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:273
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:197
FTensor::Index< 'j', 3 > j
const double T
@ row_major
Definition: Layout.hpp:13
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition: Types.hpp:108
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
std::vector< double, std::allocator< double > > DoubleAllocator
Definition: Types.hpp:70
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
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:503
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:44
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec< double, DoubleAllocator >(ublas::vector< double, DoubleAllocator > &data)
Definition: Templates.hpp:150
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:918
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr)
Get FTensor2 from array.
Definition: Templates.hpp:656
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:399
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:563
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr)
Get FTensor1 from array.
Definition: Templates.hpp:624
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:875
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1116
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1080
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:949
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
Definition: Templates.hpp:684
FTensor::Tensor1< FTensor::PackPtr< double *, 2 >, 2 > getFTensor1FromPtr< 2 >(double *ptr)
Definition: Templates.hpp:534
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
Definition: Templates.hpp:594
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
Definition: Templates.hpp:720
FTensor::Tensor1< FTensor::PackPtr< double *, DIM >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition: Templates.hpp:527
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:886
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
Definition: Templates.hpp:541
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
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.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
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:574
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
std::string toString(X x)
Definition: Templates.hpp:119
FTensor::Tensor2< FTensor::PackPtr< double *, DIM1 *DIM2 >, DIM1, DIM2 > getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
Definition: Templates.hpp:555
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:764
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
boost::shared_ptr< std::vector< T > > ShardVec
Definition: Templates.hpp:23
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:905
const double r
rate factor
const int N
Definition: speed_test.cpp:3
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:108
static FTensor::Tensor1< FTensor::PackPtr< double *, S >, 2 > get(MatrixDouble &data)
Definition: Templates.hpp:176
static FTensor::Tensor1< FTensor::PackPtr< double *, S >, 3 > get(MatrixDouble &data)
Definition: Templates.hpp:161
static FTensor::Tensor2_symmetric< FTensor::PackPtr< T *, S >, 2 > get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:297
static FTensor::Tensor2_symmetric< FTensor::PackPtr< T *, S >, 3 > get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:281
static FTensor::Ddg< FTensor::PackPtr< double *, S >, 1, 1 > get(MatrixDouble &data)
Definition: Templates.hpp:332
static FTensor::Ddg< FTensor::PackPtr< double *, S >, 2, 2 > get(MatrixDouble &data)
Definition: Templates.hpp:347
static FTensor::Ddg< FTensor::PackPtr< double *, S >, 3, 3 > get(MatrixDouble &data)
Definition: Templates.hpp:364
static FTensor::Tensor4< FTensor::PackPtr< double *, S >, 1, 1, 1, 1 > get(MatrixDouble &data)
Definition: Templates.hpp:421
static FTensor::Tensor4< FTensor::PackPtr< double *, S >, 2, 2, 2, 2 > get(MatrixDouble &data)
Definition: Templates.hpp:437
static FTensor::Tensor4< FTensor::PackPtr< double *, S >, 3, 3, 3, 3 > get(MatrixDouble &data)
Definition: Templates.hpp:456
unsigned int operator()(const id_type &value) const
Definition: Templates.hpp:114
KeyExtractor1 key1
Definition: Templates.hpp:97
result_type operator()(Arg &arg) const
Definition: Templates.hpp:92
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition: Templates.hpp:88
KeyExtractor2 key2
Definition: Templates.hpp:98
KeyExtractor1::result_type result_type
Definition: Templates.hpp:86
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:102
Do nothing, used to rebuild database.
Definition: Templates.hpp:1101
Extract entity handle form multi-index container.
Definition: Templates.hpp:1056
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:1058