v0.9.1
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 of for dimension 2
304  *
305  * Specialisation for symmetric tensor 2
306  *
307  * @tparam
308  * @param data
309  * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, 2>
310  */
311 template <>
314  if (data.size1() != 3)
316  "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
317  "of rows should be 3 but is " +
318  boost::lexical_cast<std::string>(data.size1()));
320  &data(0, 0), &data(1, 0), &data(2, 0));
321 }
322 
323 /**
324  * @brief Get symmetric tensor rank 2 form matrix
325  *
326  * Specialisation for symmetric tensor 2
327  *
328  * @tparam Tensor_Dim
329  * @param data
330  * @return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, Tensor_Dim>
331  */
332 template <int Tensor_Dim>
334  Tensor_Dim>
337  DoubleAllocator>(data);
338 }
339 
340 /**
341  * @brief Get symmetric tensor rank 4 on first two and last indices from
342  * form data matrix
343  *
344  * @tparam Tensor_Dim01 dimension of frirst two indicies
345  * @tparam Tensor_Dim23 dimension of second two indicies
346  * @tparam T the type of object stored
347  * @tparam L the storage organization
348  * @tparam A the type of Storage array
349  * @param data data container
350  * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
351  */
352 template <int Tensor_Dim01, int Tensor_Dim23, class T, class L, class A>
353 static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
354 getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
355  static_assert(
356  !std::is_same<T, T>::value,
357  "Such getFTensor4DdgFromMat specialisation is not implemented");
358 }
359 
360 /**
361  * @brief Get symmetric tensor rank 4 on first two and last indices from
362  * form data matrix
363  *
364  * @param data matrix container which has 36 rows
365  * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
366  */
367 template <>
370  if (data.size1() != 36)
372  "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
373  "of rows should be 36 but is " +
374  boost::lexical_cast<std::string>(data.size1()));
375 
377  &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
378  &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
379  &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
380  &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
381  &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
382  &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
383  &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
384  &data(35, 0));
385 }
386 
387 /**
388  * @brief Make Tensor1 from pointer
389  *
390  * @tparam DIM
391  * @param ptr
392  * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
393  */
394 template <int DIM>
396 getFTensor1FromPtr(double *ptr) {
397  static_assert(DIM != 3,
398  "Such getFTensor1FromPtr specialization is not implemented");
399 };
400 
401 template <>
403 getFTensor1FromPtr<3>(double *ptr) {
405  &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
406 };
407 
408 /**
409  * @brief Make Tensor2 from pointer
410  *
411  * @tparam DIM
412  * @param ptr
413  * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
414  */
415 template <int DIM1, int DIM2>
417 getFTensor2FromPtr(double *ptr) {
418  static_assert(DIM1 != 3, "Such getFTensor2FromPtr is not implemented");
419  static_assert(DIM2 >= 2 && DIM2 <= 3,
420  "Such getFTensor2FromPtr is not implemented");
421 };
422 
423 template <>
425  2> inline getFTensor2FromPtr<3, 2>(double *ptr) {
427  &ptr[HVEC0_0], &ptr[HVEC0_1],
428 
429  &ptr[HVEC1_0], &ptr[HVEC1_1],
430 
431  &ptr[HVEC2_0], &ptr[HVEC2_1]);
432 };
433 
434 template <>
436  3> inline getFTensor2FromPtr<3, 3>(double *ptr) {
438  &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
439 
440  &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
441 
442  &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
443 };
444 
445 /**
446  * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
447  *
448  * @tparam T
449  * @param t
450  * @return double
451  */
452 template <class T> static inline double dEterminant(T &t) {
453  return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
454  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
455  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
456 }
457 
458 /**
459  * \brief Calculate inverse of tensor rank 2 at integration points
460 
461  */
462 template <int Tensor_Dim, class T, class L, class A>
463 inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
464  ublas::vector<T, A> &det_data,
465  ublas::matrix<T, L, A> &inv_jac_data) {
467  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
468  "Specialization for this template not yet implemented");
470 }
471 
472 template <>
473 inline MoFEMErrorCode
474 invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
475  MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
476 
477 /**
478  * \brief Calculate determinant 3 by 3
479 
480  */
481 template <class T1, class T2>
482 inline MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det) {
484  det = +t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
485  t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
486  t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
488 }
489 
490 /**
491  * \brief Calculate determinant 2 by 2
492 
493  */
494 template <class T1, class T2>
495 inline MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det) {
497  det = t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
499 }
500 
501 /**
502  * \brief Calculate matrix inverse 3 by 3
503 
504  */
505 template <class T1, class T2, class T3>
506 inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
508  const auto inv_det = 1. / det;
509  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
510  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
511  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
512  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
513  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
514  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
515  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
516  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
517  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
519 }
520 
521 /**
522  * \brief Calculate matrix inverse 2 by 2
523 
524  */
525 template <class T1, class T2, class T3>
526 inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
528  const auto inv_det = 1. / det;
529  inv_t(0, 0) = t(1, 1) * inv_det;
530  inv_t(0, 1) = -t(0, 1) * inv_det;
531  inv_t(1, 0) = -t(1, 0) * inv_det;
532  inv_t(1, 1) = t(0, 0) * inv_det;
534 }
535 
536 #ifdef WITH_ADOL_C
537 
538 /**
539  * \brief Calculate matrix inverse, specialization for adouble tensor
540 
541  */
542 template <>
543 inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
548  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
549  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
550  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
551  inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
552  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
553  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
554  inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
555  inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
556  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
558 }
559 
560 #endif
561 
562 /**
563  * \brief Calculate matrix inverse, specialization for symmetric tensor
564 
565  */
566 template <>
567 inline MoFEMErrorCode
568 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
573  const auto inv_det = 1. / det;
574  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
575  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
576  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
577  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
578  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
579  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
581 }
582 
583 #ifdef WITH_ADOL_C
584 
585 /**
586  * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
587 
588  */
589 template <>
590 inline MoFEMErrorCode
591 invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
596  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
597  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
598  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
599  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
600  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
601  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
603 }
604 
605 #endif
606 
607 /**
608  * \brief Calculate matrix inverse, specialization for symmetric (pointer)
609  tensor
610 
611  */
612 template <>
613 inline MoFEMErrorCode
614 invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
619  const auto inv_det = 1. / det;
620  inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
621  inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
622  inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
623  inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
624  inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
625  inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
627 }
628 
629 /**
630  * @brief Extract entity handle form multi-index container
631  *
632  */
634  template <typename Iterator>
635  static inline EntityHandle extract(const Iterator &it) {
636  return (*it)->getRefEnt();
637  }
638 };
639 
640 /**
641  * @brief Insert ordered mofem multi-index into range
642  *
643  * \code
644  * auto hi_rit = refEntsPtr->upper_bound(start);
645  * auto hi_rit = refEntsPtr->upper_bound(end);
646  * Range to_erase;
647  * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
648  * \endcode
649  *
650  * @tparam Iterator
651  * @param r
652  * @param begin_iter
653  * @param end_iter
654  * @return moab::Range::iterator
655  */
656 template <typename Extractor, typename Iterator>
657 moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
658  Iterator end_iter) {
659  moab::Range::iterator hint = r.begin();
660  while (begin_iter != end_iter) {
661  size_t j = 0;
662  auto bi = Extractor::extract(begin_iter);
663  Iterator pj = begin_iter;
664  while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
665  ++pj;
666  ++j;
667  }
668  hint = r.insert(hint, bi, bi + (j - 1));
669  begin_iter = pj;
670  }
671  return hint;
672 };
673 
674 } // namespace MoFEM
675 
676 #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:417
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:482
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:635
FTensor::Tensor1< FTensor::PackPtr< double *, DIM >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition: Templates.hpp:396
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
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::Index< 'n', 2 > n
Definition: PlasticOps.hpp:33
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:625
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
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:657
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:495
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:463
Extract entity handle form multi-index container.
Definition: Templates.hpp:633
KeyExtractor2 key2
Definition: Templates.hpp:96
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
Definition: Templates.hpp:403
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:436
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:526
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:452
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:32
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:425
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:354
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
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:25
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