v0.14.0
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<4, 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() != 4)
163 "getFTensor1FromMat<4>: wrong size of data matrix, number of "
164 "rows should be 4 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));
169 }
170};
171
172template <int S, class T, class A>
173struct GetFTensor1FromMatImpl<6, S, T, ublas::row_major, A> {
174 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
175#ifndef NDDEBUG
176 if (data.size1() != 6)
178 "getFTensor1FromMat<6>: wrong size of data matrix, number of "
179 "rows should be 6 but is " +
180 boost::lexical_cast<std::string>(data.size1()));
181#endif
183 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
184 &data(5, 0));
185 }
186};
187
188template <int S, class T, class A>
189struct GetFTensor1FromMatImpl<9, S, T, ublas::row_major, A> {
190 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
191#ifndef NDDEBUG
192 if (data.size1() != 9)
194 "getFTensor1FromMat<9>: wrong size of data matrix, number of "
195 "rows should be 9 but is " +
196 boost::lexical_cast<std::string>(data.size1()));
197#endif
199 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
200 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
201 }
202};
203
204template <int S, class T, class A>
205struct GetFTensor1FromMatImpl<2, S, T, ublas::row_major, A> {
206 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
207#ifndef NDDEBUG
208 if (data.size1() != 2)
210 "getFTensor1FromMat<2>: wrong size of data matrix, number of "
211 "rows should be 2 but is " +
212 boost::lexical_cast<std::string>(data.size1()));
213#endif
215 &data(1, 0));
216 }
217};
218
219template <int S, class T, class A>
220struct GetFTensor1FromMatImpl<1, S, T, ublas::row_major, A> {
221 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
222#ifndef NDEBUG
223 if (data.size1() != 1)
225 "getFTensor1FromMat<1>: wrong size of data matrix, number of "
226 "rows should be 1 but is " +
227 boost::lexical_cast<std::string>(data.size1()));
228#endif
230 }
231};
232
233/**
234 * \brief Get tensor rank 1 (vector) form data matrix
235 */
236template <int Tensor_Dim, int S = 1, class T, class L, class A>
238getFTensor1FromMat(ublas::matrix<T, L, A> &data) {
239 static_assert(!std::is_same<T, T>::value, "not implemented");
240}
241
242/**
243 * \brief Get tensor rank 1 (vector) form data matrix (specialization)
244 */
245template <int Tensor_Dim, int S = 1>
247 return GetFTensor1FromMatImpl<Tensor_Dim, S, double, ublas::row_major,
248 DoubleAllocator>::get(data);
249}
250
251/**
252 * \brief Get tensor rank 2 (matrix) form data matrix
253 */
254template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
255inline FTensor::Tensor2<FTensor::PackPtr<T *, 1>, Tensor_Dim0, Tensor_Dim1>
256getFTensor2FromMat(ublas::matrix<T, L, A> &data) {
257 static_assert(!std::is_same<T, T>::value,
258 "Such getFTensor2FromMat specialisation is not implemented");
259}
260
261
262/**
263 * Template specialization for getFTensor2FromMat
264 */
265template <>
268#ifndef NDEBUG
269 if (data.size1() != 16)
270 THROW_MESSAGE("getFTensor2FromMat<4, 4>: wrong size of data matrix, numer "
271 "of rows should be 16 but is " +
272 boost::lexical_cast<std::string>(data.size1()));
273#endif
274 std::array<double *, 16> ptrs;
275 for (auto i = 0; i != 16; ++i)
276 ptrs[i] = &data(i, 0);
278}
279
280/**
281 * Template specialization for getFTensor2FromMat
282 */
283template <>
286#ifndef NDEBUG
287 if (data.size1() != 36)
288 THROW_MESSAGE("getFTensor2FromMat<6, 6>: wrong size of data matrix, numer "
289 "of rows should be 36 but is " +
290 boost::lexical_cast<std::string>(data.size1()));
291#endif
292 std::array<double *, 36> ptrs;
293 for (auto i = 0; i != 36; ++i)
294 ptrs[i] = &data(i, 0);
296}
297
298/**
299 * Template specialization for getFTensor2FromMat
300 */
301template <>
304#ifndef NDEBUG
305 if (data.size1() != 81)
306 THROW_MESSAGE("getFTensor2FromMat<9, 9>: wrong size of data matrix, numer "
307 "of rows should be 81 but is " +
308 boost::lexical_cast<std::string>(data.size1()));
309#endif
310 std::array<double *, 81> ptrs;
311 for (auto i = 0; i != 81; ++i)
312 ptrs[i] = &data(i, 0);
314}
315
316/**
317 * Template specialization for getFTensor2FromMat
318 *
319 */
320template <>
323#ifndef NDEBUG
324 if (data.size1() != 9)
325 THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix; numer "
326 "of rows should be 9 but is " +
327 boost::lexical_cast<std::string>(data.size1()));
328#endif
330 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
331 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
332}
333
334/**
335 * Template specialization for getFTensor2FromMat
336 */
337template <>
340#ifndef NDEBUG
341 if (data.size1() != 6)
342 THROW_MESSAGE("getFTensor2FromMat<3,2>: wrong size of data matrix, numer "
343 "of rows should be 6 but is " +
344 boost::lexical_cast<std::string>(data.size1()));
345#endif
347 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
348 &data(5, 0));
349}
350
351/**
352 * Template specialization for getFTensor2FromMat
353 */
354template <>
357#ifndef NDEBUG
358 if (data.size1() != 18)
359 THROW_MESSAGE("getFTensor2FromMat<6,3>: wrong size of data matrix, numer "
360 "of rows should be 18 but is " +
361 boost::lexical_cast<std::string>(data.size1()));
362#endif
364 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
365 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
366 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
367 &data(15, 0), &data(16, 0), &data(17, 0));
368}
369
370/**
371 * Template specialization for getFTensor2FromMat
372 */
373template <>
376#ifndef NDEBUG
377 if (data.size1() != 4)
378 THROW_MESSAGE("getFTensor2FromMat<2,2>: wrong size of data matrix, numer "
379 "of rows should be 4 but is " +
380 boost::lexical_cast<std::string>(data.size1()));
381#endif
383 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
384}
385
386/**
387 * Template specialization for getFTensor2FromMat
388 */
389template <>
392#ifndef NDEBUG
393 if (data.size1() != 1)
394 THROW_MESSAGE("getFTensor2FromMat<1,1>: wrong size of data matrix, numer "
395 "of rows should be 1 but is " +
396 boost::lexical_cast<std::string>(data.size1()));
397#endif
398 return FTensor::Tensor2<FTensor::PackPtr<double *, 1>, 1, 1>(&data(0, 0));
399}
400
401/**
402 * \brief Get tensor rank 2 (matrix) form data matrix (specialization)
403 */
404template <int Tensor_Dim0, int Tensor_Dim1>
405static inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim0,
406 Tensor_Dim1>
408 return getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
409 DoubleAllocator>(data);
410}
411
412template <int Tensor_Dim, int S, class T, class L, class A>
414
415template <int S, class T, class L, class A>
417 static inline auto get(ublas::matrix<T, L, A> &data) {
418#ifndef NDEBUG
419 if (data.size1() != 6)
421 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
422 "of rows should be 6 but is " +
423 boost::lexical_cast<std::string>(data.size1()));
424#endif
426 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
427 &data(5, 0));
428 }
429};
430
431template <int S, class T, class L, class A>
433 static inline auto get(ublas::matrix<T, L, A> &data) {
434#ifndef NDEBUG
435 if (data.size1() != 3)
437 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
438 "of rows should be 3 but is " +
439 boost::lexical_cast<std::string>(data.size1()));
440#endif
442 &data(0, 0), &data(1, 0), &data(2, 0));
443 }
444};
445
446/**
447 * \brief Get symmetric tensor rank 2 (matrix) form data matrix
448 */
449template <int Tensor_Dim, int S, class T, class L, class A>
450static inline auto getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
452}
453
454template <int Tensor_Dim, int S = 1>
455static inline auto getFTensor2SymmetricFromMat(MatrixDouble &data) {
456 return getFTensor2SymmetricFromMat<Tensor_Dim, S, double, ublas::row_major,
457 DoubleAllocator>(data);
458}
459
460template <int Tensor_Dim01, int Tensor_Dim23, int S, class T, class L, class A>
462
463template <int S, class T, class A>
464struct GetFTensor4DdgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
465 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
466#ifndef NDEBUG
467 if (data.size1() != 1)
469 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
470 "of rows should be 1 but is " +
471 boost::lexical_cast<std::string>(data.size1()));
472#endif
473 return FTensor::Ddg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
474 }
475};
476
477template <int S, class T, class A>
478struct GetFTensor4DdgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
479 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
480#ifndef NDEBUG
481 if (data.size1() != 9) {
483 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
484 "of rows should be 9 but is " +
485 boost::lexical_cast<std::string>(data.size1()));
486 }
487#endif
489 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
490 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
491 }
492};
493
494template <int S, class T, class A>
495struct GetFTensor4DdgFromMatImpl<3, 3, 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() != 36) {
499 cerr << data.size1() << endl;
501 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
502 "of rows should be 36 but is " +
503 boost::lexical_cast<std::string>(data.size1()));
504 }
505#endif
507 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
508 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
509 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
510 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
511 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
512 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
513 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
514 &data(35, 0)};
515 }
516};
517
518/**
519 * @brief Get symmetric tensor rank 4 on first two and last indices from
520 * form data matrix
521 *
522 * @tparam Tensor_Dim01 dimension of first two indicies
523 * @tparam Tensor_Dim23 dimension of second two indicies
524 * @tparam T the type of object stored
525 * @tparam L the storage organization
526 * @tparam A the type of Storage array
527 * @param data data container
528 * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
529 */
530template <int Tensor_Dim01, int Tensor_Dim23, int S = 1, class T, class L,
531 class A>
532static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
533getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
534 static_assert(!std::is_same<T, T>::value,
535 "Such getFTensor4DdgFromMat specialisation is not implemented");
536}
537
538template <int Tensor_Dim01, int Tensor_Dim23, int S = 1>
539static inline auto getFTensor4DdgFromMat(MatrixDouble &data) {
540 return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, double,
541 ublas::row_major,
542 DoubleAllocator>::get(data);
543}
544
545template <int Tensor_Dim01, int Tensor_Dim2, int S, class T, class L, class A>
547
548template <int S, class T, class A>
549struct GetFTensor3DgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
550 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
551#ifndef NDEBUG
552 if (data.size1() != 1)
554 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
555 "of rows should be 1 but is " +
556 boost::lexical_cast<std::string>(data.size1()));
557#endif
558 return FTensor::Dg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
559 }
560};
561
562template <int S, class T, class A>
563struct GetFTensor3DgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
564 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
565#ifndef NDEBUG
566 if (data.size1() != 6) {
568 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
569 "of rows should be 6 but is " +
570 boost::lexical_cast<std::string>(data.size1()));
571 }
572#endif
574 &data(0, 0), &data(1, 0), &data(2, 0),
575 &data(3, 0), &data(4, 0), &data(5, 0)};
576 }
577};
578
579template <int S, class T, class A>
580struct GetFTensor3DgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
581 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
582#ifndef NDEBUG
583 if (data.size1() != 18) {
584 cerr << data.size1() << endl;
586 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
587 "of rows should be 18 but is " +
588 boost::lexical_cast<std::string>(data.size1()));
589 }
590#endif
592 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
593 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
594 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
595 &data(15, 0), &data(16, 0), &data(17, 0)};
596 }
597};
598
599/**
600 * @brief Get symmetric tensor rank 3 on the first two indices from
601 * form data matrix
602 *
603 * @tparam Tensor_Dim01 dimension of first two indicies
604 * @tparam Tensor_Dim2 dimension of last index
605 * @tparam T the type of object stored
606 * @tparam L the storage organization
607 * @tparam A the type of Storage array
608 * @param data data container
609 * @return FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
610 */
611template <int Tensor_Dim01, int Tensor_Dim2, int S = 1, class T, class L,
612 class A>
613static inline FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim2>
614getFTensor3DgFromMat(ublas::matrix<T, L, A> &data) {
615 static_assert(!std::is_same<T, T>::value,
616 "Such getFTensor3DgFromMat specialisation is not implemented");
617}
618
619template <int Tensor_Dim01, int Tensor_Dim2, int S = 1>
620static inline auto getFTensor3DgFromMat(MatrixDouble &data) {
621 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, double,
622 ublas::row_major, DoubleAllocator>::get(data);
623}
624
625template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
626 int S, class T, class L, class A>
628
629template <int S, class T, class A>
630struct GetFTensor4FromMatImpl<1, 1, 1, 1, S, T, ublas::row_major, A> {
631 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
632#ifndef NDEBUG
633 if (data.size1() != 1)
635 "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
636 "of rows should be 1 but is " +
637 boost::lexical_cast<std::string>(data.size1()));
638#endif
640 &data(0, 0)};
641 }
642};
643
644template <int S, class T, class A>
645struct GetFTensor4FromMatImpl<2, 2, 2, 2, S, T, ublas::row_major, A> {
646 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
647#ifndef NDEBUG
648 if (data.size1() != 16) {
650 "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
651 "of rows should be 16 but is " +
652 boost::lexical_cast<std::string>(data.size1()));
653 }
654#endif
656 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
657 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
658 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
659 &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
660 }
661};
662
663template <int S, class T, class A>
664struct GetFTensor4FromMatImpl<3, 3, 3, 3, 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() != 81) {
668 cerr << data.size1() << endl;
670 "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
671 "of rows should be 81 but is " +
672 boost::lexical_cast<std::string>(data.size1()));
673 }
674#endif
676 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
677 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
678 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
679 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
680 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
681 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
682 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
683 &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
684 &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
685 &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
686 &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
687 &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
688 &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
689 &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
690 &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
691 &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
692 &data(80, 0)};
693 }
694};
695
696/**
697 * @brief Get tensor rank 4 (non symmetric) form data matrix
698 *
699 * @tparam Tensor_Dim0 dimension of frirst index
700 * @tparam Tensor_Dim1 dimension of second index
701 * @tparam Tensor_Dim2 dimension of third index
702 * @tparam Tensor_Dim3 dimension of fourth index
703 * @tparam T the type of object stored
704 * @tparam L the storage organization
705 * @tparam A the type of Storage array
706 * @param data data container
707 * @return FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
708 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
709 */
710template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
711 int S = 1, class T, class L, class A>
712static inline FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
713 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
714getFTensor4FromMat(ublas::matrix<T, L, A> &data) {
715 static_assert(!std::is_same<T, T>::value,
716 "Such getFTensor4FromMat specialisation is not implemented");
717}
718
719template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
720 int S = 1>
721static inline auto getFTensor4FromMat(MatrixDouble &data) {
722 return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
723 Tensor_Dim3, S, double, ublas::row_major,
724 DoubleAllocator>::get(data);
725}
726
727template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class T,
728 class L, class A>
730
731template <int S, class T, class A>
732struct GetFTensor3FromMatImpl<1, 1, 1, S, T, ublas::row_major, A> {
733 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
734#ifndef NDEBUG
735 if (data.size1() != 1)
737 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
738 "of rows should be 1 but is " +
739 boost::lexical_cast<std::string>(data.size1()));
740#endif
742 &data(0, 0)};
743 }
744};
745
746template <int S, class T, class A>
747struct GetFTensor3FromMatImpl<2, 2, 2, S, T, ublas::row_major, A> {
748 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
749#ifndef NDEBUG
750 if (data.size1() != 8)
752 "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
753 "of rows should be 8 but is " +
754 boost::lexical_cast<std::string>(data.size1()));
755#endif
757 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
758 &data(5, 0), &data(6, 0), &data(7, 0)
759
760 };
761 }
762};
763
764template <int S, class T, class A>
765struct GetFTensor3FromMatImpl<3, 2, 2, S, T, ublas::row_major, A> {
766 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
767#ifndef NDEBUG
768 if (data.size1() != 12)
770 "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
771 "of rows should be 12 but is " +
772 boost::lexical_cast<std::string>(data.size1()));
773#endif
775 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
776 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
777 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
778 }
779};
780
781template <int S, class T, class A>
782struct GetFTensor3FromMatImpl<2, 2, 3, S, T, ublas::row_major, A> {
783 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
784#ifndef NDEBUG
785 if (data.size1() != 12)
787 "getFTensor3FromMat<2, 2, 3>: wrong size of data matrix, number "
788 "of rows should be 12 but is " +
789 boost::lexical_cast<std::string>(data.size1()));
790#endif
792 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
793 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
794 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
795 }
796};
797
798template <int S, class T, class A>
799struct GetFTensor3FromMatImpl<3, 3, 3, S, T, ublas::row_major, A> {
800 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
801#ifndef NDEBUG
802 if (data.size1() != 27)
804 "getFTensor3FromMat<3, 3, 3>: wrong size of data matrix, number "
805 "of rows should be 27 but is " +
806 boost::lexical_cast<std::string>(data.size1()));
807#endif
809 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
810 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
811 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
812 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
813 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
814 &data(25, 0), &data(26, 0)};
815 }
816};
817
818template <int S, class T, class A>
819struct GetFTensor3FromMatImpl<6, 3, 3, S, T, ublas::row_major, A> {
820 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
821#ifndef NDEBUG
822 if (data.size1() != 54)
824 "getFTensor3FromMat<6, 3, 3>: wrong size of data matrix, number "
825 "of rows should be 54 but is " +
826 boost::lexical_cast<std::string>(data.size1()));
827#endif
828 std::array<double *, 54> ptrs;
829 for (auto i = 0; i != 54; ++i)
830 ptrs[i] = &data(i, 0);
832 }
833};
834
835template <int S, class T, class A>
836struct GetFTensor3FromMatImpl<3, 3, 6, S, T, ublas::row_major, A> {
837 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
838#ifndef NDEBUG
839 if (data.size1() != 54)
841 "getFTensor3FromMat<3, 3, 6>: wrong size of data matrix, number "
842 "of rows should be 54 but is " +
843 boost::lexical_cast<std::string>(data.size1()));
844#endif
845 std::array<double *, 54> ptrs;
846 for (auto i = 0; i != 54; ++i)
847 ptrs[i] = &data(i, 0);
849 }
850};
851
852/**
853 * @brief Get tensor rank 3 (non symmetries) form data matrix
854 *
855 * @tparam Tensor_Dim0 dimension of first index
856 * @tparam Tensor_Dim1 dimension of second index
857 * @tparam Tensor_Dim2 dimension of third index
858 * @tparam S shift size
859 * @tparam T the type of object stored
860 * @tparam L the storage organization
861 * @tparam A the type of Storage array
862 * @param data data container
863 * @return FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
864 Tensor_Dim1, Tensor_Dim2>
865 */
866template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1, class T,
867 class L, class A>
868static inline FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
869 Tensor_Dim1, Tensor_Dim2>
870getFTensor3FromMat(ublas::matrix<T, L, A> &data) {
871 static_assert(!std::is_same<T, T>::value,
872 "Such getFTensor3FromMat specialisation is not implemented");
873}
874
875template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1>
876static inline auto getFTensor3FromMat(MatrixDouble &data) {
877 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
878 double, ublas::row_major,
879 DoubleAllocator>::get(data);
880}
881
882template<int DIM, int S = DIM>
884
885template <int S> struct GetFTensor1FromPtrImpl<2, S> {
887 inline static auto get(double *ptr) {
889 &ptr[HVEC1]);
890 }
891};
892
893template <int S> struct GetFTensor1FromPtrImpl<3, S> {
895 inline static auto get(double *ptr) {
897 &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
898 }
899};
900
901/**
902 * @brief Make Tensor1 from pointer
903 *
904 * @tparam DIM
905 * @param ptr
906 * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
907 */
908template <int DIM, int S = DIM>
910getFTensor1FromPtr(double *ptr) {
912};
913
914/**
915 * @brief Make Tensor2 from pointer
916 *
917 * @tparam DIM
918 * @param ptr
919 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
920 */
921template <int DIM1, int DIM2>
923getFTensor2FromPtr(double *ptr) {
924 static_assert(DIM1 == DIM1 || DIM2 != DIM2,
925 "Such getFTensor2FromPtr is not implemented");
926};
927
928template <>
930 2> inline getFTensor2FromPtr<3, 2>(double *ptr) {
932 &ptr[HVEC0_0], &ptr[HVEC0_1],
933
934 &ptr[HVEC1_0], &ptr[HVEC1_1],
935
936 &ptr[HVEC2_0], &ptr[HVEC2_1]);
937};
938
939template <>
941 3> inline getFTensor2FromPtr<3, 3>(double *ptr) {
943 &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
944
945 &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
946
947 &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
948};
949
950template <>
952 2> inline getFTensor2FromPtr<2, 2>(double *ptr) {
954 &ptr[0], &ptr[1], &ptr[2], &ptr[3]);
955};
956
957/*
958 * @brief Make Tensor3 from pointer
959 *
960 * @tparam DIM
961 * @param ptr
962 * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3>, DIM1,
963 * DIM2, DIM3>
964 */
965template <int DIM1, int DIM2, int DIM3>
967 DIM2, DIM3>
968getFTensor3FromPtr(double *ptr) {
969 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
970 "Such getFTensor2FromPtr is not implemented");
971};
972
973template <>
975getFTensor3FromPtr<3, 2, 2>(double *ptr) {
977 &ptr[0], &ptr[1], &ptr[2],
978
979 &ptr[3], &ptr[4], &ptr[5],
980
981 &ptr[6], &ptr[7], &ptr[8],
982
983 &ptr[9], &ptr[10], &ptr[11]
984
985 );
986};
987
988/**
989 * @brief Make symmetric Tensor2 from pointer, taking lower triangle of matrix
990 *
991 * @tparam DIM
992 * @param ptr
993 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
994 */
995template <int DIM>
998 static_assert(DIM,
999 "Such getFTensor2SymmetricUpperFromPtr is not implemented");
1000}
1001
1002template <>
1006 &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
1007
1008 &ptr[HVEC1_0], &ptr[HVEC1_1],
1009
1010 &ptr[HVEC2_2]);
1011};
1012
1013template <>
1017 &ptr[0], &ptr[1], &ptr[3]);
1018};
1019
1020template <int DIM, int S> struct GetFTensor1FromArray;
1021
1022template <int S> struct GetFTensor1FromArray<2, S> {
1024 template <typename V> static inline auto get(V &data) {
1025 using T = typename std::remove_reference<decltype(data[0])>::type;
1026 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 2>{&data[0], &data[1]};
1027 }
1028};
1029
1030template <int S> struct GetFTensor1FromArray<3, S> {
1032 template <typename V> static inline auto get(V &data) {
1033 using T = typename std::remove_reference<decltype(data[0])>::type;
1035 &data[0], &data[1], &data[2]};
1036 }
1037};
1038
1039template <int S> struct GetFTensor1FromArray<4, S> {
1041 template <typename V> static inline auto get(V &data) {
1042 using T = typename std::remove_reference<decltype(data[0])>::type;
1043 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 4>{&data[0], &data[1],
1044 &data[2], &data[3]};
1045 }
1046};
1047
1048template <int S> struct GetFTensor1FromArray<6, S> {
1050 template <typename V> static inline auto get(V &data) {
1051 using T = typename std::remove_reference<decltype(data[0])>::type;
1053 &data[0], &data[1], &data[2], &data[3], &data[4], &data[5]};
1054 }
1055};
1056
1057template <int S> struct GetFTensor1FromArray<9, S> {
1059 template <typename V> static inline auto get(V &data) {
1060 using T = typename std::remove_reference<decltype(data[0])>::type;
1062 &data[0], &data[1], &data[2], &data[3], &data[4],
1063 &data[5], &data[6], &data[7], &data[8]};
1064 }
1065};
1066
1067/**
1068 * @brief Get FTensor1 from array
1069 *
1070 * \todo Generalise for different arrays and data types
1071 *
1072 * @tparam DIM
1073 * @param data
1074 * @return FTensor::Tensor1<FTensor::PackPtr<double *, S>, DIM>
1075 */
1076template <int DIM, int S> inline auto getFTensor1FromArray(VectorDouble &data) {
1078}
1079
1080/** @copydoc getFTensor1FromArray */
1081template <int DIM, int S = 0>
1083
1084template <> inline auto getFTensor1FromArray<3, 0>(VectorDouble3 &data) {
1086}
1087
1088template <int DIM, int S>
1090getFTensor1FromMat(MatrixDouble &data, const size_t rr);
1091
1092template <>
1094getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1095 return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2>{&data(rr + 0, 0),
1096 &data(rr + 1, 0)};
1097}
1098
1099template <>
1101getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1103 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1104}
1105
1106template <>
1108getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1110 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0), &data(rr + 3, 0)};
1111}
1112
1113template <>
1115getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1117 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0),
1118 &data(rr + 3, 0), &data(rr + 4, 0), &data(rr + 5, 0),
1119 &data(rr + 6, 0), &data(rr + 7, 0), &data(rr + 8, 0)};
1120}
1121
1122
1123/**
1124 * @brief Get FTensor1 from array
1125 *
1126 * \todo Generalise for diffrent arrays and data types
1127 *
1128 * @tparam DIM
1129 * @param data
1130 * @param rr
1131 * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
1132 */
1133template <int DIM, int S>
1136 static_assert(DIM != DIM, "not implemented");
1138}
1139
1140template <>
1143 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
1144 &data(rr + 1, 1)};
1145}
1146
1147template <>
1151 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1152}
1153
1154/**
1155 * @brief Get FTensor2 from array
1156 *
1157 * \note Generalise for other data types
1158 *
1159 * @tparam DIM1
1160 * @tparam DIM2
1161 * @tparam S
1162 * @param data
1163 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
1164 */
1165template <int DIM1, int DIM2, int S, class T, class L, class A>
1167
1168template <int DIM1, int DIM2, class T, class L, class A>
1170
1171template <int S, class T, class L, class A>
1172struct GetFTensor2FromArrayImpl<2, 2, S, T, L, A> {
1174 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1175 const size_t cc) {
1177 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1178
1179 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1180 }
1181};
1182
1183template <int S, class T, class L, class A>
1184struct GetFTensor2FromArrayImpl<3, 3, S, T, L, A> {
1186 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1187 const size_t cc) {
1189 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1190 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1191 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2)};
1192 }
1193};
1194
1195template <class T, class L, class A>
1198 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1199 const size_t cc, const int ss = 0) {
1201 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1202
1203 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
1204 }
1205};
1206
1207template <class T, class L, class A>
1210 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1211 const size_t cc, const int ss = 0) {
1213 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1214 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1215 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
1216 ss);
1217 }
1218};
1219
1220template <int DIM1, int DIM2, int S>
1222getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc = 0) {
1223 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, double, ublas::row_major,
1224 VecAllocator<double>>::get(data, rr, cc);
1225}
1226
1227template <int DIM1, int DIM2>
1229getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc,
1230 const int ss) {
1231 return GetFTensor2FromArrayRawPtrImpl<DIM1, DIM2, double, ublas::row_major,
1232 VecAllocator<double>>::get(data, rr, cc, ss);
1233}
1234
1235template <int S, typename T, typename L, typename A>
1236inline auto getFTensor2FromArray2by2(ublas::matrix<T, L, A> &data,
1237 const FTensor::Number<S> &,
1238 const size_t rr, const size_t cc = 0) {
1240}
1241
1242template <int S, typename T, typename L, typename A>
1243inline auto getFTensor2FromArray3by3(ublas::matrix<T, L, A> &data,
1244 const FTensor::Number<S> &,
1245 const size_t rr, const size_t cc = 0) {
1247}
1248
1249#ifdef WITH_ADOL_C
1250
1251template <int DIM1, int DIM2, int S>
1252inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
1253 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, adouble, ublas::row_major,
1254 VecAllocator<adouble>>::get(data, rr);
1255}
1256
1257#endif
1258
1259// list of lapack wrappers
1260/**
1261 * @brief compute matrix inverse with lapack dgetri
1262 *
1263 * @param mat input square matrix / output inverse matrix
1264 * @return MoFEMErrorCode
1265 */
1268
1269 const size_t M = mat.size1();
1270 const size_t N = mat.size2();
1271
1272 if (M != N)
1273 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1274 "The input matrix for inverse computation is not square %d != %d",
1275 M, N);
1276
1277 int *ipv = new int[N];
1278 int lwork = N * N;
1279 double *work = new double[lwork];
1280 int info;
1281 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1282 if (info != 0)
1283 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1284 "lapack error info = %d", info);
1285 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1286 if (info != 0)
1287 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1288 "lapack error info = %d", info);
1289
1290 delete[] ipv;
1291 delete[] work;
1292
1294}
1295/**
1296 * @brief solve linear system with lapack dgesv
1297 *
1298 * @param mat input lhs square matrix / output L and U from the factorization
1299 * @param f input rhs vector / output solution vector
1300 * @return MoFEMErrorCode
1301 */
1304
1305 const size_t M = mat.size1();
1306 const size_t N = mat.size2();
1307
1308 if (M == 0 || M != N)
1309 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1310 "The input matrix for inverse computation is not square %d != %d",
1311 M, N);
1312 if (f.size() != M)
1313 f.resize(M, false);
1314
1315 const int nrhs = 1;
1316 int info;
1317 int *ipiv = new int[M];
1318 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1319 &*f.data().begin(), nrhs);
1320
1321 if (info != 0) {
1322 SETERRQ1(PETSC_COMM_SELF, 1, "error lapack solve dgesv info = %d", info);
1323 }
1324
1325 delete[] ipiv;
1327}
1328
1329/**
1330 * @brief Solve linear system of equations using Lapack
1331 *
1332 * @param mat
1333 * @param f
1334 * @return MoFEMErrorCode
1335 */
1337 VectorDouble &f) {
1339 // copy matrix since on output lapack returns factorisation
1340 auto mat_copy = mat;
1341 CHKERR solveLinearSystem(mat_copy, f);
1343}
1344
1345/**
1346 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1347 *
1348 * @param mat input symmetric matrix
1349 * @param eig output eigen values sorted
1350 * @param eigen_vec output matrix of row eigen vectors
1351 * @return MoFEMErrorCode
1352 */
1354 VectorDouble &eig,
1355 MatrixDouble &eigen_vec) {
1357
1358 const size_t M = mat.size1();
1359 const size_t N = mat.size2();
1360
1361 if (M == 0 || M != N)
1362 SETERRQ2(
1363 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1364 "The input matrix for eigen value computation is not square %d != %d",
1365 M, N);
1366 if (eig.size() != M)
1367 eig.resize(M, false);
1368
1369 eigen_vec = mat;
1370 const int n = M;
1371 const int lda = M;
1372 const int size = (M + 2) * M;
1373 int lwork = size;
1374 double *work = new double[size];
1375
1376 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
1377 &*eig.data().begin(), work, lwork) > 0)
1378 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1379 "The algorithm failed to compute eigenvalues.");
1380
1381 delete[] work;
1383}
1384/**
1385 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1386 *
1387 * @tparam DIM
1388 * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
1389 * @param eig output eigen values sorted
1390 * @return MoFEMErrorCode
1391 */
1392template <int DIM, typename T1, typename T2>
1393inline MoFEMErrorCode
1397
1398 const int n = DIM;
1399 const int lda = DIM;
1400 const int lwork = (DIM + 2) * DIM;
1401 std::array<double, (DIM + 2) * DIM> work;
1402
1403 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
1404 lwork) > 0)
1405 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1406 "The algorithm failed to compute eigenvalues.");
1408}
1409
1410/**
1411 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1412 *
1413 * @tparam DIM
1414 * @param mat input tensor pointer of size DIM x DIM
1415 * @param eig output eigen values sorted
1416 * @param eigen_vec output matrix of row eigen vectors
1417 * @return MoFEMErrorCode
1418 */
1419template <int DIM, typename T1, typename T2, typename T3>
1423 FTensor::Tensor2<T3, DIM, DIM> &eigen_vec) {
1425 for (int ii = 0; ii != DIM; ii++)
1426 for (int jj = 0; jj != DIM; jj++)
1427 eigen_vec(ii, jj) = mat(ii, jj);
1428
1429 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
1430
1432}
1433
1434/**
1435 * @deprecated do not use, is kept for backward compatibility
1436*/
1437template <int DIM>
1441 return computeEigenValuesSymmetric<DIM, double, double>(eigen_vec, eig);
1442}
1443
1444/**
1445 * @deprecated do not use, is kept for backward compatibility
1446*/
1447template <int DIM>
1452 return computeEigenValuesSymmetric<DIM, double, double, double>(mat, eig,
1453 eigen_vec);
1454}
1455
1456/**
1457 * @deprecated do not use, is kept for backward compatibility
1458*/
1459template <int DIM>
1464 return computeEigenValuesSymmetric<DIM, FTensor::PackPtr<double *, 1>, double,
1465 double>(mat, eig, eigen_vec);
1466}
1467
1468/**
1469 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
1470 *
1471 * @tparam T
1472 * @param t
1473 * @return double
1474 */
1475template <typename T> static inline auto determinantTensor3by3(T &t) {
1476 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
1477 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
1478 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
1479}
1480
1481/**
1482 * \brief Calculate inverse of tensor rank 2 at integration points
1483
1484 */
1485template <int Tensor_Dim, class T, class L, class A>
1486inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
1487 ublas::vector<T, A> &det_data,
1488 ublas::matrix<T, L, A> &inv_jac_data) {
1490 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1491 "Specialization for this template not yet implemented");
1493}
1494
1495template <>
1496inline MoFEMErrorCode
1498 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
1499
1500/**
1501 * \brief Calculate determinant 3 by 3
1502
1503 */
1504template <class T1, class T2>
1507 det = determinantTensor3by3(t);
1509}
1510
1511/**
1512 * \brief Calculate determinant 2 by 2
1513
1514 */
1515template <class T1, class T2>
1518 det = t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
1520}
1521
1522/**
1523 * \brief Calculate matrix inverse 3 by 3
1524
1525 */
1526template <class T1, class T2, class T3>
1527inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
1529 const auto inv_det = 1. / det;
1530 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1531 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1532 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1533 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
1534 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1535 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1536 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
1537 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
1538 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1540}
1541
1542/**
1543 * \brief Calculate matrix inverse 2 by 2
1544
1545 */
1546template <class T1, class T2, class T3>
1547inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
1549 const auto inv_det = 1. / det;
1550 inv_t(0, 0) = t(1, 1) * inv_det;
1551 inv_t(0, 1) = -t(0, 1) * inv_det;
1552 inv_t(1, 0) = -t(1, 0) * inv_det;
1553 inv_t(1, 1) = t(0, 0) * inv_det;
1555}
1556
1557#ifdef WITH_ADOL_C
1558
1559/**
1560 * \brief Calculate matrix inverse, specialization for adouble tensor
1561
1562 */
1563template <>
1564inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
1569 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1570 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1571 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1572 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
1573 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1574 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1575 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
1576 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
1577 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1579}
1580
1581#endif
1582
1583/**
1584 * \brief Calculate matrix inverse, specialization for symmetric tensor
1585
1586 */
1587template <>
1588inline MoFEMErrorCode
1589invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1594 const auto inv_det = 1. / det;
1595 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1596 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1597 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1598 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1599 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1600 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1602}
1603
1604#ifdef WITH_ADOL_C
1605
1606/**
1607 * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1608
1609 */
1610template <>
1611inline MoFEMErrorCode
1612invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
1617 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1618 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1619 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1620 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1621 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1622 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1624}
1625
1626#endif
1627
1628/**
1629 * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1630 tensor
1631
1632 */
1633template <>
1634inline MoFEMErrorCode
1635invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1640 const auto inv_det = 1. / det;
1641 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1642 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1643 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1644 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1645 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1646 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1648}
1649
1650/**
1651 * @brief Extract entity handle form multi-index container
1652 *
1653 */
1655 template <typename Iterator>
1656 static inline EntityHandle extract(const Iterator &it) {
1657 return (*it)->getEnt();
1658 }
1659};
1660
1661/**
1662 * @brief Insert ordered mofem multi-index into range
1663 *
1664 * \note Inserted range has to be ordered.
1665 *
1666 * \code
1667 * auto hi_rit = refEntsPtr->upper_bound(start);
1668 * auto hi_rit = refEntsPtr->upper_bound(end);
1669 * Range to_erase;
1670 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1671 * \endcode
1672 *
1673 * @tparam Iterator
1674 * @param r
1675 * @param begin_iter
1676 * @param end_iter
1677 * @return moab::Range::iterator
1678 */
1679template <typename Extractor, typename Iterator>
1680moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1681 Iterator end_iter) {
1682 moab::Range::iterator hint = r.begin();
1683 while (begin_iter != end_iter) {
1684 size_t j = 0;
1685 auto bi = Extractor::extract(begin_iter);
1686 Iterator pj = begin_iter;
1687 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1688 ++pj;
1689 ++j;
1690 }
1691 hint = r.insert(hint, bi, bi + (j - 1));
1692 begin_iter = pj;
1693 }
1694 return hint;
1695};
1696
1697/**
1698 * @brief Do nothing, used to rebuild database
1699 *
1700 */
1703 template <typename T> inline void operator()(T &e) {}
1704};
1705
1706/**
1707 * @brief Template used to reconstruct multi-index
1708 *
1709 * @tparam MI multi-index
1710 * @tparam Modifier
1711 * @param mi
1712 * @param mo
1713 * @return MoFEMErrorCode
1714 */
1715template <typename MI, typename MO = Modify_change_nothing>
1717 MO &&mo = Modify_change_nothing()) {
1719 for (auto it = mi.begin(); it != mi.end(); ++it) {
1720 if (!const_cast<MI &>(mi).modify(it, mo))
1721 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1722 "Houston we have a problem");
1723 }
1725}
1726
1728 TempMeshset(moab::Interface &moab) : moab(moab) {
1729 rval = moab.create_meshset(MESHSET_SET, meshset);
1731 }
1733 operator EntityHandle() const { return meshset; }
1734 auto get_ptr() { return &meshset; }
1735
1736private:
1738 rval = moab.delete_entities(&meshset, 1);
1740 }
1742 moab::Interface &moab;
1743};
1744
1745/**
1746 * @brief Create smart pointer to temporary meshset
1747 *
1748 */
1749inline auto get_temp_meshset_ptr(moab::Interface &moab) {
1750 return boost::make_shared<TempMeshset>(moab);
1751};
1752
1753inline auto id_from_handle(const EntityHandle h) {
1754 return static_cast<EntityID>(h & MB_ID_MASK);
1755};
1756
1757/**
1758 * @brief get type from entity handle
1759 *
1760 */
1761inline auto type_from_handle(const EntityHandle h) {
1762 return static_cast<EntityType>(h >> MB_ID_WIDTH);
1763};
1764
1765/**
1766 * @brief get entity handle from type and id
1767 *
1768 */
1769inline auto ent_form_type_and_id(const EntityType type, const EntityID id) {
1770 return (static_cast<EntityHandle>(type) << MB_ID_WIDTH) | id;
1771};
1772
1773/**
1774 * @brief get entity dimension form handle
1775 *
1776 */
1778 return moab::CN::Dimension(type_from_handle(h));
1779};
1780
1781/**
1782 * @brief get entity type name from handle
1783 *
1784 */
1786 return moab::CN::EntityTypeName(type_from_handle(h));
1787};
1788
1789/**
1790 * @brief get field bit id from bit number
1791 *
1792 */
1793inline auto field_bit_from_bit_number(const int bit_number) {
1794 return BitFieldId().set(bit_number - 1);
1795};
1796
1797/**
1798 * @brief Insert ranges
1799 *
1800 * @tparam I
1801 * @param f
1802 * @param s
1803 * @param tester
1804 * @param inserter
1805 * @return auto
1806 */
1807template <typename I>
1808auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
1809 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
1811
1812 auto first = f;
1813 while (first != s)
1814 if (tester(first)) {
1815
1816 auto second = first;
1817 ++second;
1818
1819 while (second != s) {
1820 if (tester(second))
1821 ++second;
1822 else
1823 break;
1824 }
1825
1826 CHKERR inserter(first, second);
1827
1828 first = second;
1829 if (first != s)
1830 ++first;
1831
1832 } else {
1833 ++first;
1834 }
1835
1837}
1838
1839/**
1840 * @brief Create Array
1841 *
1842 * See:
1843 * <a
1844 * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
1845 * stack overflow</a>
1846 *
1847 * @tparam Dest
1848 * @tparam Arg
1849 * @param arg
1850 * @return constexpr auto
1851 */
1852template <typename Dest = void, typename... Arg>
1853constexpr auto make_array(Arg &&...arg) {
1854 if constexpr (std::is_same<void, Dest>::value)
1855 return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
1856 {std::forward<Arg>(arg)...}};
1857 else
1858 return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
1859}
1860
1861} // namespace MoFEM
1862
1863#endif //__TEMPLATES_HPP__
static Index< 'M', 3 > M
T data[Tensor_Dim]
#define MB_ID_MASK
Definition: definitions.h:234
#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 DEPRECATED
Definition: definitions.h:17
#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
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
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:870
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:614
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1761
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1516
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1222
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:256
FTensor::Tensor3< FTensor::PackPtr< double *, 12 >, 3, 2, 2 > getFTensor3FromPtr< 3, 2, 2 >(double *ptr)
Definition: Templates.hpp:975
FTensor::Tensor2< FTensor::PackPtr< double *, 4 >, 2, 2 > getFTensor2FromPtr< 2, 2 >(double *ptr)
Definition: Templates.hpp:952
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1716
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1680
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1547
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
Definition: Templates.hpp:1266
auto id_from_handle(const EntityHandle h)
Definition: Templates.hpp:1753
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
Definition: Templates.hpp:1777
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
Definition: Templates.hpp:1302
auto getFTensor1FromArray< 3, 0 >(VectorDouble3 &data)
Definition: Templates.hpp:1084
auto type_name_from_handle(const EntityHandle h)
get entity type name from handle
Definition: Templates.hpp:1785
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:1486
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:533
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:1015
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:238
auto getFTensor2FromArray2by2(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1236
constexpr auto make_array(Arg &&...arg)
Create Array.
Definition: Templates.hpp:1853
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:923
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:450
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition: Templates.hpp:910
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
Definition: Templates.hpp:1793
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:714
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:941
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr)
Get FTensor1 from array.
Definition: Templates.hpp:1135
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1243
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1749
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1353
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:997
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1475
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:930
FTensor::Tensor3< FTensor::PackPtr< double *, DIM1 *DIM2 *DIM3 >, DIM1, DIM2, DIM3 > getFTensor3FromPtr(double *ptr)
Definition: Templates.hpp:968
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:1808
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
Definition: Templates.hpp:1769
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
Definition: Templates.hpp:1076
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 9 >, 3 > getFTensor2SymmetricLowerFromPtr< 3 >(double *ptr)
Definition: Templates.hpp:1004
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:221
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:206
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(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:174
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:190
static auto get(double *ptr)
Definition: Templates.hpp:887
static auto get(double *ptr)
Definition: Templates.hpp:895
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
Definition: Templates.hpp:1174
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
Definition: Templates.hpp:1186
Get FTensor2 from array.
Definition: Templates.hpp:1166
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
Definition: Templates.hpp:1198
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
Definition: Templates.hpp:1210
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:433
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:417
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:550
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:564
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:581
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:733
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:748
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:783
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:766
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:800
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:837
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:820
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:465
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:479
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:631
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:646
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:665
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:1701
Extract entity handle form multi-index container.
Definition: Templates.hpp:1654
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:1656
virtual ~TempMeshset()
Definition: Templates.hpp:1732
EntityHandle meshset
Definition: Templates.hpp:1741
moab::Interface & moab
Definition: Templates.hpp:1742
TempMeshset(moab::Interface &moab)
Definition: Templates.hpp:1728