v0.16.0
Loading...
Searching...
No Matches
Templates.hpp
Go to the documentation of this file.
1/** \file Templates.hpp
2 * \brief Templates declarations
3 */
4
5#ifndef __TEMPLATES_HPP__
6#define __TEMPLATES_HPP__
7
8namespace MoFEM {
9
10/** @name Data layout helpers
11 * Traits describing tensor storage layout conventions.
12 * @{
13 */
14enum class DataLayout { GaussByCoeffs, CoeffsByGauss };
15template <DataLayout DL> struct DataLayoutTraits;
16template <> struct DataLayoutTraits<DataLayout::GaussByCoeffs> {
17 static constexpr bool isGaussByCoeffs = true;
18 static constexpr bool isCoeffsByGauss = false;
19};
20template <> struct DataLayoutTraits<DataLayout::CoeffsByGauss> {
21 static constexpr bool isGaussByCoeffs = false;
22 static constexpr bool isCoeffsByGauss = true;
23};
24/** @} */
25
26/** @name Adaptor functions
27 * Views that wrap raw storage as vector or matrix adaptors.
28 * @{
29 */
30/**
31 * @brief Get Vector adaptor
32 *
33 * \code
34 *
35 * double *a;
36 * CHKERR VecGetArray(v,&a);
37 *
38 * for(int n = 0; n != nodes; ++n) {
39 *
40 * auto a = getVectorAdaptor(&a[3*n], 3);
41 * double dot = inner_prod(a, a);
42 *
43 * }
44 *
45 * CHKERR VecRetsoreArray(v,&a);
46 * \endcode
47 *
48 */
49template <typename T1> inline auto getVectorAdaptor(T1 ptr, const size_t n) {
50 typedef typename std::remove_pointer<T1>::type T;
51 return VectorShallowArrayAdaptor<T>(n,
52 ublas::shallow_array_adaptor<T>(n, ptr));
53};
54
55/**
56 * @brief Get Matrix adaptor
57 *
58 * \code
59 *
60 * double *a;
61 * CHKERR VecGetArray(v,&a);
62 *
63 * for(int n = 0; n != nodes; ++n) {
64 *
65 * auto F = getMatrixAdaptor(&a[3*3*n], 3, 3);
66 * MatrixDouble C = prod(F, trans(F));
67 *
68 * }
69 *
70 * CHKERR VecRetsoreArray(v,&a);
71 * \endcode
72 *
73 */
74template <typename T1>
75inline auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m) {
76 typedef typename std::remove_pointer<T1>::type T;
77 return MatrixShallowArrayAdaptor<T>(
78 n, m, ublas::shallow_array_adaptor<T>(n * m, ptr));
79};
80/** @} */
81
82/** @name Generic utility helpers
83 * Shared template utilities used across the header.
84 * @{
85 */
86/**
87 * This small utility that cascades two key extractors will be
88 * used throughout the boost example
89 * <a
90 * href=http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp>
91 * http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp
92 * </a>
93 */
94template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
95public:
96 typedef typename KeyExtractor1::result_type result_type;
97
98 KeyFromKey(const KeyExtractor1 &key1_ = KeyExtractor1(),
99 const KeyExtractor2 &key2_ = KeyExtractor2())
100 : key1(key1_), key2(key2_) {}
101
102 template <typename Arg> result_type operator()(Arg &arg) const {
103 return key1(key2(arg));
104 }
105
106private:
107 KeyExtractor1 key1;
108 KeyExtractor2 key2;
109};
110
111template <typename id_type> struct LtBit {
112 inline bool operator()(const id_type &valueA, const id_type &valueB) const {
113 return valueA.to_ulong() < valueB.to_ulong();
114 }
115};
116
117template <typename id_type> struct EqBit {
118 inline bool operator()(const id_type &valueA, const id_type &valueB) const {
119 return valueA.to_ulong() == valueB.to_ulong();
120 }
121};
122
123template <typename id_type> struct HashBit {
124 inline unsigned int operator()(const id_type &value) const {
125 return value.to_ulong();
126 }
127};
128
129template <class X> inline std::string toString(X x) {
130 std::ostringstream buffer;
131 buffer << x;
132 return buffer.str();
133}
134/** @} */
135
136/** @name Tensor0 functions
137 * Scalar tensor access helpers.
138 * @{
139 */
140template <int S, class V> struct GetFTensor0FromVecImpl {
141 static inline auto get(V &data) {
142 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0))>>;
143 return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(data.data().data());
144 }
145};
146
147/**
148* \brief Get tensor rank 0 (scalar) form data vector
149
150Example how to use it.
151\code
152VectorDouble vec;
153vec.resize(nb_gauss_pts,false);
154vec.clear();
155auto t0 = getFTensor0FromVec<1>(vec);
156for(int gg = 0;gg!=nb_gauss_pts;gg++) {
157
158 ++t0;
159}
160\endcode
161
162*/
163template <int S = 1, class V> static inline auto getFTensor0FromVec(V &data) {
164 return GetFTensor0FromVecImpl<S, V>::get(data);
165}
166
167template <int S, class M> struct GetFTensor0FromMatImpl {
168 static inline auto get(M &data) {
169#ifndef NDEBUG
170 // We can do rows, or columns, does not matter for rank 0 tensor, but we need to be sure that one of them is 1
171 if (data.size1() != 1 && data.size2() != 1)
172 THROW_MESSAGE("getFTensor0FromMat: wrong size of data matrix, number of "
173 "rows should be 1 but is " +
174 boost::lexical_cast<std::string>(data.size1()));
175#endif
176 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
177 return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(data.data().data());
178 }
179};
180
181/**
182* \brief Get tensor rank 0 (scalar) form data vector
183
184Example how to use it.
185\code
186MatrixDouble vec;
187vec.resize(1, nb_gauss_pts, false);
188vec.clear();
189auto t0 = getFTensor0FromMat<1>(vec);
190for(int gg = 0;gg!=nb_gauss_pts;gg++) {
191
192 ++t0;
193}
194\endcode
195
196*/
197template <int S = 1, class M> inline auto getFTensor0FromMat(M &data) {
198 return GetFTensor0FromMatImpl<S, M>::get(data);
199}
200
201template <int S = 1, class T, class A>
202inline auto getFTensor0FromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
203#ifndef NDEBUG
204 if (data_ptr.expired())
205 THROW_MESSAGE("getFTensor0FromMat: data pointer expired");
206#endif
207 return GetFTensor0FromMatImpl<S, MatrixDouble>::get(*data_ptr.lock());
208}
209
210template <int S = 1>
211inline auto getFTensor0FromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
212#ifndef NDEBUG
213 if (!data_ptr)
214 THROW_MESSAGE("getFTensor0FromMat: data pointer is not available");
215#endif
216 return GetFTensor0FromMatImpl<S, MatrixDouble>::get(*data_ptr);
217}
218
219/**
220 * \brief Get tensor rank 0 (scalar) from pointer
221 */
222template <int S = 1, class T> static inline auto getFTensor0FromPtr(T *ptr) {
223 return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(ptr);
224}
225/** @} */
226
227/** @name Tensor1 functions
228 * Rank-1 tensor construction helpers.
229 * @{
230 */
231template <int Tensor_Dim, int S, class DL, class M>
232struct GetFTensor1FromMatImpl;
233
234template <int Tensor_Dim, int S, class M>
235struct GetFTensor1FromMatImpl<Tensor_Dim, S,
236 DataLayoutTraits<DataLayout::GaussByCoeffs>, M> {
237 static inline auto get(M &data, int rr, int cc) {
238 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
239 constexpr int stride = S == -1 ? Tensor_Dim : S;
240 static_assert(
241 stride % Tensor_Dim == 0,
242 "getFTensor1FromMat: stride should be a multiple of Tensor_Dim");
243
244#ifndef NDEBUG
245 if (data.size2() != Tensor_Dim)
246 THROW_MESSAGE(
247 "getFTensor1FromMat<" + boost::lexical_cast<std::string>(Tensor_Dim) +
248 ">: wrong size of data matrix, number of columns should be " +
249 boost::lexical_cast<std::string>(Tensor_Dim) + " but is " +
250 boost::lexical_cast<std::string>(data.size2()));
251 auto *first = &data(rr + 0, cc + 0);
252 auto *last = &data(rr + 0, cc + 0 + Tensor_Dim - 1);
253 if (last - first != Tensor_Dim - 1)
254 THROW_MESSAGE("getFTensor1FromMat<" +
255 boost::lexical_cast<std::string>(Tensor_Dim) +
256 ">: row slice is not contiguous");
257
258#endif
259 return FTensor::Tensor1<FTensor::CursorPtr<T *, stride>, Tensor_Dim>(
260 &data(rr + 0, cc + 0));
261 }
262};
263
264template <int S, class M>
265struct GetFTensor1FromMatImpl<3, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
266 M> {
267 static inline auto get(M &data, int rr, int cc) {
268 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
269 constexpr int stride = S == -1 ? 1 : S;
270
271#ifndef NDEBUG
272 if (data.size1() % 3 != 0)
273 THROW_MESSAGE(
274 "getFTensor1FromMat<3>: wrong size of data matrix, number of "
275 "rows modulo of 3 but is " +
276 boost::lexical_cast<std::string>(data.size1()));
277#endif
278 return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 3>(
279 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0));
280 }
281};
282
283template <int S, class M>
284struct GetFTensor1FromMatImpl<4, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
285 M> {
286 static inline auto get(M &data, int rr, int cc) {
287 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
288 constexpr int stride = S == -1 ? 1 : S;
289#ifndef NDEBUG
290 if (data.size1() % 4 != 0)
291 THROW_MESSAGE(
292 "getFTensor1FromMat<4>: wrong size of data matrix, number of "
293 "rows should be 4 but is " +
294 boost::lexical_cast<std::string>(data.size1()));
295#endif
296 return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 4>(
297 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
298 &data(rr + 3, cc + 0));
299 }
300};
301
302template <int S, class M>
303struct GetFTensor1FromMatImpl<6, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
304 M> {
305 static inline auto get(M &data, int rr, int cc) {
306 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
307 constexpr int stride = S == -1 ? 1 : S;
308#ifndef NDEBUG
309 if (data.size1() % 6 != 0)
310 THROW_MESSAGE(
311 "getFTensor1FromMat<6>: wrong size of data matrix, number of "
312 "rows should be 6 but is " +
313 boost::lexical_cast<std::string>(data.size1()));
314#endif
315 return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 6>(
316 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
317 &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0));
318 }
319};
320
321template <int S, class M>
322struct GetFTensor1FromMatImpl<9, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
323 M> {
324 static inline auto get(M &data, int rr, int cc) {
325 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
326 constexpr int stride = S == -1 ? 1 : S;
327#ifndef NDEBUG
328 if (data.size1() % 9 != 0)
329 THROW_MESSAGE(
330 "getFTensor1FromMat<9>: wrong size of data matrix, number of "
331 "rows should be 9 but is " +
332 boost::lexical_cast<std::string>(data.size1()));
333#endif
334 return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 9>(
335 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
336 &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0),
337 &data(rr + 6, cc + 0), &data(rr + 7, cc + 0), &data(rr + 8, cc + 0));
338 }
339};
340
341template <int S, class M>
342struct GetFTensor1FromMatImpl<2, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
343 M> {
344 static inline auto get(M &data, int rr, int cc) {
345 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
346 constexpr int stride = S == -1 ? 1 : S;
347#ifndef NDEBUG
348 if (data.size1() % 2 != 0)
349 THROW_MESSAGE(
350 "getFTensor1FromMat<2>: wrong size of data matrix, number of "
351 "row should be 2 but is " +
352 boost::lexical_cast<std::string>(data.size1()));
353#endif
354 return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 2>(
355 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0));
356 }
357};
358
359template <int S, class M>
360struct GetFTensor1FromMatImpl<1, S, DataLayoutTraits<DataLayout::CoeffsByGauss>,
361 M> {
362 static inline auto get(M &data, int rr, int cc) {
363 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
364 constexpr int stride = S == -1 ? 1 : S;
365 return FTensor::Tensor1<FTensor::PackPtr<T *, stride>, 1>(
366 &data(rr + 0, cc + 0));
367 }
368};
369
370/**
371 * \brief Get tensor rank 1 (vector) form data matrix
372 */
373template <int Tensor_Dim, int S = -1,
374 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
375 class M = MatrixDouble>
376inline auto getFTensor1FromMat(M &data, int rr = 0, int cc = 0) {
377 return GetFTensor1FromMatImpl<Tensor_Dim, S, DL, M>::get(data, rr, cc);
378}
379
380template <int Tensor_Dim, int S = -1,
381 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
382 class M = MatrixDouble>
383using GetFTensor1FromMatType =
384 decltype(GetFTensor1FromMatImpl<Tensor_Dim, S, DL, M>::get(
385 std::declval<M &>(), 0, 0));
386
387template <int Tensor_Dim, int S = -1,
388 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
389inline auto getFTensor1FromMat(boost::weak_ptr<MatrixDouble> data, int rr = 0,
390 int cc = 0) {
391#ifndef NDEBUG
392 if (!data.lock()) {
393 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
394 "Data matrix is not available");
395 }
396#endif
397 return GetFTensor1FromMatImpl<Tensor_Dim, S, DL, MatrixDouble>::get(
398 *data.lock(), rr, cc);
399}
400
401template <int Tensor_Dim, int S = -1,
402 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
403inline auto getFTensor1FromMat(boost::shared_ptr<MatrixDouble> data, int rr = 0,
404 int cc = 0) {
405#ifndef NDEBUG
406 if (!data) {
407 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
408 "Data matrix is not available");
409 }
410#endif
411 return GetFTensor1FromMatImpl<Tensor_Dim, S, DL, MatrixDouble>::get(*data, rr,
412 cc);
413}
414
415template <int Tensor_Dim, int S = Tensor_Dim, typename T = double>
416inline auto getFTensor1FromPtr(T *ptr) {
417 return FTensor::Tensor1<FTensor::CursorPtr<T *, S>, Tensor_Dim>(ptr);
418}
419
420template <int Tensor_Dim, int S = Tensor_Dim>
421inline auto getFTensor1FromPtr(std::complex<double> *ptr) {
422 return getFTensor1FromPtr<Tensor_Dim, S, std::complex<double>>(ptr);
423}
424
425#ifdef WITH_ADOL_C
426template <int Tensor_Dim, int S = Tensor_Dim>
427inline auto getFTensor1FromPtr(adouble *ptr) {
428 return getFTensor1FromPtr<Tensor_Dim, S, adouble>(ptr);
429}
430#endif
431
432template <int Tensor_Dim, int S = Tensor_Dim, typename T = VectorDouble>
433inline auto getFTensor1FromArray(T &data) {
434 static_assert(S % Tensor_Dim == 0,
435 "getFTensor1FromArray(VectorDouble&) requires S to be a "
436 "multiple of Tensor_Dim");
437
438#ifndef NDEBUG
439 if (data.size() % S != 0) {
440 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
441 "getFTensor1FromArray: data size should be divisible by "
442 "pack stride dimension");
443 }
444#endif
445 return getFTensor1FromPtr<Tensor_Dim, S>(data.data().data());
446}
447
448/**
449 * @brief Get FTensor1 from array
450 *
451 * \todo Generalize for different arrays and data types
452 *
453 * @tparam DIM
454 * @param data
455 * @param rr
456 * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
457 */
458template <int DIM, int S>
459inline auto getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr);
460
461template <>
462inline auto getFTensor1FromArrayDiag<2, 2>(MatrixDouble &data,
463 const size_t rr) {
464 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
465 &data(rr + 1, 1)};
466}
467
468template <>
469inline auto getFTensor1FromArrayDiag<3, 3>(MatrixDouble &data,
470 const size_t rr) {
471 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>{
472 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
473}
474/** @} */
475
476/** @name Tensor2 functions
477 * Rank-2 tensor construction and array helper routines.
478 * @{
479 */
480template <int Tensor_Dim1, int Tensor_Dim2, int S, class DL, class M>
481struct GetFTensor2FromMatImpl {
482 static inline auto get(M &data, int rr = 0, int cc = 0) {
483 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
484 constexpr int stride =
485 S == -1 ? (DL::isGaussByCoeffs ? Tensor_Dim1 * Tensor_Dim2 : 1) : S;
486 if constexpr (DL::isGaussByCoeffs) {
487 static_assert(stride % (Tensor_Dim1 * Tensor_Dim2) == 0,
488 "getFTensor2FromMat: stride should be a multiple of "
489 "Tensor_Dim1 * Tensor_Dim2");
490
491#ifndef NDEBUG
492 if (data.size2() != Tensor_Dim1 * Tensor_Dim2) {
493 THROW_MESSAGE(
494 "getFTensor2FromMat<" +
495 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
496 boost::lexical_cast<std::string>(Tensor_Dim2) +
497 ">: wrong size of columns in data matrix, should be " +
498 boost::lexical_cast<std::string>(Tensor_Dim1 * Tensor_Dim2) +
499 " but is " + boost::lexical_cast<std::string>(data.size2()));
500 }
501 auto *first = &data(rr + 0, cc + 0);
502 auto *last = &data(rr + 0, cc + 0 + Tensor_Dim1 * Tensor_Dim2 - 1);
503 if (last - first != Tensor_Dim1 * Tensor_Dim2 - 1) {
504 THROW_MESSAGE("getFTensor2FromMat<" +
505 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
506 boost::lexical_cast<std::string>(Tensor_Dim2) +
507 ">: row slice is not contiguous");
508 }
509#endif
510 return FTensor::Tensor2<FTensor::CursorPtr<T *, stride>, Tensor_Dim1,
511 Tensor_Dim2>(&data(rr + 0, cc + 0));
512 } else if constexpr (DL::isCoeffsByGauss) {
513 std::array<T *, Tensor_Dim1 * Tensor_Dim2> ptrs;
514#ifndef NDEBUG
515 if (data.size1() != Tensor_Dim1 * Tensor_Dim2) {
516 THROW_MESSAGE(
517 "getFTensor2FromMat<" +
518 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
519 boost::lexical_cast<std::string>(Tensor_Dim2) +
520 ">: wrong size of rows in data matrix, should be " +
521 boost::lexical_cast<std::string>(Tensor_Dim1 * Tensor_Dim2) +
522 " but is " + boost::lexical_cast<std::string>(data.size1()));
523 }
524#endif
525 for (auto i = 0; i != Tensor_Dim1 * Tensor_Dim2; ++i)
526 ptrs[i] = &data(rr + i, cc + 0);
527 return FTensor::Tensor2<FTensor::PackPtr<T *, stride>, Tensor_Dim1,
528 Tensor_Dim2>(ptrs);
529 } else {
530 static_assert(!std::is_same<M, M>::value,
531 "Unsupported data layout for getFTensor2FromMatImpl");
532 }
533 }
534};
535
536/**
537 * \brief Get tensor rank 2 (matrix) form data matrix
538 */
539template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
540 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
541 class M = MatrixDouble>
545
546template <int Tensor_Dim0, int Tensor_Dim1, int S = -1,
547 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
548 class M = MatrixDouble>
549using GetFTensor2FromMatType =
550 decltype(GetFTensor2FromMatImpl<Tensor_Dim0, Tensor_Dim1, S, DL, M>::get(
551 std::declval<M &>(), 0, 0));
552
553/**
554 * \brief Get tensor rank 2 (matrix) form data matrix
555 */
556template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
557 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
558inline auto getFTensor2FromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
559#ifndef NDEBUG
560 if (!data_ptr.lock()) {
561 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
562 "Data matrix is not available");
563 }
564#endif
565 return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
566 MatrixDouble>::get(*(data_ptr.lock()));
567}
568
569/**
570 * \brief Get tensor rank 2 (matrix) form data matrix
571 */
572template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
573 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
574inline auto getFTensor2FromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
575#ifndef NDEBUG
576 if (!data_ptr) {
577 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
578 "Data matrix is not available");
579 }
580#endif
581 return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
582 MatrixDouble>::get(*data_ptr);
583}
584
585template <int Tensor_Dim1, int Tensor_Dim2, int S = -1>
586inline auto getFTensor2FromVec(VectorDouble &data) {
587 auto matrix_data =
588 getMatrixAdaptor(&data[0], data.size() / Tensor_Dim1 * Tensor_Dim2,
589 Tensor_Dim1 * Tensor_Dim2);
590 return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S,
591 DataLayoutTraits<DataLayout::GaussByCoeffs>,
592 decltype(matrix_data)>::get(matrix_data);
593}
594
595template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
596 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
597inline auto getFTensor2FromPtr(double *ptr) {
598 if constexpr (DL::isGaussByCoeffs) {
599 auto matrix_data = getMatrixAdaptor(ptr, 1, Tensor_Dim1 * Tensor_Dim2);
600 return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
601 decltype(matrix_data)>::get(matrix_data);
602 } else if constexpr (DL::isCoeffsByGauss) {
603 auto matrix_data = getMatrixAdaptor(ptr, Tensor_Dim1 * Tensor_Dim2, 1);
604 return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
605 decltype(matrix_data)>::get(matrix_data);
606 }
607}
608
609template <int Tensor_Dim1, int Tensor_Dim2, int S = -1,
610 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
611inline auto getFTensor2FromPtr(std::complex<double> *ptr) {
612 if constexpr (DL::isGaussByCoeffs) {
613 auto matrix_data = getMatrixAdaptor(ptr, 1, Tensor_Dim1 * Tensor_Dim2);
614 return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
615 decltype(matrix_data)>::get(matrix_data);
616 } else if constexpr (DL::isCoeffsByGauss) {
617 auto matrix_data = getMatrixAdaptor(ptr, Tensor_Dim1 * Tensor_Dim2, 1);
618 return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, S, DL,
619 decltype(matrix_data)>::get(matrix_data);
620 }
621}
622
623/**
624 * @brief Make Tensor2 for HVec base from pointer
625 *
626 * @tparam DIM
627 * @param ptr
628 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
629 */
630template <int DIM1, int DIM2> inline auto getFTensor2HVecFromPtr(double *ptr);
631
632template <> inline auto getFTensor2HVecFromPtr<3, 2>(double *ptr) {
633 return FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2>(
634 ptr + HVEC0_0, ptr + HVEC0_1,
635
636 ptr + HVEC1_0, ptr + HVEC1_1,
637
638 ptr + HVEC2_0, ptr + HVEC2_1);
639}
640
641template <> inline auto getFTensor2HVecFromPtr<3, 3>(double *ptr) {
642 return FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>(
643 ptr + HVEC0_0, ptr + HVEC0_1, ptr + HVEC0_2,
644
645 ptr + HVEC1_0, ptr + HVEC1_1, ptr + HVEC1_2,
646
647 ptr + HVEC2_0, ptr + HVEC2_1, ptr + HVEC2_2);
648}
649/**
650 * @brief Get FTensor2 from array
651 *
652 * \note Generalize for other data types
653 *
654 * @tparam DIM1
655 * @tparam DIM2
656 * @tparam S
657 * @param data
658 * @return FTensor::Tensor2<FTensor::PackPtr<T *, S>, DIM1, DIM2>
659 */
660template <int DIM1, int DIM2, int S, class M> struct GetFTensor2FromArrayImpl;
661
662/**
663 * @brief Get FTensor2 from array
664 *
665 * \note Generalize for other data types
666 *
667 * @tparam DIM1
668 * @tparam DIM2
669 * @tparam S
670 * @param data
671 * @return FTensor::Tensor2<T *, DIM1, DIM2>
672 */
673template <int DIM1, int DIM2, class M> struct GetFTensor2FromArrayRawPtrImpl;
674
675template <int S, class M> struct GetFTensor2FromArrayImpl<2, 2, S, M> {
676 GetFTensor2FromArrayImpl() = delete;
677 inline static auto get(M &data, const size_t rr, const size_t cc) {
678 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
679 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 2, 2>{
680 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
681
682 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
683 }
684};
685
686template <int S, class M> struct GetFTensor2FromArrayImpl<3, 3, S, M> {
687 GetFTensor2FromArrayImpl() = delete;
688 inline static auto get(M &data, const size_t rr, const size_t cc) {
689 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
690 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 3, 3>{
691 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
692 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
693 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2)};
694 }
695};
696
697template <class M> struct GetFTensor2FromArrayRawPtrImpl<2, 2, M> {
698 GetFTensor2FromArrayRawPtrImpl() = delete;
699 inline static auto get(M &data, const size_t rr, const size_t cc,
700 const int ss = 0) {
701 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
702 return FTensor::Tensor2<T *, 2, 2>(
703 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
704
705 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
706 }
707};
708
709template <class M> struct GetFTensor2FromArrayRawPtrImpl<3, 3, M> {
710 GetFTensor2FromArrayRawPtrImpl() = delete;
711 inline static auto get(M &data, const size_t rr, const size_t cc,
712 const int ss = 0) {
713 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
714 return FTensor::Tensor2<T *, 3, 3>(
715 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
716 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
717 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
718 ss);
719 }
720};
721
722template <int DIM1, int DIM2, int S>
723inline auto getFTensor2FromArray(MatrixDouble &data, const size_t rr,
724 const size_t cc = 0) {
725 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, MatrixDouble>::get(data, rr,
726 cc);
727}
728
729template <int DIM1, int DIM2>
730inline auto getFTensor2FromArray(MatrixDouble &data, const size_t rr,
731 const size_t cc, const int ss) {
732 return GetFTensor2FromArrayRawPtrImpl<DIM1, DIM2, MatrixDouble>::get(data, rr,
733 cc, ss);
734}
735
736#ifdef WITH_ADOL_C
737template <int DIM1, int DIM2, int S>
738inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
739 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, MatrixADouble>::get(data, rr);
740}
741#endif
742/** @} */
743
744/** @name Tensor2_symmetric functions
745 * Symmetric rank-2 tensor construction and conversion helpers.
746 * @{
747 */
748template <int Tensor_Dim, int S, class DL, class M>
749struct GetFTensor2SymmetricFromMatImpl {};
750
751template <int Tensor_Dim, int S, class M>
752struct GetFTensor2SymmetricFromMatImpl<
753 Tensor_Dim, S, DataLayoutTraits<DataLayout::GaussByCoeffs>, M> {
754 static inline auto get(M &data, int rr = 0, int cc = 0) {
755 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
756 constexpr int symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
757 constexpr int stride = S == -1 ? symm_size : S;
758 static_assert(stride % symm_size == 0,
759 "getFTensor2SymmetricFromMat: stride should be a multiple of "
760 "symm_size");
761
762#ifndef NDEBUG
763 if (data.size2() != symm_size)
764 THROW_MESSAGE(
765 "getFTensor2SymmetricFromMat<" +
766 boost::lexical_cast<std::string>(Tensor_Dim) +
767 ">: wrong size of data matrix, number of columns should be " +
768 boost::lexical_cast<std::string>(symm_size) + " but is " +
769 boost::lexical_cast<std::string>(data.size2()));
770 auto *first = &data(rr + 0, cc + 0);
771 auto *last = &data(rr + 0, cc + 0 + symm_size - 1);
772 if (last - first != symm_size - 1)
773 THROW_MESSAGE("getFTensor2SymmetricFromMat<" +
774 boost::lexical_cast<std::string>(Tensor_Dim) +
775 ">: row slice is not contiguous");
776#endif
777 return FTensor::Tensor2_symmetric<FTensor::CursorPtr<T *, stride>,
778 Tensor_Dim>(&data(rr + 0, cc + 0));
779 }
780};
781
782template <int S, class M>
783struct GetFTensor2SymmetricFromMatImpl<
784 3, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
785 static inline auto get(M &data, int rr = 0, int cc = 0) {
786 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
787 constexpr int stride = S == -1 ? 1 : S;
788#ifndef NDEBUG
789 constexpr int symm_size = 6;
790 if (data.size1() != symm_size)
791 THROW_MESSAGE(
792 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, "
793 "number of rows should be 6 but is " +
794 boost::lexical_cast<std::string>(data.size1()));
795#endif
796 return FTensor::Tensor2_symmetric<FTensor::PackPtr<T *, stride>, 3>(
797 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
798 &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0));
799 }
800};
801
802template <int S, class M>
803struct GetFTensor2SymmetricFromMatImpl<
804 2, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
805 static inline auto get(M &data, int rr = 0, int cc = 0) {
806 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
807 constexpr int stride = S == -1 ? 1 : S;
808#ifndef NDEBUG
809 constexpr int symm_size = 3;
810 if (data.size1() != symm_size)
811 THROW_MESSAGE(
812 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, "
813 "number of rows should be 3 but is " +
814 boost::lexical_cast<std::string>(data.size1()));
815#endif
816 return FTensor::Tensor2_symmetric<FTensor::PackPtr<T *, stride>, 2>(
817 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0));
818 }
819};
820
821/**
822 * \brief Get symmetric tensor rank 2 (matrix) form data matrix
823 */
824template <int Tensor_Dim, int S = -1,
825 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
826 class M = MatrixDouble>
830
831template <int Tensor_Dim, int S = -1,
832 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
833 class M = MatrixDouble>
834using GetFTensor2SymmetricFromMatType =
835 decltype(GetFTensor2SymmetricFromMatImpl<Tensor_Dim, S, DL, M>::get(
836 std::declval<M &>(), 0, 0));
837
838template <int Tensor_Dim, int S = -1,
839 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
840static inline auto
841getFTensor2SymmetricFromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
842#ifndef NDEBUG
843 if (data_ptr.expired()) {
844 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
845 "Data matrix is not available");
846 }
847#endif
848 return GetFTensor2SymmetricFromMatImpl<Tensor_Dim, S, DL, MatrixDouble>::get(
849 *data_ptr.lock());
850}
851
852template <int Tensor_Dim, int S = -1,
853 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
854static inline auto
855getFTensor2SymmetricFromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
856#ifndef NDEBUG
857 if (!data_ptr) {
858 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
859 "Data matrix is not available");
860 }
861#endif
862 return GetFTensor2SymmetricFromMatImpl<Tensor_Dim, S, DL, MatrixDouble>::get(
863 *data_ptr);
864}
865
866template <int Tensor_Dim> inline auto getFTensor2SymmetricFromPtr(double *ptr) {
867 auto matrix_data =
868 getMatrixAdaptor(ptr, 1, Tensor_Dim * (Tensor_Dim + 1) / 2);
869 return GetFTensor2SymmetricFromMatImpl<
870 Tensor_Dim, -1, DataLayoutTraits<DataLayout::GaussByCoeffs>,
871 decltype(matrix_data)>::get(matrix_data, 0, 0);
872}
873
874#ifdef WITH_ADOL_C
875template <int Tensor_Dim>
876inline auto getFTensor2SymmetricFromPtr(adouble *ptr) {
877 auto matrix_data =
878 getMatrixAdaptor(ptr, 1, Tensor_Dim * (Tensor_Dim + 1) / 2);
879 return GetFTensor2SymmetricFromMatImpl<
880 Tensor_Dim, -1, DataLayoutTraits<DataLayout::GaussByCoeffs>,
881 decltype(matrix_data)>::get(matrix_data, 0, 0);
882}
883#endif
884/**
885 * @brief Make symmetric Tensor2 from pointer, taking lower triangle of matrix
886 *
887 * @tparam DIM
888 * @param ptr
889 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
890 */
891template <int DIM> inline auto getFTensor2SymmetricLowerFromPtr(double *ptr);
892
893template <> inline auto getFTensor2SymmetricLowerFromPtr<3>(double *ptr) {
894 return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 9>, 3>(
895 ptr + HVEC0_0, ptr + HVEC0_1, ptr + HVEC0_2,
896
897 ptr + HVEC1_0, ptr + HVEC1_1,
898
899 ptr + HVEC2_2);
900};
901
902template <> inline auto getFTensor2SymmetricLowerFromPtr<2>(double *ptr) {
903 return FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 4>, 2>(
904 ptr + 0, ptr + 1, ptr + 3);
905};
906// Template function for conversion of symmetric tensor to non-symmetric
907
908template <int DIM, class T>
909inline auto to_non_symm(const FTensor::Tensor2_symmetric<T, DIM> &symm) {
910 FTensor::Tensor2<T, DIM, DIM> full;
911 FTensor::Index<'i', DIM> i;
912 FTensor::Index<'j', DIM> j;
913 full(i, j) = symm(i, j);
914 return full;
915}
916/** @} */
917
918/** @name Tensor3 functions
919 * Rank-3 tensor construction helpers.
920 * @{
921 */
922template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class DL,
923 class M>
924struct GetFTensor3FromMatImpl;
925
926template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class M>
927struct GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
928 DataLayoutTraits<DataLayout::GaussByCoeffs>, M> {
929 static inline auto get(M &data, int rr = 0, int cc = 0) {
930 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
931 constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
932 constexpr int stride = S == -1 ? tensor_size : S;
933 static_assert(
934 stride % tensor_size == 0,
935 "getFTensor3FromMat: stride should be a multiple of tensor_size");
936
937#ifndef NDEBUG
938 if (data.size2() != tensor_size) {
939 THROW_MESSAGE("getFTensor3FromMat<" +
940 boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
941 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
942 boost::lexical_cast<std::string>(Tensor_Dim2) +
943 ">: wrong size of columns in data matrix, should be " +
944 boost::lexical_cast<std::string>(tensor_size) + " but is " +
945 boost::lexical_cast<std::string>(data.size2()));
946 }
947 auto *first = &data(rr + 0, cc + 0);
948 auto *last = &data(rr + 0, cc + 0 + tensor_size - 1);
949 if (last - first != tensor_size - 1)
950 THROW_MESSAGE("getFTensor3FromMat<" +
951 boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
952 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
953 boost::lexical_cast<std::string>(Tensor_Dim2) +
954 ">: row slice is not contiguous");
955#endif
956 return FTensor::Tensor3<FTensor::CursorPtr<T *, stride>, Tensor_Dim0,
957 Tensor_Dim1, Tensor_Dim2>{&data(rr + 0, cc + 0)};
958 }
959};
960
961template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class M>
962struct GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
963 DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
964 static inline auto get(M &data, int rr = 0, int cc = 0) {
965 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
966 constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
967 constexpr int stride = S == -1 ? 1 : S;
968 std::array<T *, tensor_size> ptrs;
969#ifndef NDEBUG
970 if (data.size1() != tensor_size) {
971 THROW_MESSAGE("getFTensor3FromMat<" +
972 boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
973 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
974 boost::lexical_cast<std::string>(Tensor_Dim2) +
975 ">: wrong size of rows in data matrix, should be " +
976 boost::lexical_cast<std::string>(tensor_size) + " but is " +
977 boost::lexical_cast<std::string>(data.size1()));
978 }
979#endif
980 for (auto i = 0; i != tensor_size; ++i)
981 ptrs[i] = &data(rr + i, cc + 0);
982 return FTensor::Tensor3<FTensor::PackPtr<T *, stride>, Tensor_Dim0,
983 Tensor_Dim1, Tensor_Dim2>(ptrs);
984 }
985};
986
987/**
988 * @brief Get tensor rank 3 (non symmetries) form data matrix
989 *
990 * @tparam Tensor_Dim0 dimension of first index
991 * @tparam Tensor_Dim1 dimension of second index
992 * @tparam Tensor_Dim2 dimension of third index
993 * @tparam S shift size
994 * @tparam DL data layout
995 * @tparam M matrix type
996 * @param data data container
997 * @return FTensor::Tensor3<FTensor::PackPtr<T *, S>, Tensor_Dim0,
998 Tensor_Dim1, Tensor_Dim2>
999 */
1000template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = -1,
1001 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
1002 class M = MatrixDouble>
1003static inline auto getFTensor3FromMat(M &data, int rr = 0, int cc = 0) {
1004 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL,
1005 M>::get(data, rr, cc);
1006}
1007
1008template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = -1,
1009 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
1010 class M = MatrixDouble>
1011using GetFTensor3FromMatType =
1012 decltype(GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
1013 DL, M>::get(std::declval<M &>(), 0, 0));
1014
1015template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = -1,
1016 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
1017static inline auto getFTensor3FromMat(boost::weak_ptr<MatrixDouble> data_ptr,
1018 int rr = 0, int cc = 0) {
1019#ifndef NDEBUG
1020 if (!data_ptr.lock()) {
1021 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
1022 "Data matrix is not available");
1023 }
1024#endif
1025 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL,
1026 MatrixDouble>::get(*data_ptr.lock(), rr, cc);
1027}
1028
1029template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = -1,
1030 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
1031static inline auto getFTensor3FromMat(boost::shared_ptr<MatrixDouble> data_ptr,
1032 int rr = 0, int cc = 0) {
1033#ifndef NDEBUG
1034 if (!data_ptr) {
1035 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
1036 "Data matrix is not available");
1037 }
1038#endif
1039 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL,
1040 MatrixDouble>::get(*data_ptr, rr, cc);
1041}
1042/*
1043 * @brief Make Tensor3 from pointer
1044 *
1045 * @tparam DIM
1046 * @param ptr
1047 * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3>, DIM1,
1048 * DIM2, DIM3>
1049 */
1050template <int DIM1, int DIM2, int DIM3>
1051inline FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2 * DIM3>, DIM1,
1052 DIM2, DIM3>
1053getFTensor3FromPtr(double *ptr) {
1054 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
1055 "Such getFTensor3FromPtr is not implemented");
1056}
1057
1058template <>
1059inline FTensor::Tensor3<FTensor::PackPtr<double *, 12>, 3, 2, 2>
1060getFTensor3FromPtr<3, 2, 2>(double *ptr) {
1061 return FTensor::Tensor3<FTensor::PackPtr<double *, 12>, 3, 2, 2>(
1062 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
1063 ptr + 8, ptr + 9, ptr + 10, ptr + 11);
1064}
1065
1066template <>
1067inline FTensor::Tensor3<FTensor::PackPtr<double *, 27>, 3, 3, 3>
1068getFTensor3FromPtr<3, 3, 3>(double *ptr) {
1069 return FTensor::Tensor3<FTensor::PackPtr<double *, 27>, 3, 3, 3>(
1070 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
1071 ptr + 8, ptr + 9, ptr + 10, ptr + 11, ptr + 12, ptr + 13, ptr + 14,
1072 ptr + 15, ptr + 16, ptr + 17, ptr + 18, ptr + 19, ptr + 20, ptr + 21,
1073 ptr + 22, ptr + 23, ptr + 24, ptr + 25, ptr + 26);
1074}
1075/** @} */
1076
1077/** @name Dg functions
1078 * Helpers for third-order tensors symmetric in the first two indices.
1079 * @{
1080 */
1081template <int Tensor_Dim01, int Tensor_Dim2, int S, class DL, class M>
1082struct GetFTensor3DgFromMatImpl {};
1083
1084template <int Tensor_Dim01, int Tensor_Dim2, int S, class M>
1085struct GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S,
1086 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1087 M> {
1088 static inline auto get(M &data, int rr = 0, int cc = 0) {
1089 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1090 constexpr int symm_size_01 = (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2;
1091 constexpr int dg_size = symm_size_01 * Tensor_Dim2;
1092 constexpr int stride = S == -1 ? dg_size : S;
1093 static_assert(
1094 stride % dg_size == 0,
1095 "getFTensor3DgFromMat: stride should be a multiple of dg_size");
1096
1097#ifndef NDEBUG
1098 if (data.size2() != dg_size) {
1099 THROW_MESSAGE(
1100 "getFTensor3DgFromMat<" +
1101 boost::lexical_cast<std::string>(Tensor_Dim01) + ", " +
1102 boost::lexical_cast<std::string>(Tensor_Dim2) +
1103 ">: wrong size of data matrix, number of columns should be " +
1104 boost::lexical_cast<std::string>(dg_size) + " but is " +
1105 boost::lexical_cast<std::string>(data.size2()));
1106 }
1107 auto *first = &data(rr + 0, cc + 0);
1108 auto *last = &data(rr + 0, cc + 0 + dg_size - 1);
1109 if (last - first != dg_size - 1)
1110 THROW_MESSAGE("getFTensor3DgFromMat<" +
1111 boost::lexical_cast<std::string>(Tensor_Dim01) + ", " +
1112 boost::lexical_cast<std::string>(Tensor_Dim2) +
1113 ">: row slice is not contiguous");
1114#endif
1115 return FTensor::Dg<FTensor::CursorPtr<T *, stride>, Tensor_Dim01,
1116 Tensor_Dim2>{&data(rr + 0, cc + 0)};
1117 }
1118};
1119
1120template <int S, class M>
1121struct GetFTensor3DgFromMatImpl<
1122 1, 1, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
1123 static inline auto get(M &data, int rr = 0, int cc = 0) {
1124 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1125 constexpr int stride = S == -1 ? 1 : S;
1126
1127#ifndef NDEBUG
1128 constexpr int dg_size = 1;
1129 if (data.size1() != dg_size)
1130 THROW_MESSAGE(
1131 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
1132 "of rows should be 1 but is " +
1133 boost::lexical_cast<std::string>(data.size1()));
1134#endif
1135 return FTensor::Dg<FTensor::PackPtr<T *, stride>, 1, 1>{
1136 &data(rr + 0, cc + 0)};
1137 }
1138};
1139
1140template <int S, class M>
1141struct GetFTensor3DgFromMatImpl<
1142 2, 2, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
1143 static inline auto get(M &data, int rr = 0, int cc = 0) {
1144 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1145 constexpr int stride = S == -1 ? 1 : S;
1146#ifndef NDEBUG
1147 constexpr int dg_size = 6;
1148 if (data.size1() != dg_size) {
1149 THROW_MESSAGE(
1150 "getFTensor3DgFromMat<2, 2>: wrong size of data matrix, number "
1151 "of rows should be 6 but is " +
1152 boost::lexical_cast<std::string>(data.size1()));
1153 }
1154#endif
1155 return FTensor::Dg<FTensor::PackPtr<T *, stride>, 2, 2>{
1156 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
1157 &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0)};
1158 }
1159};
1160
1161template <int S, class M>
1162struct GetFTensor3DgFromMatImpl<
1163 3, 3, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
1164 static inline auto get(M &data, int rr = 0, int cc = 0) {
1165 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1166 constexpr int stride = S == -1 ? 1 : S;
1167#ifndef NDEBUG
1168 constexpr int dg_size = 18;
1169 if (data.size1() != dg_size) {
1170 THROW_MESSAGE(
1171 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
1172 "of rows should be 18 but is " +
1173 boost::lexical_cast<std::string>(data.size1()));
1174 }
1175#endif
1176 return FTensor::Dg<FTensor::PackPtr<T *, stride>, 3, 3>{
1177 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
1178 &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0),
1179 &data(rr + 6, cc + 0), &data(rr + 7, cc + 0), &data(rr + 8, cc + 0),
1180 &data(rr + 9, cc + 0), &data(rr + 10, cc + 0), &data(rr + 11, cc + 0),
1181 &data(rr + 12, cc + 0), &data(rr + 13, cc + 0), &data(rr + 14, cc + 0),
1182 &data(rr + 15, cc + 0), &data(rr + 16, cc + 0), &data(rr + 17, cc + 0)};
1183 }
1184};
1185
1186/**
1187 * @brief Get symmetric tensor rank 3 on the first two indices from
1188 * form data matrix
1189 *
1190 * @tparam Tensor_Dim01 dimension of first two indicies
1191 * @tparam Tensor_Dim2 dimension of last index
1192 * @tparam S Memory shift
1193 * @tparam DL data layout
1194 * @tparam M matrix type
1195 * @param data data container
1196 * @return FTensor::Dg<FTensor::PackPtr<T *, S>, Tensor_Dim01, Tensor_Dim2>
1197 */
1198template <int Tensor_Dim01, int Tensor_Dim2, int S = -1,
1199 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
1200 class M = MatrixDouble>
1201static inline auto getFTensor3DgFromMat(M &data) {
1202 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, DL, M>::get(
1203 data);
1204}
1205
1206template <int Tensor_Dim01, int Tensor_Dim2, int S = -1,
1207 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
1208 class M = MatrixDouble>
1209using GetFTensor3DgFromMatType =
1210 decltype(GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, DL, M>::get(
1211 std::declval<M &>(), 0, 0));
1212
1213template <int Tensor_Dim01, int Tensor_Dim2, int S = -1,
1214 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
1215static inline auto
1216getFTensor3DgFromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
1217#ifndef NDEBUG
1218 if (!data_ptr.lock()) {
1219 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
1220 "Data matrix is not available");
1221 }
1222#endif
1223 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, DL,
1224 MatrixDouble>::get(*(data_ptr.lock()));
1225}
1226
1227template <int Tensor_Dim01, int Tensor_Dim2, int S = -1,
1228 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
1229static inline auto
1230getFTensor3DgFromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
1231#ifndef NDEBUG
1232 if (!data_ptr) {
1233 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
1234 "Data matrix is not available");
1235 }
1236#endif
1237 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, DL,
1238 MatrixDouble>::get(*data_ptr);
1239}
1240
1241template <int Tensor_Dim01, int Tensor_Dim2, int S, class T = double>
1242static inline auto getFTensor3DgFromPtr(T *ptr) {
1243 constexpr auto dg_size =
1244 ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) * Tensor_Dim2;
1245 auto matrix_data = getMatrixAdaptor(ptr, 1, dg_size);
1246 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S,
1247 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1248 decltype(matrix_data)>::get(matrix_data);
1249}
1250/** @} */
1251
1252/** @name Tensor4 functions
1253 * Rank-4 tensor construction helpers.
1254 * @{
1255 */
1256template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
1257 int Tensor_Dim3, int S, std::size_t... Is>
1258static inline auto makeFTensor4FromPtrArray(
1259 const std::array<T *, Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3>
1260 &ptrs,
1261 std::index_sequence<Is...>) {
1262 return FTensor::Tensor4<FTensor::PackPtr<T *, S>, Tensor_Dim0, Tensor_Dim1,
1263 Tensor_Dim2, Tensor_Dim3>(ptrs[Is]...);
1264}
1265
1266template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
1267 int S, class DL, class M>
1268struct GetFTensor4FromMatImpl {
1269 static inline auto get(M &data, int rr = 0, int cc = 0) {
1270 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1271 constexpr int tensor_size =
1272 Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3;
1273 constexpr int stride =
1274 S == -1 ? (DL::isGaussByCoeffs ? tensor_size : 1) : S;
1275 if constexpr (DL::isGaussByCoeffs) {
1276 static_assert(
1277 stride % tensor_size == 0,
1278 "getFTensor4FromMat: stride should be a multiple of tensor_size");
1279
1280#ifndef NDEBUG
1281 if (data.size2() != tensor_size) {
1282 THROW_MESSAGE("getFTensor4FromMat<" +
1283 boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
1284 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
1285 boost::lexical_cast<std::string>(Tensor_Dim2) + "," +
1286 boost::lexical_cast<std::string>(Tensor_Dim3) +
1287 ">: wrong size of columns in data matrix, should be " +
1288 boost::lexical_cast<std::string>(tensor_size) +
1289 " but is " +
1290 boost::lexical_cast<std::string>(data.size2()));
1291 }
1292 auto *first = &data(rr + 0, cc + 0);
1293 auto *last = &data(rr + 0, cc + 0 + tensor_size - 1);
1294 if (last - first != tensor_size - 1) {
1295 THROW_MESSAGE("getFTensor4FromMat<" +
1296 boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
1297 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
1298 boost::lexical_cast<std::string>(Tensor_Dim2) + "," +
1299 boost::lexical_cast<std::string>(Tensor_Dim3) +
1300 ">: row slice is not contiguous");
1301 }
1302#endif
1303 return FTensor::Tensor4<FTensor::CursorPtr<T *, stride>, Tensor_Dim0,
1304 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>{
1305 &data(rr + 0, cc + 0)};
1306 } else if constexpr (DL::isCoeffsByGauss) {
1307 std::array<T *, tensor_size> ptrs;
1308#ifndef NDEBUG
1309 if (data.size1() != tensor_size) {
1310 THROW_MESSAGE("getFTensor4FromMat<" +
1311 boost::lexical_cast<std::string>(Tensor_Dim0) + "," +
1312 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
1313 boost::lexical_cast<std::string>(Tensor_Dim2) + "," +
1314 boost::lexical_cast<std::string>(Tensor_Dim3) +
1315 ">: wrong size of rows in data matrix, should be " +
1316 boost::lexical_cast<std::string>(tensor_size) +
1317 " but is " +
1318 boost::lexical_cast<std::string>(data.size1()));
1319 }
1320#endif
1321 for (auto i = 0; i != tensor_size; ++i)
1322 ptrs[i] = &data(rr + i, cc + 0);
1323 return makeFTensor4FromPtrArray<T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
1324 Tensor_Dim3, stride>(
1325 ptrs, std::make_index_sequence<tensor_size>{});
1326 } else {
1327 static_assert(!std::is_same<M, M>::value,
1328 "Unsupported data layout for getFTensor4FromMatImpl");
1329 }
1330 }
1331};
1332
1333/**
1334 * @brief Get tensor rank 4 (non symmetric) form data matrix
1335 *
1336 * @tparam Tensor_Dim0 dimension of first index
1337 * @tparam Tensor_Dim1 dimension of second index
1338 * @tparam Tensor_Dim2 dimension of third index
1339 * @tparam Tensor_Dim3 dimension of fourth index
1340 * @tparam S Memory shift
1341 * @tparam DL Data layout
1342 * @tparam M Matrix type
1343 * @param data data container
1344 * @return FTensor tensor rank 4 view matching the selected data layout
1345 */
1346template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
1347 int S = -1, class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
1348 class M = MatrixDouble>
1349static inline auto getFTensor4FromMat(M &data) {
1350 return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
1351 Tensor_Dim3, S, DL, M>::get(data);
1352}
1353
1354template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
1355 int S = -1, class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
1356 class M = MatrixDouble>
1357using GetFTensor4FromMatType =
1358 decltype(GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
1359 Tensor_Dim3, S, DL, M>::get(
1360 std::declval<M &>(), 0, 0));
1361
1362template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
1363 int S = -1, class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
1364static inline auto getFTensor4FromMat(boost::weak_ptr<MatrixDouble> data_ptr) {
1365#ifndef NDEBUG
1366 if (!data_ptr.lock()) {
1367 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
1368 "Data matrix is not available");
1369 }
1370#endif
1371 return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
1372 Tensor_Dim3, S, DL,
1373 MatrixDouble>::get(*data_ptr.lock());
1374}
1375
1376template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
1377 int S = -1, class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
1378static inline auto
1379getFTensor4FromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
1380#ifndef NDEBUG
1381 if (!data_ptr) {
1382 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
1383 "Data matrix is not available");
1384 }
1385#endif
1386 return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
1387 Tensor_Dim3, S, DL,
1388 MatrixDouble>::get(*data_ptr);
1389}
1390
1391template <int DIM1, int DIM2, int DIM3, int DIM4, int S = -1, class T = double>
1392inline auto getFTensor4FromPtr(T *ptr) {
1393 auto matrix_data = getMatrixAdaptor(ptr, 1, DIM1 * DIM2 * DIM3 * DIM4);
1394 return GetFTensor4FromMatImpl<DIM1, DIM2, DIM3, DIM4, S,
1395 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1396 decltype(matrix_data)>::get(matrix_data);
1397}
1398/** @} */
1399
1400/** @name Ddg functions
1401 * Helpers for fourth-order tensors symmetric in index pairs.
1402 * @{
1403 */
1404template <int Tensor_Dim01, int Tensor_Dim23, int S, class DL, class M>
1405struct GetFTensor4DdgFromMatImpl {};
1406
1407template <int Tensor_Dim01, int Tensor_Dim23, int S, class M>
1408struct GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S,
1409 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1410 M> {
1411 static inline auto get(M &data, int rr = 0, int cc = 0) {
1412 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1413 constexpr int symm_size_01 = (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2;
1414 constexpr int symm_size_23 = (Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2;
1415 constexpr int ddg_size = symm_size_01 * symm_size_23;
1416 constexpr int stride = S == -1 ? ddg_size : S;
1417 static_assert(
1418 stride % ddg_size == 0,
1419 "Stride should be a multiple of the size of the symmetric tensor");
1420
1421#ifndef NDEBUG
1422 if (data.size2() != ddg_size) {
1423 THROW_MESSAGE(
1424 "getFTensor4DdgFromMat<" +
1425 boost::lexical_cast<std::string>(Tensor_Dim01) + ", " +
1426 boost::lexical_cast<std::string>(Tensor_Dim23) +
1427 ">: wrong size of data matrix, number of columns should be " +
1428 boost::lexical_cast<std::string>(ddg_size) + " but is " +
1429 boost::lexical_cast<std::string>(data.size2()));
1430 }
1431 auto *first = &data(rr + 0, cc + 0);
1432 auto *last = &data(rr + 0, cc + 0 + ddg_size - 1);
1433 if (last - first != ddg_size - 1)
1434 THROW_MESSAGE("getFTensor4DdgFromMat<" +
1435 boost::lexical_cast<std::string>(Tensor_Dim01) + ", " +
1436 boost::lexical_cast<std::string>(Tensor_Dim23) +
1437 ">: row slice is not contiguous");
1438#endif
1439 return FTensor::Ddg<FTensor::CursorPtr<T *, stride>, Tensor_Dim01,
1440 Tensor_Dim23>{&data(rr + 0, cc + 0)};
1441 }
1442};
1443
1444template <int S, class M>
1445struct GetFTensor4DdgFromMatImpl<
1446 1, 1, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
1447 static inline auto get(M &data, int rr = 0, int cc = 0) {
1448 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1449 constexpr int stride = S == -1 ? 1 : S;
1450#ifndef NDEBUG
1451 constexpr int ddg_size = 1;
1452 if (data.size1() != ddg_size)
1453 THROW_MESSAGE(
1454 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
1455 "of rows should be 1 but is " +
1456 boost::lexical_cast<std::string>(data.size1()));
1457#endif
1458 return FTensor::Ddg<FTensor::PackPtr<T *, stride>, 1, 1>{
1459 &data(rr + 0, cc + 0)};
1460 }
1461};
1462
1463template <int S, class M>
1464struct GetFTensor4DdgFromMatImpl<
1465 2, 2, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
1466 static inline auto get(M &data, int rr = 0, int cc = 0) {
1467 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1468 constexpr int stride = S == -1 ? 1 : S;
1469#ifndef NDEBUG
1470 constexpr int ddg_size = 9;
1471 if (data.size1() != ddg_size) {
1472 THROW_MESSAGE(
1473 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
1474 "of rows should be 9 but is " +
1475 boost::lexical_cast<std::string>(data.size1()));
1476 }
1477#endif
1478 return FTensor::Ddg<FTensor::PackPtr<T *, stride>, 2, 2>{
1479 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
1480 &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0),
1481 &data(rr + 6, cc + 0), &data(rr + 7, cc + 0), &data(rr + 8, cc + 0)};
1482 }
1483};
1484
1485template <int S, class M>
1486struct GetFTensor4DdgFromMatImpl<
1487 3, 3, S, DataLayoutTraits<DataLayout::CoeffsByGauss>, M> {
1488 static inline auto get(M &data, int rr = 0, int cc = 0) {
1489 using T = std::remove_cv_t<std::remove_reference_t<decltype(data(0, 0))>>;
1490 constexpr int stride = S == -1 ? 1 : S;
1491#ifndef NDEBUG
1492 constexpr int ddg_size = 36;
1493 if (data.size1() != ddg_size) {
1494 THROW_MESSAGE(
1495 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
1496 "of rows should be 36 but is " +
1497 boost::lexical_cast<std::string>(data.size1()));
1498 }
1499#endif
1500 return FTensor::Ddg<FTensor::PackPtr<T *, stride>, 3, 3>{
1501 &data(rr + 0, cc + 0), &data(rr + 1, cc + 0), &data(rr + 2, cc + 0),
1502 &data(rr + 3, cc + 0), &data(rr + 4, cc + 0), &data(rr + 5, cc + 0),
1503 &data(rr + 6, cc + 0), &data(rr + 7, cc + 0), &data(rr + 8, cc + 0),
1504 &data(rr + 9, cc + 0), &data(rr + 10, cc + 0), &data(rr + 11, cc + 0),
1505 &data(rr + 12, cc + 0), &data(rr + 13, cc + 0), &data(rr + 14, cc + 0),
1506 &data(rr + 15, cc + 0), &data(rr + 16, cc + 0), &data(rr + 17, cc + 0),
1507 &data(rr + 18, cc + 0), &data(rr + 19, cc + 0), &data(rr + 20, cc + 0),
1508 &data(rr + 21, cc + 0), &data(rr + 22, cc + 0), &data(rr + 23, cc + 0),
1509 &data(rr + 24, cc + 0), &data(rr + 25, cc + 0), &data(rr + 26, cc + 0),
1510 &data(rr + 27, cc + 0), &data(rr + 28, cc + 0), &data(rr + 29, cc + 0),
1511 &data(rr + 30, cc + 0), &data(rr + 31, cc + 0), &data(rr + 32, cc + 0),
1512 &data(rr + 33, cc + 0), &data(rr + 34, cc + 0), &data(rr + 35, cc + 0)};
1513 }
1514};
1515
1516/**
1517 * @brief Get symmetric tensor rank 4 on first two and last indices from
1518 * form data matrix
1519 *
1520 * @tparam Tensor_Dim01 dimension of first two indicies
1521 * @tparam Tensor_Dim23 dimension of second two indicies
1522 * @tparam Memory shift
1523 * @tparam M matrix type
1524 * @param data data container
1525 * @return FTensor::Ddg<FTensor::PackPtr<T *, S>, Tensor_Dim01, TensorDim23>
1526 */
1527template <int Tensor_Dim01, int Tensor_Dim23, int S = -1,
1528 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
1529 class M = MatrixDouble>
1534
1535template <int Tensor_Dim01, int Tensor_Dim23, int S = -1,
1536 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>,
1537 class M = MatrixDouble>
1538using GetFTensor4DdgFromMatType =
1539 decltype(GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, DL,
1540 M>::get(std::declval<M &>(), 0, 0));
1541
1542template <int Tensor_Dim01, int Tensor_Dim23, int S = -1,
1543 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
1544static inline auto getFTensor4DdgFromMat(boost::weak_ptr<MatrixDouble> data) {
1545#ifndef NDEBUG
1546 if (!data.lock()) {
1547 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
1548 "Data matrix is not available");
1549 }
1550#endif
1551 return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, DL,
1552 MatrixDouble>::get(*data.lock());
1553}
1554
1555template <int Tensor_Dim01, int Tensor_Dim23, int S = -1,
1556 class DL = DataLayoutTraits<DataLayout::GaussByCoeffs>>
1557static inline auto
1558getFTensor4DdgFromMat(boost::shared_ptr<MatrixDouble> data_ptr) {
1559#ifndef NDEBUG
1560 if (!data_ptr) {
1561 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL,
1562 "Data matrix is not available");
1563 }
1564#endif
1565 return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, DL,
1566 MatrixDouble>::get(*data_ptr);
1567}
1568
1569template <int Tensor_Dim01, int Tensor_Dim23, int S, class T = double>
1570static inline auto getFTensor4DdgFromPtr(T *ptr) {
1571 constexpr auto symm_size01 = (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2;
1572 constexpr auto symm_size23 = (Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2;
1573 constexpr auto symm_size = symm_size01 * symm_size23;
1574 auto matrix_data = getMatrixAdaptor(ptr, 1, symm_size);
1575 return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S,
1576 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1577 decltype(matrix_data)>::get(matrix_data);
1578}
1579/** @} */
1580
1581/** @name Matrix sizeing helpers
1582 * Matrix size helpers
1583 * @{
1584*/
1585
1586template <class Tensor, class DL> struct MatrixSizeHelper;
1587
1588template <class T, int S, int Tensor_Dim>
1589struct MatrixSizeHelper<
1590 FTensor::Tensor1<FTensor::CursorPtr<T *, S>, Tensor_Dim>,
1591 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1592 template <class M>
1593 static inline auto size(M &m, const size_t nb_integration_pts) {
1594 if (m.size1() != nb_integration_pts || m.size2() != Tensor_Dim) {
1595 m.resize(nb_integration_pts, Tensor_Dim, false);
1596 }
1597 return [&m]() {
1598 return GetFTensor1FromMatImpl<Tensor_Dim, S,
1599 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1600 M>::get(m, 0, 0);
1601 };
1602 }
1603 template <class M>
1604 static inline auto get(M &m, const size_t nb_integration_pts) {
1605#ifndef NDEBUG
1606 if (S == -1 || S == Tensor_Dim) {
1607 if (m.size1() != nb_integration_pts) {
1608 CHK_THROW_MESSAGE(
1609 MOFEM_DATA_INCONSISTENCY,
1610 "Matrix size is not compatible with the number of "
1611 "integration points: expected " +
1612 boost::lexical_cast<std::string>(nb_integration_pts) +
1613 " but got " + boost::lexical_cast<std::string>(m.size1()));
1614 }
1615 }
1616#endif
1617 return [&m]() {
1618 return GetFTensor1FromMatImpl<Tensor_Dim, S,
1619 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1620 M>::get(m, 0, 0);
1621 };
1622 }
1623};
1624
1625template <class T, int S, int Tensor_Dim0, int Tensor_Dim1>
1626struct MatrixSizeHelper<
1627 FTensor::Tensor2<FTensor::CursorPtr<T *, S>, Tensor_Dim0, Tensor_Dim1>,
1628 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1629 template <class M>
1630 static inline auto size(M &m, const size_t nb_integration_pts) {
1631 if (m.size1() != nb_integration_pts ||
1632 m.size2() != Tensor_Dim0 * Tensor_Dim1) {
1633 m.resize(nb_integration_pts, Tensor_Dim0 * Tensor_Dim1, false);
1634 }
1635 return [&m]() {
1636 return GetFTensor2FromMatImpl<
1637 Tensor_Dim0, Tensor_Dim1, S,
1638 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1639 };
1640 }
1641 template <class M>
1642 static inline auto get(M &m, const size_t nb_integration_pts) {
1643#ifndef NDEBUG
1644 if (S == -1 || S == Tensor_Dim0 * Tensor_Dim1) {
1645 if (m.size1() != nb_integration_pts) {
1646 CHK_THROW_MESSAGE(
1647 MOFEM_DATA_INCONSISTENCY,
1648 "Matrix size is not compatible with the number of "
1649 "integration points: expected " +
1650 boost::lexical_cast<std::string>(nb_integration_pts) +
1651 " but got " + boost::lexical_cast<std::string>(m.size1()));
1652 }
1653 }
1654#endif
1655 return [&m]() {
1656 return GetFTensor2FromMatImpl<
1657 Tensor_Dim0, Tensor_Dim1, S,
1658 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1659 };
1660 }
1661};
1662
1663template <class T, int S, int Tensor_Dim>
1664struct MatrixSizeHelper<
1665 FTensor::Tensor2_symmetric<FTensor::CursorPtr<T *, S>, Tensor_Dim>,
1666 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1667 template <class M>
1668 static inline auto size(M &m, const size_t nb_integration_pts) {
1669 constexpr int symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1670 if (m.size1() != nb_integration_pts || m.size2() != symm_size) {
1671 m.resize(nb_integration_pts, symm_size, false);
1672 }
1673 return [&m]() {
1674 return GetFTensor2SymmetricFromMatImpl<
1675 Tensor_Dim, S, DataLayoutTraits<DataLayout::GaussByCoeffs>,
1676 M>::get(m, 0, 0);
1677 };
1678 }
1679 template <class M>
1680 static inline auto get(M &m, const size_t nb_integration_pts) {
1681#ifndef NDEBUG
1682 constexpr int symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1683 if (S == -1 || S == symm_size) {
1684 if (m.size1() != nb_integration_pts) {
1685 CHK_THROW_MESSAGE(
1686 MOFEM_DATA_INCONSISTENCY,
1687 "Matrix size is not compatible with the number of "
1688 "integration points: expected " +
1689 boost::lexical_cast<std::string>(nb_integration_pts) +
1690 " but got " + boost::lexical_cast<std::string>(m.size1()));
1691 }
1692 }
1693#endif
1694 return [&m]() {
1695 return GetFTensor2SymmetricFromMatImpl<
1696 Tensor_Dim, S, DataLayoutTraits<DataLayout::GaussByCoeffs>,
1697 M>::get(m, 0, 0);
1698 };
1699 }
1700};
1701
1702template <class T, int S, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
1703struct MatrixSizeHelper<
1704 FTensor::Tensor3<FTensor::CursorPtr<T *, S>, Tensor_Dim0, Tensor_Dim1,
1705 Tensor_Dim2>,
1706 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1707 template <class M>
1708 static inline auto size(M &m, const size_t nb_integration_pts) {
1709 constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
1710 if (m.size1() != nb_integration_pts || m.size2() != tensor_size) {
1711 m.resize(nb_integration_pts, tensor_size, false);
1712 }
1713 return [&m]() {
1714 return GetFTensor3FromMatImpl<
1715 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
1716 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1717 };
1718 }
1719 template <class M>
1720 static inline auto get(M &m, const size_t nb_integration_pts) {
1721#ifndef NDEBUG
1722 constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
1723 if (S == -1 || S == tensor_size) {
1724 if (m.size1() != nb_integration_pts) {
1725 CHK_THROW_MESSAGE(
1726 MOFEM_DATA_INCONSISTENCY,
1727 "Matrix size is not compatible with the number of "
1728 "integration points: expected " +
1729 boost::lexical_cast<std::string>(nb_integration_pts) +
1730 " but got " + boost::lexical_cast<std::string>(m.size1()));
1731 }
1732 }
1733#endif
1734 return [&m]() {
1735 return GetFTensor3FromMatImpl<
1736 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
1737 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1738 };
1739 }
1740};
1741
1742template <class T, int S, int Tensor_Dim01, int Tensor_Dim2>
1743struct MatrixSizeHelper<
1744 FTensor::Dg<FTensor::CursorPtr<T *, S>, Tensor_Dim01, Tensor_Dim2>,
1745 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1746 template <class M>
1747 static inline auto size(M &m, const size_t nb_integration_pts) {
1748 constexpr int dg_size =
1749 ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) * Tensor_Dim2;
1750 if (m.size1() != nb_integration_pts || m.size2() != dg_size) {
1751 m.resize(nb_integration_pts, dg_size, false);
1752 }
1753 return [&m]() {
1754 return GetFTensor3DgFromMatImpl<
1755 Tensor_Dim01, Tensor_Dim2, S,
1756 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1757 };
1758 }
1759 template <class M>
1760 static inline auto get(M &m, const size_t nb_integration_pts) {
1761#ifndef NDEBUG
1762 constexpr int dg_size =
1763 ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) * Tensor_Dim2;
1764 if (S == -1 || S == dg_size) {
1765 if (m.size1() != nb_integration_pts) {
1766 CHK_THROW_MESSAGE(
1767 MOFEM_DATA_INCONSISTENCY,
1768 "Matrix size is not compatible with the number of "
1769 "integration points: expected " +
1770 boost::lexical_cast<std::string>(nb_integration_pts) +
1771 " but got " + boost::lexical_cast<std::string>(m.size1()));
1772 }
1773 }
1774#endif
1775 return [&m]() {
1776 return GetFTensor3DgFromMatImpl<
1777 Tensor_Dim01, Tensor_Dim2, S,
1778 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1779 };
1780 }
1781};
1782
1783template <class T, int S, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
1784 int Tensor_Dim3>
1785struct MatrixSizeHelper<
1786 FTensor::Tensor4<FTensor::CursorPtr<T *, S>, Tensor_Dim0, Tensor_Dim1,
1787 Tensor_Dim2, Tensor_Dim3>,
1788 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1789 template <class M>
1790 static inline auto size(M &m, const size_t nb_integration_pts) {
1791 constexpr int tensor_size =
1792 Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3;
1793 if (m.size1() != nb_integration_pts || m.size2() != tensor_size) {
1794 m.resize(nb_integration_pts, tensor_size, false);
1795 }
1796 return [&m]() {
1797 return GetFTensor4FromMatImpl<
1798 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S,
1799 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1800 };
1801 }
1802 template <class M>
1803 static inline auto get(M &m, const size_t nb_integration_pts) {
1804#ifndef NDEBUG
1805 if (S == -1 || S == Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3) {
1806 if (m.size1() != nb_integration_pts) {
1807 CHK_THROW_MESSAGE(
1808 MOFEM_DATA_INCONSISTENCY,
1809 "Matrix size is not compatible with the number of "
1810 "integration points: expected " +
1811 boost::lexical_cast<std::string>(nb_integration_pts) +
1812 " but got " + boost::lexical_cast<std::string>(m.size1()));
1813 }
1814 }
1815#endif
1816 return [&m]() {
1817 return GetFTensor4FromMatImpl<
1818 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S,
1819 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1820 };
1821 }
1822};
1823
1824template <class T, int S, int Tensor_Dim01, int Tensor_Dim23>
1825struct MatrixSizeHelper<
1826 FTensor::Ddg<FTensor::CursorPtr<T *, S>, Tensor_Dim01, Tensor_Dim23>,
1827 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1828 template <class M>
1829 static inline auto size(M &m, const size_t nb_integration_pts) {
1830 constexpr int ddg_size =
1831 ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) *
1832 ((Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2);
1833 if (m.size1() != nb_integration_pts || m.size2() != ddg_size) {
1834 m.resize(nb_integration_pts, ddg_size, false);
1835 }
1836 return [&m]() {
1837 return GetFTensor4DdgFromMatImpl<
1838 Tensor_Dim01, Tensor_Dim23, S,
1839 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1840 };
1841 }
1842 template <class M>
1843 static inline auto get(M &m, const size_t nb_integration_pts) {
1844#ifndef NDEBUG
1845 constexpr int ddg_size =
1846 ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) *
1847 ((Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2);
1848 if (S == -1 || S == ddg_size) {
1849
1850 if (m.size1() != nb_integration_pts) {
1851 CHK_THROW_MESSAGE(
1852 MOFEM_DATA_INCONSISTENCY,
1853 "Matrix size is not compatible with the number of "
1854 "integration points: expected " +
1855 boost::lexical_cast<std::string>(nb_integration_pts) +
1856 " but got " + boost::lexical_cast<std::string>(m.size1()));
1857 }
1858 }
1859#endif
1860 return [&m]() {
1861 return GetFTensor4DdgFromMatImpl<
1862 Tensor_Dim01, Tensor_Dim23, S,
1863 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1864 };
1865 }
1866};
1867
1868/** @} */
1869
1870/** @name Voigt notation helpers
1871 * Conversions between tensor views and Voigt-like storage.
1872 * @{
1873 */
1874// Template functions for conversion to Voigt notation for 2nd order tensors
1875
1876template <int DIM, typename T>
1877inline auto getVoigtVec(T &t_mat) {
1878 if constexpr (DIM == 3) {
1879 return std::array<double, 9>{t_mat(0, 0), t_mat(1, 1), t_mat(2, 2),
1880 t_mat(0, 1), t_mat(1, 0), t_mat(0, 2),
1881 t_mat(2, 0), t_mat(1, 2), t_mat(2, 1)};
1882 } else {
1883 return std::array<double, 5>{t_mat(0, 0), t_mat(1, 1), 1.0, t_mat(0, 1),
1884 t_mat(1, 0)};
1885 }
1886}
1887
1888template <int DIM, typename T> inline auto getVoigtVecSymm(T &t_mat) {
1889 constexpr double sqr2 = boost::math::constants::root_two<double>();
1890
1891 if constexpr (DIM == 3) {
1892 return std::array<double, 9>{t_mat(0, 0),
1893 t_mat(1, 1),
1894 t_mat(2, 2),
1895 sqr2 * t_mat(0, 1),
1896 sqr2 * t_mat(0, 2),
1897 sqr2 * t_mat(1, 2),
1898 0.0,
1899 0.0,
1900 0.0};
1901 } else {
1902 return std::array<double, 4>{t_mat(0, 0), t_mat(1, 1), 0.0,
1903 sqr2 * t_mat(0, 1)};
1904 }
1905}
1906
1907template <typename T>
1908inline auto getVoigtVecAxisymm(T &t_mat, const double hoop_term) {
1909 array<double, 5> vec{t_mat(0, 0), t_mat(1, 1), 1. + hoop_term, t_mat(0, 1),
1910 t_mat(1, 0)};
1911
1912 return vec;
1913}
1914
1915template <typename T>
1916inline auto getVoigtVecSymmAxisymm(T &t_mat, const double hoop_term) {
1917 constexpr double sqr2 = boost::math::constants::root_two<double>();
1918
1919 array<double, 4> vec_sym{t_mat(0, 0), t_mat(1, 1), hoop_term,
1920 sqr2 * t_mat(0, 1)};
1921 return vec_sym;
1922}
1923/** @} */
1924
1925/** @name Linear algebra helpers
1926 * LAPACK-backed matrix and eigenvalue utilities.
1927 * @{
1928 */
1929// list of lapack wrappers
1930/**
1931 * @brief compute matrix inverse with lapack dgetri
1932 *
1933 * @param mat input square matrix / output inverse matrix
1934 * @return MoFEMErrorCode
1935 */
1936inline MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat) {
1937 MoFEMFunctionBegin;
1938
1939 const size_t M = mat.size1();
1940 const size_t N = mat.size2();
1941
1942 if (M != N)
1943 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1944 "The input matrix for inverse computation is not square %zu != %zu",
1945 M, N);
1946
1947 int *ipv = new int[N];
1948 int lwork = N * N;
1949 double *work = new double[lwork];
1950 int info;
1951 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1952 if (info != 0)
1953 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "lapack error info = %d",
1954 info);
1955 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1956 if (info != 0)
1957 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1958 "lapack error info = %d", info);
1959
1960 delete[] ipv;
1961 delete[] work;
1962
1963 MoFEMFunctionReturn(0);
1964}
1965/**
1966 * @brief solve linear system with lapack dgesv
1967 *
1968 * @param mat input lhs square matrix / output L and U from the factorization
1969 * @param f input rhs vector / output solution vector
1970 * @return MoFEMErrorCode
1971 */
1972inline MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f) {
1973 MoFEMFunctionBegin;
1974
1975 const size_t M = mat.size1();
1976 const size_t N = mat.size2();
1977
1978 if (M == 0 || M != N)
1979 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1980 "The input matrix for inverse computation is not square %zu != %zu",
1981 M, N);
1982 if (f.size() != M)
1983 f.resize(M, false);
1984
1985 const int nrhs = 1;
1986 int info;
1987 int *ipiv = new int[M];
1988 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1989 &*f.data().begin(), M);
1990 if (info != 0) {
1991 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1992 "error lapack solve dgesv info = %d", info);
1993 }
1994
1995 delete[] ipiv;
1996 MoFEMFunctionReturn(0);
1997}
1998
1999/**
2000 * @brief Solve linear system of equations using Lapack
2001 *
2002 * @param mat
2003 * @param f
2004 * @return MoFEMErrorCode
2005 */
2006inline MoFEMErrorCode solveLinearSystem(const MatrixDouble &mat,
2007 VectorDouble &f) {
2008 MoFEMFunctionBegin;
2009 // copy matrix since on output lapack returns factorisation
2010 auto mat_copy = mat;
2011 CHKERR solveLinearSystem(mat_copy, f);
2012 MoFEMFunctionReturn(0);
2013}
2014
2015/**
2016 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
2017 *
2018 * @param mat input symmetric matrix
2019 * @param eig output eigen values sorted
2020 * @param eigen_vec output matrix of row eigen vectors
2021 * @return MoFEMErrorCode
2022 */
2023inline MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat,
2024 VectorDouble &eig,
2025 MatrixDouble &eigen_vec) {
2026 MoFEMFunctionBegin;
2027
2028 const size_t M = mat.size1();
2029 const size_t N = mat.size2();
2030
2031 if (M == 0 || M != N)
2032 SETERRQ(
2033 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2034 "The input matrix for eigen value computation is not square %zu != %zu",
2035 M, N);
2036 if (eig.size() != M)
2037 eig.resize(M, false);
2038
2039 eigen_vec = mat;
2040 const int n = M;
2041 const int lda = M;
2042 const int size = (M + 2) * M;
2043 int lwork = size;
2044 double *work = new double[size];
2045
2046 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
2047 &*eig.data().begin(), work, lwork) > 0)
2048 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
2049 "The algorithm failed to compute eigenvalues.");
2050
2051 delete[] work;
2052 MoFEMFunctionReturn(0);
2053}
2054/**
2055 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
2056 *
2057 * @tparam DIM
2058 * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
2059 * @param eig output eigen values sorted
2060 * @return MoFEMErrorCode
2061 */
2062template <int DIM, typename T1, typename T2>
2063inline MoFEMErrorCode
2064computeEigenValuesSymmetric(FTensor::Tensor2<T1, DIM, DIM> &eigen_vec,
2065 FTensor::Tensor1<T2, DIM> &eig) {
2066 MoFEMFunctionBegin;
2067
2068 const int n = DIM;
2069 const int lda = DIM;
2070 const int lwork = (DIM + 2) * DIM;
2071 std::array<double, (DIM + 2) * DIM> work;
2072
2073 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
2074 lwork) > 0)
2075 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
2076 "The algorithm failed to compute eigenvalues.");
2077 MoFEMFunctionReturn(0);
2078}
2079
2080/**
2081 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
2082 *
2083 * @tparam DIM
2084 * @param mat input tensor pointer of size DIM x DIM
2085 * @param eig output eigen values sorted
2086 * @param eigen_vec output matrix of row eigen vectors
2087 * @return MoFEMErrorCode
2088 */
2089template <int DIM, typename T1, typename T2, typename T3>
2090inline MoFEMErrorCode
2091computeEigenValuesSymmetric(const FTensor::Tensor2_symmetric<T1, DIM> &mat,
2092 FTensor::Tensor1<T2, DIM> &eig,
2093 FTensor::Tensor2<T3, DIM, DIM> &eigen_vec) {
2094 MoFEMFunctionBegin;
2095 for (int ii = 0; ii != DIM; ii++)
2096 for (int jj = 0; jj != DIM; jj++)
2097 eigen_vec(ii, jj) = mat(ii, jj);
2098
2099 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
2100
2101 MoFEMFunctionReturn(0);
2102}
2103/** @} */
2104
2105/** @name Tensor algebra helpers
2106 * Determinant and inverse routines for tensor-valued objects.
2107 * @{
2108 */
2109/**
2110 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
2111 *
2112 * @tparam T
2113 * @param t
2114 * @return double
2115 */
2116template <typename T> static inline auto determinantTensor3by3(T &t) {
2117 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
2118 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
2119 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
2120}
2121
2122/**
2123 * @brief Calculate the determinant of a 2x2 matrix or a tensor of rank 2
2124 *
2125 * @tparam T
2126 * @param t
2127 * @return double
2128 */
2129template <typename T> static inline auto determinantTensor2by2(T &t) {
2130 return t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
2131}
2132
2133template <typename T, int DIM> struct DeterminantTensorImpl;
2134
2135template <typename T> struct DeterminantTensorImpl<T, 3> {
2136 static inline auto get(T &t) { return determinantTensor3by3(t); }
2137};
2138
2139template <typename T> struct DeterminantTensorImpl<T, 2> {
2140 static auto get(T &t) { return determinantTensor2by2(t); }
2141};
2142
2143/**
2144 * @brief Calculate the determinant of a tensor of rank DIM
2145 */
2146template <typename T, int DIM>
2150
2151/**
2152 * @brief Calculate the determinant of a tensor of rank DIM
2153 */
2154template <typename T, int DIM>
2158
2159/**
2160 * \brief Calculate inverse of tensor rank 2 at integration points
2161
2162 */
2163template <int Tensor_Dim, class T, class L, class A>
2164inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
2165 ublas::vector<T, A> &det_data,
2166 ublas::matrix<T, L, A> &inv_jac_data) {
2167 MoFEMFunctionBeginHot;
2168 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2169 "Specialization for this template not yet implemented");
2170 MoFEMFunctionReturnHot(0);
2171}
2172
2173template <>
2174inline MoFEMErrorCode
2175invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
2176 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
2177
2178/**
2179 * \brief Calculate determinant 3 by 3
2180
2181 */
2182template <class T1, class T2>
2188
2189/**
2190 * \brief Calculate determinant 2 by 2
2191
2192 */
2193template <class T1, class T2>
2199
2200/**
2201 * \brief Calculate matrix inverse 3 by 3
2202
2203 */
2204template <class T1, class T2, class T3>
2205inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
2206 MoFEMFunctionBeginHot;
2207 const auto a = t(0, 0), b = t(0, 1), c = t(0, 2);
2208 const auto d = t(1, 0), e = t(1, 1), f = t(1, 2);
2209 const auto g = t(2, 0), h = t(2, 1), i = t(2, 2);
2210 const auto inv_det = 1.0 / det;
2211 inv_t(0, 0) = (e * i - f * h) * inv_det;
2212 inv_t(0, 1) = (c * h - b * i) * inv_det;
2213 inv_t(0, 2) = (b * f - c * e) * inv_det;
2214 inv_t(1, 0) = (f * g - d * i) * inv_det;
2215 inv_t(1, 1) = (a * i - c * g) * inv_det;
2216 inv_t(1, 2) = (c * d - a * f) * inv_det;
2217 inv_t(2, 0) = (d * h - e * g) * inv_det;
2218 inv_t(2, 1) = (b * g - a * h) * inv_det;
2219 inv_t(2, 2) = (a * e - b * d) * inv_det;
2220 MoFEMFunctionReturnHot(0);
2221}
2222
2223/**
2224 * \brief Calculate matrix inverse 2 by 2
2225
2226 */
2227template <class T1, class T2, class T3>
2228inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
2229 MoFEMFunctionBeginHot;
2230 const auto a = t(0, 0);
2231 const auto b = t(0, 1);
2232 const auto c = t(1, 0);
2233 const auto d = t(1, 1);
2234 inv_t(0, 0) = d / det;
2235 inv_t(0, 1) = -b / det;
2236 inv_t(1, 0) = -c / det;
2237 inv_t(1, 1) = a / det;
2238 MoFEMFunctionReturnHot(0);
2239}
2240
2241/**
2242 * \brief Calculate matrix inverse, overload for symmetric tensor
2243
2244 */
2245template <typename T1, typename T2, typename T3>
2246inline MoFEMErrorCode
2247invertTensor3by3(FTensor::Tensor2_symmetric<T1, 3> &t, T2 &det,
2248 FTensor::Tensor2_symmetric<T3, 3> &inv_t) {
2249 MoFEMFunctionBeginHot;
2250 const auto a = t(0, 0);
2251 const auto b = t(0, 1);
2252 const auto c = t(0, 2);
2253 const auto d = t(1, 1);
2254 const auto e = t(1, 2);
2255 const auto f = t(2, 2);
2256 const auto inv_det = 1.0 / det;
2257 inv_t(0, 0) = (d * f - e * e) * inv_det;
2258 inv_t(0, 1) = (c * e - b * f) * inv_det;
2259 inv_t(0, 2) = (b * e - c * d) * inv_det;
2260 inv_t(1, 1) = (a * f - c * c) * inv_det;
2261 inv_t(1, 2) = (b * c - a * e) * inv_det;
2262 inv_t(2, 2) = (a * d - b * b) * inv_det;
2263 MoFEMFunctionReturnHot(0);
2264}
2265
2266/**
2267 * \brief Calculate matrix inverse, overload for symmetric tensor
2268
2269 */
2270template <typename T1, typename T2, typename T3>
2271inline MoFEMErrorCode
2272invertTensor2by2(FTensor::Tensor2_symmetric<T1, 2> &t, T2 &det,
2273 FTensor::Tensor2_symmetric<T3, 2> &inv_t) {
2274 MoFEMFunctionBeginHot;
2275 const auto a = t(0, 0);
2276 const auto b = t(0, 1);
2277 const auto d = t(1, 1);
2278 const auto inv_det = 1.0 / det;
2279 inv_t(0, 0) = d * inv_det;
2280 inv_t(0, 1) = -b * inv_det;
2281 inv_t(1, 1) = a * inv_det;
2282 MoFEMFunctionReturnHot(0);
2283}
2284
2285template <typename T1, typename T2, typename T3, int DIM>
2286struct InvertTensorImpl;
2287
2288template <typename T1, typename T2, typename T3>
2289struct InvertTensorImpl<T1, T2, T3, 3> {
2290 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
2291 return invertTensor3by3(t, det, inv_t);
2292 }
2293};
2294
2295template <typename T1, typename T2, typename T3>
2296struct InvertTensorImpl<T1, T2, T3, 2> {
2297 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
2298 return invertTensor2by2(t, det, inv_t);
2299 }
2300};
2301
2302template <typename T1, typename T2, typename T3, int DIM>
2303static inline MoFEMErrorCode
2310
2311template <typename T1, typename T2, typename T3, int DIM>
2312static inline MoFEMErrorCode
2319/** @} */
2320
2321/** @name Mesh and range helpers
2322 * Utility types and helpers for mesh entities, ranges, and handles.
2323 * @{
2324 */
2325/**
2326 * @brief Extract entity handle form multi-index container
2327 *
2328 */
2329struct RefEntExtractor {
2330 template <typename Iterator>
2331 static inline EntityHandle extract(const Iterator &it) {
2332 return (*it)->getEnt();
2333 }
2334};
2335
2336/**
2337 * @brief Insert ordered mofem multi-index into range
2338 *
2339 * \note Inserted range has to be ordered.
2340 *
2341 * \code
2342 * auto hi_rit = refEntsPtr->upper_bound(start);
2343 * auto hi_rit = refEntsPtr->upper_bound(end);
2344 * Range to_erase;
2345 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
2346 * \endcode
2347 *
2348 * @tparam Iterator
2349 * @param r
2350 * @param begin_iter
2351 * @param end_iter
2352 * @return moab::Range::iterator
2353 */
2354template <typename Extractor, typename Iterator>
2355moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
2356 Iterator end_iter) {
2357 moab::Range::iterator hint = r.begin();
2358 while (begin_iter != end_iter) {
2359 size_t j = 0;
2360 auto bi = Extractor::extract(begin_iter);
2361 Iterator pj = begin_iter;
2362 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
2363 ++pj;
2364 ++j;
2365 }
2366 hint = r.insert(hint, bi, bi + (j - 1));
2367 begin_iter = pj;
2368 }
2369 return hint;
2370};
2371
2372/**
2373 * @brief Do nothing, used to rebuild database
2374 *
2375 */
2376struct Modify_change_nothing {
2377 Modify_change_nothing() = default;
2378 template <typename T> inline void operator()(T &e) {}
2379};
2380
2381/**
2382 * @brief Template used to reconstruct multi-index
2383 *
2384 * @tparam MI multi-index
2385 * @tparam Modifier
2386 * @param mi
2387 * @param mo
2388 * @return MoFEMErrorCode
2389 */
2390template <typename MI, typename MO = Modify_change_nothing>
2391inline MoFEMErrorCode reconstructMultiIndex(const MI &mi,
2392 MO &&mo = Modify_change_nothing()) {
2393 MoFEMFunctionBegin;
2394 for (auto it = mi.begin(); it != mi.end(); ++it) {
2395 if (!const_cast<MI &>(mi).modify(it, mo))
2396 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2397 "Houston we have a problem");
2398 }
2399 MoFEMFunctionReturn(0);
2400}
2401
2402struct TempMeshset {
2403 TempMeshset(moab::Interface &moab) : moab(moab) {
2404 rval = moab.create_meshset(MESHSET_SET, meshset);
2405 MOAB_THROW(rval);
2406 }
2407 virtual ~TempMeshset() { delete_meshset(); }
2408 operator EntityHandle() const { return meshset; }
2409 auto get_ptr() { return &meshset; }
2410
2411private:
2412 void delete_meshset() {
2413 rval = moab.delete_entities(&meshset, 1);
2414 MOAB_THROW(rval);
2415 }
2416 EntityHandle meshset;
2417 moab::Interface &moab;
2418};
2419
2420/**
2421 * @brief Create smart pointer to temporary meshset
2422 *
2423 */
2424inline auto get_temp_meshset_ptr(moab::Interface &moab) {
2425 return boost::make_shared<TempMeshset>(moab);
2426};
2427
2428inline auto id_from_handle(const EntityHandle h) {
2429 return static_cast<EntityID>(h & MB_ID_MASK);
2430};
2431
2432/**
2433 * @brief get type from entity handle
2434 *
2435 */
2436inline auto type_from_handle(const EntityHandle h) {
2437 return static_cast<EntityType>(h >> MB_ID_WIDTH);
2438};
2439
2440/**
2441 * @brief get entity handle from type and id
2442 *
2443 */
2444inline auto ent_form_type_and_id(const EntityType type, const EntityID id) {
2445 return (static_cast<EntityHandle>(type) << MB_ID_WIDTH) | id;
2446};
2447
2448/**
2449 * @brief get entity dimension form handle
2450 *
2451 */
2452inline auto dimension_from_handle(const EntityHandle h) {
2453 return moab::CN::Dimension(type_from_handle(h));
2454};
2455
2456/**
2457 * @brief get entity type name from handle
2458 *
2459 */
2460inline auto type_name_from_handle(const EntityHandle h) {
2461 return moab::CN::EntityTypeName(type_from_handle(h));
2462};
2463
2464/**
2465 * @brief get field bit id from bit number
2466 *
2467 */
2468inline auto field_bit_from_bit_number(const int bit_number) {
2469 return BitFieldId().set(bit_number - 1);
2470};
2471
2472/**
2473 * @brief Insert ranges
2474 *
2475 * @tparam I
2476 * @param f
2477 * @param s
2478 * @param tester
2479 * @param inserter
2480 * @return auto
2481 */
2482template <typename I>
2483auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
2484 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
2485 MoFEMFunctionBegin;
2486
2487 auto first = f;
2488 while (first != s)
2489 if (tester(first)) {
2490
2491 auto second = first;
2492 ++second;
2493
2494 while (second != s) {
2495 if (tester(second))
2496 ++second;
2497 else
2498 break;
2499 }
2500
2501 CHKERR inserter(first, second);
2502
2503 first = second;
2504 if (first != s)
2505 ++first;
2506
2507 } else {
2508 ++first;
2509 }
2510
2511 MoFEMFunctionReturn(0);
2512}
2513
2514/**
2515 * @brief Create Array
2516 *
2517 * See:
2518 * <a
2519 * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
2520 * stack overflow</a>
2521 *
2522 * @tparam Dest
2523 * @tparam Arg
2524 * @param arg
2525 * @return constexpr auto
2526 */
2527template <typename Dest = void, typename... Arg>
2528constexpr auto make_array(Arg &&...arg) {
2529 if constexpr (std::is_same<void, Dest>::value)
2530 return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
2531 {std::forward<Arg>(arg)...}};
2532 else
2533 return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
2534}
2535/** @} */
2536
2537/** @name FTensor index aliases
2538 * Convenient aliases for commonly used FTensor indices.
2539 * @{
2540 */
2541template <int DIM> using i_FTIndex = FTensor::Index<'i', DIM>;
2542template <int DIM> using j_FTIndex = FTensor::Index<'j', DIM>;
2543template <int DIM> using k_FTIndex = FTensor::Index<'k', DIM>;
2544template <int DIM> using l_FTIndex = FTensor::Index<'l', DIM>;
2545template <int DIM> using m_FTIndex = FTensor::Index<'m', DIM>;
2546template <int DIM> using n_FTIndex = FTensor::Index<'n', DIM>;
2547template <int DIM> using o_FTIndex = FTensor::Index<'o', DIM>;
2548template <int DIM> using p_FTIndex = FTensor::Index<'p', DIM>;
2549template <int DIM> using x_FTIndex = FTensor::Index<'x', DIM>;
2550template <int DIM> using y_FTIndex = FTensor::Index<'y', DIM>;
2551template <int DIM> using I_FTIndex = FTensor::Index<'I', DIM>;
2552template <int DIM> using J_FTIndex = FTensor::Index<'J', DIM>;
2553template <int DIM> using K_FTIndex = FTensor::Index<'K', DIM>;
2554template <int DIM> using L_FTIndex = FTensor::Index<'L', DIM>;
2555template <int DIM> using M_FTIndex = FTensor::Index<'M', DIM>;
2556template <int DIM> using N_FTIndex = FTensor::Index<'N', DIM>;
2557template <int DIM> using O_FTIndex = FTensor::Index<'O', DIM>;
2558template <int DIM> using Z_FTIndex = FTensor::Index<'Z', DIM>;
2559template <int DIM> using Y_FTIndex = FTensor::Index<'Y', DIM>;
2560
2561#define FTENSOR_INDEX(DIM, I) I##_FTIndex<DIM> I;
2562
2563#define FTENSOR_INDEX_DECL(r, DIM, I) FTENSOR_INDEX(DIM, I)
2564
2565#define FTENSOR_INDEXES(DIM, ...) \
2566 BOOST_PP_SEQ_FOR_EACH(FTENSOR_INDEX_DECL, DIM, \
2567 BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__))
2568
2569template <typename T> struct parent_of;
2570/** @} */
2571} // namespace MoFEM
2572
2573#endif //__TEMPLATES_HPP__
type
std::string type
Definition JsonConfigManager.cpp:28
a
constexpr double a
Definition approx_sphere.cpp:30
FTensor::Tensor1::data
T data[Tensor_Dim]
Definition Tensor1_value.hpp:10
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
MB_ID_MASK
#define MB_ID_MASK
Definition definitions.h:247
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition definitions.h:557
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MB_ID_WIDTH
#define MB_ID_WIDTH
Definition definitions.h:240
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
HVEC1_1
@ HVEC1_1
Definition definitions.h:209
HVEC0_1
@ HVEC0_1
Definition definitions.h:208
HVEC1_0
@ HVEC1_0
Definition definitions.h:206
HVEC2_1
@ HVEC2_1
Definition definitions.h:210
HVEC1_2
@ HVEC1_2
Definition definitions.h:212
HVEC2_2
@ HVEC2_2
Definition definitions.h:213
HVEC2_0
@ HVEC2_0
Definition definitions.h:207
HVEC0_2
@ HVEC0_2
Definition definitions.h:211
HVEC0_0
@ HVEC0_0
Definition definitions.h:205
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition definitions.h:577
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:39
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
lapack_dgesv
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition lapack_wrap.h:188
lapack_dgetrf
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition lapack_wrap.h:169
lapack_dsyev
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition lapack_wrap.h:273
lapack_dgetri
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition lapack_wrap.h:197
DIM2
constexpr int DIM2
Definition level_set.cpp:22
DIM1
constexpr int DIM1
Definition level_set.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
FTensor
Tensors class implemented by Walter Landry.
Definition FTensor.hpp:51
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition Exceptions.hpp:74
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::BitFieldId
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition Types.hpp:42
MoFEM::Types::VectorShallowArrayAdaptor
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition Types.hpp:114
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM::Types::MatrixADouble
UBlasMatrix< adouble > MatrixADouble
Definition Types.hpp:80
MoFEM::Types::MatrixShallowArrayAdaptor
ublas::matrix< T, ublas::row_major, ublas::shallow_array_adaptor< T > > MatrixShallowArrayAdaptor
Definition Types.hpp:120
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition Types.hpp:68
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::getVectorAdaptor
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition Templates.hpp:49
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition Templates.hpp:2436
MoFEM::getFTensor4FromMat
static auto getFTensor4FromMat(M &data)
Get tensor rank 4 (non symmetric) form data matrix.
Definition Templates.hpp:1349
MoFEM::GetFTensor2SymmetricFromMatType
decltype(GetFTensor2SymmetricFromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2SymmetricFromMatType
Definition Templates.hpp:836
MoFEM::getVoigtVecAxisymm
auto getVoigtVecAxisymm(T &t_mat, const double hoop_term)
Definition Templates.hpp:1908
MoFEM::getFTensor2FromMat
auto getFTensor2FromMat(M &data)
Get tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:542
MoFEM::to_non_symm
auto to_non_symm(const FTensor::Tensor2_symmetric< T, DIM > &symm)
Definition Templates.hpp:909
MoFEM::getFTensor3FromMat
static auto getFTensor3FromMat(M &data, int rr=0, int cc=0)
Get tensor rank 3 (non symmetries) form data matrix.
Definition Templates.hpp:1003
MoFEM::getFTensor2SymmetricLowerFromPtr< 3 >
auto getFTensor2SymmetricLowerFromPtr< 3 >(double *ptr)
Definition Templates.hpp:893
MoFEM::getFTensor2SymmetricFromPtr
auto getFTensor2SymmetricFromPtr(double *ptr)
Definition Templates.hpp:866
MoFEM::GetFTensor4FromMatType
decltype(GetFTensor4FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4FromMatType
Definition Templates.hpp:1360
MoFEM::getFTensor2SymmetricLowerFromPtr< 2 >
auto getFTensor2SymmetricLowerFromPtr< 2 >(double *ptr)
Definition Templates.hpp:902
MoFEM::getFTensor3FromPtr< 3, 2, 2 >
FTensor::Tensor3< FTensor::PackPtr< double *, 12 >, 3, 2, 2 > getFTensor3FromPtr< 3, 2, 2 >(double *ptr)
Definition Templates.hpp:1060
MoFEM::getVoigtVec
auto getVoigtVec(T &t_mat)
Definition Templates.hpp:1877
MoFEM::reconstructMultiIndex
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition Templates.hpp:2391
MoFEM::getFTensor2FromVec
auto getFTensor2FromVec(VectorDouble &data)
Definition Templates.hpp:586
MoFEM::insertOrdered
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition Templates.hpp:2355
MoFEM::invertTensor2by2
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition Templates.hpp:2228
MoFEM::getFTensor2SymmetricFromMat
static auto getFTensor2SymmetricFromMat(M &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:827
MoFEM::computeMatrixInverse
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
Definition Templates.hpp:1936
MoFEM::getFTensor1FromArrayDiag
auto getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr)
Get FTensor1 from array.
MoFEM::getFTensor2HVecFromPtr< 3, 2 >
auto getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition Templates.hpp:632
MoFEM::id_from_handle
auto id_from_handle(const EntityHandle h)
Definition Templates.hpp:2428
MoFEM::getFTensor3DgFromPtr
static auto getFTensor3DgFromPtr(T *ptr)
Definition Templates.hpp:1242
MoFEM::getFTensor1FromArray
auto getFTensor1FromArray(T &data)
Definition Templates.hpp:433
MoFEM::dimension_from_handle
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
Definition Templates.hpp:2452
MoFEM::getFTensor1FromArrayDiag< 2, 2 >
auto getFTensor1FromArrayDiag< 2, 2 >(MatrixDouble &data, const size_t rr)
Definition Templates.hpp:462
MoFEM::solveLinearSystem
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
Definition Templates.hpp:1972
MoFEM::determinantTensor2by2
static auto determinantTensor2by2(T &t)
Calculate the determinant of a 2x2 matrix or a tensor of rank 2.
Definition Templates.hpp:2129
MoFEM::getFTensor2HVecFromPtr
auto getFTensor2HVecFromPtr(double *ptr)
Make Tensor2 for HVec base from pointer.
MoFEM::type_name_from_handle
auto type_name_from_handle(const EntityHandle h)
get entity type name from handle
Definition Templates.hpp:2460
MoFEM::getFTensor3DgFromMat
static auto getFTensor3DgFromMat(M &data)
Get symmetric tensor rank 3 on the first two indices from form data matrix.
Definition Templates.hpp:1201
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
auto getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition Templates.hpp:641
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition Templates.hpp:2164
MoFEM::GetFTensor4DdgFromMatType
decltype(GetFTensor4DdgFromMatImpl< Tensor_Dim01, Tensor_Dim23, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4DdgFromMatType
Definition Templates.hpp:1540
MoFEM::invertTensor3by3< 3, double, ublas::row_major, DoubleAllocator >
MoFEMErrorCode invertTensor3by3< 3, double, ublas::row_major, DoubleAllocator >(MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data)
Definition DataOperators.cpp:149
MoFEM::getFTensor1FromPtr
auto getFTensor1FromPtr(T *ptr)
Definition Templates.hpp:416
MoFEM::invertTensor
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
Definition Templates.hpp:2304
MoFEM::getFTensor4DdgFromPtr
static auto getFTensor4DdgFromPtr(T *ptr)
Definition Templates.hpp:1570
MoFEM::make_array
constexpr auto make_array(Arg &&...arg)
Create Array.
Definition Templates.hpp:2528
MoFEM::getMatrixAdaptor
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition Templates.hpp:75
MoFEM::getFTensor1FromMat
auto getFTensor1FromMat(M &data, int rr=0, int cc=0)
Get tensor rank 1 (vector) form data matrix.
Definition Templates.hpp:376
MoFEM::getFTensor4FromPtr
auto getFTensor4FromPtr(T *ptr)
Definition Templates.hpp:1392
MoFEM::getFTensor2FromArray
auto getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc=0)
Definition Templates.hpp:723
MoFEM::determinantTensor
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
Definition Templates.hpp:2147
MoFEM::field_bit_from_bit_number
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
Definition Templates.hpp:2468
MoFEM::getFTensor0FromPtr
static auto getFTensor0FromPtr(T *ptr)
Get tensor rank 0 (scalar) from pointer.
Definition Templates.hpp:222
MoFEM::toString
std::string toString(X x)
Definition Templates.hpp:129
MoFEM::GetFTensor1FromMatType
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
Definition Templates.hpp:385
MoFEM::getFTensor2SymmetricLowerFromPtr
auto getFTensor2SymmetricLowerFromPtr(double *ptr)
Make symmetric Tensor2 from pointer, taking lower triangle of matrix.
MoFEM::getFTensor4DdgFromMat
static auto getFTensor4DdgFromMat(M &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
Definition Templates.hpp:1530
MoFEM::getFTensor2FromPtr
auto getFTensor2FromPtr(double *ptr)
Definition Templates.hpp:597
MoFEM::getVoigtVecSymm
auto getVoigtVecSymm(T &t_mat)
Definition Templates.hpp:1888
MoFEM::makeFTensor4FromPtrArray
static auto makeFTensor4FromPtrArray(const std::array< T *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 *Tensor_Dim3 > &ptrs, std::index_sequence< Is... >)
Definition Templates.hpp:1258
MoFEM::getFTensor0FromMat
auto getFTensor0FromMat(M &data)
Get tensor rank 0 (scalar) form data vector.
Definition Templates.hpp:197
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition Templates.hpp:2424
MoFEM::GetFTensor3DgFromMatType
decltype(GetFTensor3DgFromMatImpl< Tensor_Dim01, Tensor_Dim2, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor3DgFromMatType
Definition Templates.hpp:1211
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition Templates.hpp:2023
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(V &data)
Get tensor rank 0 (scalar) form data vector.
Definition Templates.hpp:163
MoFEM::getFTensor3FromPtr< 3, 3, 3 >
FTensor::Tensor3< FTensor::PackPtr< double *, 27 >, 3, 3, 3 > getFTensor3FromPtr< 3, 3, 3 >(double *ptr)
Definition Templates.hpp:1068
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition Templates.hpp:2116
MoFEM::GetFTensor3FromMatType
decltype(GetFTensor3FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor3FromMatType
Definition Templates.hpp:1013
MoFEM::getVoigtVecSymmAxisymm
auto getVoigtVecSymmAxisymm(T &t_mat, const double hoop_term)
Definition Templates.hpp:1916
MoFEM::getFTensor3FromPtr
FTensor::Tensor3< FTensor::PackPtr< double *, DIM1 *DIM2 *DIM3 >, DIM1, DIM2, DIM3 > getFTensor3FromPtr(double *ptr)
Definition Templates.hpp:1053
MoFEM::GetFTensor2FromMatType
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
Definition Templates.hpp:551
MoFEM::rangeInserter
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
Definition Templates.hpp:2483
MoFEM::getFTensor1FromArrayDiag< 3, 3 >
auto getFTensor1FromArrayDiag< 3, 3 >(MatrixDouble &data, const size_t rr)
Definition Templates.hpp:469
MoFEM::ent_form_type_and_id
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
Definition Templates.hpp:2444
I
constexpr IntegrationType I
Definition operators_tests.cpp:31
h
double h
Definition photon_diffusion.cpp:60
t
constexpr double t
plate stiffness
Definition plate.cpp:58
g
constexpr double g
Definition shallow_wave.cpp:63
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
N
const int N
Definition speed_test.cpp:3
MoFEM::EqBit::operator()
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition Templates.hpp:118
MoFEM::GetFTensor0FromMatImpl::get
static auto get(M &data)
Definition Templates.hpp:168
MoFEM::GetFTensor0FromVecImpl::get
static auto get(V &data)
Definition Templates.hpp:141
MoFEM::GetFTensor1FromMatImpl< Tensor_Dim, S, DataLayoutTraits< DataLayout::GaussByCoeffs >, M >::get
static auto get(M &data, int rr, int cc)
Definition Templates.hpp:237
MoFEM::GetFTensor2FromArrayImpl< 2, 2, S, M >::get
static auto get(M &data, const size_t rr, const size_t cc)
Definition Templates.hpp:677
MoFEM::GetFTensor2FromArrayImpl< 3, 3, S, M >::get
static auto get(M &data, const size_t rr, const size_t cc)
Definition Templates.hpp:688
MoFEM::GetFTensor2FromArrayImpl
Get FTensor2 from array.
Definition Templates.hpp:660
MoFEM::GetFTensor2FromArrayRawPtrImpl< 2, 2, M >::get
static auto get(M &data, const size_t rr, const size_t cc, const int ss=0)
Definition Templates.hpp:699
MoFEM::GetFTensor2FromArrayRawPtrImpl< 3, 3, M >::get
static auto get(M &data, const size_t rr, const size_t cc, const int ss=0)
Definition Templates.hpp:711
MoFEM::GetFTensor2FromMatImpl::get
static auto get(M &data, int rr=0, int cc=0)
Definition Templates.hpp:482
MoFEM::GetFTensor2SymmetricFromMatImpl< Tensor_Dim, S, DataLayoutTraits< DataLayout::GaussByCoeffs >, M >::get
static auto get(M &data, int rr=0, int cc=0)
Definition Templates.hpp:754
MoFEM::GetFTensor3DgFromMatImpl< Tensor_Dim01, Tensor_Dim2, S, DataLayoutTraits< DataLayout::GaussByCoeffs >, M >::get
static auto get(M &data, int rr=0, int cc=0)
Definition Templates.hpp:1088
MoFEM::GetFTensor3FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DataLayoutTraits< DataLayout::CoeffsByGauss >, M >::get
static auto get(M &data, int rr=0, int cc=0)
Definition Templates.hpp:964
MoFEM::GetFTensor3FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DataLayoutTraits< DataLayout::GaussByCoeffs >, M >::get
static auto get(M &data, int rr=0, int cc=0)
Definition Templates.hpp:929
MoFEM::GetFTensor4DdgFromMatImpl< Tensor_Dim01, Tensor_Dim23, S, DataLayoutTraits< DataLayout::GaussByCoeffs >, M >::get
static auto get(M &data, int rr=0, int cc=0)
Definition Templates.hpp:1411
MoFEM::GetFTensor4FromMatImpl::get
static auto get(M &data, int rr=0, int cc=0)
Definition Templates.hpp:1269
MoFEM::HashBit::operator()
unsigned int operator()(const id_type &value) const
Definition Templates.hpp:124
MoFEM::InvertTensorImpl< T1, T2, T3, 2 >::invert
static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t)
Definition Templates.hpp:2297
MoFEM::InvertTensorImpl< T1, T2, T3, 3 >::invert
static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t)
Definition Templates.hpp:2290
MoFEM::KeyFromKey::key1
KeyExtractor1 key1
Definition Templates.hpp:107
MoFEM::KeyFromKey::operator()
result_type operator()(Arg &arg) const
Definition Templates.hpp:102
MoFEM::KeyFromKey::KeyFromKey
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition Templates.hpp:98
MoFEM::KeyFromKey::key2
KeyExtractor2 key2
Definition Templates.hpp:108
MoFEM::KeyFromKey::result_type
KeyExtractor1::result_type result_type
Definition Templates.hpp:96
MoFEM::LtBit::operator()
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition Templates.hpp:112
MoFEM::MatrixSizeHelper< FTensor::Tensor2_symmetric< FTensor::CursorPtr< T *, S >, Tensor_Dim >, DataLayoutTraits< DataLayout::GaussByCoeffs > >::get
static auto get(M &m, const size_t nb_integration_pts)
Definition Templates.hpp:1680
MoFEM::MatrixSizeHelper< FTensor::Tensor2_symmetric< FTensor::CursorPtr< T *, S >, Tensor_Dim >, DataLayoutTraits< DataLayout::GaussByCoeffs > >::size
static auto size(M &m, const size_t nb_integration_pts)
Definition Templates.hpp:1668
MoFEM::Modify_change_nothing
Do nothing, used to rebuild database.
Definition Templates.hpp:2376
MoFEM::Modify_change_nothing::operator()
void operator()(T &e)
Definition Templates.hpp:2378
MoFEM::Modify_change_nothing::Modify_change_nothing
Modify_change_nothing()=default
MoFEM::RefEntExtractor
Extract entity handle form multi-index container.
Definition Templates.hpp:2329
MoFEM::RefEntExtractor::extract
static EntityHandle extract(const Iterator &it)
Definition Templates.hpp:2331
MoFEM::TempMeshset::meshset
EntityHandle meshset
Definition Templates.hpp:2416
MoFEM::TempMeshset::moab
moab::Interface & moab
Definition Templates.hpp:2417
MoFEM::TempMeshset::TempMeshset
TempMeshset(moab::Interface &moab)
Definition Templates.hpp:2403
MoFEM::parent_of
Definition Templates.hpp:2569
Generated by Doxygen 1.9.8 and hosted at University of Glasgow