v0.15.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
10template <typename T> using ShardVec = boost::shared_ptr<std::vector<T>>;
11
12/**
13 * @brief Get Vector adaptor
14 *
15 * \code
16 *
17 * double *a;
18 * CHKERR VecGetArray(v,&a);
19 *
20 * for(int n = 0; n != nodes; ++n) {
21 *
22 * auto a = getVectorAdaptor(&a[3*n], 3);
23 * double dot = inner_prod(a, a);
24 *
25 * }
26 *
27 * CHKERR VecRetsoreArray(v,&a);
28 * \endcode
29 *
30 */
31template <typename T1> inline auto getVectorAdaptor(T1 ptr, const size_t n) {
32 typedef typename std::remove_pointer<T1>::type T;
34 ublas::shallow_array_adaptor<T>(n, ptr));
35};
36
37/**
38 * @brief Get Matrix adaptor
39 *
40 * \code
41 *
42 * double *a;
43 * CHKERR VecGetArray(v,&a);
44 *
45 * for(int n = 0; n != nodes; ++n) {
46 *
47 * auto F = getMatrixAdaptor(&a[3*3*n], 3, 3);
48 * MatrixDouble C = prod(F, trans(F));
49 *
50 * }
51 *
52 * CHKERR VecRetsoreArray(v,&a);
53 * \endcode
54 *
55 */
56template <typename T1>
57inline auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m) {
58 typedef typename std::remove_pointer<T1>::type T;
60 n, m, ublas::shallow_array_adaptor<T>(n * m, ptr));
61};
62
63/**
64 * This small utility that cascades two key extractors will be
65 * used throughout the boost example
66 * <a
67 * href=http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp>
68 * http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp
69 * </a>
70 */
71template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
72public:
73 typedef typename KeyExtractor1::result_type result_type;
74
75 KeyFromKey(const KeyExtractor1 &key1_ = KeyExtractor1(),
76 const KeyExtractor2 &key2_ = KeyExtractor2())
77 : key1(key1_), key2(key2_) {}
78
79 template <typename Arg> result_type operator()(Arg &arg) const {
80 return key1(key2(arg));
81 }
82
83private:
84 KeyExtractor1 key1;
85 KeyExtractor2 key2;
86};
87
88template <typename id_type> struct LtBit {
89 inline bool operator()(const id_type &valueA, const id_type &valueB) const {
90 return valueA.to_ulong() < valueB.to_ulong();
91 }
92};
93
94template <typename id_type> struct EqBit {
95 inline bool operator()(const id_type &valueA, const id_type &valueB) const {
96 return valueA.to_ulong() == valueB.to_ulong();
97 }
98};
99
100template <typename id_type> struct HashBit {
101 inline unsigned int operator()(const id_type &value) const {
102 return value.to_ulong();
103 }
104};
105
106template <class X> inline std::string toString(X x) {
107 std::ostringstream buffer;
108 buffer << x;
109 return buffer.str();
110}
111
112template <int S, class T, class A> struct GetFTensor0FromVecImpl {
113 static inline auto get(ublas::vector<T, A> &data) {
114 return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(&*data.data().begin());
115 }
116};
117
118/**
119* \brief Get tensor rank 0 (scalar) form data vector
120
121Example how to use it.
122\code
123VectorDouble vec;
124vec.resize(nb_gauss_pts,false);
125vec.clear();
126auto t0 = getFTensor0FromData(data);
127for(int gg = 0;gg!=nb_gauss_pts;gg++) {
128
129 ++t0;
130}
131\endcode
132
133*/
134template <int S = 1, class T, class A>
135static inline auto getFTensor0FromVec(ublas::vector<T, A> &data) {
137}
138
139template <int Tensor_Dim, int S, class T, class L, class A>
141
142template <int S, class T, class A>
143struct GetFTensor1FromMatImpl<3, S, T, ublas::row_major, A> {
144 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
145#ifndef NDDEBUG
146 if (data.size1() != 3)
148 "getFTensor1FromMat<3>: wrong size of data matrix, number of "
149 "rows should be 3 but is " +
150 boost::lexical_cast<std::string>(data.size1()));
151#endif
153 &data(0, 0), &data(1, 0), &data(2, 0));
154 }
155};
156
157template <int S, class T, class A>
158struct GetFTensor1FromMatImpl<4, S, T, ublas::row_major, A> {
159 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
160#ifndef NDDEBUG
161 if (data.size1() != 4)
163 "getFTensor1FromMat<4>: wrong size of data matrix, number of "
164 "rows should be 4 but is " +
165 boost::lexical_cast<std::string>(data.size1()));
166#endif
168 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
169 }
170};
171
172template <int S, class T, class A>
173struct GetFTensor1FromMatImpl<6, S, T, ublas::row_major, A> {
174 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
175#ifndef NDDEBUG
176 if (data.size1() != 6)
178 "getFTensor1FromMat<6>: wrong size of data matrix, number of "
179 "rows should be 6 but is " +
180 boost::lexical_cast<std::string>(data.size1()));
181#endif
183 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
184 &data(5, 0));
185 }
186};
187
188template <int S, class T, class A>
189struct GetFTensor1FromMatImpl<9, S, T, ublas::row_major, A> {
190 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
191#ifndef NDDEBUG
192 if (data.size1() != 9)
194 "getFTensor1FromMat<9>: wrong size of data matrix, number of "
195 "rows should be 9 but is " +
196 boost::lexical_cast<std::string>(data.size1()));
197#endif
199 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
200 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
201 }
202};
203
204template <int S, class T, class A>
205struct GetFTensor1FromMatImpl<2, S, T, ublas::row_major, A> {
206 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
207#ifndef NDDEBUG
208 if (data.size1() != 2)
210 "getFTensor1FromMat<2>: wrong size of data matrix, number of "
211 "rows should be 2 but is " +
212 boost::lexical_cast<std::string>(data.size1()));
213#endif
215 &data(1, 0));
216 }
217};
218
219template <int S, class T, class A>
220struct GetFTensor1FromMatImpl<1, S, T, ublas::row_major, A> {
221 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
222#ifndef NDEBUG
223 if (data.size1() != 1)
225 "getFTensor1FromMat<1>: wrong size of data matrix, number of "
226 "rows should be 1 but is " +
227 boost::lexical_cast<std::string>(data.size1()));
228#endif
230 }
231};
232
233/**
234 * \brief Get tensor rank 1 (vector) form data matrix
235 */
236template <int Tensor_Dim, int S = 1, class T, class L, class A>
238getFTensor1FromMat(ublas::matrix<T, L, A> &data) {
239 static_assert(!std::is_same<T, T>::value, "not implemented");
240}
241
242/**
243 * \brief Get tensor rank 1 (vector) form data matrix (specialization)
244 */
245template <int Tensor_Dim, int S = 1>
247 return GetFTensor1FromMatImpl<Tensor_Dim, S, double, ublas::row_major,
248 DoubleAllocator>::get(data);
249}
250
251template <int Tensor_Dim1, int Tensor_Dim2, int S, class T, class L, class A>
253 static inline auto get(ublas::matrix<T, L, A> &data) {
254#ifndef NDEBUG
255 if (data.size1() != Tensor_Dim1 * Tensor_Dim2) {
257 "getFTensor2FromMat<" +
258 boost::lexical_cast<std::string>(Tensor_Dim1) + "," +
259 boost::lexical_cast<std::string>(Tensor_Dim2) +
260 ">: wrong size of rows in data matrix, should be " +
261 boost::lexical_cast<std::string>(Tensor_Dim1 * Tensor_Dim2) +
262 " but is " + boost::lexical_cast<std::string>(data.size1()));
263 }
264#endif
265 std::array<double *, Tensor_Dim1 * Tensor_Dim2> ptrs;
266 for (auto i = 0; i != Tensor_Dim1 * Tensor_Dim2; ++i)
267 ptrs[i] = &data(i, 0);
269 Tensor_Dim2>(ptrs);
270 }
271};
272
273/**
274 * \brief Get tensor rank 2 (matrix) form data matrix
275 */
276template <int Tensor_Dim1, int Tensor_Dim2>
277inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim1, Tensor_Dim2>
279 return GetFTensor2FromMatImpl<Tensor_Dim1, Tensor_Dim2, 1, double,
280 ublas::row_major, DoubleAllocator>::get(data);
281}
282
283template <int Tensor_Dim1, int Tensor_Dim2>
284inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim1, Tensor_Dim2>
286
287/**
288 * Template specialization for getFTensor2FromMat
289 */
290template <>
295
296template <int Tensor_Dim, int S, class T, class L, class A>
298
299template <int S, class T, class L, class A>
301 static inline auto get(ublas::matrix<T, L, A> &data) {
302#ifndef NDEBUG
303 if (data.size1() != 6)
305 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
306 "of rows should be 6 but is " +
307 boost::lexical_cast<std::string>(data.size1()));
308#endif
310 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
311 &data(5, 0));
312 }
313};
314
315template <int S, class T, class L, class A>
317 static inline auto get(ublas::matrix<T, L, A> &data) {
318#ifndef NDEBUG
319 if (data.size1() != 3)
321 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, number "
322 "of rows should be 3 but is " +
323 boost::lexical_cast<std::string>(data.size1()));
324#endif
326 &data(0, 0), &data(1, 0), &data(2, 0));
327 }
328};
329
330/**
331 * \brief Get symmetric tensor rank 2 (matrix) form data matrix
332 */
333template <int Tensor_Dim, int S, class T, class L, class A>
334static inline auto getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
336}
337
338template <int Tensor_Dim, int S = 1>
339static inline auto getFTensor2SymmetricFromMat(MatrixDouble &data) {
340 return getFTensor2SymmetricFromMat<Tensor_Dim, S, double, ublas::row_major,
341 DoubleAllocator>(data);
342}
343
344template <int Tensor_Dim01, int Tensor_Dim23, int S, class T, class L, class A>
346
347template <int S, class T, class A>
348struct GetFTensor4DdgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
349 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
350#ifndef NDEBUG
351 if (data.size1() != 1)
353 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
354 "of rows should be 1 but is " +
355 boost::lexical_cast<std::string>(data.size1()));
356#endif
357 return FTensor::Ddg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
358 }
359};
360
361template <int S, class T, class A>
362struct GetFTensor4DdgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
363 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
364#ifndef NDEBUG
365 if (data.size1() != 9) {
367 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
368 "of rows should be 9 but is " +
369 boost::lexical_cast<std::string>(data.size1()));
370 }
371#endif
373 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
374 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
375 }
376};
377
378template <int S, class T, class A>
379struct GetFTensor4DdgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
380 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
381#ifndef NDEBUG
382 if (data.size1() != 36) {
383 cerr << data.size1() << endl;
385 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
386 "of rows should be 36 but is " +
387 boost::lexical_cast<std::string>(data.size1()));
388 }
389#endif
391 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
392 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
393 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
394 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
395 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
396 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
397 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
398 &data(35, 0)};
399 }
400};
401
402/**
403 * @brief Get symmetric tensor rank 4 on first two and last indices from
404 * form data matrix
405 *
406 * @tparam Tensor_Dim01 dimension of first two indicies
407 * @tparam Tensor_Dim23 dimension of second two indicies
408 * @tparam Memory shift
409 * @tparam T the type of object stored
410 * @tparam L the storage organization
411 * @tparam A the type of Storage array
412 * @param data data container
413 * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
414 */
415template <int Tensor_Dim01, int Tensor_Dim23, int S = 1, class T, class L,
416 class A>
417static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
418getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
419 static_assert(!std::is_same<T, T>::value,
420 "Such getFTensor4DdgFromMat specialisation is not implemented");
421}
422
423template <int Tensor_Dim01, int Tensor_Dim23, int S = 1>
424static inline auto getFTensor4DdgFromMat(MatrixDouble &data) {
425 return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, double,
426 ublas::row_major,
427 DoubleAllocator>::get(data);
428}
429
430template <int Tensor_Dim01, int Tensor_Dim23, int S, class T>
432
433template <int S, class T> struct GetFTensor4DdgFromPtrImpl<3, 3, S, T> {
434 static inline auto get(T *ptr) {
436
437 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5,
438 ptr + 6, ptr + 7, ptr + 8, ptr + 9, ptr + 10, ptr + 11,
439 ptr + 12, ptr + 13, ptr + 14, ptr + 15, ptr + 16, ptr + 17,
440 ptr + 18, ptr + 19, ptr + 20, ptr + 21, ptr + 22, ptr + 23,
441 ptr + 24, ptr + 25, ptr + 26, ptr + 27, ptr + 28, ptr + 29,
442 ptr + 30, ptr + 31, ptr + 32, ptr + 33, ptr + 34, ptr + 35};
443 }
444};
445
446template <int Tensor_Dim01, int Tensor_Dim23, int S = 1, class T = double>
450
451template <int Tensor_Dim01, int Tensor_Dim2, int S, class T, class L, class A>
453
454template <int S, class T, class A>
455struct GetFTensor3DgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
456 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
457#ifndef NDEBUG
458 if (data.size1() != 1)
460 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
461 "of rows should be 1 but is " +
462 boost::lexical_cast<std::string>(data.size1()));
463#endif
464 return FTensor::Dg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
465 }
466};
467
468template <int S, class T, class A>
469struct GetFTensor3DgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
470 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
471#ifndef NDEBUG
472 if (data.size1() != 6) {
474 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
475 "of rows should be 6 but is " +
476 boost::lexical_cast<std::string>(data.size1()));
477 }
478#endif
480 &data(0, 0), &data(1, 0), &data(2, 0),
481 &data(3, 0), &data(4, 0), &data(5, 0)};
482 }
483};
484
485template <int S, class T, class A>
486struct GetFTensor3DgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
487 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
488#ifndef NDEBUG
489 if (data.size1() != 18) {
490 cerr << data.size1() << endl;
492 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
493 "of rows should be 18 but is " +
494 boost::lexical_cast<std::string>(data.size1()));
495 }
496#endif
498 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
499 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
500 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
501 &data(15, 0), &data(16, 0), &data(17, 0)};
502 }
503};
504
505/**
506 * @brief Get symmetric tensor rank 3 on the first two indices from
507 * form data matrix
508 *
509 * @tparam Tensor_Dim01 dimension of first two indicies
510 * @tparam Tensor_Dim2 dimension of last index
511 * @tparam T the type of object stored
512 * @tparam L the storage organization
513 * @tparam A the type of Storage array
514 * @param data data container
515 * @return FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
516 */
517template <int Tensor_Dim01, int Tensor_Dim2, int S = 1, class T, class L,
518 class A>
519static inline FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim2>
520getFTensor3DgFromMat(ublas::matrix<T, L, A> &data) {
521 static_assert(!std::is_same<T, T>::value,
522 "Such getFTensor3DgFromMat specialisation is not implemented");
523}
524
525template <int Tensor_Dim01, int Tensor_Dim2, int S = 1>
526static inline auto getFTensor3DgFromMat(MatrixDouble &data) {
527 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, double,
528 ublas::row_major, DoubleAllocator>::get(data);
529}
530
531template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
532 int S, class T, class L, class A>
534
535template <int S, class T, class A>
536struct GetFTensor4FromMatImpl<1, 1, 1, 1, S, T, ublas::row_major, A> {
537 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
538#ifndef NDEBUG
539 if (data.size1() != 1)
541 "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
542 "of rows should be 1 but is " +
543 boost::lexical_cast<std::string>(data.size1()));
544#endif
546 &data(0, 0)};
547 }
548};
549
550template <int S, class T, class A>
551struct GetFTensor4FromMatImpl<2, 2, 2, 2, S, T, ublas::row_major, A> {
552 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
553#ifndef NDEBUG
554 if (data.size1() != 16) {
556 "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
557 "of rows should be 16 but is " +
558 boost::lexical_cast<std::string>(data.size1()));
559 }
560#endif
562 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
563 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
564 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
565 &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
566 }
567};
568
569template <int S, class T, class A>
570struct GetFTensor4FromMatImpl<3, 3, 3, 3, S, T, ublas::row_major, A> {
571 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
572#ifndef NDEBUG
573 if (data.size1() != 81) {
574 cerr << data.size1() << endl;
576 "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
577 "of rows should be 81 but is " +
578 boost::lexical_cast<std::string>(data.size1()));
579 }
580#endif
582 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
583 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
584 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
585 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
586 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
587 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
588 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
589 &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
590 &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
591 &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
592 &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
593 &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
594 &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
595 &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
596 &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
597 &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
598 &data(80, 0)};
599 }
600};
601
602/**
603 * @brief Get tensor rank 4 (non symmetric) form data matrix
604 *
605 * @tparam Tensor_Dim0 dimension of frirst index
606 * @tparam Tensor_Dim1 dimension of second index
607 * @tparam Tensor_Dim2 dimension of third index
608 * @tparam Tensor_Dim3 dimension of fourth index
609 * @tparam T the type of object stored
610 * @tparam L the storage organization
611 * @tparam A the type of Storage array
612 * @param data data container
613 * @return FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
614 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
615 */
616template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
617 int S = 1, class T, class L, class A>
618static inline FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
619 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
620getFTensor4FromMat(ublas::matrix<T, L, A> &data) {
621 static_assert(!std::is_same<T, T>::value,
622 "Such getFTensor4FromMat specialisation is not implemented");
623}
624
625template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
626 int S = 1>
627static inline auto getFTensor4FromMat(MatrixDouble &data) {
628 return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
629 Tensor_Dim3, S, double, ublas::row_major,
630 DoubleAllocator>::get(data);
631}
632
633template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class T,
634 class L, class A>
636
637template <int S, class T, class A>
638struct GetFTensor3FromMatImpl<1, 1, 1, S, T, ublas::row_major, A> {
639 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
640#ifndef NDEBUG
641 if (data.size1() != 1)
643 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
644 "of rows should be 1 but is " +
645 boost::lexical_cast<std::string>(data.size1()));
646#endif
648 &data(0, 0)};
649 }
650};
651
652template <int S, class T, class A>
653struct GetFTensor3FromMatImpl<2, 2, 2, S, T, ublas::row_major, A> {
654 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
655#ifndef NDEBUG
656 if (data.size1() != 8)
658 "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
659 "of rows should be 8 but is " +
660 boost::lexical_cast<std::string>(data.size1()));
661#endif
663 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
664 &data(5, 0), &data(6, 0), &data(7, 0)
665
666 };
667 }
668};
669
670template <int S, class T, class A>
671struct GetFTensor3FromMatImpl<3, 2, 2, S, T, ublas::row_major, A> {
672 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
673#ifndef NDEBUG
674 if (data.size1() != 12)
676 "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
677 "of rows should be 12 but is " +
678 boost::lexical_cast<std::string>(data.size1()));
679#endif
681 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
682 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
683 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
684 }
685};
686
687template <int S, class T, class A>
688struct GetFTensor3FromMatImpl<2, 2, 3, S, T, ublas::row_major, A> {
689 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
690#ifndef NDEBUG
691 if (data.size1() != 12)
693 "getFTensor3FromMat<2, 2, 3>: wrong size of data matrix, number "
694 "of rows should be 12 but is " +
695 boost::lexical_cast<std::string>(data.size1()));
696#endif
698 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
699 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
700 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
701 }
702};
703
704template <int S, class T, class A>
705struct GetFTensor3FromMatImpl<3, 3, 3, S, T, ublas::row_major, A> {
706 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
707#ifndef NDEBUG
708 if (data.size1() != 27)
710 "getFTensor3FromMat<3, 3, 3>: wrong size of data matrix, number "
711 "of rows should be 27 but is " +
712 boost::lexical_cast<std::string>(data.size1()));
713#endif
715 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
716 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
717 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
718 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
719 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
720 &data(25, 0), &data(26, 0)};
721 }
722};
723
724template <int S, class T, class A>
725struct GetFTensor3FromMatImpl<6, 3, 3, S, T, ublas::row_major, A> {
726 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
727#ifndef NDEBUG
728 if (data.size1() != 54)
730 "getFTensor3FromMat<6, 3, 3>: wrong size of data matrix, number "
731 "of rows should be 54 but is " +
732 boost::lexical_cast<std::string>(data.size1()));
733#endif
734 std::array<double *, 54> ptrs;
735 for (auto i = 0; i != 54; ++i)
736 ptrs[i] = &data(i, 0);
738 }
739};
740
741template <int S, class T, class A>
742struct GetFTensor3FromMatImpl<3, 3, 6, S, T, ublas::row_major, A> {
743 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
744#ifndef NDEBUG
745 if (data.size1() != 54)
747 "getFTensor3FromMat<3, 3, 6>: wrong size of data matrix, number "
748 "of rows should be 54 but is " +
749 boost::lexical_cast<std::string>(data.size1()));
750#endif
751 std::array<double *, 54> ptrs;
752 for (auto i = 0; i != 54; ++i)
753 ptrs[i] = &data(i, 0);
755 }
756};
757
758/**
759 * @brief Get tensor rank 3 (non symmetries) form data matrix
760 *
761 * @tparam Tensor_Dim0 dimension of first index
762 * @tparam Tensor_Dim1 dimension of second index
763 * @tparam Tensor_Dim2 dimension of third index
764 * @tparam S shift size
765 * @tparam T the type of object stored
766 * @tparam L the storage organization
767 * @tparam A the type of Storage array
768 * @param data data container
769 * @return FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
770 Tensor_Dim1, Tensor_Dim2>
771 */
772template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1, class T,
773 class L, class A>
774static inline FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
775 Tensor_Dim1, Tensor_Dim2>
776getFTensor3FromMat(ublas::matrix<T, L, A> &data) {
777 static_assert(!std::is_same<T, T>::value,
778 "Such getFTensor3FromMat specialisation is not implemented");
779}
780
781template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1>
782static inline auto getFTensor3FromMat(MatrixDouble &data) {
783 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
784 double, ublas::row_major,
785 DoubleAllocator>::get(data);
786}
787
788template <int DIM, int S, typename T> struct GetFTensor1FromPtrImpl;
789
790template <int S, typename T> struct GetFTensor1FromPtrImpl<1, S, T> {
792 inline static auto get(T *ptr) {
794 }
795};
796
797template <int S, typename T> struct GetFTensor1FromPtrImpl<2, S, T> {
799 inline static auto get(T *ptr) {
801 ptr + HVEC1);
802 }
803};
804
805template <int S, typename T> struct GetFTensor1FromPtrImpl<3, S, T> {
807 inline static auto get(T *ptr) {
809 ptr + HVEC0, ptr + HVEC1, ptr + HVEC2);
810 }
811};
812
813template <int S, typename T> struct GetFTensor1FromPtrImpl<4, S, T> {
815 inline static auto get(T *ptr) {
816 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 4>(ptr + 0, ptr + 1,
817 ptr + 2, ptr + 3);
818 }
819};
820
821template <int S, typename T> struct GetFTensor1FromPtrImpl<6, S, T> {
823 inline static auto get(T *ptr) {
825 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5);
826 }
827};
828
829/**
830 * @brief Make Tensor1 from pointer
831 *
832 * @tparam DIM
833 * @param ptr
834 * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
835 */
836template <int DIM, int S = DIM>
841
842#ifdef WITH_ADOL_C
843template <int DIM, int S = DIM>
848#endif
849
850template <int DIM, int S = DIM>
852getFTensor1FromPtr(std::complex<double> *ptr) {
854};
855
856template <int DIM1, int DIM2, int S, typename T> struct GetFTensor2FromPtr;
857
858template <int S, typename T> struct GetFTensor2FromPtr<3, 2, S, T> {
860 inline static auto get(T *ptr) {
862 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5);
863 }
864};
865
866template <int S, typename T> struct GetFTensor2FromPtr<6, 6, S, T> {
868 inline static auto get(T *ptr) {
870 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
871 ptr + 8, ptr + 9, ptr + 10, ptr + 11, ptr + 12, ptr + 13, ptr + 14,
872 ptr + 15, ptr + 16, ptr + 17, ptr + 18, ptr + 19, ptr + 20, ptr + 21,
873 ptr + 22, ptr + 23, ptr + 24, ptr + 25, ptr + 26, ptr + 27, ptr + 28,
874 ptr + 29, ptr + 30, ptr + 31, ptr + 32, ptr + 33, ptr + 34, ptr + 35);
875 }
876};
877
878template <int S, typename T> struct GetFTensor2FromPtr<3, 3, S, T> {
880 inline static auto get(T *ptr) {
882 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
883 ptr + 8);
884 }
885};
886
887template <int S, typename T> struct GetFTensor2FromPtr<2, 2, S, T> {
889 inline static auto get(T *ptr) {
890 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 2, 2>(&ptr[0], &ptr[1],
891 &ptr[2], &ptr[3]);
892 }
893};
894
895template <int S, typename T> struct GetFTensor2FromPtr<1, 3, S, T> {
897 inline static auto get(T *ptr) {
898 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 3>(&ptr[0], &ptr[1],
899 &ptr[2]);
900 }
901};
902
903template <int S, typename T> struct GetFTensor2FromPtr<1, 2, S, T> {
905 inline static auto get(T *ptr) {
906 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 2>(&ptr[0], &ptr[1]);
907 }
908};
909
910template <int S, typename T> struct GetFTensor2FromPtr<1, 1, S, T> {
912 inline static auto get(T *ptr) {
913 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 1>(&ptr[0]);
914 }
915};
916
917/**
918 * @brief Make Tensor2 from pointer
919 *
920 * @tparam DIM
921 * @param ptr
922 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
923 */
924template <int DIM1, int DIM2, int S = DIM1 * DIM2>
925inline auto getFTensor2FromPtr(double *ptr) {
927};
928
929/**
930 * @brief Make Tensor2 from pointer
931 *
932 * @tparam DIM
933 * @param ptr
934 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
935 */
936template <int DIM1, int DIM2, int S = DIM1 * DIM2>
937inline auto getFTensor2FromPtr(std::complex<double> *ptr) {
939};
940
941/**
942 * @brief Make Tensor2 for HVec base from pointer
943 *
944 * @tparam DIM
945 * @param ptr
946 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
947 */
948template <int DIM1, int DIM2>
951 static_assert(DIM1 == DIM1 || DIM2 != DIM2,
952 "Such getFTensor2HVecFromPtr is not implemented");
953};
954
955template <>
957 2> inline getFTensor2HVecFromPtr<3, 2>(double *ptr) {
959 ptr + HVEC0_0, ptr + HVEC0_1,
960
961 ptr + HVEC1_0, ptr + HVEC1_1,
962
963 ptr + HVEC2_0, ptr + HVEC2_1);
964};
965
966template <>
968 3> inline getFTensor2HVecFromPtr<3, 3>(double *ptr) {
970 ptr + HVEC0_0, ptr + HVEC0_1, ptr + HVEC0_2,
971
972 ptr + HVEC1_0, ptr + HVEC1_1, ptr + HVEC1_2,
973
974 ptr + HVEC2_0, ptr + HVEC2_1, ptr + HVEC2_2);
975};
976
977/*
978 * @brief Make Tensor3 from pointer
979 *
980 * @tparam DIM
981 * @param ptr
982 * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3>, DIM1,
983 * DIM2, DIM3>
984 */
985template <int DIM1, int DIM2, int DIM3>
987 DIM2, DIM3>
988getFTensor3FromPtr(double *ptr) {
989 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
990 "Such getFTensor3FromPtr is not implemented");
991};
992
993template <>
995getFTensor3FromPtr<3, 2, 2>(double *ptr) {
997 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
998 ptr + 8, ptr + 9, ptr + 10, ptr + 11);
999};
1000
1001template <>
1003getFTensor3FromPtr<3, 3, 3>(double *ptr) {
1005 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
1006 ptr + 8, ptr + 9, ptr + 10, ptr + 11, ptr + 12, ptr + 13, ptr + 14,
1007 ptr + 15, ptr + 16, ptr + 17, ptr + 18, ptr + 19, ptr + 20, ptr + 21,
1008 ptr + 22, ptr + 23, ptr + 24, ptr + 25, ptr + 26);
1009};
1010
1011/*
1012 * @brief Make Tensor3 from pointer
1013 *
1014 * @tparam DIM
1015 * @param ptr
1016 * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3 *DIM4>,
1017 * DIM1, DIM2, DIM3, DIM4>
1018 */
1019template <int DIM1, int DIM2, int DIM3, int DIM4>
1021 DIM1, DIM2, DIM3, DIM4>
1023 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3 || DIM4 != DIM4,
1024 "Such getFTensor4FromPtr is not implemented");
1025};
1026
1027template <>
1029getFTensor4FromPtr<3, 3, 3, 3>(double *ptr) {
1031 ptr + 0, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
1032 ptr + 8, ptr + 9, ptr + 10, ptr + 11, ptr + 12, ptr + 13, ptr + 14,
1033 ptr + 15, ptr + 16, ptr + 17, ptr + 18, ptr + 19, ptr + 20, ptr + 21,
1034 ptr + 22, ptr + 23, ptr + 24, ptr + 25, ptr + 26, ptr + 27, ptr + 28,
1035 ptr + 29, ptr + 30, ptr + 31, ptr + 32, ptr + 33, ptr + 34, ptr + 35,
1036 ptr + 36, ptr + 37, ptr + 38, ptr + 39, ptr + 40, ptr + 41, ptr + 42,
1037 ptr + 43, ptr + 44, ptr + 45, ptr + 46, ptr + 47, ptr + 48, ptr + 49,
1038 ptr + 50, ptr + 51, ptr + 52, ptr + 53, ptr + 54, ptr + 55, ptr + 56,
1039 ptr + 57, ptr + 58, ptr + 59, ptr + 60, ptr + 61, ptr + 62, ptr + 63,
1040 ptr + 64, ptr + 65, ptr + 66, ptr + 67, ptr + 68, ptr + 69, ptr + 70,
1041 ptr + 71, ptr + 72, ptr + 73, ptr + 74, ptr + 75, ptr + 76, ptr + 77,
1042 ptr + 78, ptr + 79, ptr + 80);
1043};
1044
1045/**
1046 * @brief Make symmetric Tensor2 from pointer
1047 *
1048 * @tparam DIM
1049 * @param ptr
1050 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
1051 */
1052template <int DIM>
1054 FTensor::PackPtr<double *, (DIM * (DIM + 1)) / 2>, DIM>
1056 static_assert(DIM, "Such getFTensor2SymmetricFromPtr is not implemented");
1057}
1058
1059template <>
1061getFTensor2SymmetricFromPtr<3>(double *ptr) {
1063 ptr + 0, ptr + 1, ptr + 2,
1064
1065 ptr + 3, ptr + 4,
1066
1067 ptr + 5);
1068};
1069
1070template <>
1072getFTensor2SymmetricFromPtr<2>(double *ptr) {
1074 &ptr[0], &ptr[1], &ptr[2]);
1075};
1076
1077#ifdef WITH_ADOL_C
1078
1079/**
1080 * @brief Make symmetric Tensor2 from pointer
1081 *
1082 * @tparam DIM
1083 * @param ptr
1084 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
1085 */
1086template <int DIM>
1088 FTensor::PackPtr<adouble *, (DIM * (DIM + 1)) / 2>, DIM>
1090 static_assert(DIM, "Such getFTensor2SymmetricFromPtr is not implemented");
1091}
1092
1093template <>
1097 ptr + 0, ptr + 1, ptr + 2,
1098
1099 ptr + 3, ptr + 4,
1100
1101 ptr + 5);
1102};
1103
1104template <>
1108 ptr + 0, ptr + 1, ptr + 2);
1109};
1110
1111#endif
1112
1113/**
1114 * @brief Make symmetric Tensor2 from pointer, taking lower triangle of matrix
1115 *
1116 * @tparam DIM
1117 * @param ptr
1118 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
1119 */
1120template <int DIM>
1123 static_assert(DIM,
1124 "Such getFTensor2SymmetricLowerFromPtr is not implemented");
1125}
1126
1127template <>
1131 ptr + HVEC0_0, ptr + HVEC0_1, ptr + HVEC0_2,
1132
1133 ptr + HVEC1_0, ptr + HVEC1_1,
1134
1135 ptr + HVEC2_2);
1136};
1137
1138template <>
1142 ptr + 0, ptr + 1, ptr + 3);
1143};
1144
1145template <int DIM, int S> struct GetFTensor1FromArray;
1146
1147template <int S> struct GetFTensor1FromArray<1, S> {
1149 template <typename V> static inline auto get(V &data) {
1150 using T = typename std::remove_reference<decltype(data[0])>::type;
1151 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 1>{&data[0]};
1152 }
1153};
1154
1155template <int S> struct GetFTensor1FromArray<2, S> {
1157 template <typename V> static inline auto get(V &data) {
1158 using T = typename std::remove_reference<decltype(data[0])>::type;
1159 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 2>{&data[0], &data[1]};
1160 }
1161};
1162
1163template <int S> struct GetFTensor1FromArray<3, S> {
1165 template <typename V> static inline auto get(V &data) {
1166 using T = typename std::remove_reference<decltype(data[0])>::type;
1167 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 3>{&data[0], &data[1],
1168 &data[2]};
1169 }
1170};
1171
1172template <int S> struct GetFTensor1FromArray<4, S> {
1174 template <typename V> static inline auto get(V &data) {
1175 using T = typename std::remove_reference<decltype(data[0])>::type;
1176 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 4>{&data[0], &data[1],
1177 &data[2], &data[3]};
1178 }
1179};
1180
1181template <int S> struct GetFTensor1FromArray<6, S> {
1183 template <typename V> static inline auto get(V &data) {
1184 using T = typename std::remove_reference<decltype(data[0])>::type;
1186 &data[0], &data[1], &data[2], &data[3], &data[4], &data[5]};
1187 }
1188};
1189
1190template <int S> struct GetFTensor1FromArray<9, S> {
1192 template <typename V> static inline auto get(V &data) {
1193 using T = typename std::remove_reference<decltype(data[0])>::type;
1195 &data[0], &data[1], &data[2], &data[3], &data[4],
1196 &data[5], &data[6], &data[7], &data[8]};
1197 }
1198};
1199
1200/**
1201 * @brief Get FTensor1 from array
1202 *
1203 * \todo Generalise for different arrays and data types
1204 *
1205 * @tparam DIM
1206 * @param data
1207 * @return FTensor::Tensor1<FTensor::PackPtr<double *, S>, DIM>
1208 */
1209template <int DIM, int S> inline auto getFTensor1FromArray(VectorDouble &data) {
1211}
1212
1213/** @copydoc getFTensor1FromArray */
1214template <int DIM, int S = 0>
1216
1217template <> inline auto getFTensor1FromArray<3, 0>(VectorDouble3 &data) {
1219}
1220
1221template <int DIM, int S>
1223getFTensor1FromMat(MatrixDouble &data, const size_t rr);
1224
1225template <>
1227getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1228 return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2>{&data(rr + 0, 0),
1229 &data(rr + 1, 0)};
1230}
1231
1232template <>
1234getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1236 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1237}
1238
1239template <>
1241getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1243 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0), &data(rr + 3, 0)};
1244}
1245
1246template <>
1248getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1250 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0),
1251 &data(rr + 3, 0), &data(rr + 4, 0), &data(rr + 5, 0),
1252 &data(rr + 6, 0), &data(rr + 7, 0), &data(rr + 8, 0)};
1253}
1254
1255/**
1256 * @brief Get FTensor1 from array
1257 *
1258 * \todo Generalise for diffrent arrays and data types
1259 *
1260 * @tparam DIM
1261 * @param data
1262 * @param rr
1263 * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
1264 */
1265template <int DIM, int S>
1268 static_assert(DIM != DIM, "not implemented");
1270}
1271
1272template <>
1275 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
1276 &data(rr + 1, 1)};
1277}
1278
1279template <>
1283 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1284}
1285
1286/**
1287 * @brief Get FTensor2 from array
1288 *
1289 * \note Generalise for other data types
1290 *
1291 * @tparam DIM1
1292 * @tparam DIM2
1293 * @tparam S
1294 * @param data
1295 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
1296 */
1297template <int DIM1, int DIM2, int S, class T, class L, class A>
1299
1300template <int DIM1, int DIM2, class T, class L, class A>
1302
1303template <int S, class T, class L, class A>
1304struct GetFTensor2FromArrayImpl<2, 2, S, T, L, A> {
1306 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1307 const size_t cc) {
1309 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1310
1311 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1312 }
1313};
1314
1315template <int S, class T, class L, class A>
1316struct GetFTensor2FromArrayImpl<3, 3, S, T, L, A> {
1318 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1319 const size_t cc) {
1321 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1322 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1323 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2)};
1324 }
1325};
1326
1327template <class T, class L, class A>
1330 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1331 const size_t cc, const int ss = 0) {
1333 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1334
1335 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
1336 }
1337};
1338
1339template <class T, class L, class A>
1342 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1343 const size_t cc, const int ss = 0) {
1345 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1346 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1347 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
1348 ss);
1349 }
1350};
1351
1352template <int DIM1, int DIM2, int S>
1354getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc = 0) {
1355 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, double, ublas::row_major,
1356 VecAllocator<double>>::get(data, rr, cc);
1357}
1358
1359template <int DIM1, int DIM2>
1361getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc,
1362 const int ss) {
1363 return GetFTensor2FromArrayRawPtrImpl<DIM1, DIM2, double, ublas::row_major,
1364 VecAllocator<double>>::get(data, rr, cc,
1365 ss);
1366}
1367
1368template <int S, typename T, typename L, typename A>
1369inline auto getFTensor2FromArray2by2(ublas::matrix<T, L, A> &data,
1370 const FTensor::Number<S> &,
1371 const size_t rr, const size_t cc = 0) {
1373}
1374
1375template <int S, typename T, typename L, typename A>
1376inline auto getFTensor2FromArray3by3(ublas::matrix<T, L, A> &data,
1377 const FTensor::Number<S> &,
1378 const size_t rr, const size_t cc = 0) {
1380}
1381
1382#ifdef WITH_ADOL_C
1383
1384template <int DIM1, int DIM2, int S>
1385inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
1386 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, adouble, ublas::row_major,
1387 VecAllocator<adouble>>::get(data, rr);
1388}
1389
1390#endif
1391
1392// Template functions for conversion to Voigt notation for 2nd order tensors
1393
1394template <int DIM, typename T> inline auto getVoigtVec(T &t_mat) {
1395 if constexpr (DIM == 3) {
1396 return std::array<double, 9>{t_mat(0, 0), t_mat(1, 1), t_mat(2, 2),
1397 t_mat(0, 1), t_mat(1, 0), t_mat(0, 2),
1398 t_mat(2, 0), t_mat(1, 2), t_mat(2, 1)};
1399 } else {
1400 return std::array<double, 5>{t_mat(0, 0), t_mat(1, 1), 1.0, t_mat(0, 1),
1401 t_mat(1, 0)};
1402 }
1403};
1404
1405template <int DIM, typename T> inline auto getVoigtVecSymm(T &t_mat) {
1406 constexpr double sqr2 = boost::math::constants::root_two<double>();
1407
1408 if constexpr (DIM == 3) {
1409 return std::array<double, 9>{t_mat(0, 0),
1410 t_mat(1, 1),
1411 t_mat(2, 2),
1412 sqr2 * t_mat(0, 1),
1413 sqr2 * t_mat(0, 2),
1414 sqr2 * t_mat(1, 2),
1415 0.0,
1416 0.0,
1417 0.0};
1418 } else {
1419 return std::array<double, 4>{t_mat(0, 0), t_mat(1, 1), 0.0,
1420 sqr2 * t_mat(0, 1)};
1421 }
1422};
1423
1424template <typename T>
1425inline auto getVoigtVecAxisymm(T &t_mat, const double hoop_term) {
1426 array<double, 5> vec{t_mat(0, 0), t_mat(1, 1), 1. + hoop_term, t_mat(0, 1),
1427 t_mat(1, 0)};
1428
1429 return vec;
1430};
1431
1432template <typename T>
1433inline auto getVoigtVecSymmAxisymm(T &t_mat, const double hoop_term) {
1434 constexpr double sqr2 = boost::math::constants::root_two<double>();
1435
1436 array<double, 4> vec_sym{t_mat(0, 0), t_mat(1, 1), hoop_term,
1437 sqr2 * t_mat(0, 1)};
1438 return vec_sym;
1439};
1440
1441// Template function for conversion of symmetric tensor to non-symmetric
1442
1443template <int DIM, class T>
1446 FTensor::Index<'i', DIM> i;
1447 FTensor::Index<'j', DIM> j;
1448 full(i, j) = symm(i, j);
1449 return full;
1450}
1451
1452// list of lapack wrappers
1453/**
1454 * @brief compute matrix inverse with lapack dgetri
1455 *
1456 * @param mat input square matrix / output inverse matrix
1457 * @return MoFEMErrorCode
1458 */
1461
1462 const size_t M = mat.size1();
1463 const size_t N = mat.size2();
1464
1465 if (M != N)
1466 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1467 "The input matrix for inverse computation is not square %ld != %ld",
1468 M, N);
1469
1470 int *ipv = new int[N];
1471 int lwork = N * N;
1472 double *work = new double[lwork];
1473 int info;
1474 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1475 if (info != 0)
1476 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "lapack error info = %d",
1477 info);
1478 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1479 if (info != 0)
1480 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1481 "lapack error info = %d", info);
1482
1483 delete[] ipv;
1484 delete[] work;
1485
1487}
1488/**
1489 * @brief solve linear system with lapack dgesv
1490 *
1491 * @param mat input lhs square matrix / output L and U from the factorization
1492 * @param f input rhs vector / output solution vector
1493 * @return MoFEMErrorCode
1494 */
1497
1498 const size_t M = mat.size1();
1499 const size_t N = mat.size2();
1500
1501 if (M == 0 || M != N)
1502 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1503 "The input matrix for inverse computation is not square %ld != %ld",
1504 M, N);
1505 if (f.size() != M)
1506 f.resize(M, false);
1507
1508 const int nrhs = 1;
1509 int info;
1510 int *ipiv = new int[M];
1511 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1512 &*f.data().begin(), M);
1513 if (info != 0) {
1514 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1515 "error lapack solve dgesv info = %d", info);
1516 }
1517
1518 delete[] ipiv;
1520}
1521
1522/**
1523 * @brief Solve linear system of equations using Lapack
1524 *
1525 * @param mat
1526 * @param f
1527 * @return MoFEMErrorCode
1528 */
1530 VectorDouble &f) {
1532 // copy matrix since on output lapack returns factorisation
1533 auto mat_copy = mat;
1534 CHKERR solveLinearSystem(mat_copy, f);
1536}
1537
1538/**
1539 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1540 *
1541 * @param mat input symmetric matrix
1542 * @param eig output eigen values sorted
1543 * @param eigen_vec output matrix of row eigen vectors
1544 * @return MoFEMErrorCode
1545 */
1547 VectorDouble &eig,
1548 MatrixDouble &eigen_vec) {
1550
1551 const size_t M = mat.size1();
1552 const size_t N = mat.size2();
1553
1554 if (M == 0 || M != N)
1555 SETERRQ(
1556 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1557 "The input matrix for eigen value computation is not square %ld != %ld",
1558 M, N);
1559 if (eig.size() != M)
1560 eig.resize(M, false);
1561
1562 eigen_vec = mat;
1563 const int n = M;
1564 const int lda = M;
1565 const int size = (M + 2) * M;
1566 int lwork = size;
1567 double *work = new double[size];
1568
1569 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
1570 &*eig.data().begin(), work, lwork) > 0)
1571 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1572 "The algorithm failed to compute eigenvalues.");
1573
1574 delete[] work;
1576}
1577/**
1578 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1579 *
1580 * @tparam DIM
1581 * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
1582 * @param eig output eigen values sorted
1583 * @return MoFEMErrorCode
1584 */
1585template <int DIM, typename T1, typename T2>
1586inline MoFEMErrorCode
1590
1591 const int n = DIM;
1592 const int lda = DIM;
1593 const int lwork = (DIM + 2) * DIM;
1594 std::array<double, (DIM + 2) * DIM> work;
1595
1596 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
1597 lwork) > 0)
1598 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1599 "The algorithm failed to compute eigenvalues.");
1601}
1602
1603/**
1604 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1605 *
1606 * @tparam DIM
1607 * @param mat input tensor pointer of size DIM x DIM
1608 * @param eig output eigen values sorted
1609 * @param eigen_vec output matrix of row eigen vectors
1610 * @return MoFEMErrorCode
1611 */
1612template <int DIM, typename T1, typename T2, typename T3>
1613inline MoFEMErrorCode
1616 FTensor::Tensor2<T3, DIM, DIM> &eigen_vec) {
1618 for (int ii = 0; ii != DIM; ii++)
1619 for (int jj = 0; jj != DIM; jj++)
1620 eigen_vec(ii, jj) = mat(ii, jj);
1621
1622 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
1623
1625}
1626
1627/**
1628 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
1629 *
1630 * @tparam T
1631 * @param t
1632 * @return double
1633 */
1634template <typename T> static inline auto determinantTensor3by3(T &t) {
1635 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
1636 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
1637 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
1638}
1639
1640/**
1641 * @brief Calculate the determinant of a 2x2 matrix or a tensor of rank 2
1642 *
1643 * @tparam T
1644 * @param t
1645 * @return double
1646 */
1647template <typename T> static inline auto determinantTensor2by2(T &t) {
1648 return t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
1649}
1650
1651template <typename T, int DIM> struct DeterminantTensorImpl;
1652
1653template <typename T> struct DeterminantTensorImpl<T, 3> {
1654 static inline auto get(T &t) { return determinantTensor3by3(t); }
1655};
1656
1657template <typename T> struct DeterminantTensorImpl<T, 2> {
1658 static auto get(T &t) { return determinantTensor2by2(t); }
1659};
1660
1661/**
1662 * @brief Calculate the determinant of a tensor of rank DIM
1663 */
1664template <typename T, int DIM>
1668
1669/**
1670 * @brief Calculate the determinant of a tensor of rank DIM
1671 */
1672template <typename T, int DIM>
1676
1677/**
1678 * \brief Calculate inverse of tensor rank 2 at integration points
1679
1680 */
1681template <int Tensor_Dim, class T, class L, class A>
1682inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
1683 ublas::vector<T, A> &det_data,
1684 ublas::matrix<T, L, A> &inv_jac_data) {
1686 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1687 "Specialization for this template not yet implemented");
1689}
1690
1691template <>
1692inline MoFEMErrorCode
1694 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
1695
1696/**
1697 * \brief Calculate determinant 3 by 3
1698
1699 */
1700template <class T1, class T2>
1706
1707/**
1708 * \brief Calculate determinant 2 by 2
1709
1710 */
1711template <class T1, class T2>
1717
1718/**
1719 * \brief Calculate matrix inverse 3 by 3
1720
1721 */
1722template <class T1, class T2, class T3>
1723inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
1725 const auto inv_det = 1. / det;
1726 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1727 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1728 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1729 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
1730 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1731 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1732 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
1733 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
1734 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1736}
1737
1738/**
1739 * \brief Calculate matrix inverse 2 by 2
1740
1741 */
1742template <class T1, class T2, class T3>
1743inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
1745 const auto inv_det = 1. / det;
1746 inv_t(0, 0) = t(1, 1) * inv_det;
1747 inv_t(0, 1) = -t(0, 1) * inv_det;
1748 inv_t(1, 0) = -t(1, 0) * inv_det;
1749 inv_t(1, 1) = t(0, 0) * inv_det;
1751}
1752
1753#ifdef WITH_ADOL_C
1754
1755/**
1756 * \brief Calculate matrix inverse, specialization for adouble tensor
1757
1758 */
1759template <>
1760inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
1765 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1766 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1767 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1768 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
1769 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1770 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1771 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
1772 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
1773 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1775}
1776
1777#endif
1778
1779/**
1780 * \brief Calculate matrix inverse, specialization for symmetric tensor
1781
1782 */
1783template <>
1784inline MoFEMErrorCode
1785invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1790 const auto inv_det = 1. / det;
1791 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1792 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1793 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1794 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1795 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1796 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1798}
1799
1800#ifdef WITH_ADOL_C
1801
1802/**
1803 * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1804
1805 */
1806template <>
1807inline MoFEMErrorCode
1808invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
1813 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1814 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1815 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1816 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1817 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1818 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1820}
1821
1822#endif
1823
1824/**
1825 * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1826 tensor
1827
1828 */
1829template <>
1830inline MoFEMErrorCode
1831invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1836 const auto inv_det = 1. / det;
1837 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1838 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1839 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1840 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1841 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1842 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1844}
1845
1846template <typename T1, typename T2, typename T3, int DIM>
1848
1849template <typename T1, typename T2, typename T3>
1850struct InvertTensorImpl<T1, T2, T3, 3> {
1851 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
1852 return invertTensor3by3(t, det, inv_t);
1853 }
1854};
1855
1856template <typename T1, typename T2, typename T3>
1857struct InvertTensorImpl<T1, T2, T3, 2> {
1858 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
1859 return invertTensor2by2(t, det, inv_t);
1860 }
1861};
1862
1863template <typename T1, typename T2, typename T3, int DIM>
1864static inline MoFEMErrorCode
1871
1872template <typename T1, typename T2, typename T3, int DIM>
1873static inline MoFEMErrorCode
1880
1881/**
1882 * @brief Extract entity handle form multi-index container
1883 *
1884 */
1886 template <typename Iterator>
1887 static inline EntityHandle extract(const Iterator &it) {
1888 return (*it)->getEnt();
1889 }
1890};
1891
1892/**
1893 * @brief Insert ordered mofem multi-index into range
1894 *
1895 * \note Inserted range has to be ordered.
1896 *
1897 * \code
1898 * auto hi_rit = refEntsPtr->upper_bound(start);
1899 * auto hi_rit = refEntsPtr->upper_bound(end);
1900 * Range to_erase;
1901 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1902 * \endcode
1903 *
1904 * @tparam Iterator
1905 * @param r
1906 * @param begin_iter
1907 * @param end_iter
1908 * @return moab::Range::iterator
1909 */
1910template <typename Extractor, typename Iterator>
1911moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1912 Iterator end_iter) {
1913 moab::Range::iterator hint = r.begin();
1914 while (begin_iter != end_iter) {
1915 size_t j = 0;
1916 auto bi = Extractor::extract(begin_iter);
1917 Iterator pj = begin_iter;
1918 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1919 ++pj;
1920 ++j;
1921 }
1922 hint = r.insert(hint, bi, bi + (j - 1));
1923 begin_iter = pj;
1924 }
1925 return hint;
1926};
1927
1928/**
1929 * @brief Do nothing, used to rebuild database
1930 *
1931 */
1934 template <typename T> inline void operator()(T &e) {}
1935};
1936
1937/**
1938 * @brief Template used to reconstruct multi-index
1939 *
1940 * @tparam MI multi-index
1941 * @tparam Modifier
1942 * @param mi
1943 * @param mo
1944 * @return MoFEMErrorCode
1945 */
1946template <typename MI, typename MO = Modify_change_nothing>
1948 MO &&mo = Modify_change_nothing()) {
1950 for (auto it = mi.begin(); it != mi.end(); ++it) {
1951 if (!const_cast<MI &>(mi).modify(it, mo))
1952 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1953 "Houston we have a problem");
1954 }
1956}
1957
1959 TempMeshset(moab::Interface &moab) : moab(moab) {
1960 rval = moab.create_meshset(MESHSET_SET, meshset);
1962 }
1964 operator EntityHandle() const { return meshset; }
1965 auto get_ptr() { return &meshset; }
1966
1967private:
1969 rval = moab.delete_entities(&meshset, 1);
1971 }
1973 moab::Interface &moab;
1974};
1975
1976/**
1977 * @brief Create smart pointer to temporary meshset
1978 *
1979 */
1980inline auto get_temp_meshset_ptr(moab::Interface &moab) {
1981 return boost::make_shared<TempMeshset>(moab);
1982};
1983
1984inline auto id_from_handle(const EntityHandle h) {
1985 return static_cast<EntityID>(h & MB_ID_MASK);
1986};
1987
1988/**
1989 * @brief get type from entity handle
1990 *
1991 */
1992inline auto type_from_handle(const EntityHandle h) {
1993 return static_cast<EntityType>(h >> MB_ID_WIDTH);
1994};
1995
1996/**
1997 * @brief get entity handle from type and id
1998 *
1999 */
2000inline auto ent_form_type_and_id(const EntityType type, const EntityID id) {
2001 return (static_cast<EntityHandle>(type) << MB_ID_WIDTH) | id;
2002};
2003
2004/**
2005 * @brief get entity dimension form handle
2006 *
2007 */
2009 return moab::CN::Dimension(type_from_handle(h));
2010};
2011
2012/**
2013 * @brief get entity type name from handle
2014 *
2015 */
2017 return moab::CN::EntityTypeName(type_from_handle(h));
2018};
2019
2020/**
2021 * @brief get field bit id from bit number
2022 *
2023 */
2024inline auto field_bit_from_bit_number(const int bit_number) {
2025 return BitFieldId().set(bit_number - 1);
2026};
2027
2028/**
2029 * @brief Insert ranges
2030 *
2031 * @tparam I
2032 * @param f
2033 * @param s
2034 * @param tester
2035 * @param inserter
2036 * @return auto
2037 */
2038template <typename I>
2039auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
2040 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
2042
2043 auto first = f;
2044 while (first != s)
2045 if (tester(first)) {
2046
2047 auto second = first;
2048 ++second;
2049
2050 while (second != s) {
2051 if (tester(second))
2052 ++second;
2053 else
2054 break;
2055 }
2056
2057 CHKERR inserter(first, second);
2058
2059 first = second;
2060 if (first != s)
2061 ++first;
2062
2063 } else {
2064 ++first;
2065 }
2066
2068}
2069
2070/**
2071 * @brief Create Array
2072 *
2073 * See:
2074 * <a
2075 * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
2076 * stack overflow</a>
2077 *
2078 * @tparam Dest
2079 * @tparam Arg
2080 * @param arg
2081 * @return constexpr auto
2082 */
2083template <typename Dest = void, typename... Arg>
2084constexpr auto make_array(Arg &&...arg) {
2085 if constexpr (std::is_same<void, Dest>::value)
2086 return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
2087 {std::forward<Arg>(arg)...}};
2088 else
2089 return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
2090}
2091
2092template <int DIM> using i_FTIndex = FTensor::Index<'i', DIM>;
2093template <int DIM> using j_FTIndex = FTensor::Index<'j', DIM>;
2094template <int DIM> using k_FTIndex = FTensor::Index<'k', DIM>;
2095template <int DIM> using l_FTIndex = FTensor::Index<'l', DIM>;
2096template <int DIM> using m_FTIndex = FTensor::Index<'m', DIM>;
2097template <int DIM> using n_FTIndex = FTensor::Index<'n', DIM>;
2098template <int DIM> using o_FTIndex = FTensor::Index<'o', DIM>;
2099template <int DIM> using p_FTIndex = FTensor::Index<'p', DIM>;
2100template <int DIM> using x_FTIndex = FTensor::Index<'x', DIM>;
2101template <int DIM> using y_FTIndex = FTensor::Index<'y', DIM>;
2102template <int DIM> using I_FTIndex = FTensor::Index<'I', DIM>;
2103template <int DIM> using J_FTIndex = FTensor::Index<'J', DIM>;
2104template <int DIM> using K_FTIndex = FTensor::Index<'K', DIM>;
2105template <int DIM> using L_FTIndex = FTensor::Index<'L', DIM>;
2106template <int DIM> using M_FTIndex = FTensor::Index<'M', DIM>;
2107template <int DIM> using N_FTIndex = FTensor::Index<'N', DIM>;
2108
2109#define FTENSOR_INDEX(DIM, I) I##_FTIndex<DIM> I;
2110
2111} // namespace MoFEM
2112
2113#endif //__TEMPLATES_HPP__
T data[Tensor_Dim]
#define MB_ID_MASK
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MB_ID_WIDTH
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ HVEC0
@ HVEC1
@ HVEC2
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ HVEC1_1
@ HVEC0_1
@ HVEC1_0
@ HVEC2_1
@ HVEC1_2
@ HVEC2_2
@ HVEC2_0
@ HVEC0_2
@ HVEC0_0
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
constexpr int DIM2
Definition level_set.cpp:22
constexpr int DIM1
Definition level_set.cpp:21
FTensor::Index< 'j', 3 > j
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::vector< T, std::allocator< T > > VecAllocator
Definition Types.hpp:59
VectorBoundedArray< double, 3 > VectorDouble3
Definition Types.hpp:92
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition Types.hpp:42
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition Types.hpp:114
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
UBlasMatrix< adouble > MatrixADouble
Definition Types.hpp:80
VecAllocator< double > DoubleAllocator
Definition Types.hpp:62
UBlasVector< double > VectorDouble
Definition Types.hpp:68
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition Types.hpp:121
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static FTensor::Tensor3< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 3 (non symmetries) form data matrix.
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition Templates.hpp:31
static FTensor::Dg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim2 > getFTensor3DgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 3 on the first two indices from form data matrix.
auto type_from_handle(const EntityHandle h)
get type from entity handle
auto getVoigtVecAxisymm(T &t_mat, const double hoop_term)
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc=0)
auto to_non_symm(const FTensor::Tensor2_symmetric< T, DIM > &symm)
FTensor::Tensor3< FTensor::PackPtr< double *, 12 >, 3, 2, 2 > getFTensor3FromPtr< 3, 2, 2 >(double *ptr)
FTensor::Tensor2< FTensor::PackPtr< double *, DIM1 *DIM2 >, DIM1, DIM2 > getFTensor2HVecFromPtr(double *ptr)
Make Tensor2 for HVec base from pointer.
auto getVoigtVec(T &t_mat)
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
auto id_from_handle(const EntityHandle h)
auto getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 3 >, 2 > getFTensor2SymmetricFromPtr< 2 >(double *ptr)
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
auto getFTensor1FromArray< 3, 0 >(VectorDouble3 &data)
static auto determinantTensor2by2(T &t)
Calculate the determinant of a 2x2 matrix or a tensor of rank 2.
auto type_name_from_handle(const EntityHandle h)
get entity type name from handle
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
MoFEMErrorCode invertTensor3by3< 3, double, ublas::row_major, DoubleAllocator >(MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data)
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 4 >, 2 > getFTensor2SymmetricLowerFromPtr< 2 >(double *ptr)
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
auto getFTensor2FromArray2by2(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
static auto getFTensor4DdgFromPtr(T *ptr)
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
constexpr auto make_array(Arg &&...arg)
Create Array.
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition Templates.hpp:57
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
std::string toString(X x)
static FTensor::Tensor4< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3 > getFTensor4FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 4 (non symmetric) form data matrix.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 6 >, 3 > getFTensor2SymmetricFromPtr< 3 >(double *ptr)
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr)
Get FTensor1 from array.
FTensor::Tensor4< FTensor::PackPtr< double *, DIM1 *DIM2 *DIM3 *DIM4 >, DIM1, DIM2, DIM3, DIM4 > getFTensor4FromPtr(double *ptr)
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
auto getVoigtVecSymm(T &t_mat)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, DIM *DIM >, DIM > getFTensor2SymmetricLowerFromPtr(double *ptr)
Make symmetric Tensor2 from pointer, taking lower triangle of matrix.
FTensor::Tensor4< FTensor::PackPtr< double *, 81 >, 3, 3, 3, 3 > getFTensor4FromPtr< 3, 3, 3, 3 >(double *ptr)
FTensor::Tensor3< FTensor::PackPtr< double *, 27 >, 3, 3, 3 > getFTensor3FromPtr< 3, 3, 3 >(double *ptr)
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
auto getVoigtVecSymmAxisymm(T &t_mat, const double hoop_term)
FTensor::Tensor3< FTensor::PackPtr< double *, DIM1 *DIM2 *DIM3 >, DIM1, DIM2, DIM3 > getFTensor3FromPtr(double *ptr)
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(DIM *(DIM+1))/2 >, DIM > getFTensor2SymmetricFromPtr(double *ptr)
Make symmetric Tensor2 from pointer.
boost::shared_ptr< std::vector< T > > ShardVec
Definition Templates.hpp:10
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromVec(VectorDouble &data)
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 9 >, 3 > getFTensor2SymmetricLowerFromPtr< 3 >(double *ptr)
constexpr IntegrationType I
constexpr AssemblyType A
double h
constexpr double t
plate stiffness
Definition plate.cpp:58
FTensor::Index< 'm', 3 > m
const int N
Definition speed_test.cpp:3
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition Templates.hpp:95
static auto get(ublas::vector< T, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc)
Get FTensor2 from array.
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
static auto get(ublas::matrix< T, L, A > &data, const size_t rr, const size_t cc, const int ss=0)
static auto get(ublas::matrix< T, L, A > &data)
static auto get(ublas::matrix< T, L, A > &data)
static auto get(ublas::matrix< T, L, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
unsigned int operator()(const id_type &value) const
static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t)
static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t)
KeyExtractor1 key1
Definition Templates.hpp:84
result_type operator()(Arg &arg) const
Definition Templates.hpp:79
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition Templates.hpp:75
KeyExtractor2 key2
Definition Templates.hpp:85
KeyExtractor1::result_type result_type
Definition Templates.hpp:73
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition Templates.hpp:89
Do nothing, used to rebuild database.
Extract entity handle form multi-index container.
static EntityHandle extract(const Iterator &it)
EntityHandle meshset
moab::Interface & moab
TempMeshset(moab::Interface &moab)