v0.13.1
Loading...
Searching...
No Matches
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
8namespace MoFEM {
9
10template <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 */
31template <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 */
56template <typename T1>
57inline 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 */
71template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
72public:
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
83private:
84 KeyExtractor1 key1;
85 KeyExtractor2 key2;
86};
87
88template <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
94template <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
100template <typename id_type> struct HashBit {
101 inline unsigned int operator()(const id_type &value) const {
102 return value.to_ulong();
103 }
104};
105
106template <class X> inline std::string toString(X x) {
107 std::ostringstream buffer;
108 buffer << x;
109 return buffer.str();
110}
111
112template <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
121Example how to use it.
122\code
123VectorDouble vec;
124vec.resize(nb_gauss_pts,false);
125vec.clear();
126auto t0 = getFTensor0FromData(data);
127for(int gg = 0;gg!=nb_gauss_pts;gg++) {
128
129 ++t0;
130}
131\endcode
132
133*/
134template <int S = 1, class T, class A>
135static inline auto getFTensor0FromVec(ublas::vector<T, A> &data) {
137}
138
139template <int Tensor_Dim, int S, class T, class L, class A>
141
142template <int S, class T, class A>
143struct 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
157template <int S, class T, class A>
158struct GetFTensor1FromMatImpl<6, 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() != 6)
163 "getFTensor1FromMat<6>: wrong size of data matrix, number of "
164 "rows should be 6 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), &data(4, 0),
169 &data(5, 0));
170 }
171};
172
173template <int S, class T, class A>
174struct GetFTensor1FromMatImpl<2, S, T, ublas::row_major, A> {
175 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
176#ifndef NDDEBUG
177 if (data.size1() != 2)
179 "getFTensor1FromMat<2>: wrong size of data matrix, number of "
180 "rows should be 2 but is " +
181 boost::lexical_cast<std::string>(data.size1()));
182#endif
184 &data(1, 0));
185 }
186};
187
188template <int S, class T, class A>
189struct GetFTensor1FromMatImpl<1, S, T, ublas::row_major, A> {
190 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
191#ifndef NDEBUG
192 if (data.size1() != 1)
194 "getFTensor1FromMat<1>: wrong size of data matrix, number of "
195 "rows should be 1 but is " +
196 boost::lexical_cast<std::string>(data.size1()));
197#endif
199 }
200};
201
202/**
203 * \brief Get tensor rank 1 (vector) form data matrix
204 */
205template <int Tensor_Dim, int S = 1, class T, class L, class A>
207getFTensor1FromMat(ublas::matrix<T, L, A> &data) {
208 static_assert(!std::is_same<T, T>::value, "not implemented");
209}
210
211/**
212 * \brief Get tensor rank 1 (vector) form data matrix (specialization)
213 */
214template <int Tensor_Dim, int S = 1>
216 return GetFTensor1FromMatImpl<Tensor_Dim, S, double, ublas::row_major,
217 DoubleAllocator>::get(data);
218}
219
220/**
221 * \brief Get tensor rank 2 (matrix) form data matrix
222 */
223template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
224inline FTensor::Tensor2<FTensor::PackPtr<T *, 1>, Tensor_Dim0, Tensor_Dim1>
225getFTensor2FromMat(ublas::matrix<T, L, A> &data) {
226 static_assert(!std::is_same<T, T>::value,
227 "Such getFTensor2FromMat specialisation is not implemented");
228}
229
230/**
231 * Template specialization for getFTensor2FromMat
232 */
233template <>
236#ifndef NDEBUG
237 if (data.size1() != 36)
238 THROW_MESSAGE("getFTensor2FromMat<6, 6>: wrong size of data matrix, numer "
239 "of rows should be 36 but is " +
240 boost::lexical_cast<std::string>(data.size1()));
241#endif
242 std::array<double *, 36> ptrs;
243 for (auto i = 0; i != 36; ++i)
244 ptrs[i] = &data(i, 0);
246}
247
248/**
249 * Template specialization for getFTensor2FromMat
250 *
251 */
252template <>
255#ifndef NDEBUG
256 if (data.size1() != 9)
257 THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix; numer "
258 "of rows should be 9 but is " +
259 boost::lexical_cast<std::string>(data.size1()));
260#endif
262 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
263 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
264}
265
266/**
267 * Template specialization for getFTensor2FromMat
268 */
269template <>
272#ifndef NDEBUG
273 if (data.size1() != 6)
274 THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix, numer "
275 "of rows should be 6 but is " +
276 boost::lexical_cast<std::string>(data.size1()));
277#endif
279 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
280 &data(5, 0));
281}
282
283/**
284 * Template specialization for getFTensor2FromMat
285 */
286template <>
289#ifndef NDEBUG
290 if (data.size1() != 18)
291 THROW_MESSAGE("getFTensor2FromMat<6,3>: wrong size of data matrix, numer "
292 "of rows should be 18 but is " +
293 boost::lexical_cast<std::string>(data.size1()));
294#endif
296 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
297 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
298 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
299 &data(15, 0), &data(16, 0), &data(17, 0));
300}
301
302/**
303 * Template specialization for getFTensor2FromMat
304 */
305template <>
308#ifndef NDEBUG
309 if (data.size1() != 4)
310 THROW_MESSAGE("getFTensor2FromMat<2,2>: wrong size of data matrix, numer "
311 "of rows should be 4 but is " +
312 boost::lexical_cast<std::string>(data.size1()));
313#endif
315 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
316}
317
318/**
319 * Template specialization for getFTensor2FromMat
320 */
321template <>
324#ifndef NDEBUG
325 if (data.size1() != 1)
326 THROW_MESSAGE("getFTensor2FromMat<1,1>: wrong size of data matrix, numer "
327 "of rows should be 1 but is " +
328 boost::lexical_cast<std::string>(data.size1()));
329#endif
330 return FTensor::Tensor2<FTensor::PackPtr<double *, 1>, 1, 1>(&data(0, 0));
331}
332
333/**
334 * \brief Get tensor rank 2 (matrix) form data matrix (specialization)
335 */
336template <int Tensor_Dim0, int Tensor_Dim1>
337static inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim0,
338 Tensor_Dim1>
340 return getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
341 DoubleAllocator>(data);
342}
343
344template <int Tensor_Dim, int S, class T, class L, class A>
346
347template <int S, class T, class L, class A>
349 static inline auto get(ublas::matrix<T, L, A> &data) {
350#ifndef NDEBUG
351 if (data.size1() != 6)
353 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
354 "of rows should be 6 but is " +
355 boost::lexical_cast<std::string>(data.size1()));
356#endif
358 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
359 &data(5, 0));
360 }
361};
362
363template <int S, class T, class L, class A>
365 static inline auto get(ublas::matrix<T, L, A> &data) {
366#ifndef NDEBUG
367 if (data.size1() != 3)
369 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
370 "of rows should be 3 but is " +
371 boost::lexical_cast<std::string>(data.size1()));
372#endif
374 &data(0, 0), &data(1, 0), &data(2, 0));
375 }
376};
377
378/**
379 * \brief Get symmetric tensor rank 2 (matrix) form data matrix
380 */
381template <int Tensor_Dim, int S, class T, class L, class A>
382static inline auto getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
384}
385
386template <int Tensor_Dim, int S = 1>
387static inline auto getFTensor2SymmetricFromMat(MatrixDouble &data) {
388 return getFTensor2SymmetricFromMat<Tensor_Dim, S, double, ublas::row_major,
389 DoubleAllocator>(data);
390}
391
392template <int Tensor_Dim01, int Tensor_Dim23, int S, class T, class L, class A>
394
395template <int S, class T, class A>
396struct GetFTensor4DdgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
397 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
398#ifndef NDEBUG
399 if (data.size1() != 1)
401 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
402 "of rows should be 1 but is " +
403 boost::lexical_cast<std::string>(data.size1()));
404#endif
405 return FTensor::Ddg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
406 }
407};
408
409template <int S, class T, class A>
410struct GetFTensor4DdgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
411 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
412#ifndef NDEBUG
413 if (data.size1() != 9) {
415 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
416 "of rows should be 9 but is " +
417 boost::lexical_cast<std::string>(data.size1()));
418 }
419#endif
421 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
422 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
423 }
424};
425
426template <int S, class T, class A>
427struct GetFTensor4DdgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
428 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
429#ifndef NDEBUG
430 if (data.size1() != 36) {
431 cerr << data.size1() << endl;
433 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
434 "of rows should be 36 but is " +
435 boost::lexical_cast<std::string>(data.size1()));
436 }
437#endif
439 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
440 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
441 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
442 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
443 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
444 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
445 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
446 &data(35, 0)};
447 }
448};
449
450/**
451 * @brief Get symmetric tensor rank 4 on first two and last indices from
452 * form data matrix
453 *
454 * @tparam Tensor_Dim01 dimension of first two indicies
455 * @tparam Tensor_Dim23 dimension of second two indicies
456 * @tparam T the type of object stored
457 * @tparam L the storage organization
458 * @tparam A the type of Storage array
459 * @param data data container
460 * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
461 */
462template <int Tensor_Dim01, int Tensor_Dim23, int S = 1, class T, class L,
463 class A>
464static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
465getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
466 static_assert(!std::is_same<T, T>::value,
467 "Such getFTensor4DdgFromMat specialisation is not implemented");
468}
469
470template <int Tensor_Dim01, int Tensor_Dim23, int S = 1>
471static inline auto getFTensor4DdgFromMat(MatrixDouble &data) {
472 return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, double,
473 ublas::row_major,
474 DoubleAllocator>::get(data);
475}
476
477template <int Tensor_Dim01, int Tensor_Dim2, int S, class T, class L, class A>
479
480template <int S, class T, class A>
481struct GetFTensor3DgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
482 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
483#ifndef NDEBUG
484 if (data.size1() != 1)
486 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
487 "of rows should be 1 but is " +
488 boost::lexical_cast<std::string>(data.size1()));
489#endif
490 return FTensor::Dg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
491 }
492};
493
494template <int S, class T, class A>
495struct GetFTensor3DgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
496 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
497#ifndef NDEBUG
498 if (data.size1() != 6) {
500 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
501 "of rows should be 6 but is " +
502 boost::lexical_cast<std::string>(data.size1()));
503 }
504#endif
506 &data(0, 0), &data(1, 0), &data(2, 0),
507 &data(3, 0), &data(4, 0), &data(5, 0)};
508 }
509};
510
511template <int S, class T, class A>
512struct GetFTensor3DgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
513 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
514#ifndef NDEBUG
515 if (data.size1() != 18) {
516 cerr << data.size1() << endl;
518 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
519 "of rows should be 18 but is " +
520 boost::lexical_cast<std::string>(data.size1()));
521 }
522#endif
524 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
525 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
526 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
527 &data(15, 0), &data(16, 0), &data(17, 0)};
528 }
529};
530
531/**
532 * @brief Get symmetric tensor rank 3 on the first two indices from
533 * form data matrix
534 *
535 * @tparam Tensor_Dim01 dimension of first two indicies
536 * @tparam Tensor_Dim2 dimension of last index
537 * @tparam T the type of object stored
538 * @tparam L the storage organization
539 * @tparam A the type of Storage array
540 * @param data data container
541 * @return FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
542 */
543template <int Tensor_Dim01, int Tensor_Dim2, int S = 1, class T, class L,
544 class A>
545static inline FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim2>
546getFTensor3DgFromMat(ublas::matrix<T, L, A> &data) {
547 static_assert(!std::is_same<T, T>::value,
548 "Such getFTensor3DgFromMat specialisation is not implemented");
549}
550
551template <int Tensor_Dim01, int Tensor_Dim2, int S = 1>
552static inline auto getFTensor3DgFromMat(MatrixDouble &data) {
553 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, double,
554 ublas::row_major, DoubleAllocator>::get(data);
555}
556
557template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
558 int S, class T, class L, class A>
560
561template <int S, class T, class A>
562struct GetFTensor4FromMatImpl<1, 1, 1, 1, S, T, ublas::row_major, A> {
563 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
564#ifndef NDEBUG
565 if (data.size1() != 1)
567 "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
568 "of rows should be 1 but is " +
569 boost::lexical_cast<std::string>(data.size1()));
570#endif
572 &data(0, 0)};
573 }
574};
575
576template <int S, class T, class A>
577struct GetFTensor4FromMatImpl<2, 2, 2, 2, S, T, ublas::row_major, A> {
578 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
579#ifndef NDEBUG
580 if (data.size1() != 16) {
582 "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
583 "of rows should be 16 but is " +
584 boost::lexical_cast<std::string>(data.size1()));
585 }
586#endif
588 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
589 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
590 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
591 &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
592 }
593};
594
595template <int S, class T, class A>
596struct GetFTensor4FromMatImpl<3, 3, 3, 3, S, T, ublas::row_major, A> {
597 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
598#ifndef NDEBUG
599 if (data.size1() != 81) {
600 cerr << data.size1() << endl;
602 "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
603 "of rows should be 81 but is " +
604 boost::lexical_cast<std::string>(data.size1()));
605 }
606#endif
608 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
609 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
610 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
611 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
612 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
613 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
614 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
615 &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
616 &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
617 &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
618 &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
619 &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
620 &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
621 &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
622 &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
623 &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
624 &data(80, 0)};
625 }
626};
627
628/**
629 * @brief Get tensor rank 4 (non symmetric) form data matrix
630 *
631 * @tparam Tensor_Dim0 dimension of frirst index
632 * @tparam Tensor_Dim1 dimension of second index
633 * @tparam Tensor_Dim2 dimension of third index
634 * @tparam Tensor_Dim3 dimension of fourth index
635 * @tparam T the type of object stored
636 * @tparam L the storage organization
637 * @tparam A the type of Storage array
638 * @param data data container
639 * @return FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
640 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
641 */
642template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
643 int S = 1, class T, class L, class A>
644static inline FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
645 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
646getFTensor4FromMat(ublas::matrix<T, L, A> &data) {
647 static_assert(!std::is_same<T, T>::value,
648 "Such getFTensor4FromMat specialisation is not implemented");
649}
650
651template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
652 int S = 1>
653static inline auto getFTensor4FromMat(MatrixDouble &data) {
654 return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
655 Tensor_Dim3, S, double, ublas::row_major,
656 DoubleAllocator>::get(data);
657}
658
659template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class T,
660 class L, class A>
662
663template <int S, class T, class A>
664struct GetFTensor3FromMatImpl<1, 1, 1, S, T, ublas::row_major, A> {
665 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
666#ifndef NDEBUG
667 if (data.size1() != 1)
669 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
670 "of rows should be 1 but is " +
671 boost::lexical_cast<std::string>(data.size1()));
672#endif
674 &data(0, 0)};
675 }
676};
677
678template <int S, class T, class A>
679struct GetFTensor3FromMatImpl<2, 2, 2, S, T, ublas::row_major, A> {
680 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
681#ifndef NDEBUG
682 if (data.size1() != 8)
684 "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
685 "of rows should be 8 but is " +
686 boost::lexical_cast<std::string>(data.size1()));
687#endif
689 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
690 &data(5, 0), &data(6, 0), &data(7, 0)
691
692 };
693 }
694};
695
696template <int S, class T, class A>
697struct GetFTensor3FromMatImpl<3, 2, 2, S, T, ublas::row_major, A> {
698 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
699#ifndef NDEBUG
700 if (data.size1() != 12)
702 "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
703 "of rows should be 12 but is " +
704 boost::lexical_cast<std::string>(data.size1()));
705#endif
707 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
708 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
709 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
710 }
711};
712
713template <int S, class T, class A>
714struct GetFTensor3FromMatImpl<2, 2, 3, S, T, ublas::row_major, A> {
715 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
716#ifndef NDEBUG
717 if (data.size1() != 12)
719 "getFTensor3FromMat<2, 2, 3>: wrong size of data matrix, number "
720 "of rows should be 12 but is " +
721 boost::lexical_cast<std::string>(data.size1()));
722#endif
724 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
725 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
726 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
727 }
728};
729
730template <int S, class T, class A>
731struct GetFTensor3FromMatImpl<3, 3, 3, S, T, ublas::row_major, A> {
732 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
733#ifndef NDEBUG
734 if (data.size1() != 27)
736 "getFTensor3FromMat<3, 3, 3>: wrong size of data matrix, number "
737 "of rows should be 27 but is " +
738 boost::lexical_cast<std::string>(data.size1()));
739#endif
741 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
742 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
743 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
744 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
745 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
746 &data(25, 0), &data(26, 0)};
747 }
748};
749
750template <int S, class T, class A>
751struct GetFTensor3FromMatImpl<6, 3, 3, S, T, ublas::row_major, A> {
752 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
753#ifndef NDEBUG
754 if (data.size1() != 54)
756 "getFTensor3FromMat<6, 3, 3>: wrong size of data matrix, number "
757 "of rows should be 54 but is " +
758 boost::lexical_cast<std::string>(data.size1()));
759#endif
760 std::array<double *, 54> ptrs;
761 for (auto i = 0; i != 54; ++i)
762 ptrs[i] = &data(i, 0);
764 }
765};
766
767template <int S, class T, class A>
768struct GetFTensor3FromMatImpl<3, 3, 6, S, T, ublas::row_major, A> {
769 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
770#ifndef NDEBUG
771 if (data.size1() != 54)
773 "getFTensor3FromMat<3, 3, 6>: wrong size of data matrix, number "
774 "of rows should be 54 but is " +
775 boost::lexical_cast<std::string>(data.size1()));
776#endif
777 std::array<double *, 54> ptrs;
778 for (auto i = 0; i != 54; ++i)
779 ptrs[i] = &data(i, 0);
781 }
782};
783
784/**
785 * @brief Get tensor rank 3 (non symmetries) form data matrix
786 *
787 * @tparam Tensor_Dim0 dimension of first index
788 * @tparam Tensor_Dim1 dimension of second index
789 * @tparam Tensor_Dim2 dimension of third index
790 * @tparam S shift size
791 * @tparam T the type of object stored
792 * @tparam L the storage organization
793 * @tparam A the type of Storage array
794 * @param data data container
795 * @return FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
796 Tensor_Dim1, Tensor_Dim2>
797 */
798template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1, class T,
799 class L, class A>
800static inline FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
801 Tensor_Dim1, Tensor_Dim2>
802getFTensor3FromMat(ublas::matrix<T, L, A> &data) {
803 static_assert(!std::is_same<T, T>::value,
804 "Such getFTensor3FromMat specialisation is not implemented");
805}
806
807template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1>
808static inline auto getFTensor3FromMat(MatrixDouble &data) {
809 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
810 double, ublas::row_major,
811 DoubleAllocator>::get(data);
812}
813
814template<int DIM, int S = DIM>
816
817template <int S> struct GetFTensor1FromPtrImpl<2, S> {
819 inline static auto get(double *ptr) {
821 &ptr[HVEC1]);
822 }
823};
824
825template <int S> struct GetFTensor1FromPtrImpl<3, S> {
827 inline static auto get(double *ptr) {
829 &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
830 }
831};
832
833/**
834 * @brief Make Tensor1 from pointer
835 *
836 * @tparam DIM
837 * @param ptr
838 * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
839 */
840template <int DIM, int S = DIM>
842getFTensor1FromPtr(double *ptr) {
844};
845
846/**
847 * @brief Make Tensor2 from pointer
848 *
849 * @tparam DIM
850 * @param ptr
851 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
852 */
853template <int DIM1, int DIM2>
855getFTensor2FromPtr(double *ptr) {
856 static_assert(DIM1 == DIM1 || DIM2 != DIM2,
857 "Such getFTensor2FromPtr is not implemented");
858};
859
860template <>
862 2> inline getFTensor2FromPtr<3, 2>(double *ptr) {
864 &ptr[HVEC0_0], &ptr[HVEC0_1],
865
866 &ptr[HVEC1_0], &ptr[HVEC1_1],
867
868 &ptr[HVEC2_0], &ptr[HVEC2_1]);
869};
870
871template <>
873 3> inline getFTensor2FromPtr<3, 3>(double *ptr) {
875 &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
876
877 &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
878
879 &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
880};
881
882template <>
884 2> inline getFTensor2FromPtr<2, 2>(double *ptr) {
886 &ptr[0], &ptr[1], &ptr[2], &ptr[3]);
887};
888
889/*
890 * @brief Make Tensor3 from pointer
891 *
892 * @tparam DIM
893 * @param ptr
894 * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3>, DIM1,
895 * DIM2, DIM3>
896 */
897template <int DIM1, int DIM2, int DIM3>
899 DIM2, DIM3>
900getFTensor3FromPtr(double *ptr) {
901 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
902 "Such getFTensor2FromPtr is not implemented");
903};
904
905template <>
907getFTensor3FromPtr<3, 2, 2>(double *ptr) {
909 &ptr[0], &ptr[1], &ptr[2],
910
911 &ptr[3], &ptr[4], &ptr[5],
912
913 &ptr[6], &ptr[7], &ptr[8],
914
915 &ptr[9], &ptr[10], &ptr[11]
916
917 );
918};
919
920/**
921 * @brief Make symmetric Tensor2 from pointer, taking lower triangle of matrix
922 *
923 * @tparam DIM
924 * @param ptr
925 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
926 */
927template <int DIM>
930 static_assert(DIM,
931 "Such getFTensor2SymmetricUpperFromPtr is not implemented");
932}
933
934template <>
938 &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
939
940 &ptr[HVEC1_0], &ptr[HVEC1_1],
941
942 &ptr[HVEC2_2]);
943};
944
945template <>
949 &ptr[0], &ptr[1], &ptr[3]);
950};
951
952template <int DIM, int S> struct GetFTensor1FromArray;
953
954template <int S> struct GetFTensor1FromArray<2, S> {
956 template <typename V> static inline auto get(V &data) {
957 using T = typename std::remove_reference<decltype(data[0])>::type;
958 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 2>{&data[0], &data[1]};
959 }
960};
961
962template <int S> struct GetFTensor1FromArray<3, S> {
964 template <typename V> static inline auto get(V &data) {
965 using T = typename std::remove_reference<decltype(data[0])>::type;
967 &data[0], &data[1], &data[2]};
968 }
969};
970
971template <int S> struct GetFTensor1FromArray<6, S> {
973 template <typename V> static inline auto get(V &data) {
974 using T = typename std::remove_reference<decltype(data[0])>::type;
976 &data[0], &data[1], &data[2], &data[3], &data[4], &data[5]};
977 }
978};
979
980/**
981 * @brief Get FTensor1 from array
982 *
983 * \todo Generalise for different arrays and data types
984 *
985 * @tparam DIM
986 * @param data
987 * @return FTensor::Tensor1<FTensor::PackPtr<double *, S>, DIM>
988 */
989template <int DIM, int S> inline auto getFTensor1FromArray(VectorDouble &data) {
991}
992
993template <int DIM, int S>
995getFTensor1FromMat(MatrixDouble &data, const size_t rr);
996
997template <>
999getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1000 return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2>{&data(rr + 0, 0),
1001 &data(rr + 1, 0)};
1002}
1003
1004template <>
1006getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1008 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1009}
1010
1011/**
1012 * @brief Get FTensor1 from array
1013 *
1014 * \todo Generalise for diffrent arrays and data types
1015 *
1016 * @tparam DIM
1017 * @param data
1018 * @param rr
1019 * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
1020 */
1021template <int DIM, int S>
1024 static_assert(DIM != DIM, "not implemented");
1026}
1027
1028template <>
1031 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
1032 &data(rr + 1, 1)};
1033}
1034
1035template <>
1039 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1040}
1041
1042/**
1043 * @brief Get FTensor2 from array
1044 *
1045 * \note Generalise for other data types
1046 *
1047 * @tparam DIM1
1048 * @tparam DIM2
1049 * @tparam S
1050 * @param data
1051 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
1052 */
1053template <int DIM1, int DIM2, int S, class T, class L, class A>
1055
1056template <int S, class T, class L, class A>
1057struct GetFTensor2FromArrayImpl<2, 2, S, T, L, A> {
1059 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr) {
1061 &data(rr + 0, 0), &data(rr + 0, 1), &data(rr + 1, 0), &data(rr + 1, 1)};
1062 }
1063};
1064
1065template <int S, class T, class L, class A>
1066struct GetFTensor2FromArrayImpl<3, 3, S, T, L, A> {
1068 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr) {
1070 &data(rr + 0, 0), &data(rr + 0, 1), &data(rr + 0, 2),
1071 &data(rr + 1, 0), &data(rr + 1, 1), &data(rr + 1, 2),
1072 &data(rr + 2, 0), &data(rr + 2, 1), &data(rr + 2, 2)};
1073 }
1074};
1075
1076template <int DIM1, int DIM2, int S>
1078getFTensor2FromArray(MatrixDouble &data, const size_t rr) {
1079 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, double, ublas::row_major,
1080 VecAllocator<double>>::get(data, rr);
1081}
1082
1083template <int S, typename T, typename L, typename A>
1084inline auto getFTensor2FromArray2by2(ublas::matrix<T, L, A> &data,
1085 const FTensor::Number<S> &,
1086 const size_t rr) {
1088}
1089
1090template <int S, typename T, typename L, typename A>
1091inline auto getFTensor2FromArray3by3(ublas::matrix<T, L, A> &data,
1092 const FTensor::Number<S> &,
1093 const size_t rr) {
1095}
1096
1097#ifdef WITH_ADOL_C
1098
1099template <int DIM1, int DIM2, int S>
1100inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
1101 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, adouble, ublas::row_major,
1102 VecAllocator<adouble>>::get(data, rr);
1103}
1104
1105#endif
1106
1107// list of lapack wrappers
1108/**
1109 * @brief compute matrix inverse with lapack dgetri
1110 *
1111 * @param mat input square matrix / output inverse matrix
1112 * @return MoFEMErrorCode
1113 */
1116
1117 const size_t M = mat.size1();
1118 const size_t N = mat.size2();
1119
1120 if (M != N)
1121 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1122 "The input matrix for inverse computation is not square %d != %d",
1123 M, N);
1124
1125 int *ipv = new int[N];
1126 int lwork = N * N;
1127 double *work = new double[lwork];
1128 int info;
1129 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1130 if (info != 0)
1131 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1132 "lapack error info = %d", info);
1133 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1134 if (info != 0)
1135 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1136 "lapack error info = %d", info);
1137
1138 delete[] ipv;
1139 delete[] work;
1140
1142}
1143/**
1144 * @brief solve linear system with lapack dgesv
1145 *
1146 * @param mat input lhs square matrix / output L and U from the factorization
1147 * @param f input rhs vector / output solution vector
1148 * @return MoFEMErrorCode
1149 */
1152
1153 const size_t M = mat.size1();
1154 const size_t N = mat.size2();
1155
1156 if (M == 0 || M != N)
1157 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1158 "The input matrix for inverse computation is not square %d != %d",
1159 M, N);
1160 if (f.size() != M)
1161 f.resize(M, false);
1162
1163 const int nrhs = 1;
1164 int info;
1165 int *ipiv = new int[M];
1166 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1167 &*f.data().begin(), nrhs);
1168
1169 if (info != 0) {
1170 SETERRQ1(PETSC_COMM_SELF, 1, "error lapack solve dgesv info = %d", info);
1171 }
1172
1173 delete[] ipiv;
1175}
1176
1177/**
1178 * @brief Solve linear system of equations using Lapack
1179 *
1180 * @param mat
1181 * @param f
1182 * @return MoFEMErrorCode
1183 */
1185 VectorDouble &f) {
1187 // copy matrix since on output lapack returns factorisation
1188 auto mat_copy = mat;
1189 CHKERR solveLinearSystem(mat_copy, f);
1191}
1192
1193/**
1194 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1195 *
1196 * @param mat input symmetric matrix
1197 * @param eig output eigen values sorted
1198 * @param eigen_vec output matrix of row eigen vectors
1199 * @return MoFEMErrorCode
1200 */
1202 VectorDouble &eig,
1203 MatrixDouble &eigen_vec) {
1205
1206 const size_t M = mat.size1();
1207 const size_t N = mat.size2();
1208
1209 if (M == 0 || M != N)
1210 SETERRQ2(
1211 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1212 "The input matrix for eigen value computation is not square %d != %d",
1213 M, N);
1214 if (eig.size() != M)
1215 eig.resize(M, false);
1216
1217 eigen_vec = mat;
1218 const int n = M;
1219 const int lda = M;
1220 const int size = (M + 2) * M;
1221 int lwork = size;
1222 double *work = new double[size];
1223
1224 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
1225 &*eig.data().begin(), work, lwork) > 0)
1226 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1227 "The algorithm failed to compute eigenvalues.");
1228
1229 delete[] work;
1231}
1232/**
1233 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1234 *
1235 * @tparam DIM
1236 * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
1237 * @param eig output eigen values sorted
1238 * @return MoFEMErrorCode
1239 */
1240template <int DIM>
1241inline MoFEMErrorCode
1245
1246 const int n = DIM;
1247 const int lda = DIM;
1248 const int lwork = (DIM + 2) * DIM;
1249 std::array<double, (DIM + 2) * DIM> work;
1250
1251 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
1252 lwork) > 0)
1253 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1254 "The algorithm failed to compute eigenvalues.");
1256}
1257/**
1258 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1259 *
1260 * @tparam DIM
1261 * @param mat input tensor pointer of size DIM x DIM
1262 * @param eig output eigen values sorted
1263 * @param eigen_vec output matrix of row eigen vectors
1264 * @return MoFEMErrorCode
1265 */
1266template <int DIM>
1272 for (int ii = 0; ii != DIM; ii++)
1273 for (int jj = 0; jj != DIM; jj++)
1274 eigen_vec(ii, jj) = mat(ii, jj);
1275
1276 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1277
1279}
1280
1281/**
1282 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1283 *
1284 * @tparam DIM
1285 * @param mat input tensor of size DIM x DIM
1286 * @param eig output eigen values sorted
1287 * @param eigen_vec output matrix of row eigen vectors
1288 * @return MoFEMErrorCode
1289 */
1290template <int DIM>
1291inline MoFEMErrorCode
1296 for (int ii = 0; ii != DIM; ii++)
1297 for (int jj = 0; jj != DIM; jj++)
1298 eigen_vec(ii, jj) = mat(ii, jj);
1299
1300 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1301
1303}
1304
1305/**
1306 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
1307 *
1308 * @tparam T
1309 * @param t
1310 * @return double
1311 */
1312template <typename T> static inline auto determinantTensor3by3(T &t) {
1313 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
1314 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
1315 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
1316}
1317
1318/**
1319 * \brief Calculate inverse of tensor rank 2 at integration points
1320
1321 */
1322template <int Tensor_Dim, class T, class L, class A>
1323inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
1324 ublas::vector<T, A> &det_data,
1325 ublas::matrix<T, L, A> &inv_jac_data) {
1327 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1328 "Specialization for this template not yet implemented");
1330}
1331
1332template <>
1333inline MoFEMErrorCode
1335 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
1336
1337/**
1338 * \brief Calculate determinant 3 by 3
1339
1340 */
1341template <class T1, class T2>
1344 det = determinantTensor3by3(t);
1346}
1347
1348/**
1349 * \brief Calculate determinant 2 by 2
1350
1351 */
1352template <class T1, class T2>
1355 det = t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
1357}
1358
1359/**
1360 * \brief Calculate matrix inverse 3 by 3
1361
1362 */
1363template <class T1, class T2, class T3>
1364inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
1366 const auto inv_det = 1. / det;
1367 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1368 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1369 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1370 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
1371 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1372 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1373 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
1374 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
1375 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1377}
1378
1379/**
1380 * \brief Calculate matrix inverse 2 by 2
1381
1382 */
1383template <class T1, class T2, class T3>
1384inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
1386 const auto inv_det = 1. / det;
1387 inv_t(0, 0) = t(1, 1) * inv_det;
1388 inv_t(0, 1) = -t(0, 1) * inv_det;
1389 inv_t(1, 0) = -t(1, 0) * inv_det;
1390 inv_t(1, 1) = t(0, 0) * inv_det;
1392}
1393
1394#ifdef WITH_ADOL_C
1395
1396/**
1397 * \brief Calculate matrix inverse, specialization for adouble tensor
1398
1399 */
1400template <>
1401inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
1406 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1407 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1408 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1409 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
1410 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1411 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1412 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
1413 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
1414 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1416}
1417
1418#endif
1419
1420/**
1421 * \brief Calculate matrix inverse, specialization for symmetric tensor
1422
1423 */
1424template <>
1425inline MoFEMErrorCode
1426invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1431 const auto inv_det = 1. / det;
1432 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1433 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1434 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1435 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1436 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1437 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1439}
1440
1441#ifdef WITH_ADOL_C
1442
1443/**
1444 * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1445
1446 */
1447template <>
1448inline MoFEMErrorCode
1449invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
1454 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1455 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1456 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1457 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1458 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1459 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1461}
1462
1463#endif
1464
1465/**
1466 * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1467 tensor
1468
1469 */
1470template <>
1471inline MoFEMErrorCode
1472invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1477 const auto inv_det = 1. / det;
1478 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1479 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1480 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1481 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1482 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1483 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1485}
1486
1487/**
1488 * @brief Extract entity handle form multi-index container
1489 *
1490 */
1492 template <typename Iterator>
1493 static inline EntityHandle extract(const Iterator &it) {
1494 return (*it)->getEnt();
1495 }
1496};
1497
1498/**
1499 * @brief Insert ordered mofem multi-index into range
1500 *
1501 * \code
1502 * auto hi_rit = refEntsPtr->upper_bound(start);
1503 * auto hi_rit = refEntsPtr->upper_bound(end);
1504 * Range to_erase;
1505 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1506 * \endcode
1507 *
1508 * @tparam Iterator
1509 * @param r
1510 * @param begin_iter
1511 * @param end_iter
1512 * @return moab::Range::iterator
1513 */
1514template <typename Extractor, typename Iterator>
1515moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1516 Iterator end_iter) {
1517 moab::Range::iterator hint = r.begin();
1518 while (begin_iter != end_iter) {
1519 size_t j = 0;
1520 auto bi = Extractor::extract(begin_iter);
1521 Iterator pj = begin_iter;
1522 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1523 ++pj;
1524 ++j;
1525 }
1526 hint = r.insert(hint, bi, bi + (j - 1));
1527 begin_iter = pj;
1528 }
1529 return hint;
1530};
1531
1532/**
1533 * @brief Do nothing, used to rebuild database
1534 *
1535 */
1538 template <typename T> inline void operator()(T &e) {}
1539};
1540
1541/**
1542 * @brief Template used to reconstruct multi-index
1543 *
1544 * @tparam MI multi-index
1545 * @tparam Modifier
1546 * @param mi
1547 * @param mo
1548 * @return MoFEMErrorCode
1549 */
1550template <typename MI, typename MO = Modify_change_nothing>
1552 MO &&mo = Modify_change_nothing()) {
1554 for (auto it = mi.begin(); it != mi.end(); ++it) {
1555 if (!const_cast<MI &>(mi).modify(it, mo))
1556 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1557 "Houston we have a problem");
1558 }
1560}
1561
1563 TempMeshset(moab::Interface &moab) : moab(moab) {
1564 rval = moab.create_meshset(MESHSET_SET, meshset);
1566 }
1568 operator EntityHandle() const { return meshset; }
1569 auto get_ptr() { return &meshset; }
1570
1571private:
1573 rval = moab.delete_entities(&meshset, 1);
1575 }
1577 moab::Interface &moab;
1578};
1579
1580/**
1581 * @brief Create smart pointer to temprary meshset
1582 *
1583 */
1584inline auto get_temp_meshset_ptr(moab::Interface &moab) {
1585 return boost::make_shared<TempMeshset>(moab);
1586};
1587
1588/**
1589 * @brief get type from entity handle
1590 *
1591 */
1592inline auto type_from_handle(const EntityHandle h) {
1593 return static_cast<EntityType>(h >> MB_ID_WIDTH);
1594};
1595
1596/**
1597 * @brief get entity dimension form handle
1598 *
1599 */
1601 return moab::CN::Dimension(type_from_handle(h));
1602};
1603
1604/**
1605 * @brief get field bit id from bit number
1606 *
1607 */
1608inline auto field_bit_from_bit_number(const int bit_number) {
1609 return BitFieldId().set(bit_number - 1);
1610};
1611
1612/**
1613 * @brief Insert ranges
1614 *
1615 * @tparam I
1616 * @param f
1617 * @param s
1618 * @param tester
1619 * @param inserter
1620 * @return auto
1621 */
1622template <typename I>
1623auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
1624 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
1626
1627 auto first = f;
1628 while (first != s)
1629 if (tester(first)) {
1630
1631 auto second = first;
1632 ++second;
1633
1634 while (second != s) {
1635 if (tester(second))
1636 ++second;
1637 else
1638 break;
1639 }
1640
1641 CHKERR inserter(first, second);
1642
1643 first = second;
1644 if (first != s)
1645 ++first;
1646
1647 } else {
1648 ++first;
1649 }
1650
1652}
1653
1654// template <typename T, typename I>
1655// std::vector<std::pair<T, T>> getPairRange(I s, I e, boost::function<T(I)>){
1656
1657// };
1658
1659} // namespace MoFEM
1660
1661#endif //__TEMPLATES_HPP__
static Index< 'M', 3 > M
static Index< 'I', 3 > I
T data[Tensor_Dim]
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MB_ID_WIDTH
Definition: definitions.h:227
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ HVEC0
Definition: definitions.h:186
@ HVEC1
Definition: definitions.h:186
@ HVEC2
Definition: definitions.h:186
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
@ HVEC1_1
Definition: definitions.h:196
@ HVEC0_1
Definition: definitions.h:195
@ HVEC1_0
Definition: definitions.h:193
@ HVEC2_1
Definition: definitions.h:197
@ HVEC1_2
Definition: definitions.h:199
@ HVEC2_2
Definition: definitions.h:200
@ HVEC2_0
Definition: definitions.h:194
@ HVEC0_2
Definition: definitions.h:198
@ HVEC0_0
Definition: definitions.h:192
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
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
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
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
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
FTensor::Index< 'j', 3 > j
const double T
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::vector< T, std::allocator< T > > VecAllocator
Definition: Types.hpp:59
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:42
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition: Types.hpp:114
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasMatrix< adouble > MatrixADouble
Definition: Types.hpp:80
VecAllocator< double > DoubleAllocator
Definition: Types.hpp:62
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition: Types.hpp:121
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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:802
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:31
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:546
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1592
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1353
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:225
FTensor::Tensor3< FTensor::PackPtr< double *, 12 >, 3, 2, 2 > getFTensor3FromPtr< 3, 2, 2 >(double *ptr)
Definition: Templates.hpp:907
FTensor::Tensor2< FTensor::PackPtr< double *, 4 >, 2, 2 > getFTensor2FromPtr< 2, 2 >(double *ptr)
Definition: Templates.hpp:884
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1551
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1515
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1384
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:1091
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
Definition: Templates.hpp:1114
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
Definition: Templates.hpp:1600
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
Definition: Templates.hpp:1150
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:1323
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:465
MoFEMErrorCode invertTensor3by3< 3, double, ublas::row_major, DoubleAllocator >(MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data)
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 4 >, 2 > getFTensor2SymmetricLowerFromPtr< 2 >(double *ptr)
Definition: Templates.hpp:947
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:207
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:57
FTensor::Tensor2< FTensor::PackPtr< double *, DIM1 *DIM2 >, DIM1, DIM2 > getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
Definition: Templates.hpp:855
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:382
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition: Templates.hpp:842
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
Definition: Templates.hpp:1608
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
std::string toString(X x)
Definition: Templates.hpp:106
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:646
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:873
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr)
Get FTensor1 from array.
Definition: Templates.hpp:1023
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr)
Definition: Templates.hpp:1078
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temprary meshset.
Definition: Templates.hpp:1584
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1201
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:929
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1312
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:862
FTensor::Tensor3< FTensor::PackPtr< double *, DIM1 *DIM2 *DIM3 >, DIM1, DIM2, DIM3 > getFTensor3FromPtr(double *ptr)
Definition: Templates.hpp:900
boost::shared_ptr< std::vector< T > > ShardVec
Definition: Templates.hpp:10
auto getFTensor2FromArray2by2(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:1084
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:1623
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
Definition: Templates.hpp:989
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 9 >, 3 > getFTensor2SymmetricLowerFromPtr< 3 >(double *ptr)
Definition: Templates.hpp:936
double h
constexpr AssemblyType A
Definition: plastic.cpp:46
constexpr double t
plate stiffness
Definition: plate.cpp:59
const int N
Definition: speed_test.cpp:3
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:95
static auto get(ublas::vector< T, A > &data)
Definition: Templates.hpp:113
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:190
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:175
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:144
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:159
static auto get(double *ptr)
Definition: Templates.hpp:819
static auto get(double *ptr)
Definition: Templates.hpp:827
static auto get(ublas::matrix< T, L, A > &data, const size_t rr)
Definition: Templates.hpp:1059
static auto get(ublas::matrix< T, L, A > &data, const size_t rr)
Definition: Templates.hpp:1068
Get FTensor2 from array.
Definition: Templates.hpp:1054
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:365
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:349
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:482
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:496
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:513
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:665
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:680
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:715
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:698
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:732
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:769
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:752
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:397
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:411
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:428
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:563
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:578
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:597
unsigned int operator()(const id_type &value) const
Definition: Templates.hpp:101
KeyExtractor1 key1
Definition: Templates.hpp:84
result_type operator()(Arg &arg) const
Definition: Templates.hpp:79
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition: Templates.hpp:75
KeyExtractor2 key2
Definition: Templates.hpp:85
KeyExtractor1::result_type result_type
Definition: Templates.hpp:73
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:89
Do nothing, used to rebuild database.
Definition: Templates.hpp:1536
Extract entity handle form multi-index container.
Definition: Templates.hpp:1491
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:1493
virtual ~TempMeshset()
Definition: Templates.hpp:1567
EntityHandle meshset
Definition: Templates.hpp:1576
moab::Interface & moab
Definition: Templates.hpp:1577
TempMeshset(moab::Interface &moab)
Definition: Templates.hpp:1563