5#ifndef __TEMPLATES_HPP__
6#define __TEMPLATES_HPP__
10template <
typename T>
using ShardVec = boost::shared_ptr<std::vector<T>>;
32 typedef typename std::remove_pointer<T1>::type
T;
34 ublas::shallow_array_adaptor<T>(
n, ptr));
58 typedef typename std::remove_pointer<T1>::type
T;
60 n,
m, ublas::shallow_array_adaptor<T>(
n *
m, ptr));
71template <
class KeyExtractor1,
class KeyExtractor2>
struct KeyFromKey {
76 const KeyExtractor2 &key2_ = KeyExtractor2())
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();
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();
102 return value.to_ulong();
106template <
class X>
inline std::string
toString(X x) {
107 std::ostringstream buffer;
113 static inline auto get(ublas::vector<T, A> &data) {
134template <
int S = 1,
class T,
class A>
139template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
142template <
int S,
class T,
class A>
144 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
153 &data(0, 0), &data(1, 0), &data(2, 0));
157template <
int S,
class T,
class A>
159 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
168 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
173template <
int S,
class T,
class A>
175 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
188template <
int S,
class T,
class A>
190 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
205template <
int Tensor_Dim,
int S = 1,
class T,
class L,
class A>
208 static_assert(!std::is_same<T, T>::value,
"not implemented");
214template <
int Tensor_Dim,
int S = 1>
223template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
226 static_assert(!std::is_same<T, T>::value,
227 "Such getFTensor2FromMat specialisation is not implemented");
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()));
242 std::array<double *, 36> ptrs;
243 for (
auto i = 0;
i != 36; ++
i)
244 ptrs[
i] = &data(
i, 0);
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()));
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));
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()));
279 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
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()));
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));
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()));
315 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
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()));
336template <
int Tensor_Dim0,
int Tensor_Dim1>
344template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
347template <
int S,
class T,
class L,
class A>
349 static inline auto get(ublas::matrix<T, L, A> &data) {
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()));
358 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
363template <
int S,
class T,
class L,
class A>
365 static inline auto get(ublas::matrix<T, L, A> &data) {
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()));
374 &data(0, 0), &data(1, 0), &data(2, 0));
381template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
386template <
int Tensor_Dim,
int S = 1>
392template <
int Tensor_Dim01,
int Tensor_Dim23,
int S,
class T,
class L,
class A>
395template <
int S,
class T,
class A>
397 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
409template <
int S,
class T,
class A>
411 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
426template <
int S,
class T,
class A>
428 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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),
462template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1,
class T,
class L,
466 static_assert(!std::is_same<T, T>::value,
467 "Such getFTensor4DdgFromMat specialisation is not implemented");
470template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1>
477template <
int Tensor_Dim01,
int Tensor_Dim2,
int S,
class T,
class L,
class A>
480template <
int S,
class T,
class A>
482 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
494template <
int S,
class T,
class A>
496 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
506 &data(0, 0), &data(1, 0), &data(2, 0),
507 &data(3, 0), &data(4, 0), &data(5, 0)};
511template <
int S,
class T,
class A>
513 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
543template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1,
class T,
class L,
547 static_assert(!std::is_same<T, T>::value,
548 "Such getFTensor3DgFromMat specialisation is not implemented");
551template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1>
557template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
558 int S,
class T,
class L,
class A>
561template <
int S,
class T,
class A>
563 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
576template <
int S,
class T,
class A>
578 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
595template <
int S,
class T,
class A>
597 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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),
642template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
643 int S = 1,
class T,
class L,
class A>
645 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
647 static_assert(!std::is_same<T, T>::value,
648 "Such getFTensor4FromMat specialisation is not implemented");
651template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
655 Tensor_Dim3, S,
double, ublas::row_major,
659template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S,
class T,
663template <
int S,
class T,
class A>
665 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
678template <
int S,
class T,
class A>
680 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)
696template <
int S,
class T,
class A>
698 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
713template <
int S,
class T,
class A>
715 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
730template <
int S,
class T,
class A>
732 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
750template <
int S,
class T,
class A>
752 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
760 std::array<double *, 54> ptrs;
761 for (
auto i = 0;
i != 54; ++
i)
762 ptrs[
i] = &data(
i, 0);
767template <
int S,
class T,
class A>
769 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
777 std::array<double *, 54> ptrs;
778 for (
auto i = 0;
i != 54; ++
i)
779 ptrs[
i] = &data(
i, 0);
798template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1,
class T,
801 Tensor_Dim1, Tensor_Dim2>
803 static_assert(!std::is_same<T, T>::value,
804 "Such getFTensor3FromMat specialisation is not implemented");
807template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1>
814template<
int DIM,
int S = DIM>
819 inline static auto get(
double *ptr) {
827 inline static auto get(
double *ptr) {
840template <
int DIM,
int S = DIM>
853template <
int DIM1,
int DIM2>
856 static_assert(DIM1 == DIM1 || DIM2 != DIM2,
857 "Such getFTensor2FromPtr is not implemented");
886 &ptr[0], &ptr[1], &ptr[2], &ptr[3]);
897template <
int DIM1,
int DIM2,
int DIM3>
901 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
902 "Such getFTensor2FromPtr is not implemented");
909 &ptr[0], &ptr[1], &ptr[2],
911 &ptr[3], &ptr[4], &ptr[5],
913 &ptr[6], &ptr[7], &ptr[8],
915 &ptr[9], &ptr[10], &ptr[11]
931 "Such getFTensor2SymmetricUpperFromPtr is not implemented");
949 &ptr[0], &ptr[1], &ptr[3]);
956 template <
typename V>
static inline auto get(V &data) {
957 using T =
typename std::remove_reference<
decltype(data[0])>::type;
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]};
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]};
993template <
int DIM,
int S>
1008 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1021template <
int DIM,
int S>
1024 static_assert(DIM != DIM,
"not implemented");
1039 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1053template <
int DIM1,
int DIM2,
int S,
class T,
class L,
class A>
1056template <
int DIM1,
int DIM2,
class T,
class L,
class A>
1059template <
int S,
class T,
class L,
class A>
1062 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1065 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1067 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1071template <
int S,
class T,
class L,
class A>
1074 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
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)};
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),
1091 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
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),
1108template <
int DIM1,
int DIM2,
int S>
1115template <
int DIM1,
int DIM2>
1123template <
int S,
typename T,
typename L,
typename A>
1126 const size_t rr,
const size_t cc = 0) {
1130template <
int S,
typename T,
typename L,
typename A>
1133 const size_t rr,
const size_t cc = 0) {
1139template <
int DIM1,
int DIM2,
int S>
1157 const size_t M = mat.size1();
1158 const size_t N = mat.size2();
1162 "The input matrix for inverse computation is not square %d != %d",
1165 int *ipv =
new int[
N];
1167 double *work =
new double[lwork];
1172 "lapack error info = %d", info);
1176 "lapack error info = %d", info);
1193 const size_t M = mat.size1();
1194 const size_t N = mat.size2();
1196 if (
M == 0 ||
M !=
N)
1198 "The input matrix for inverse computation is not square %d != %d",
1205 int *ipiv =
new int[
M];
1207 &*f.data().begin(), nrhs);
1210 SETERRQ1(PETSC_COMM_SELF, 1,
"error lapack solve dgesv info = %d", info);
1228 auto mat_copy = mat;
1246 const size_t M = mat.size1();
1247 const size_t N = mat.size2();
1249 if (
M == 0 ||
M !=
N)
1252 "The input matrix for eigen value computation is not square %d != %d",
1254 if (eig.size() !=
M)
1255 eig.resize(
M,
false);
1260 const int size = (
M + 2) *
M;
1262 double *work =
new double[size];
1264 if (
lapack_dsyev(
'V',
'U',
n, &*eigen_vec.data().begin(), lda,
1265 &*eig.data().begin(), work, lwork) > 0)
1267 "The algorithm failed to compute eigenvalues.");
1287 const int lda = DIM;
1288 const int lwork = (DIM + 2) * DIM;
1289 std::array<
double, (DIM + 2) * DIM> work;
1294 "The algorithm failed to compute eigenvalues.");
1312 for (
int ii = 0; ii != DIM; ii++)
1313 for (
int jj = 0; jj != DIM; jj++)
1314 eigen_vec(ii, jj) = mat(ii, jj);
1316 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1336 for (
int ii = 0; ii != DIM; ii++)
1337 for (
int jj = 0; jj != DIM; jj++)
1338 eigen_vec(ii, jj) = mat(ii, jj);
1340 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
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);
1362template <
int Tensor_Dim,
class T,
class L,
class A>
1364 ublas::vector<T, A> &det_data,
1365 ublas::matrix<T, L, A> &inv_jac_data) {
1368 "Specialization for this template not yet implemented");
1381template <
class T1,
class T2>
1392template <
class T1,
class T2>
1395 det =
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0);
1403template <
class T1,
class T2,
class T3>
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;
1423template <
class T1,
class T2,
class T3>
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;
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;
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;
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;
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;
1532 template <
typename Iterator>
1534 return (*it)->getEnt();
1556template <
typename Extractor,
typename Iterator>
1558 Iterator end_iter) {
1559 moab::Range::iterator hint = r.begin();
1560 while (begin_iter != end_iter) {
1562 auto bi = Extractor::extract(begin_iter);
1563 Iterator pj = begin_iter;
1564 while (pj != end_iter && (bi +
j) == Extractor::extract(pj)) {
1568 hint = r.insert(hint, bi, bi + (
j - 1));
1592template <
typename MI,
typename MO = Modify_change_nothing>
1596 for (
auto it = mi.begin(); it != mi.end(); ++it) {
1597 if (!
const_cast<MI &
>(mi).modify(it, mo))
1599 "Houston we have a problem");
1627 return boost::make_shared<TempMeshset>(moab);
1664template <
typename I>
1671 if (tester(first)) {
1673 auto second = first;
1676 while (second != s) {
1683 CHKERR inserter(first, second);
1709template <
typename Dest = void,
typename... 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)...}};
1715 return std::array<Dest,
sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
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)
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
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)
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
FTensor::Index< 'j', 3 > j
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::vector< T, std::allocator< T > > VecAllocator
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
UBlasMatrix< double > MatrixDouble
UBlasMatrix< adouble > MatrixADouble
VecAllocator< double > DoubleAllocator
UBlasVector< double > VectorDouble
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
implementation of Data Operators for Forces and Sources
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.
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
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.
auto type_from_handle(const EntityHandle h)
get type from entity handle
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc=0)
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.
FTensor::Tensor3< FTensor::PackPtr< double *, 12 >, 3, 2, 2 > getFTensor3FromPtr< 3, 2, 2 >(double *ptr)
FTensor::Tensor2< FTensor::PackPtr< double *, 4 >, 2, 2 > getFTensor2FromPtr< 2, 2 >(double *ptr)
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
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.
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.
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)
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
auto getFTensor2FromArray2by2(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
constexpr auto make_array(Arg &&...arg)
Create Array.
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
FTensor::Tensor2< FTensor::PackPtr< double *, DIM1 *DIM2 >, DIM1, DIM2 > getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
std::string toString(X x)
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.
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr)
Get FTensor1 from array.
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, DIM *DIM >, DIM > getFTensor2SymmetricLowerFromPtr(double *ptr)
Make symmetric Tensor2 from pointer, taking lower triangle of matrix.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
FTensor::Tensor3< FTensor::PackPtr< double *, DIM1 *DIM2 *DIM3 >, DIM1, DIM2, DIM3 > getFTensor3FromPtr(double *ptr)
boost::shared_ptr< std::vector< T > > ShardVec
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 9 >, 3 > getFTensor2SymmetricLowerFromPtr< 3 >(double *ptr)
constexpr IntegrationType I
constexpr double t
plate stiffness
bool operator()(const id_type &valueA, const id_type &valueB) const
static auto get(ublas::vector< T, A > &data)
GetFTensor1FromArray()=delete
GetFTensor1FromArray()=delete
GetFTensor1FromArray()=delete
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
GetFTensor1FromPtrImpl()=delete
static auto get(double *ptr)
GetFTensor1FromPtrImpl()=delete
static auto get(double *ptr)
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
GetFTensor2FromArrayImpl()=delete
GetFTensor2FromArrayImpl()=delete
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
GetFTensor2FromArrayRawPtrImpl()=delete
GetFTensor2FromArrayRawPtrImpl()=delete
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
static auto get(ublas::matrix< T, L, A > &data)
static auto get(ublas::matrix< T, L, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
unsigned int operator()(const id_type &value) const
result_type operator()(Arg &arg) const
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
KeyExtractor1::result_type result_type
bool operator()(const id_type &valueA, const id_type &valueB) const
Do nothing, used to rebuild database.
Modify_change_nothing()=default
TempMeshset(moab::Interface &moab)