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() != 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()));
168 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
172template <
int S,
class T,
class A>
174 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
183 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
188template <
int S,
class T,
class A>
190 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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));
204template <
int S,
class T,
class A>
206 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
219template <
int S,
class T,
class A>
221 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
236template <
int Tensor_Dim,
int S = 1,
class T,
class L,
class A>
239 static_assert(!std::is_same<T, T>::value,
"not implemented");
245template <
int Tensor_Dim,
int S = 1>
254template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
257 static_assert(!std::is_same<T, T>::value,
258 "Such getFTensor2FromMat specialisation is not implemented");
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()));
274 std::array<double *, 16> ptrs;
275 for (
auto i = 0;
i != 16; ++
i)
276 ptrs[
i] = &data(
i, 0);
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()));
292 std::array<double *, 36> ptrs;
293 for (
auto i = 0;
i != 36; ++
i)
294 ptrs[
i] = &data(
i, 0);
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()));
310 std::array<double *, 81> ptrs;
311 for (
auto i = 0;
i != 81; ++
i)
312 ptrs[
i] = &data(
i, 0);
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()));
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));
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()));
347 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
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()));
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));
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()));
383 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
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()));
404template <
int Tensor_Dim0,
int Tensor_Dim1>
412template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
415template <
int S,
class T,
class L,
class A>
417 static inline auto get(ublas::matrix<T, L, A> &data) {
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()));
426 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
431template <
int S,
class T,
class L,
class A>
433 static inline auto get(ublas::matrix<T, L, A> &data) {
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()));
442 &data(0, 0), &data(1, 0), &data(2, 0));
449template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
454template <
int Tensor_Dim,
int S = 1>
460template <
int Tensor_Dim01,
int Tensor_Dim23,
int S,
class T,
class L,
class A>
463template <
int S,
class T,
class A>
465 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
477template <
int S,
class T,
class A>
479 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
494template <
int S,
class T,
class A>
496 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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),
530template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1,
class T,
class L,
534 static_assert(!std::is_same<T, T>::value,
535 "Such getFTensor4DdgFromMat specialisation is not implemented");
538template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1>
545template <
int Tensor_Dim01,
int Tensor_Dim2,
int S,
class T,
class L,
class A>
548template <
int S,
class T,
class A>
550 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
562template <
int S,
class T,
class A>
564 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
574 &data(0, 0), &data(1, 0), &data(2, 0),
575 &data(3, 0), &data(4, 0), &data(5, 0)};
579template <
int S,
class T,
class A>
581 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
611template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1,
class T,
class L,
615 static_assert(!std::is_same<T, T>::value,
616 "Such getFTensor3DgFromMat specialisation is not implemented");
619template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1>
625template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
626 int S,
class T,
class L,
class A>
629template <
int S,
class T,
class A>
631 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
644template <
int S,
class T,
class A>
646 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
663template <
int S,
class T,
class A>
665 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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),
710template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
711 int S = 1,
class T,
class L,
class A>
713 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
715 static_assert(!std::is_same<T, T>::value,
716 "Such getFTensor4FromMat specialisation is not implemented");
719template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
723 Tensor_Dim3, S,
double, ublas::row_major,
727template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S,
class T,
731template <
int S,
class T,
class A>
733 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
746template <
int S,
class T,
class A>
748 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)
764template <
int S,
class T,
class A>
766 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
781template <
int S,
class T,
class A>
783 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
798template <
int S,
class T,
class A>
800 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
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)};
818template <
int S,
class T,
class A>
820 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
828 std::array<double *, 54> ptrs;
829 for (
auto i = 0;
i != 54; ++
i)
830 ptrs[
i] = &data(
i, 0);
835template <
int S,
class T,
class A>
837 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
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()));
845 std::array<double *, 54> ptrs;
846 for (
auto i = 0;
i != 54; ++
i)
847 ptrs[
i] = &data(
i, 0);
866template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1,
class T,
869 Tensor_Dim1, Tensor_Dim2>
871 static_assert(!std::is_same<T, T>::value,
872 "Such getFTensor3FromMat specialisation is not implemented");
875template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1>
882template<
int DIM,
int S = DIM>
887 inline static auto get(
double *ptr) {
895 inline static auto get(
double *ptr) {
908template <
int DIM,
int S = DIM>
921template <
int DIM1,
int DIM2>
924 static_assert(DIM1 == DIM1 || DIM2 != DIM2,
925 "Such getFTensor2FromPtr is not implemented");
954 &ptr[0], &ptr[1], &ptr[2], &ptr[3]);
965template <
int DIM1,
int DIM2,
int DIM3>
969 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
970 "Such getFTensor2FromPtr is not implemented");
977 &ptr[0], &ptr[1], &ptr[2],
979 &ptr[3], &ptr[4], &ptr[5],
981 &ptr[6], &ptr[7], &ptr[8],
983 &ptr[9], &ptr[10], &ptr[11]
999 "Such getFTensor2SymmetricUpperFromPtr is not implemented");
1017 &ptr[0], &ptr[1], &ptr[3]);
1024 template <
typename V>
static inline auto get(V &data) {
1025 using T =
typename std::remove_reference<
decltype(data[0])>::type;
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]};
1041 template <
typename V>
static inline auto get(V &data) {
1042 using T =
typename std::remove_reference<
decltype(data[0])>::type;
1044 &data[2], &data[3]};
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]};
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]};
1081template <
int DIM,
int S = 0>
1088template <
int DIM,
int S>
1103 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1110 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0), &data(rr + 3, 0)};
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)};
1133template <
int DIM,
int S>
1136 static_assert(DIM != DIM,
"not implemented");
1151 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1165template <
int DIM1,
int DIM2,
int S,
class T,
class L,
class A>
1168template <
int DIM1,
int DIM2,
class T,
class L,
class A>
1171template <
int S,
class T,
class L,
class A>
1174 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
1177 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1179 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1183template <
int S,
class T,
class L,
class A>
1186 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr,
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)};
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),
1203 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
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),
1220template <
int DIM1,
int DIM2,
int S>
1227template <
int DIM1,
int DIM2>
1235template <
int S,
typename T,
typename L,
typename A>
1238 const size_t rr,
const size_t cc = 0) {
1242template <
int S,
typename T,
typename L,
typename A>
1245 const size_t rr,
const size_t cc = 0) {
1251template <
int DIM1,
int DIM2,
int S>
1269 const size_t M = mat.size1();
1270 const size_t N = mat.size2();
1274 "The input matrix for inverse computation is not square %d != %d",
1277 int *ipv =
new int[
N];
1279 double *work =
new double[lwork];
1284 "lapack error info = %d", info);
1288 "lapack error info = %d", info);
1305 const size_t M = mat.size1();
1306 const size_t N = mat.size2();
1308 if (
M == 0 ||
M !=
N)
1310 "The input matrix for inverse computation is not square %d != %d",
1317 int *ipiv =
new int[
M];
1319 &*f.data().begin(), nrhs);
1322 SETERRQ1(PETSC_COMM_SELF, 1,
"error lapack solve dgesv info = %d", info);
1340 auto mat_copy = mat;
1358 const size_t M = mat.size1();
1359 const size_t N = mat.size2();
1361 if (
M == 0 ||
M !=
N)
1364 "The input matrix for eigen value computation is not square %d != %d",
1366 if (eig.size() !=
M)
1367 eig.resize(
M,
false);
1372 const int size = (
M + 2) *
M;
1374 double *work =
new double[size];
1376 if (
lapack_dsyev(
'V',
'U',
n, &*eigen_vec.data().begin(), lda,
1377 &*eig.data().begin(), work, lwork) > 0)
1379 "The algorithm failed to compute eigenvalues.");
1392template <
int DIM,
typename T1,
typename T2>
1399 const int lda = DIM;
1400 const int lwork = (DIM + 2) * DIM;
1401 std::array<
double, (DIM + 2) * DIM> work;
1406 "The algorithm failed to compute eigenvalues.");
1419template <
int DIM,
typename T1,
typename T2,
typename T3>
1425 for (
int ii = 0; ii != DIM; ii++)
1426 for (
int jj = 0; jj != DIM; jj++)
1427 eigen_vec(ii, jj) = mat(ii, jj);
1441 return computeEigenValuesSymmetric<DIM, double, double>(eigen_vec, eig);
1452 return computeEigenValuesSymmetric<DIM, double, double, double>(mat, eig,
1464 return computeEigenValuesSymmetric<DIM, FTensor::PackPtr<double *, 1>,
double,
1465 double>(mat, eig, eigen_vec);
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);
1485template <
int Tensor_Dim,
class T,
class L,
class A>
1487 ublas::vector<T, A> &det_data,
1488 ublas::matrix<T, L, A> &inv_jac_data) {
1491 "Specialization for this template not yet implemented");
1504template <
class T1,
class T2>
1515template <
class T1,
class T2>
1518 det =
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0);
1526template <
class T1,
class T2,
class T3>
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;
1546template <
class T1,
class T2,
class T3>
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;
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;
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;
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;
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;
1655 template <
typename Iterator>
1657 return (*it)->getEnt();
1679template <
typename Extractor,
typename Iterator>
1681 Iterator end_iter) {
1682 moab::Range::iterator hint = r.begin();
1683 while (begin_iter != end_iter) {
1685 auto bi = Extractor::extract(begin_iter);
1686 Iterator pj = begin_iter;
1687 while (pj != end_iter && (bi +
j) == Extractor::extract(pj)) {
1691 hint = r.insert(hint, bi, bi + (
j - 1));
1715template <
typename MI,
typename MO = Modify_change_nothing>
1719 for (
auto it = mi.begin(); it != mi.end(); ++it) {
1720 if (!
const_cast<MI &
>(mi).modify(it, mo))
1722 "Houston we have a problem");
1750 return boost::make_shared<TempMeshset>(moab);
1807template <
typename I>
1814 if (tester(first)) {
1816 auto second = first;
1819 while (second != s) {
1826 CHKERR inserter(first, second);
1852template <
typename Dest = void,
typename... 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)...}};
1858 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
VectorBoundedArray< double, 3 > VectorDouble3
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 id_from_handle(const EntityHandle h)
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
auto getFTensor1FromArray< 3, 0 >(VectorDouble3 &data)
auto type_name_from_handle(const EntityHandle h)
get entity type name from handle
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 ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
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
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)
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)