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,
1360 M>::get(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,
1605 [[maybe_unused]] const size_t nb_integration_pts) {
1606#ifndef NDEBUG
1607 if (S == -1 || S == Tensor_Dim) {
1608 if (m.size1() != nb_integration_pts) {
1609 CHK_THROW_MESSAGE(
1610 MOFEM_DATA_INCONSISTENCY,
1611 "Matrix size is not compatible with the number of "
1612 "integration points: expected " +
1613 boost::lexical_cast<std::string>(nb_integration_pts) +
1614 " but got " + boost::lexical_cast<std::string>(m.size1()));
1615 }
1616 }
1617#endif
1618 return [&m]() {
1619 return GetFTensor1FromMatImpl<Tensor_Dim, S,
1620 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1621 M>::get(m, 0, 0);
1622 };
1623 }
1624};
1625
1626template <class T, int S, int Tensor_Dim0, int Tensor_Dim1>
1627struct MatrixSizeHelper<
1628 FTensor::Tensor2<FTensor::CursorPtr<T *, S>, Tensor_Dim0, Tensor_Dim1>,
1629 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1630 template <class M>
1631 static inline auto size(M &m, const size_t nb_integration_pts) {
1632 if (m.size1() != nb_integration_pts ||
1633 m.size2() != Tensor_Dim0 * Tensor_Dim1) {
1634 m.resize(nb_integration_pts, Tensor_Dim0 * Tensor_Dim1, false);
1635 }
1636 return [&m]() {
1637 return GetFTensor2FromMatImpl<Tensor_Dim0, Tensor_Dim1, S,
1638 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1639 M>::get(m, 0, 0);
1640 };
1641 }
1642 template <class M>
1643 static inline auto get(M &m,
1644 [[maybe_unused]] const size_t nb_integration_pts) {
1645#ifndef NDEBUG
1646 if (S == -1 || S == Tensor_Dim0 * Tensor_Dim1) {
1647 if (m.size1() != nb_integration_pts) {
1648 CHK_THROW_MESSAGE(
1649 MOFEM_DATA_INCONSISTENCY,
1650 "Matrix size is not compatible with the number of "
1651 "integration points: expected " +
1652 boost::lexical_cast<std::string>(nb_integration_pts) +
1653 " but got " + boost::lexical_cast<std::string>(m.size1()));
1654 }
1655 }
1656#endif
1657 return [&m]() {
1658 return GetFTensor2FromMatImpl<Tensor_Dim0, Tensor_Dim1, S,
1659 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1660 M>::get(m, 0, 0);
1661 };
1662 }
1663};
1664
1665template <class T, int S, int Tensor_Dim>
1666struct MatrixSizeHelper<
1667 FTensor::Tensor2_symmetric<FTensor::CursorPtr<T *, S>, Tensor_Dim>,
1668 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1669 template <class M>
1670 static inline auto size(M &m, const size_t nb_integration_pts) {
1671 constexpr int symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1672 if (m.size1() != nb_integration_pts || m.size2() != symm_size) {
1673 m.resize(nb_integration_pts, symm_size, false);
1674 }
1675 return [&m]() {
1676 return GetFTensor2SymmetricFromMatImpl<
1677 Tensor_Dim, S, DataLayoutTraits<DataLayout::GaussByCoeffs>,
1678 M>::get(m, 0, 0);
1679 };
1680 }
1681 template <class M>
1682 static inline auto get(M &m,
1683 [[maybe_unused]] const size_t nb_integration_pts) {
1684#ifndef NDEBUG
1685 constexpr int symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1686 if (S == -1 || S == symm_size) {
1687 if (m.size1() != nb_integration_pts) {
1688 CHK_THROW_MESSAGE(
1689 MOFEM_DATA_INCONSISTENCY,
1690 "Matrix size is not compatible with the number of "
1691 "integration points: expected " +
1692 boost::lexical_cast<std::string>(nb_integration_pts) +
1693 " but got " + boost::lexical_cast<std::string>(m.size1()));
1694 }
1695 }
1696#endif
1697 return [&m]() {
1698 return GetFTensor2SymmetricFromMatImpl<
1699 Tensor_Dim, S, DataLayoutTraits<DataLayout::GaussByCoeffs>,
1700 M>::get(m, 0, 0);
1701 };
1702 }
1703};
1704
1705template <class T, int S, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2>
1706struct MatrixSizeHelper<FTensor::Tensor3<FTensor::CursorPtr<T *, S>,
1707 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2>,
1708 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1709 template <class M>
1710 static inline auto size(M &m, const size_t nb_integration_pts) {
1711 constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
1712 if (m.size1() != nb_integration_pts || m.size2() != tensor_size) {
1713 m.resize(nb_integration_pts, tensor_size, false);
1714 }
1715 return [&m]() {
1716 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
1717 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1718 M>::get(m, 0, 0);
1719 };
1720 }
1721 template <class M>
1722 static inline auto get(M &m,
1723 [[maybe_unused]] const size_t nb_integration_pts) {
1724#ifndef NDEBUG
1725 constexpr int tensor_size = Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2;
1726 if (S == -1 || S == tensor_size) {
1727 if (m.size1() != nb_integration_pts) {
1728 CHK_THROW_MESSAGE(
1729 MOFEM_DATA_INCONSISTENCY,
1730 "Matrix size is not compatible with the number of "
1731 "integration points: expected " +
1732 boost::lexical_cast<std::string>(nb_integration_pts) +
1733 " but got " + boost::lexical_cast<std::string>(m.size1()));
1734 }
1735 }
1736#endif
1737 return [&m]() {
1738 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
1739 DataLayoutTraits<DataLayout::GaussByCoeffs>,
1740 M>::get(m, 0, 0);
1741 };
1742 }
1743};
1744
1745template <class T, int S, int Tensor_Dim01, int Tensor_Dim2>
1746struct MatrixSizeHelper<
1747 FTensor::Dg<FTensor::CursorPtr<T *, S>, Tensor_Dim01, Tensor_Dim2>,
1748 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1749 template <class M>
1750 static inline auto size(M &m, const size_t nb_integration_pts) {
1751 constexpr int dg_size =
1752 ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) * Tensor_Dim2;
1753 if (m.size1() != nb_integration_pts || m.size2() != dg_size) {
1754 m.resize(nb_integration_pts, dg_size, false);
1755 }
1756 return [&m]() {
1757 return GetFTensor3DgFromMatImpl<
1758 Tensor_Dim01, Tensor_Dim2, S,
1759 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1760 };
1761 }
1762 template <class M>
1763 static inline auto get(M &m,
1764 [[maybe_unused]] const size_t nb_integration_pts) {
1765#ifndef NDEBUG
1766 constexpr int dg_size =
1767 ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) * Tensor_Dim2;
1768 if (S == -1 || S == dg_size) {
1769 if (m.size1() != nb_integration_pts) {
1770 CHK_THROW_MESSAGE(
1771 MOFEM_DATA_INCONSISTENCY,
1772 "Matrix size is not compatible with the number of "
1773 "integration points: expected " +
1774 boost::lexical_cast<std::string>(nb_integration_pts) +
1775 " but got " + boost::lexical_cast<std::string>(m.size1()));
1776 }
1777 }
1778#endif
1779 return [&m]() {
1780 return GetFTensor3DgFromMatImpl<
1781 Tensor_Dim01, Tensor_Dim2, S,
1782 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1783 };
1784 }
1785};
1786
1787template <class T, int S, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
1788 int Tensor_Dim3>
1789struct MatrixSizeHelper<
1790 FTensor::Tensor4<FTensor::CursorPtr<T *, S>, Tensor_Dim0, Tensor_Dim1,
1791 Tensor_Dim2, Tensor_Dim3>,
1792 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1793 template <class M>
1794 static inline auto size(M &m, const size_t nb_integration_pts) {
1795 constexpr int tensor_size =
1796 Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3;
1797 if (m.size1() != nb_integration_pts || m.size2() != tensor_size) {
1798 m.resize(nb_integration_pts, tensor_size, false);
1799 }
1800 return [&m]() {
1801 return GetFTensor4FromMatImpl<
1802 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S,
1803 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1804 };
1805 }
1806 template <class M>
1807 static inline auto get(M &m,
1808 [[maybe_unused]] const size_t nb_integration_pts) {
1809#ifndef NDEBUG
1810 if (S == -1 || S == Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim2 * Tensor_Dim3) {
1811 if (m.size1() != nb_integration_pts) {
1812 CHK_THROW_MESSAGE(
1813 MOFEM_DATA_INCONSISTENCY,
1814 "Matrix size is not compatible with the number of "
1815 "integration points: expected " +
1816 boost::lexical_cast<std::string>(nb_integration_pts) +
1817 " but got " + boost::lexical_cast<std::string>(m.size1()));
1818 }
1819 }
1820#endif
1821 return [&m]() {
1822 return GetFTensor4FromMatImpl<
1823 Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S,
1824 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1825 };
1826 }
1827};
1828
1829template <class T, int S, int Tensor_Dim01, int Tensor_Dim23>
1830struct MatrixSizeHelper<
1831 FTensor::Ddg<FTensor::CursorPtr<T *, S>, Tensor_Dim01, Tensor_Dim23>,
1832 DataLayoutTraits<DataLayout::GaussByCoeffs>> {
1833 template <class M>
1834 static inline auto size(M &m, const size_t nb_integration_pts) {
1835 constexpr int ddg_size = ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) *
1836 ((Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2);
1837 if (m.size1() != nb_integration_pts || m.size2() != ddg_size) {
1838 m.resize(nb_integration_pts, ddg_size, false);
1839 }
1840 return [&m]() {
1841 return GetFTensor4DdgFromMatImpl<
1842 Tensor_Dim01, Tensor_Dim23, S,
1843 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1844 };
1845 }
1846 template <class M>
1847 static inline auto get(M &m,
1848 [[maybe_unused]] const size_t nb_integration_pts) {
1849#ifndef NDEBUG
1850 constexpr int ddg_size = ((Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2) *
1851 ((Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2);
1852 if (S == -1 || S == ddg_size) {
1853
1854 if (m.size1() != nb_integration_pts) {
1855 CHK_THROW_MESSAGE(
1856 MOFEM_DATA_INCONSISTENCY,
1857 "Matrix size is not compatible with the number of "
1858 "integration points: expected " +
1859 boost::lexical_cast<std::string>(nb_integration_pts) +
1860 " but got " + boost::lexical_cast<std::string>(m.size1()));
1861 }
1862 }
1863#endif
1864 return [&m]() {
1865 return GetFTensor4DdgFromMatImpl<
1866 Tensor_Dim01, Tensor_Dim23, S,
1867 DataLayoutTraits<DataLayout::GaussByCoeffs>, M>::get(m, 0, 0);
1868 };
1869 }
1870};
1871
1872/** @} */
1873
1874/** @name Voigt notation helpers
1875 * Conversions between tensor views and Voigt-like storage.
1876 * @{
1877 */
1878// Template functions for conversion to Voigt notation for 2nd order tensors
1879
1880template <int DIM, typename T> inline auto getVoigtVec(T &t_mat) {
1881 if constexpr (DIM == 3) {
1882 return std::array<double, 9>{t_mat(0, 0), t_mat(1, 1), t_mat(2, 2),
1883 t_mat(0, 1), t_mat(1, 0), t_mat(0, 2),
1884 t_mat(2, 0), t_mat(1, 2), t_mat(2, 1)};
1885 } else {
1886 return std::array<double, 5>{t_mat(0, 0), t_mat(1, 1), 1.0, t_mat(0, 1),
1887 t_mat(1, 0)};
1888 }
1889}
1890
1891template <int DIM, typename T> inline auto getVoigtVecSymm(T &t_mat) {
1892 constexpr double sqr2 = boost::math::constants::root_two<double>();
1893
1894 if constexpr (DIM == 3) {
1895 return std::array<double, 9>{t_mat(0, 0),
1896 t_mat(1, 1),
1897 t_mat(2, 2),
1898 sqr2 * t_mat(0, 1),
1899 sqr2 * t_mat(0, 2),
1900 sqr2 * t_mat(1, 2),
1901 0.0,
1902 0.0,
1903 0.0};
1904 } else {
1905 return std::array<double, 4>{t_mat(0, 0), t_mat(1, 1), 0.0,
1906 sqr2 * t_mat(0, 1)};
1907 }
1908}
1909
1910template <typename T>
1911inline auto getVoigtVecAxisymm(T &t_mat, const double hoop_term) {
1912 array<double, 5> vec{t_mat(0, 0), t_mat(1, 1), 1. + hoop_term, t_mat(0, 1),
1913 t_mat(1, 0)};
1914
1915 return vec;
1916}
1917
1918template <typename T>
1919inline auto getVoigtVecSymmAxisymm(T &t_mat, const double hoop_term) {
1920 constexpr double sqr2 = boost::math::constants::root_two<double>();
1921
1922 array<double, 4> vec_sym{t_mat(0, 0), t_mat(1, 1), hoop_term,
1923 sqr2 * t_mat(0, 1)};
1924 return vec_sym;
1925}
1926/** @} */
1927
1928/** @name Linear algebra helpers
1929 * LAPACK-backed matrix and eigenvalue utilities.
1930 * @{
1931 */
1932// list of lapack wrappers
1933/**
1934 * @brief compute matrix inverse with lapack dgetri
1935 *
1936 * @param mat input square matrix / output inverse matrix
1937 * @return MoFEMErrorCode
1938 */
1939inline MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat) {
1940 MoFEMFunctionBegin;
1941
1942 const size_t M = mat.size1();
1943 const size_t N = mat.size2();
1944
1945 if (M != N)
1946 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1947 "The input matrix for inverse computation is not square %zu != %zu",
1948 M, N);
1949
1950 int *ipv = new int[N];
1951 int lwork = N * N;
1952 double *work = new double[lwork];
1953 int info;
1954 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1955 if (info != 0)
1956 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "lapack error info = %d",
1957 info);
1958 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1959 if (info != 0)
1960 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1961 "lapack error info = %d", info);
1962
1963 delete[] ipv;
1964 delete[] work;
1965
1966 MoFEMFunctionReturn(0);
1967}
1968/**
1969 * @brief solve linear system with lapack dgesv
1970 *
1971 * @param mat input lhs square matrix / output L and U from the factorization
1972 * @param f input rhs vector / output solution vector
1973 * @return MoFEMErrorCode
1974 */
1975inline MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f) {
1976 MoFEMFunctionBegin;
1977
1978 const size_t M = mat.size1();
1979 const size_t N = mat.size2();
1980
1981 if (M == 0 || M != N)
1982 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1983 "The input matrix for inverse computation is not square %zu != %zu",
1984 M, N);
1985 if (f.size() != M)
1986 f.resize(M, false);
1987
1988 const int nrhs = 1;
1989 int info;
1990 int *ipiv = new int[M];
1991 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1992 &*f.data().begin(), M);
1993 if (info != 0) {
1994 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1995 "error lapack solve dgesv info = %d", info);
1996 }
1997
1998 delete[] ipiv;
1999 MoFEMFunctionReturn(0);
2000}
2001
2002/**
2003 * @brief Solve linear system of equations using Lapack
2004 *
2005 * @param mat
2006 * @param f
2007 * @return MoFEMErrorCode
2008 */
2009inline MoFEMErrorCode solveLinearSystem(const MatrixDouble &mat,
2010 VectorDouble &f) {
2011 MoFEMFunctionBegin;
2012 // copy matrix since on output lapack returns factorisation
2013 auto mat_copy = mat;
2014 CHKERR solveLinearSystem(mat_copy, f);
2015 MoFEMFunctionReturn(0);
2016}
2017
2018/** check if two doubles are approximately equal */
2019template <typename T = double> inline bool isEq(const T a, const T b) {
2020 const double eps = 100 * std::numeric_limits<T>::epsilon();
2021 const auto abs_max =
2022 std::max(static_cast<T>(1), std::max(std::abs(a), std::abs(b)));
2023 return std::abs(a - b) <= eps * abs_max;
2024}
2025
2026template <int DIM, typename T = double>
2027inline auto convertToShallow(FTensor::Tensor1<T, DIM> &vec) {
2028 return getVectorAdaptor(&vec(0), DIM);
2029}
2030
2031template <int DIM1, int DIM2, typename T = double>
2032inline auto convertToShallow(FTensor::Tensor2<T, DIM1, DIM2> &mat) {
2033 return getMatrixAdaptor(&mat(0, 0), DIM1, DIM2);
2034}
2035
2036/**
2037 * @brief Get the Uniq Nb objects
2038 *
2039 * @tparam DIM
2040 * @param ptr
2041 * @return auto
2042 */
2043template <int DIM, typename T = double>
2044inline auto getUniqNb(VectorShallowArrayAdaptor<T> vec, FTensor::Number<DIM>) {
2045 std::array<double, DIM> tmp;
2046 std::copy(vec.data().begin(), vec.data().begin() + DIM, tmp.begin());
2047 std::sort(tmp.begin(), tmp.end());
2048 std::size_t n_unique = 1;
2049 for (int i = 1; i < DIM; ++i) {
2050 if (!isEq(tmp[i], tmp[i - 1])) {
2051 ++n_unique;
2052 }
2053 }
2054 return n_unique;
2055}
2056
2057template <int DIM, typename T = double>
2058inline auto sortEigenVals(VectorShallowArrayAdaptor<T> eig,
2059 MatrixShallowArrayAdaptor<T> eigen_vec,
2060 FTensor::Number<DIM>) {
2061 if constexpr (DIM == 2) {
2062 MoFEMFunctionReturnHot(0);
2063 }
2064 static_assert(DIM == 2 || DIM == 3, "That works only for DIM = 3");
2065 MoFEMFunctionBegin;
2066 int i = 0, j = 1, k = 2;
2067 if (isEq<T>(eig(0), eig(1))) {
2068 i = 0;
2069 j = 2;
2070 k = 1;
2071 } else if (isEq<T>(eig(0), eig(2))) {
2072 i = 0;
2073 j = 1;
2074 k = 2;
2075 } else if (isEq<T>(eig(1), eig(2))) {
2076 i = 1;
2077 j = 0;
2078 k = 2;
2079 }
2080 const std::array<T, 9> eigen_vec_c{
2081 eigen_vec(i, 0), eigen_vec(i, 1), eigen_vec(i, 2),
2082 eigen_vec(j, 0), eigen_vec(j, 1), eigen_vec(j, 2),
2083 eigen_vec(k, 0), eigen_vec(k, 1), eigen_vec(k, 2)};
2084 const std::array<T, 3> eig_c{eig(i), eig(j), eig(k)};
2085 std::copy(eigen_vec_c.begin(), eigen_vec_c.end(), eigen_vec.data().begin());
2086 std::copy(eig_c.begin(), eig_c.end(), eig.data().begin());
2087 MoFEMFunctionReturn(0);
2088}
2089
2090template <int DIM, typename T = double>
2091inline auto getUniqNb(FTensor::Tensor1<T, DIM> vec) {
2092 return getUniqNb<DIM, T>(convertToShallow(vec), FTensor::Number<DIM>{});
2093}
2094
2095template <int DIM, typename T = double>
2096inline auto sortEigenVals(FTensor::Tensor1<T, DIM> &eig,
2097 FTensor::Tensor2<T, DIM, DIM> &eigen_vec) {
2098 return sortEigenVals<DIM, T>(convertToShallow(eig),
2099 convertToShallow(eigen_vec),
2100 FTensor::Number<DIM>{});
2101}
2102
2103/**
2104 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
2105 *
2106 * LAPACK `dsyev` returns eigenvectors as columns of the column-major matrix
2107 * `Q`. Since MoFEM passes the raw buffer of `eigen_vec` directly to LAPACK
2108 * and then accesses the same buffer with C++ indexing, `eigen_vec(i, j)`
2109 * stores `Q(j, i)`. Therefore, in MoFEM convention, row `i` of `eigen_vec`
2110 * stores eigenvector `i`, and `Q = eigen_vec^T` in the usual linear algebra
2111 * notation with eigenvectors in columns.
2112 *
2113 * @param mat input symmetric matrix
2114 * @param eig output eigen values sorted
2115 * @param eigen_vec output matrix of eigen vectors in MoFEM convention
2116 * @return MoFEMErrorCode
2117 */
2118inline MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat,
2119 VectorDouble &eig,
2120 MatrixDouble &eigen_vec) {
2121 MoFEMFunctionBegin;
2122
2123 const size_t M = mat.size1();
2124 const size_t N = mat.size2();
2125
2126 if (M == 0 || M != N)
2127 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2128 "The input matrix for eigen value computation is not square %zu "
2129 "!= %zu",
2130 M, N);
2131 if (eig.size() != M)
2132 eig.resize(M, false);
2133
2134 eigen_vec = mat;
2135 const int n = M;
2136 const int lda = M;
2137 const int size = (M + 2) * M;
2138 int lwork = size;
2139 double *work = new double[size];
2140
2141 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
2142 &*eig.data().begin(), work, lwork) > 0)
2143 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
2144 "The algorithm failed to compute eigenvalues.");
2145
2146 delete[] work;
2147 MoFEMFunctionReturn(0);
2148}
2149
2150/**
2151 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
2152 *
2153 * LAPACK `dsyev` returns eigenvectors as columns of the column-major matrix
2154 * `Q`. Since MoFEM passes the raw buffer of `eigen_vec` directly to LAPACK
2155 * and then accesses the same buffer with C++ indexing, `eigen_vec(i, j)`
2156 * stores `Q(j, i)`. Therefore, in MoFEM convention, row `i` of `eigen_vec`
2157 * stores eigenvector `i`, and `Q = eigen_vec^T` in the usual linear algebra
2158 * notation with eigenvectors in columns.
2159 *
2160 * @tparam DIM
2161 * @param eigen_vec input / output DIM x DIM matrix of eigen vectors in MoFEM
2162 * convention
2163 * @param eig output eigen values sorted
2164 * @return MoFEMErrorCode
2165 */
2166template <int DIM, typename T1, typename T2>
2167inline MoFEMErrorCode
2168computeEigenValuesSymmetric(FTensor::Tensor2<T1, DIM, DIM> &eigen_vec,
2169 FTensor::Tensor1<T2, DIM> &eig) {
2170 MoFEMFunctionBegin;
2171
2172 const int n = DIM;
2173 const int lda = DIM;
2174 const int lwork = (DIM + 2) * DIM;
2175 std::array<double, (DIM + 2) * DIM> work;
2176
2177 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
2178 lwork) > 0)
2179 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
2180 "The algorithm failed to compute eigenvalues.");
2181 MoFEMFunctionReturn(0);
2182}
2183
2184/**
2185 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
2186 *
2187 * LAPACK `dsyev` returns eigenvectors as columns of the column-major matrix
2188 * `Q`. Since MoFEM passes the raw buffer of `eigen_vec` directly to LAPACK
2189 * and then accesses the same buffer with C++ indexing, `eigen_vec(i, j)`
2190 * stores `Q(j, i)`. Therefore, in MoFEM convention, row `i` of `eigen_vec`
2191 * stores eigenvector `i`, and `Q = eigen_vec^T` in the usual linear algebra
2192 * notation with eigenvectors in columns.
2193 *
2194 * @tparam DIM
2195 * @param mat input tensor pointer of size DIM x DIM
2196 * @param eig output eigen values sorted
2197 * @param eigen_vec output matrix of eigen vectors in MoFEM convention
2198 * @return MoFEMErrorCode
2199 */
2200template <int DIM, typename T1, typename T2, typename T3>
2201inline MoFEMErrorCode
2202computeEigenValuesSymmetric(const FTensor::Tensor2_symmetric<T1, DIM> &mat,
2203 FTensor::Tensor1<T2, DIM> &eig,
2204 FTensor::Tensor2<T3, DIM, DIM> &eigen_vec) {
2205 MoFEMFunctionBegin;
2206 for (int ii = 0; ii != DIM; ii++)
2207 for (int jj = 0; jj != DIM; jj++)
2208 eigen_vec(ii, jj) = mat(ii, jj);
2209
2210 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
2211
2212 MoFEMFunctionReturn(0);
2213}
2214
2215/** @} */
2216
2217/** @name Tensor algebra helpers
2218 * Determinant and inverse routines for tensor-valued objects.
2219 * @{
2220 */
2221/**
2222 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
2223 *
2224 * @tparam T
2225 * @param t
2226 * @return double
2227 */
2228template <typename T> static inline auto determinantTensor3by3(T &t) {
2229 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
2230 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
2231 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
2232}
2233
2234/**
2235 * @brief Calculate the determinant of a 2x2 matrix or a tensor of rank 2
2236 *
2237 * @tparam T
2238 * @param t
2239 * @return double
2240 */
2241template <typename T> static inline auto determinantTensor2by2(T &t) {
2242 return t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
2243}
2244
2245template <typename T, int DIM> struct DeterminantTensorImpl;
2246
2247template <typename T> struct DeterminantTensorImpl<T, 3> {
2248 static inline auto get(T &t) { return determinantTensor3by3(t); }
2249};
2250
2251template <typename T> struct DeterminantTensorImpl<T, 2> {
2252 static auto get(T &t) { return determinantTensor2by2(t); }
2253};
2254
2255/**
2256 * @brief Calculate the determinant of a tensor of rank DIM
2257 */
2258template <typename T, int DIM>
2262
2263/**
2264 * @brief Calculate the determinant of a tensor of rank DIM
2265 */
2266template <typename T, int DIM>
2270
2271/**
2272 * \brief Calculate inverse of tensor rank 2 at integration points
2273
2274 */
2275template <int Tensor_Dim, class T, class L, class A>
2276inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
2277 ublas::vector<T, A> &det_data,
2278 ublas::matrix<T, L, A> &inv_jac_data) {
2279 MoFEMFunctionBeginHot;
2280 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2281 "Specialization for this template not yet implemented");
2282 MoFEMFunctionReturnHot(0);
2283}
2284
2285template <>
2286inline MoFEMErrorCode
2287invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
2288 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
2289
2290/**
2291 * \brief Calculate determinant 3 by 3
2292
2293 */
2294template <class T1, class T2>
2300
2301/**
2302 * \brief Calculate determinant 2 by 2
2303
2304 */
2305template <class T1, class T2>
2311
2312/**
2313 * \brief Calculate matrix inverse 3 by 3
2314
2315 */
2316template <class T1, class T2, class T3>
2317inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
2318 MoFEMFunctionBeginHot;
2319 const auto a = t(0, 0), b = t(0, 1), c = t(0, 2);
2320 const auto d = t(1, 0), e = t(1, 1), f = t(1, 2);
2321 const auto g = t(2, 0), h = t(2, 1), i = t(2, 2);
2322 const auto inv_det = 1.0 / det;
2323 inv_t(0, 0) = (e * i - f * h) * inv_det;
2324 inv_t(0, 1) = (c * h - b * i) * inv_det;
2325 inv_t(0, 2) = (b * f - c * e) * inv_det;
2326 inv_t(1, 0) = (f * g - d * i) * inv_det;
2327 inv_t(1, 1) = (a * i - c * g) * inv_det;
2328 inv_t(1, 2) = (c * d - a * f) * inv_det;
2329 inv_t(2, 0) = (d * h - e * g) * inv_det;
2330 inv_t(2, 1) = (b * g - a * h) * inv_det;
2331 inv_t(2, 2) = (a * e - b * d) * inv_det;
2332 MoFEMFunctionReturnHot(0);
2333}
2334
2335/**
2336 * \brief Calculate matrix inverse 2 by 2
2337
2338 */
2339template <class T1, class T2, class T3>
2340inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
2341 MoFEMFunctionBeginHot;
2342 const auto a = t(0, 0);
2343 const auto b = t(0, 1);
2344 const auto c = t(1, 0);
2345 const auto d = t(1, 1);
2346 inv_t(0, 0) = d / det;
2347 inv_t(0, 1) = -b / det;
2348 inv_t(1, 0) = -c / det;
2349 inv_t(1, 1) = a / det;
2350 MoFEMFunctionReturnHot(0);
2351}
2352
2353/**
2354 * \brief Calculate matrix inverse, overload for symmetric tensor
2355
2356 */
2357template <typename T1, typename T2, typename T3>
2358inline MoFEMErrorCode
2359invertTensor3by3(FTensor::Tensor2_symmetric<T1, 3> &t, T2 &det,
2360 FTensor::Tensor2_symmetric<T3, 3> &inv_t) {
2361 MoFEMFunctionBeginHot;
2362 const auto a = t(0, 0);
2363 const auto b = t(0, 1);
2364 const auto c = t(0, 2);
2365 const auto d = t(1, 1);
2366 const auto e = t(1, 2);
2367 const auto f = t(2, 2);
2368 const auto inv_det = 1.0 / det;
2369 inv_t(0, 0) = (d * f - e * e) * inv_det;
2370 inv_t(0, 1) = (c * e - b * f) * inv_det;
2371 inv_t(0, 2) = (b * e - c * d) * inv_det;
2372 inv_t(1, 1) = (a * f - c * c) * inv_det;
2373 inv_t(1, 2) = (b * c - a * e) * inv_det;
2374 inv_t(2, 2) = (a * d - b * b) * inv_det;
2375 MoFEMFunctionReturnHot(0);
2376}
2377
2378/**
2379 * \brief Calculate matrix inverse, overload for symmetric tensor
2380
2381 */
2382template <typename T1, typename T2, typename T3>
2383inline MoFEMErrorCode
2384invertTensor2by2(FTensor::Tensor2_symmetric<T1, 2> &t, T2 &det,
2385 FTensor::Tensor2_symmetric<T3, 2> &inv_t) {
2386 MoFEMFunctionBeginHot;
2387 const auto a = t(0, 0);
2388 const auto b = t(0, 1);
2389 const auto d = t(1, 1);
2390 const auto inv_det = 1.0 / det;
2391 inv_t(0, 0) = d * inv_det;
2392 inv_t(0, 1) = -b * inv_det;
2393 inv_t(1, 1) = a * inv_det;
2394 MoFEMFunctionReturnHot(0);
2395}
2396
2397template <typename T1, typename T2, typename T3, int DIM>
2398struct InvertTensorImpl;
2399
2400template <typename T1, typename T2, typename T3>
2401struct InvertTensorImpl<T1, T2, T3, 3> {
2402 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
2403 return invertTensor3by3(t, det, inv_t);
2404 }
2405};
2406
2407template <typename T1, typename T2, typename T3>
2408struct InvertTensorImpl<T1, T2, T3, 2> {
2409 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
2410 return invertTensor2by2(t, det, inv_t);
2411 }
2412};
2413
2414template <typename T1, typename T2, typename T3, int DIM>
2415static inline MoFEMErrorCode
2422
2423template <typename T1, typename T2, typename T3, int DIM>
2424static inline MoFEMErrorCode
2431/** @} */
2432
2433/** @name Mesh and range helpers
2434 * Utility types and helpers for mesh entities, ranges, and handles.
2435 * @{
2436 */
2437/**
2438 * @brief Extract entity handle form multi-index container
2439 *
2440 */
2441struct RefEntExtractor {
2442 template <typename Iterator>
2443 static inline EntityHandle extract(const Iterator &it) {
2444 return (*it)->getEnt();
2445 }
2446};
2447
2448/**
2449 * @brief Insert ordered mofem multi-index into range
2450 *
2451 * \note Inserted range has to be ordered.
2452 *
2453 * \code
2454 * auto hi_rit = refEntsPtr->upper_bound(start);
2455 * auto hi_rit = refEntsPtr->upper_bound(end);
2456 * Range to_erase;
2457 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
2458 * \endcode
2459 *
2460 * @tparam Iterator
2461 * @param r
2462 * @param begin_iter
2463 * @param end_iter
2464 * @return moab::Range::iterator
2465 */
2466template <typename Extractor, typename Iterator>
2467moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
2468 Iterator end_iter) {
2469 moab::Range::iterator hint = r.begin();
2470 while (begin_iter != end_iter) {
2471 size_t j = 0;
2472 auto bi = Extractor::extract(begin_iter);
2473 Iterator pj = begin_iter;
2474 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
2475 ++pj;
2476 ++j;
2477 }
2478 hint = r.insert(hint, bi, bi + (j - 1));
2479 begin_iter = pj;
2480 }
2481 return hint;
2482};
2483
2484/**
2485 * @brief Do nothing, used to rebuild database
2486 *
2487 */
2488struct Modify_change_nothing {
2489 Modify_change_nothing() = default;
2490 template <typename T> inline void operator()(T &e) {}
2491};
2492
2493/**
2494 * @brief Template used to reconstruct multi-index
2495 *
2496 * @tparam MI multi-index
2497 * @tparam Modifier
2498 * @param mi
2499 * @param mo
2500 * @return MoFEMErrorCode
2501 */
2502template <typename MI, typename MO = Modify_change_nothing>
2503inline MoFEMErrorCode reconstructMultiIndex(const MI &mi,
2504 MO &&mo = Modify_change_nothing()) {
2505 MoFEMFunctionBegin;
2506 for (auto it = mi.begin(); it != mi.end(); ++it) {
2507 if (!const_cast<MI &>(mi).modify(it, mo))
2508 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2509 "Houston we have a problem");
2510 }
2511 MoFEMFunctionReturn(0);
2512}
2513
2514struct TempMeshset {
2515 TempMeshset(moab::Interface &moab) : moab(moab) {
2516 rval = moab.create_meshset(MESHSET_SET, meshset);
2517 MOAB_THROW(rval);
2518 }
2519 virtual ~TempMeshset() { delete_meshset(); }
2520 operator EntityHandle() const { return meshset; }
2521 auto get_ptr() { return &meshset; }
2522
2523private:
2524 void delete_meshset() {
2525 rval = moab.delete_entities(&meshset, 1);
2526 MOAB_THROW(rval);
2527 }
2528 EntityHandle meshset;
2529 moab::Interface &moab;
2530};
2531
2532/**
2533 * @brief Create smart pointer to temporary meshset
2534 *
2535 */
2536inline auto get_temp_meshset_ptr(moab::Interface &moab) {
2537 return boost::make_shared<TempMeshset>(moab);
2538};
2539
2540inline auto id_from_handle(const EntityHandle h) {
2541 return static_cast<EntityID>(h & MB_ID_MASK);
2542};
2543
2544/**
2545 * @brief get type from entity handle
2546 *
2547 */
2548inline auto type_from_handle(const EntityHandle h) {
2549 return static_cast<EntityType>(h >> MB_ID_WIDTH);
2550};
2551
2552/**
2553 * @brief get entity handle from type and id
2554 *
2555 */
2556inline auto ent_form_type_and_id(const EntityType type, const EntityID id) {
2557 return (static_cast<EntityHandle>(type) << MB_ID_WIDTH) | id;
2558};
2559
2560/**
2561 * @brief get entity dimension form handle
2562 *
2563 */
2564inline auto dimension_from_handle(const EntityHandle h) {
2565 return moab::CN::Dimension(type_from_handle(h));
2566};
2567
2568/**
2569 * @brief get entity type name from handle
2570 *
2571 */
2572inline auto type_name_from_handle(const EntityHandle h) {
2573 return moab::CN::EntityTypeName(type_from_handle(h));
2574};
2575
2576/**
2577 * @brief get field bit id from bit number
2578 *
2579 */
2580inline auto field_bit_from_bit_number(const int bit_number) {
2581 return BitFieldId().set(bit_number - 1);
2582};
2583
2584/**
2585 * @brief Insert ranges
2586 *
2587 * @tparam I
2588 * @param f
2589 * @param s
2590 * @param tester
2591 * @param inserter
2592 * @return auto
2593 */
2594template <typename I>
2595auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
2596 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
2597 MoFEMFunctionBegin;
2598
2599 auto first = f;
2600 while (first != s)
2601 if (tester(first)) {
2602
2603 auto second = first;
2604 ++second;
2605
2606 while (second != s) {
2607 if (tester(second))
2608 ++second;
2609 else
2610 break;
2611 }
2612
2613 CHKERR inserter(first, second);
2614
2615 first = second;
2616 if (first != s)
2617 ++first;
2618
2619 } else {
2620 ++first;
2621 }
2622
2623 MoFEMFunctionReturn(0);
2624}
2625
2626/**
2627 * @brief Create Array
2628 *
2629 * See:
2630 * <a
2631 * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
2632 * stack overflow</a>
2633 *
2634 * @tparam Dest
2635 * @tparam Arg
2636 * @param arg
2637 * @return constexpr auto
2638 */
2639template <typename Dest = void, typename... Arg>
2640constexpr auto make_array(Arg &&...arg) {
2641 if constexpr (std::is_same<void, Dest>::value)
2642 return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
2643 {std::forward<Arg>(arg)...}};
2644 else
2645 return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
2646}
2647/** @} */
2648
2649/** @name FTensor index aliases
2650 * Convenient aliases for commonly used FTensor indices.
2651 * @{
2652 */
2653template <int DIM> using i_FTIndex = FTensor::Index<'i', DIM>;
2654template <int DIM> using j_FTIndex = FTensor::Index<'j', DIM>;
2655template <int DIM> using k_FTIndex = FTensor::Index<'k', DIM>;
2656template <int DIM> using l_FTIndex = FTensor::Index<'l', DIM>;
2657template <int DIM> using m_FTIndex = FTensor::Index<'m', DIM>;
2658template <int DIM> using n_FTIndex = FTensor::Index<'n', DIM>;
2659template <int DIM> using o_FTIndex = FTensor::Index<'o', DIM>;
2660template <int DIM> using p_FTIndex = FTensor::Index<'p', DIM>;
2661template <int DIM> using x_FTIndex = FTensor::Index<'x', DIM>;
2662template <int DIM> using y_FTIndex = FTensor::Index<'y', DIM>;
2663template <int DIM> using I_FTIndex = FTensor::Index<'I', DIM>;
2664template <int DIM> using J_FTIndex = FTensor::Index<'J', DIM>;
2665template <int DIM> using K_FTIndex = FTensor::Index<'K', DIM>;
2666template <int DIM> using L_FTIndex = FTensor::Index<'L', DIM>;
2667template <int DIM> using M_FTIndex = FTensor::Index<'M', DIM>;
2668template <int DIM> using N_FTIndex = FTensor::Index<'N', DIM>;
2669template <int DIM> using O_FTIndex = FTensor::Index<'O', DIM>;
2670template <int DIM> using Z_FTIndex = FTensor::Index<'Z', DIM>;
2671template <int DIM> using Y_FTIndex = FTensor::Index<'Y', DIM>;
2672
2673#define FTENSOR_INDEX(DIM, I) I##_FTIndex<DIM> I;
2674
2675#define FTENSOR_INDEX_DECL(r, DIM, I) FTENSOR_INDEX(DIM, I)
2676
2677#define FTENSOR_INDEXES(DIM, ...) \
2678 BOOST_PP_SEQ_FOR_EACH(FTENSOR_INDEX_DECL, DIM, \
2679 BOOST_PP_VARIADIC_TO_SEQ(__VA_ARGS__))
2680
2681template <typename T> struct parent_of;
2682/** @} */
2683
2684} // namespace MoFEM
2685
2686#endif //__TEMPLATES_HPP__
type
std::string type
Definition JsonConfigManager.cpp:28
a
constexpr double a
Definition approx_sphere.cpp:30
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
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:40
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:39
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
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
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:2548
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::sortEigenVals
auto sortEigenVals(VectorShallowArrayAdaptor< T > eig, MatrixShallowArrayAdaptor< T > eigen_vec, FTensor::Number< DIM >)
Definition Templates.hpp:2058
MoFEM::getVoigtVecAxisymm
auto getVoigtVecAxisymm(T &t_mat, const double hoop_term)
Definition Templates.hpp:1911
MoFEM::getFTensor2FromMat
auto getFTensor2FromMat(M &data)
Get tensor rank 2 (matrix) form data matrix.
Definition Templates.hpp:542
MoFEM::isEq
bool isEq(const T a, const T b)
Definition Templates.hpp:2019
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::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:1880
MoFEM::reconstructMultiIndex
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition Templates.hpp:2503
MoFEM::getUniqNb
auto getUniqNb(VectorShallowArrayAdaptor< T > vec, FTensor::Number< DIM >)
Get the Uniq Nb objects.
Definition Templates.hpp:2044
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:2467
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::invertTensor2by2
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition Templates.hpp:2340
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:1939
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:2540
MoFEM::getFTensor3DgFromPtr
static auto getFTensor3DgFromPtr(T *ptr)
Definition Templates.hpp:1242
MoFEM::getFTensor1FromArray
auto getFTensor1FromArray(T &data)
Definition Templates.hpp:433
MoFEM::convertToShallow
auto convertToShallow(FTensor::Tensor1< T, DIM > &vec)
Definition Templates.hpp:2027
MoFEM::dimension_from_handle
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
Definition Templates.hpp:2564
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:1975
MoFEM::determinantTensor2by2
static auto determinantTensor2by2(T &t)
Calculate the determinant of a 2x2 matrix or a tensor of rank 2.
Definition Templates.hpp:2241
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:2572
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:2276
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:2416
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:2640
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:2259
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:2580
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:1891
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:2536
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:2118
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:2228
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:1919
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:2595
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:2556
I
constexpr IntegrationType I
Definition operators_tests.cpp:31
h
double h
Definition photon_diffusion.cpp:71
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:2409
MoFEM::InvertTensorImpl< T1, T2, T3, 3 >::invert
static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t)
Definition Templates.hpp:2402
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:1682
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:1670
MoFEM::Modify_change_nothing
Do nothing, used to rebuild database.
Definition Templates.hpp:2488
MoFEM::Modify_change_nothing::operator()
void operator()(T &e)
Definition Templates.hpp:2490
MoFEM::Modify_change_nothing::Modify_change_nothing
Modify_change_nothing()=default
MoFEM::RefEntExtractor
Extract entity handle form multi-index container.
Definition Templates.hpp:2441
MoFEM::RefEntExtractor::extract
static EntityHandle extract(const Iterator &it)
Definition Templates.hpp:2443
MoFEM::TempMeshset::meshset
EntityHandle meshset
Definition Templates.hpp:2528
MoFEM::TempMeshset::moab
moab::Interface & moab
Definition Templates.hpp:2529
MoFEM::TempMeshset::TempMeshset
TempMeshset(moab::Interface &moab)
Definition Templates.hpp:2515
MoFEM::parent_of
Definition Templates.hpp:2681
Generated by Doxygen 1.9.8 and hosted at University of Glasgow