v0.13.2
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 DIM1, int DIM2, class T, class L, class A>
1058
1059template <int S, class T, class L, class A>
1060struct GetFTensor2FromArrayImpl<2, 2, S, T, L, A> {
1062 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1063 const size_t cc) {
1065 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1066
1067 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1068 }
1069};
1070
1071template <int S, class T, class L, class A>
1072struct GetFTensor2FromArrayImpl<3, 3, S, T, L, A> {
1074 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1075 const size_t cc) {
1077 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1078 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1079 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2)};
1080 }
1081};
1082
1083template <class T, class L, class A>
1086 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1087 const size_t cc, const int ss = 0) {
1089 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1090
1091 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
1092 }
1093};
1094
1095template <class T, class L, class A>
1098 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1099 const size_t cc, const int ss = 0) {
1101 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1102 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1103 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
1104 ss);
1105 }
1106};
1107
1108template <int DIM1, int DIM2, int S>
1110getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc = 0) {
1111 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, double, ublas::row_major,
1112 VecAllocator<double>>::get(data, rr, cc);
1113}
1114
1115template <int DIM1, int DIM2>
1117getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc,
1118 const int ss) {
1119 return GetFTensor2FromArrayRawPtrImpl<DIM1, DIM2, double, ublas::row_major,
1120 VecAllocator<double>>::get(data, rr, cc, ss);
1121}
1122
1123template <int S, typename T, typename L, typename A>
1124inline auto getFTensor2FromArray2by2(ublas::matrix<T, L, A> &data,
1125 const FTensor::Number<S> &,
1126 const size_t rr, const size_t cc = 0) {
1128}
1129
1130template <int S, typename T, typename L, typename A>
1131inline auto getFTensor2FromArray3by3(ublas::matrix<T, L, A> &data,
1132 const FTensor::Number<S> &,
1133 const size_t rr, const size_t cc = 0) {
1135}
1136
1137#ifdef WITH_ADOL_C
1138
1139template <int DIM1, int DIM2, int S>
1140inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
1141 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, adouble, ublas::row_major,
1142 VecAllocator<adouble>>::get(data, rr);
1143}
1144
1145#endif
1146
1147// list of lapack wrappers
1148/**
1149 * @brief compute matrix inverse with lapack dgetri
1150 *
1151 * @param mat input square matrix / output inverse matrix
1152 * @return MoFEMErrorCode
1153 */
1156
1157 const size_t M = mat.size1();
1158 const size_t N = mat.size2();
1159
1160 if (M != N)
1161 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1162 "The input matrix for inverse computation is not square %d != %d",
1163 M, N);
1164
1165 int *ipv = new int[N];
1166 int lwork = N * N;
1167 double *work = new double[lwork];
1168 int info;
1169 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1170 if (info != 0)
1171 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1172 "lapack error info = %d", info);
1173 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1174 if (info != 0)
1175 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1176 "lapack error info = %d", info);
1177
1178 delete[] ipv;
1179 delete[] work;
1180
1182}
1183/**
1184 * @brief solve linear system with lapack dgesv
1185 *
1186 * @param mat input lhs square matrix / output L and U from the factorization
1187 * @param f input rhs vector / output solution vector
1188 * @return MoFEMErrorCode
1189 */
1192
1193 const size_t M = mat.size1();
1194 const size_t N = mat.size2();
1195
1196 if (M == 0 || M != N)
1197 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1198 "The input matrix for inverse computation is not square %d != %d",
1199 M, N);
1200 if (f.size() != M)
1201 f.resize(M, false);
1202
1203 const int nrhs = 1;
1204 int info;
1205 int *ipiv = new int[M];
1206 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1207 &*f.data().begin(), nrhs);
1208
1209 if (info != 0) {
1210 SETERRQ1(PETSC_COMM_SELF, 1, "error lapack solve dgesv info = %d", info);
1211 }
1212
1213 delete[] ipiv;
1215}
1216
1217/**
1218 * @brief Solve linear system of equations using Lapack
1219 *
1220 * @param mat
1221 * @param f
1222 * @return MoFEMErrorCode
1223 */
1225 VectorDouble &f) {
1227 // copy matrix since on output lapack returns factorisation
1228 auto mat_copy = mat;
1229 CHKERR solveLinearSystem(mat_copy, f);
1231}
1232
1233/**
1234 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1235 *
1236 * @param mat input symmetric matrix
1237 * @param eig output eigen values sorted
1238 * @param eigen_vec output matrix of row eigen vectors
1239 * @return MoFEMErrorCode
1240 */
1242 VectorDouble &eig,
1243 MatrixDouble &eigen_vec) {
1245
1246 const size_t M = mat.size1();
1247 const size_t N = mat.size2();
1248
1249 if (M == 0 || M != N)
1250 SETERRQ2(
1251 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1252 "The input matrix for eigen value computation is not square %d != %d",
1253 M, N);
1254 if (eig.size() != M)
1255 eig.resize(M, false);
1256
1257 eigen_vec = mat;
1258 const int n = M;
1259 const int lda = M;
1260 const int size = (M + 2) * M;
1261 int lwork = size;
1262 double *work = new double[size];
1263
1264 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
1265 &*eig.data().begin(), work, lwork) > 0)
1266 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1267 "The algorithm failed to compute eigenvalues.");
1268
1269 delete[] work;
1271}
1272/**
1273 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1274 *
1275 * @tparam DIM
1276 * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
1277 * @param eig output eigen values sorted
1278 * @return MoFEMErrorCode
1279 */
1280template <int DIM>
1281inline MoFEMErrorCode
1285
1286 const int n = DIM;
1287 const int lda = DIM;
1288 const int lwork = (DIM + 2) * DIM;
1289 std::array<double, (DIM + 2) * DIM> work;
1290
1291 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
1292 lwork) > 0)
1293 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1294 "The algorithm failed to compute eigenvalues.");
1296}
1297/**
1298 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1299 *
1300 * @tparam DIM
1301 * @param mat input tensor pointer of size DIM x DIM
1302 * @param eig output eigen values sorted
1303 * @param eigen_vec output matrix of row eigen vectors
1304 * @return MoFEMErrorCode
1305 */
1306template <int DIM>
1312 for (int ii = 0; ii != DIM; ii++)
1313 for (int jj = 0; jj != DIM; jj++)
1314 eigen_vec(ii, jj) = mat(ii, jj);
1315
1316 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1317
1319}
1320
1321/**
1322 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1323 *
1324 * @tparam DIM
1325 * @param mat input tensor of size DIM x DIM
1326 * @param eig output eigen values sorted
1327 * @param eigen_vec output matrix of row eigen vectors
1328 * @return MoFEMErrorCode
1329 */
1330template <int DIM>
1331inline MoFEMErrorCode
1336 for (int ii = 0; ii != DIM; ii++)
1337 for (int jj = 0; jj != DIM; jj++)
1338 eigen_vec(ii, jj) = mat(ii, jj);
1339
1340 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1341
1343}
1344
1345/**
1346 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
1347 *
1348 * @tparam T
1349 * @param t
1350 * @return double
1351 */
1352template <typename T> static inline auto determinantTensor3by3(T &t) {
1353 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
1354 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
1355 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
1356}
1357
1358/**
1359 * \brief Calculate inverse of tensor rank 2 at integration points
1360
1361 */
1362template <int Tensor_Dim, class T, class L, class A>
1363inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
1364 ublas::vector<T, A> &det_data,
1365 ublas::matrix<T, L, A> &inv_jac_data) {
1367 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1368 "Specialization for this template not yet implemented");
1370}
1371
1372template <>
1373inline MoFEMErrorCode
1375 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
1376
1377/**
1378 * \brief Calculate determinant 3 by 3
1379
1380 */
1381template <class T1, class T2>
1384 det = determinantTensor3by3(t);
1386}
1387
1388/**
1389 * \brief Calculate determinant 2 by 2
1390
1391 */
1392template <class T1, class T2>
1395 det = t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
1397}
1398
1399/**
1400 * \brief Calculate matrix inverse 3 by 3
1401
1402 */
1403template <class T1, class T2, class T3>
1404inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
1406 const auto inv_det = 1. / det;
1407 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1408 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1409 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1410 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
1411 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1412 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1413 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
1414 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
1415 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1417}
1418
1419/**
1420 * \brief Calculate matrix inverse 2 by 2
1421
1422 */
1423template <class T1, class T2, class T3>
1424inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
1426 const auto inv_det = 1. / det;
1427 inv_t(0, 0) = t(1, 1) * inv_det;
1428 inv_t(0, 1) = -t(0, 1) * inv_det;
1429 inv_t(1, 0) = -t(1, 0) * inv_det;
1430 inv_t(1, 1) = t(0, 0) * inv_det;
1432}
1433
1434#ifdef WITH_ADOL_C
1435
1436/**
1437 * \brief Calculate matrix inverse, specialization for adouble tensor
1438
1439 */
1440template <>
1441inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
1446 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1447 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1448 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1449 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
1450 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1451 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1452 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
1453 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
1454 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1456}
1457
1458#endif
1459
1460/**
1461 * \brief Calculate matrix inverse, specialization for symmetric tensor
1462
1463 */
1464template <>
1465inline MoFEMErrorCode
1466invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1471 const auto inv_det = 1. / det;
1472 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1473 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1474 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1475 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1476 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1477 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1479}
1480
1481#ifdef WITH_ADOL_C
1482
1483/**
1484 * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1485
1486 */
1487template <>
1488inline MoFEMErrorCode
1489invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
1494 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1495 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1496 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1497 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1498 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1499 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1501}
1502
1503#endif
1504
1505/**
1506 * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1507 tensor
1508
1509 */
1510template <>
1511inline MoFEMErrorCode
1512invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1517 const auto inv_det = 1. / det;
1518 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1519 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1520 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1521 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1522 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1523 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1525}
1526
1527/**
1528 * @brief Extract entity handle form multi-index container
1529 *
1530 */
1532 template <typename Iterator>
1533 static inline EntityHandle extract(const Iterator &it) {
1534 return (*it)->getEnt();
1535 }
1536};
1537
1538/**
1539 * @brief Insert ordered mofem multi-index into range
1540 *
1541 * \note Inserted range has to be ordered.
1542 *
1543 * \code
1544 * auto hi_rit = refEntsPtr->upper_bound(start);
1545 * auto hi_rit = refEntsPtr->upper_bound(end);
1546 * Range to_erase;
1547 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1548 * \endcode
1549 *
1550 * @tparam Iterator
1551 * @param r
1552 * @param begin_iter
1553 * @param end_iter
1554 * @return moab::Range::iterator
1555 */
1556template <typename Extractor, typename Iterator>
1557moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1558 Iterator end_iter) {
1559 moab::Range::iterator hint = r.begin();
1560 while (begin_iter != end_iter) {
1561 size_t j = 0;
1562 auto bi = Extractor::extract(begin_iter);
1563 Iterator pj = begin_iter;
1564 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1565 ++pj;
1566 ++j;
1567 }
1568 hint = r.insert(hint, bi, bi + (j - 1));
1569 begin_iter = pj;
1570 }
1571 return hint;
1572};
1573
1574/**
1575 * @brief Do nothing, used to rebuild database
1576 *
1577 */
1580 template <typename T> inline void operator()(T &e) {}
1581};
1582
1583/**
1584 * @brief Template used to reconstruct multi-index
1585 *
1586 * @tparam MI multi-index
1587 * @tparam Modifier
1588 * @param mi
1589 * @param mo
1590 * @return MoFEMErrorCode
1591 */
1592template <typename MI, typename MO = Modify_change_nothing>
1594 MO &&mo = Modify_change_nothing()) {
1596 for (auto it = mi.begin(); it != mi.end(); ++it) {
1597 if (!const_cast<MI &>(mi).modify(it, mo))
1598 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1599 "Houston we have a problem");
1600 }
1602}
1603
1605 TempMeshset(moab::Interface &moab) : moab(moab) {
1606 rval = moab.create_meshset(MESHSET_SET, meshset);
1608 }
1610 operator EntityHandle() const { return meshset; }
1611 auto get_ptr() { return &meshset; }
1612
1613private:
1615 rval = moab.delete_entities(&meshset, 1);
1617 }
1619 moab::Interface &moab;
1620};
1621
1622/**
1623 * @brief Create smart pointer to temporary meshset
1624 *
1625 */
1626inline auto get_temp_meshset_ptr(moab::Interface &moab) {
1627 return boost::make_shared<TempMeshset>(moab);
1628};
1629
1630/**
1631 * @brief get type from entity handle
1632 *
1633 */
1634inline auto type_from_handle(const EntityHandle h) {
1635 return static_cast<EntityType>(h >> MB_ID_WIDTH);
1636};
1637
1638/**
1639 * @brief get entity dimension form handle
1640 *
1641 */
1643 return moab::CN::Dimension(type_from_handle(h));
1644};
1645
1646/**
1647 * @brief get field bit id from bit number
1648 *
1649 */
1650inline auto field_bit_from_bit_number(const int bit_number) {
1651 return BitFieldId().set(bit_number - 1);
1652};
1653
1654/**
1655 * @brief Insert ranges
1656 *
1657 * @tparam I
1658 * @param f
1659 * @param s
1660 * @param tester
1661 * @param inserter
1662 * @return auto
1663 */
1664template <typename I>
1665auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
1666 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
1668
1669 auto first = f;
1670 while (first != s)
1671 if (tester(first)) {
1672
1673 auto second = first;
1674 ++second;
1675
1676 while (second != s) {
1677 if (tester(second))
1678 ++second;
1679 else
1680 break;
1681 }
1682
1683 CHKERR inserter(first, second);
1684
1685 first = second;
1686 if (first != s)
1687 ++first;
1688
1689 } else {
1690 ++first;
1691 }
1692
1694}
1695
1696/**
1697 * @brief Create Array
1698 *
1699 * See:
1700 * <a
1701 * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
1702 * stack overflow</a>
1703 *
1704 * @tparam Dest
1705 * @tparam Arg
1706 * @param arg
1707 * @return constexpr auto
1708 */
1709template <typename Dest = void, typename... Arg>
1710constexpr auto make_array(Arg &&...arg) {
1711 if constexpr (std::is_same<void, Dest>::value)
1712 return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
1713 {std::forward<Arg>(arg)...}};
1714 else
1715 return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
1716}
1717
1718} // namespace MoFEM
1719
1720#endif //__TEMPLATES_HPP__
static Index< 'M', 3 > M
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:1634
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1393
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1110
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:1593
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1557
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1424
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
Definition: Templates.hpp:1154
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
Definition: Templates.hpp:1642
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
Definition: Templates.hpp:1190
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:1363
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 getFTensor2FromArray2by2(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1124
constexpr auto make_array(Arg &&...arg)
Create Array.
Definition: Templates.hpp:1710
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:1650
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
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1131
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1626
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1241
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:1352
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 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:1665
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
constexpr IntegrationType I
constexpr AssemblyType A
double h
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, const size_t cc)
Definition: Templates.hpp:1062
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
Definition: Templates.hpp:1074
Get FTensor2 from array.
Definition: Templates.hpp:1054
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
Definition: Templates.hpp:1086
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
Definition: Templates.hpp:1098
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:1578
Extract entity handle form multi-index container.
Definition: Templates.hpp:1531
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:1533
virtual ~TempMeshset()
Definition: Templates.hpp:1609
EntityHandle meshset
Definition: Templates.hpp:1618
moab::Interface & moab
Definition: Templates.hpp:1619
TempMeshset(moab::Interface &moab)
Definition: Templates.hpp:1605