v0.9.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 /**
24  * @brief Get Vector adaptor
25  *
26  * \code
27  *
28  * double *a;
29  * CHKERR VecGetArray(v,&a);
30  *
31  * for(int n = 0; n != nodes; ++n) {
32  *
33  * auto a = getVectorAdaptor(&a[3*n], 3);
34  * double dot = inner_prod(a, a);
35  *
36  * }
37  *
38  * CHKERR VecRetsoreArray(v,&a);
39  * \endcode
40  *
41  */
42 template <typename T1> inline auto getVectorAdaptor(T1 ptr, const size_t n) {
43  typedef typename std::remove_pointer<T1>::type T;
45  ublas::shallow_array_adaptor<T>(n, ptr));
46 };
47 
48 /**
49  * @brief Get Matrix adaptor
50  *
51  * \code
52  *
53  * double *a;
54  * CHKERR VecGetArray(v,&a);
55  *
56  * for(int n = 0; n != nodes; ++n) {
57  *
58  * auto F = getMatrixAdaptor(&a[3*3*n], 3, 3);
59  * MatrixDouble C = prod(F, trans(F));
60  *
61  * }
62  *
63  * CHKERR VecRetsoreArray(v,&a);
64  * \endcode
65  *
66  */
67 template <typename T1>
68 inline auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m) {
69  typedef typename std::remove_pointer<T1>::type T;
71  n, m, ublas::shallow_array_adaptor<T>(n * m, ptr));
72 };
73 
74 /**
75  * This small utility that cascades two key extractors will be
76  * used throughout the boost example
77  * <a
78  * href=http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp>
79  * http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp
80  * </a>
81  */
82 template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
83 public:
84  typedef typename KeyExtractor1::result_type result_type;
85 
86  KeyFromKey(const KeyExtractor1 &key1_ = KeyExtractor1(),
87  const KeyExtractor2 &key2_ = KeyExtractor2())
88  : key1(key1_), key2(key2_) {}
89 
90  template <typename Arg> result_type operator()(Arg &arg) const {
91  return key1(key2(arg));
92  }
93 
94 private:
95  KeyExtractor1 key1;
96  KeyExtractor2 key2;
97 };
98 
99 template <typename id_type> struct LtBit {
100  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
101  return valueA.to_ulong() < valueB.to_ulong();
102  }
103 };
104 
105 template <typename id_type> struct EqBit {
106  inline bool operator()(const id_type &valueA, const id_type &valueB) const {
107  return valueA.to_ulong() == valueB.to_ulong();
108  }
109 };
110 
111 template <typename id_type> struct HashBit {
112  inline unsigned int operator()(const id_type &value) const {
113  return value.to_ulong();
114  }
115 };
116 
117 template <class X> inline std::string toString(X x) {
118  std::ostringstream buffer;
119  buffer << x;
120  return buffer.str();
121 }
122 
123 /**
124 * \brief Get tensor rank 0 (scalar) form data vector
125 
126 Example how to use it.
127 \code
128 VectorDouble vec;
129 vec.resize(nb_gauss_pts,false);
130 vec.clear();
131 auto t0 = getFTensor0FromData(data);
132 for(int gg = 0;gg!=nb_gauss_pts;gg++) {
133 
134  ++t0;
135 }
136 \endcode
137 
138 */
139 template <class T, class A>
141 getFTensor0FromVec(ublas::vector<T, A> &data) {
142  static_assert(!std::is_same<T, T>::value, "not implemented");
144 }
145 
146 template <>
149  ublas::vector<double, DoubleAllocator> &data) {
150  return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&*data.data().begin());
151 }
152 
153 /**
154  * \brief Get tensor rank 1 (vector) form data matrix
155  */
156 template <int Tensor_Dim, class T, class L, class A>
157 static inline FTensor::Tensor1<FTensor::PackPtr<T *, 1>, Tensor_Dim>
158 getFTensor1FromMat(ublas::matrix<T, L, A> &data) {
159  static_assert(!std::is_same<T, T>::value, "not implemented");
160 }
161 
162 /**
163  * \brief Get tensor rank 1 (vector) form data matrix (specialization)
164  */
165 template <int Tensor_Dim>
166 static inline FTensor::Tensor1<FTensor::PackPtr<double *, 1>, Tensor_Dim>
168  return getFTensor1FromMat<Tensor_Dim, double, ublas::row_major,
169  DoubleAllocator>(data);
170 }
171 
172 template <>
174 getFTensor1FromMat<3, double, ublas::row_major, DoubleAllocator>(
175  MatrixDouble &data) {
176  if (data.size1() != 3)
177  THROW_MESSAGE("getFTensor1FromMat<3>: wrong size of data matrix, number of "
178  "rows should be 3 but is %d" +
179  boost::lexical_cast<std::string>(data.size1()));
180 
182  &data(0, 0), &data(1, 0), &data(2, 0));
183 }
184 
185 template <>
187 getFTensor1FromMat<2, double, ublas::row_major, DoubleAllocator>(
188  MatrixDouble &data) {
189  if (data.size1() != 2)
190  THROW_MESSAGE("getFTensor1FromMat<2>: wrong size of data matrix, number of "
191  "rows should be 2 but is %d" +
192  boost::lexical_cast<std::string>(data.size1()));
193 
194  return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2>(&data(0, 0),
195  &data(1, 0));
196 }
197 
198 /**
199  * \brief Get tensor rank 2 (matrix) form data matrix
200  */
201 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
202 static inline FTensor::Tensor2<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
203  Tensor_Dim1>
204 getFTensor2FromMat(ublas::matrix<T, L, A> &data) {
205  static_assert(!std::is_same<T, T>::value,
206  "Such getFTensor2FromMat specialisation is not implemented");
207 }
208 
209 /**
210  * Template specialization for getFTensor2FromMat
211  *
212  */
213 template <>
216  if (data.size1() != 9)
217  THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix; numer "
218  "of rows should be 9 but is " +
219  boost::lexical_cast<std::string>(data.size1()));
220 
222  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
223  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
224 }
225 
226 /**
227  * Template specialization for getFTensor2FromMat
228  */
229 template <>
232  if (data.size1() != 6)
233  THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix, numer "
234  "of rows should be 6 but is " +
235  boost::lexical_cast<std::string>(data.size1()));
236 
238  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
239  &data(5, 0));
240 }
241 
242 /**
243  * Template specialization for getFTensor2FromMat
244  */
245 template <>
248  if (data.size1() != 4)
249  THROW_MESSAGE("getFTensor2FromMat<2,2>: wrong size of data matrix, numer "
250  "of rows should be 4 but is " +
251  boost::lexical_cast<std::string>(data.size1()));
252 
254  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
255 }
256 
257 /**
258  * \brief Get tensor rank 2 (matrix) form data matrix (specialization)
259  */
260 template <int Tensor_Dim0, int Tensor_Dim1>
261 static inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim0,
262  Tensor_Dim1>
264  return getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
265  DoubleAllocator>(data);
266 }
267 
268 /**
269  * \brief Get symmetric tensor rank 2 (matrix) form data matrix
270  */
271 template <int Tensor_Dim, class T, class L, class A>
273 getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
274  static_assert(
275  !std::is_same<T, T>::value,
276  "Such getFTensor2SymmetricFromMat specialisation is not implemented");
277 }
278 
279 /**
280  * @brief Get symmetric tensor rank 2 form matrix of for dimension 3
281  *
282  * Specialisation for symmetric tensor 2
283  *
284  * @tparam
285  * @param data
286  * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, 3>
287  */
288 template <>
291  if (data.size1() != 6)
293  "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
294  "of rows should be 6 but is " +
295  boost::lexical_cast<std::string>(data.size1()));
296 
298  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
299  &data(5, 0));
300 }
301 
302 /**
303  * @brief Get symmetric tensor rank 2 form matrix
304  *
305  * Specialisation for symmetric tensor 2
306  *
307  * @tparam Tensor_Dim
308  * @param data
309  * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, Tensor_Dim>
310  */
311 template <int Tensor_Dim>
313  Tensor_Dim>
316  DoubleAllocator>(data);
317 }
318 
319 /**
320  * @brief Get symmetric tensor rank 4 on first two and last indices from
321  * form data matrix
322  *
323  * @tparam Tensor_Dim01 dimension of frirst two indicies
324  * @tparam Tensor_Dim23 dimension of second two indicies
325  * @tparam T the type of object stored
326  * @tparam L the storage organization
327  * @tparam A the type of Storage array
328  * @param data data container
329  * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
330  */
331 template <int Tensor_Dim01, int Tensor_Dim23, class T, class L, class A>
332 static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
333 getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
334  static_assert(
335  !std::is_same<T, T>::value,
336  "Such getFTensor2SymmetricFromMat specialisation is not implemented");
337 }
338 
339 /**
340  * @brief Get symmetric tensor rank 4 on first two and last indices from
341  * form data matrix
342  *
343  * @param data matrix container which has 36 rows
344  * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
345  */
346 template <>
349  if (data.size1() != 36)
351  "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
352  "of rows should be 36 but is " +
353  boost::lexical_cast<std::string>(data.size1()));
354 
356  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
357  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
358  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
359  &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
360  &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
361  &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
362  &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
363  &data(35, 0));
364 }
365 
366 /**
367  * @brief Make Tensor1 from pointer
368  *
369  * @tparam DIM
370  * @param ptr
371  * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
372  */
373 template <int DIM>
375 getFTensor1FromPtr(double *ptr) {
376  static_assert(DIM != 3,
377  "Such getFTensor1FromPtr specialization is not implemented");
378 };
379 
380 template <>
382 getFTensor1FromPtr<3>(double *ptr) {
384  &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
385 };
386 
387 /**
388  * @brief Make Tensor2 from pointer
389  *
390  * @tparam DIM
391  * @param ptr
392  * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
393  */
394 template <int DIM1, int DIM2>
396 getFTensor2FromPtr(double *ptr) {
397  static_assert(DIM1 != 3, "Such getFTensor2FromPtr is not implemented");
398  static_assert(DIM2 >= 2 && DIM2 <= 3,
399  "Such getFTensor2FromPtr is not implemented");
400 };
401 
402 template <>
404  2> inline getFTensor2FromPtr<3, 2>(double *ptr) {
406  &ptr[HVEC0_0], &ptr[HVEC0_1],
407 
408  &ptr[HVEC1_0], &ptr[HVEC1_1],
409 
410  &ptr[HVEC2_0], &ptr[HVEC2_1]);
411 };
412 
413 template <>
415  3> inline getFTensor2FromPtr<3, 3>(double *ptr) {
417  &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
418 
419  &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
420 
421  &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
422 };
423 
424 /**
425  * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
426  *
427  * @tparam T
428  * @param t
429  * @return double
430  */
431 template <class T> static inline double dEterminant(T &t) {
432  return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
433  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
434  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
435 }
436 
437 /**
438  * \brief Calculate inverse of tensor rank 2 at integration points
439 
440  */
441 template <int Tensor_Dim, class T, class L, class A>
442 inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
443  ublas::vector<T, A> &det_data,
444  ublas::matrix<T, L, A> &inv_jac_data) {
446  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
447  "Specialization for this template not yet implemented");
449 }
450 
451 template <>
452 inline MoFEMErrorCode
453 invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
454  MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
455 
456 /**
457  * \brief Calculate determinant 3 by 3
458 
459  */
460 template <class T1, class T2>
461 inline MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det) {
463  det = +t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
464  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
465  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
467 }
468 
469 /**
470  * \brief Calculate determinant 2 by 2
471 
472  */
473 template <class T1, class T2>
474 inline MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det) {
476  det = t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
478 }
479 
480 /**
481  * \brief Calculate matrix inverse 3 by 3
482 
483  */
484 template <class T1, class T2, class T3>
485 inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
487  const auto inv_det = 1. / det;
488  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
489  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
490  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
491  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
492  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
493  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
494  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
495  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
496  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
498 }
499 
500 /**
501  * \brief Calculate matrix inverse 2 by 2
502 
503  */
504 template <class T1, class T2, class T3>
505 inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
507  const auto inv_det = 1. / det;
508  inv_t(0, 0) = t(1, 1) * inv_det;
509  inv_t(0, 1) = -t(0, 1) * inv_det;
510  inv_t(1, 0) = -t(1, 0) * inv_det;
511  inv_t(1, 1) = t(0, 0) * inv_det;
513 }
514 
515 #ifdef WITH_ADOL_C
516 
517 /**
518  * \brief Calculate matrix inverse, specialization for adouble tensor
519 
520  */
521 template <>
522 inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
527  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
528  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
529  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
530  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
531  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
532  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
533  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
534  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
535  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
537 }
538 
539 #endif
540 
541 /**
542  * \brief Calculate matrix inverse, specialization for symmetric tensor
543 
544  */
545 template <>
546 inline MoFEMErrorCode
547 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
552  const auto inv_det = 1. / det;
553  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
554  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
555  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
556  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
557  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
558  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
560 }
561 
562 #ifdef WITH_ADOL_C
563 
564 /**
565  * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
566 
567  */
568 template <>
569 inline MoFEMErrorCode
570 invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
575  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
576  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
577  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
578  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
579  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
580  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
582 }
583 
584 #endif
585 
586 /**
587  * \brief Calculate matrix inverse, specialization for symmetric (pointer)
588  tensor
589 
590  */
591 template <>
592 inline MoFEMErrorCode
593 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
598  const auto inv_det = 1. / det;
599  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
600  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
601  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
602  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
603  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
604  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
606 }
607 
608 /**
609  * @brief Extract entity handle form multi-index container
610  *
611  */
613  template <typename Iterator>
614  static inline EntityHandle extract(const Iterator &it) {
615  return (*it)->getRefEnt();
616  }
617 };
618 
619 /**
620  * @brief Insert ordered mofem multi-index into range
621  *
622  * \code
623  * auto hi_rit = refEntsPtr->upper_bound(start);
624  * auto hi_rit = refEntsPtr->upper_bound(end);
625  * Range to_erase;
626  * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
627  * \endcode
628  *
629  * @tparam Iterator
630  * @param r
631  * @param begin_iter
632  * @param end_iter
633  * @return moab::Range::iterator
634  */
635 template <typename Extractor, typename Iterator>
636 moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
637  Iterator end_iter) {
638  moab::Range::iterator hint = r.begin();
639  while (begin_iter != end_iter) {
640  size_t j = 0;
641  auto bi = Extractor::extract(begin_iter);
642  Iterator pj = begin_iter;
643  while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
644  ++pj;
645  ++j;
646  }
647  hint = r.insert(hint, bi, bi + (j - 1));
648  begin_iter = pj;
649  }
650  return hint;
651 };
652 
653 } // namespace MoFEM
654 
655 #endif //__TEMPLATES_HPP__
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
FTensor::Tensor2< FTensor::PackPtr< double *, DIM1 *DIM2 >, DIM1, DIM2 > getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
Definition: Templates.hpp:396
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:461
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:614
FTensor::Tensor1< FTensor::PackPtr< double *, DIM >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition: Templates.hpp:375
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition: Types.hpp:110
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec< double, DoubleAllocator >(ublas::vector< double, DoubleAllocator > &data)
Definition: Templates.hpp:148
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition: Templates.hpp:86
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:620
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
KeyExtractor1 key1
Definition: Templates.hpp:95
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:106
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:636
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:474
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:442
Extract entity handle form multi-index container.
Definition: Templates.hpp:612
KeyExtractor2 key2
Definition: Templates.hpp:96
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
Definition: Templates.hpp:382
unsigned int operator()(const id_type &value) const
Definition: Templates.hpp:112
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:415
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:505
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:431
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
static FTensor::Tensor1< FTensor::PackPtr< T *, 1 >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
Definition: Templates.hpp:158
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:42
std::string toString(X x)
Definition: Templates.hpp:117
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:100
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:404
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition: Types.hpp:103
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:333
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:68
result_type operator()(Arg &arg) const
Definition: Templates.hpp:90
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
static FTensor::Tensor2< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:204
KeyExtractor1::result_type result_type
Definition: Templates.hpp:84
std::vector< double, std::allocator< double > > DoubleAllocator
Definition: Types.hpp:70
static FTensor::Tensor2_symmetric< FTensor::PackPtr< T *, 1 >, Tensor_Dim > getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:273