19#ifndef __TEMPLATES_HPP__
20#define __TEMPLATES_HPP__
24template <
typename T>
using ShardVec = boost::shared_ptr<std::vector<T>>;
48 ublas::shallow_array_adaptor<T>(
n, ptr));
74 n,
m, ublas::shallow_array_adaptor<T>(
n *
m, ptr));
85template <
class KeyExtractor1,
class KeyExtractor2>
struct KeyFromKey {
90 const KeyExtractor2 &key2_ = KeyExtractor2())
102template <
typename id_type>
struct LtBit {
103 inline bool operator()(
const id_type &valueA,
const id_type &valueB)
const {
104 return valueA.to_ulong() < valueB.to_ulong();
108template <
typename id_type>
struct EqBit {
109 inline bool operator()(
const id_type &valueA,
const id_type &valueB)
const {
110 return valueA.to_ulong() == valueB.to_ulong();
116 return value.to_ulong();
120template <
class X>
inline std::string
toString(X x) {
121 std::ostringstream buffer;
127 static inline auto get(ublas::vector<T, A> &data) {
148template <
int S = 1,
class T,
class A>
153template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
156template <
int S,
class T,
class A>
158 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
160 if (data.size1() != 3)
162 "getFTensor1FromMat<3>: wrong size of data matrix, number of "
163 "rows should be 3 but is %d" +
164 boost::lexical_cast<std::string>(data.size1()));
167 &data(0, 0), &data(1, 0), &data(2, 0));
171template <
int S,
class T,
class A>
173 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
175 if (data.size1() != 2)
177 "getFTensor1FromMat<2>: wrong size of data matrix, number of "
178 "rows should be 2 but is %d" +
179 boost::lexical_cast<std::string>(data.size1()));
186template <
int S,
class T,
class A>
188 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
190 if (data.size1() != 1)
192 "getFTensor1FromMat<1>: wrong size of data matrix, number of "
193 "rows should be 1 but is %d" +
194 boost::lexical_cast<std::string>(data.size1()));
203template <
int Tensor_Dim,
int S = 1,
class T,
class L,
class A>
206 static_assert(!std::is_same<T, T>::value,
"not implemented");
212template <
int Tensor_Dim,
int S = 1>
221template <
int Tensor_Dim0,
int Tensor_Dim1,
class T,
class L,
class A>
224 static_assert(!std::is_same<T, T>::value,
225 "Such getFTensor2FromMat specialisation is not implemented");
236 if (data.size1() != 9)
237 THROW_MESSAGE(
"getFTensor2FromMat<3,3>: wrong size of data matrix; numer "
238 "of rows should be 9 but is " +
239 boost::lexical_cast<std::string>(data.size1()));
242 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
243 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
253 if (data.size1() != 6)
254 THROW_MESSAGE(
"getFTensor2FromMat<3,3>: wrong size of data matrix, numer "
255 "of rows should be 6 but is " +
256 boost::lexical_cast<std::string>(data.size1()));
259 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
270 if (data.size1() != 4)
271 THROW_MESSAGE(
"getFTensor2FromMat<2,2>: wrong size of data matrix, numer "
272 "of rows should be 4 but is " +
273 boost::lexical_cast<std::string>(data.size1()));
276 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
286 if (data.size1() != 1)
287 THROW_MESSAGE(
"getFTensor2FromMat<1,1>: wrong size of data matrix, numer "
288 "of rows should be 1 but is " +
289 boost::lexical_cast<std::string>(data.size1()));
297template <
int Tensor_Dim0,
int Tensor_Dim1>
305template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
308template <
int S,
class T,
class L,
class A>
310 static inline auto get(ublas::matrix<T, L, A> &data) {
312 if (data.size1() != 6)
314 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
315 "of rows should be 6 but is " +
316 boost::lexical_cast<std::string>(data.size1()));
319 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
324template <
int S,
class T,
class L,
class A>
326 static inline auto get(ublas::matrix<T, L, A> &data) {
328 if (data.size1() != 3)
330 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
331 "of rows should be 3 but is " +
332 boost::lexical_cast<std::string>(data.size1()));
335 &data(0, 0), &data(1, 0), &data(2, 0));
342template <
int Tensor_Dim,
int S,
class T,
class L,
class A>
347template <
int Tensor_Dim,
int S = 1>
353template <
int Tensor_Dim01,
int Tensor_Dim23,
int S,
class T,
class L,
class A>
356template <
int S,
class T,
class A>
358 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
360 if (data.size1() != 1)
362 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
363 "of rows should be 1 but is " +
364 boost::lexical_cast<std::string>(data.size1()));
370template <
int S,
class T,
class A>
372 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
374 if (data.size1() != 9) {
376 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
377 "of rows should be 9 but is " +
378 boost::lexical_cast<std::string>(data.size1()));
382 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
383 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
387template <
int S,
class T,
class A>
389 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
391 if (data.size1() != 36) {
392 cerr << data.size1() << endl;
394 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
395 "of rows should be 36 but is " +
396 boost::lexical_cast<std::string>(data.size1()));
400 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
401 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
402 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
403 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
404 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
405 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
406 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
423template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1,
class T,
class L,
427 static_assert(!std::is_same<T, T>::value,
428 "Such getFTensor4DdgFromMat specialisation is not implemented");
431template <
int Tensor_Dim01,
int Tensor_Dim23,
int S = 1>
438template <
int Tensor_Dim01,
int Tensor_Dim2,
int S,
class T,
class L,
class A>
441template <
int S,
class T,
class A>
443 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
445 if (data.size1() != 1)
447 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
448 "of rows should be 1 but is " +
449 boost::lexical_cast<std::string>(data.size1()));
455template <
int S,
class T,
class A>
457 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
459 if (data.size1() != 6) {
461 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
462 "of rows should be 6 but is " +
463 boost::lexical_cast<std::string>(data.size1()));
467 &data(0, 0), &data(1, 0), &data(2, 0),
468 &data(3, 0), &data(4, 0), &data(5, 0)};
472template <
int S,
class T,
class A>
474 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
476 if (data.size1() != 18) {
477 cerr << data.size1() << endl;
479 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
480 "of rows should be 18 but is " +
481 boost::lexical_cast<std::string>(data.size1()));
485 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
486 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
487 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
488 &data(15, 0), &data(16, 0), &data(17, 0)};
504template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1,
class T,
class L,
508 static_assert(!std::is_same<T, T>::value,
509 "Such getFTensor3DgFromMat specialisation is not implemented");
512template <
int Tensor_Dim01,
int Tensor_Dim2,
int S = 1>
518template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
519 int S,
class T,
class L,
class A>
522template <
int S,
class T,
class A>
524 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
526 if (data.size1() != 1)
528 "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
529 "of rows should be 1 but is " +
530 boost::lexical_cast<std::string>(data.size1()));
537template <
int S,
class T,
class A>
539 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
541 if (data.size1() != 16) {
543 "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
544 "of rows should be 16 but is " +
545 boost::lexical_cast<std::string>(data.size1()));
549 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
550 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
551 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
552 &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
556template <
int S,
class T,
class A>
558 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
560 if (data.size1() != 81) {
561 cerr << data.size1() << endl;
563 "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
564 "of rows should be 81 but is " +
565 boost::lexical_cast<std::string>(data.size1()));
569 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
570 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
571 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
572 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
573 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
574 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
575 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
576 &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
577 &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
578 &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
579 &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
580 &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
581 &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
582 &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
583 &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
584 &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
603template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
604 int S = 1,
class T,
class L,
class A>
606 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
608 static_assert(!std::is_same<T, T>::value,
609 "Such getFTensor4FromMat specialisation is not implemented");
612template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int Tensor_Dim3,
620template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S,
class T,
624template <
int S,
class T,
class A>
626 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
628 if (data.size1() != 1)
630 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
631 "of rows should be 1 but is " +
632 boost::lexical_cast<std::string>(data.size1()));
639template <
int S,
class T,
class A>
641 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
643 if (data.size1() != 8)
645 "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
646 "of rows should be 8 but is " +
647 boost::lexical_cast<std::string>(data.size1()));
650 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
651 &data(5, 0), &data(6, 0), &data(7, 0)
657template <
int S,
class T,
class A>
659 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
661 if (data.size1() != 12)
663 "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
664 "of rows should be 12 but is " +
665 boost::lexical_cast<std::string>(data.size1()));
668 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
669 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
670 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
674template <
int S,
class T,
class A>
676 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
678 if (data.size1() != 26)
680 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
681 "of rows should be 8 but is " +
682 boost::lexical_cast<std::string>(data.size1()));
685 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
686 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
687 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
688 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
689 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
690 &data(25, 0), &data(26, 0)};
708template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1,
class T,
711 Tensor_Dim1, Tensor_Dim2>
713 static_assert(!std::is_same<T, T>::value,
714 "Such getFTensor3FromMat specialisation is not implemented");
717template <
int Tensor_Dim0,
int Tensor_Dim1,
int Tensor_Dim2,
int S = 1>
734 static_assert(DIM != 3 && DIM != 2,
735 "Such getFTensor1FromPtr specialization is not implemented");
759template <
int DIM1,
int DIM2>
762 static_assert(DIM1 == DIM1 || DIM2 != DIM2,
763 "Such getFTensor2FromPtr is not implemented");
792 &ptr[0], &ptr[1], &ptr[2], &ptr[3]);
803template <
int DIM1,
int DIM2,
int DIM3>
807 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
808 "Such getFTensor2FromPtr is not implemented");
815 &ptr[0], &ptr[1], &ptr[2],
817 &ptr[3], &ptr[4], &ptr[5],
819 &ptr[6], &ptr[7], &ptr[8],
821 &ptr[9], &ptr[10], &ptr[11]
837 "Such getFTensor2SymmetricUpperFromPtr is not implemented");
855 &ptr[0], &ptr[1], &ptr[3]);
867template <
int DIM,
int S>
870 static_assert(DIM != DIM,
"not implemented");
887template <
int DIM,
int S>
890 static_assert(DIM != DIM,
"not implemented");
905 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
918template <
int DIM,
int S>
921 static_assert(DIM != DIM,
"not implemented");
936 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
950template <
int DIM1,
int DIM2,
int S,
class T,
class L,
class A>
953template <
int S,
class T,
class L,
class A>
955 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr) {
957 &data(rr + 0, 0), &data(rr + 0, 1), &data(rr + 1, 0), &data(rr + 1, 1)};
961template <
int S,
class T,
class L,
class A>
963 inline static auto get(ublas::matrix<T, L, A> &data,
const size_t rr) {
965 &data(rr + 0, 0), &data(rr + 0, 1), &data(rr + 0, 2),
966 &data(rr + 1, 0), &data(rr + 1, 1), &data(rr + 1, 2),
967 &data(rr + 2, 0), &data(rr + 2, 1), &data(rr + 2, 2)};
971template <
int DIM1,
int DIM2,
int S>
978template <
int S,
typename T,
typename L,
typename A>
985template <
int S,
typename T,
typename L,
typename A>
994template <
int DIM1,
int DIM2,
int S>
1012 const size_t M = mat.size1();
1013 const size_t N = mat.size2();
1017 "The input matrix for inverse computation is not square %d != %d",
1020 int *ipv =
new int[
N];
1022 double *work =
new double[lwork];
1027 "lapack error info = %d", info);
1031 "lapack error info = %d", info);
1048 const size_t M = mat.size1();
1049 const size_t N = mat.size2();
1051 if (
M == 0 ||
M !=
N)
1053 "The input matrix for inverse computation is not square %d != %d",
1060 int *ipiv =
new int[
M];
1062 &*
f.data().begin(), nrhs);
1065 SETERRQ1(PETSC_COMM_SELF, 1,
"error lapack solve dgesv info = %d", info);
1083 auto mat_copy = mat;
1101 const size_t M = mat.size1();
1102 const size_t N = mat.size2();
1104 if (
M == 0 ||
M !=
N)
1107 "The input matrix for eigen value computation is not square %d != %d",
1109 if (eig.size() !=
M)
1110 eig.resize(
M,
false);
1115 const int size = (
M + 2) *
M;
1117 double *work =
new double[size];
1119 if (
lapack_dsyev(
'V',
'U',
n, &*eigen_vec.data().begin(), lda,
1120 &*eig.data().begin(), work, lwork) > 0)
1122 "The algorithm failed to compute eigenvalues.");
1142 const int lda = DIM;
1143 const int lwork = (DIM + 2) * DIM;
1144 std::array<
double, (DIM + 2) * DIM> work;
1149 "The algorithm failed to compute eigenvalues.");
1167 for (
int ii = 0; ii != DIM; ii++)
1168 for (
int jj = 0; jj != DIM; jj++)
1169 eigen_vec(ii, jj) = mat(ii, jj);
1171 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1191 for (
int ii = 0; ii != DIM; ii++)
1192 for (
int jj = 0; jj != DIM; jj++)
1193 eigen_vec(ii, jj) = mat(ii, jj);
1195 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1208 return t(0, 0) *
t(1, 1) *
t(2, 2) +
t(1, 0) *
t(2, 1) *
t(0, 2) +
1209 t(2, 0) *
t(0, 1) *
t(1, 2) -
t(0, 0) *
t(2, 1) *
t(1, 2) -
1210 t(2, 0) *
t(1, 1) *
t(0, 2) -
t(1, 0) *
t(0, 1) *
t(2, 2);
1217template <
int Tensor_Dim,
class T,
class L,
class A>
1219 ublas::vector<T, A> &det_data,
1220 ublas::matrix<T, L, A> &inv_jac_data) {
1223 "Specialization for this template not yet implemented");
1229invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
1236template <
class T1,
class T2>
1247template <
class T1,
class T2>
1250 det =
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0);
1258template <
class T1,
class T2,
class T3>
1261 const auto inv_det = 1. / det;
1262 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1263 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1264 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1265 inv_t(1, 0) = (
t(1, 2) *
t(2, 0) -
t(1, 0) *
t(2, 2)) * inv_det;
1266 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1267 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1268 inv_t(2, 0) = (
t(1, 0) *
t(2, 1) -
t(1, 1) *
t(2, 0)) * inv_det;
1269 inv_t(2, 1) = (
t(0, 1) *
t(2, 0) -
t(0, 0) *
t(2, 1)) * inv_det;
1270 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1278template <
class T1,
class T2,
class T3>
1281 const auto inv_det = 1. / det;
1282 inv_t(0, 0) =
t(1, 1) * inv_det;
1283 inv_t(0, 1) = -
t(0, 1) * inv_det;
1284 inv_t(1, 0) = -
t(1, 0) * inv_det;
1285 inv_t(1, 1) =
t(0, 0) * inv_det;
1301 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) / det;
1302 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) / det;
1303 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) / det;
1304 inv_t(1, 0) = (
t(1, 2) *
t(2, 0) -
t(1, 0) *
t(2, 2)) / det;
1305 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) / det;
1306 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) / det;
1307 inv_t(2, 0) = (
t(1, 0) *
t(2, 1) -
t(1, 1) *
t(2, 0)) / det;
1308 inv_t(2, 1) = (
t(0, 1) *
t(2, 0) -
t(0, 0) *
t(2, 1)) / det;
1309 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) / det;
1321invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>,
double,
1326 const auto inv_det = 1. / det;
1327 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1328 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1329 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1330 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1331 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1332 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1344invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>,
adouble,
1349 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) / det;
1350 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) / det;
1351 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) / det;
1352 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) / det;
1353 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) / det;
1354 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) / det;
1367invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>,
double,
1372 const auto inv_det = 1. / det;
1373 inv_t(0, 0) = (
t(1, 1) *
t(2, 2) -
t(1, 2) *
t(2, 1)) * inv_det;
1374 inv_t(0, 1) = (
t(0, 2) *
t(2, 1) -
t(0, 1) *
t(2, 2)) * inv_det;
1375 inv_t(0, 2) = (
t(0, 1) *
t(1, 2) -
t(0, 2) *
t(1, 1)) * inv_det;
1376 inv_t(1, 1) = (
t(0, 0) *
t(2, 2) -
t(0, 2) *
t(2, 0)) * inv_det;
1377 inv_t(1, 2) = (
t(0, 2) *
t(1, 0) -
t(0, 0) *
t(1, 2)) * inv_det;
1378 inv_t(2, 2) = (
t(0, 0) *
t(1, 1) -
t(0, 1) *
t(1, 0)) * inv_det;
1387 template <
typename Iterator>
1388 static inline EntityHandle
extract(
const Iterator &it) {
1389 return (*it)->getEnt();
1409template <
typename Extractor,
typename Iterator>
1411 Iterator end_iter) {
1412 moab::Range::iterator hint =
r.begin();
1413 while (begin_iter != end_iter) {
1415 auto bi = Extractor::extract(begin_iter);
1416 Iterator pj = begin_iter;
1417 while (pj != end_iter && (bi +
j) == Extractor::extract(pj)) {
1421 hint =
r.insert(hint, bi, bi + (
j - 1));
1445template <
typename MI,
typename MO = Modify_change_nothing>
1449 for (
auto it = mi.begin(); it != mi.end(); ++it) {
1450 if (!
const_cast<MI &
>(mi).modify(it, mo))
1452 "Houston we have a problem");
1480 return boost::make_shared<TempMeshset>(moab);
1517template <
typename I>
1524 if (tester(first)) {
1526 auto second = first;
1529 while (second != s) {
1536 CHKERR inserter(first, second);
#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.
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< 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.
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
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
FTensor::Tensor1< FTensor::PackPtr< double *, DIM >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
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.
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 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.
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.
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temprary meshset.
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
DeprecatedCoreInterface Interface
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.
FTensor::Tensor1< FTensor::PackPtr< double *, 2 >, 2 > getFTensor1FromPtr< 2 >(double *ptr)
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
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 getFTensor2FromArray2by2(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 9 >, 3 > getFTensor2SymmetricLowerFromPtr< 3 >(double *ptr)
const double r
rate factor
double h
convective heat coefficient
constexpr double t
plate stiffness
FTensor::Index< 'm', 3 > m
bool operator()(const id_type &valueA, const id_type &valueB) const
static auto get(ublas::vector< T, 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, L, A > &data, const size_t rr)
static auto get(ublas::matrix< T, L, A > &data, const size_t rr)
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)
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)