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// list of lapack wrappers
1393/**
1394 * @brief compute matrix inverse with lapack dgetri
1395 *
1396 * @param mat input square matrix / output inverse matrix
1397 * @return MoFEMErrorCode
1398 */
1401
1402 const size_t M = mat.size1();
1403 const size_t N = mat.size2();
1404
1405 if (M != N)
1406 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1407 "The input matrix for inverse computation is not square %ld != %ld",
1408 M, N);
1409
1410 int *ipv = new int[N];
1411 int lwork = N * N;
1412 double *work = new double[lwork];
1413 int info;
1414 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1415 if (info != 0)
1416 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1417 "lapack error info = %d", info);
1418 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1419 if (info != 0)
1420 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1421 "lapack error info = %d", info);
1422
1423 delete[] ipv;
1424 delete[] work;
1425
1427}
1428/**
1429 * @brief solve linear system with lapack dgesv
1430 *
1431 * @param mat input lhs square matrix / output L and U from the factorization
1432 * @param f input rhs vector / output solution vector
1433 * @return MoFEMErrorCode
1434 */
1437
1438 const size_t M = mat.size1();
1439 const size_t N = mat.size2();
1440
1441 if (M == 0 || M != N)
1442 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1443 "The input matrix for inverse computation is not square %ld != %ld",
1444 M, N);
1445 if (f.size() != M)
1446 f.resize(M, false);
1447
1448 const int nrhs = 1;
1449 int info;
1450 int *ipiv = new int[M];
1451 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1452 &*f.data().begin(), M);
1453 if (info != 0) {
1454 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1455 "error lapack solve dgesv info = %d", info);
1456 }
1457
1458 delete[] ipiv;
1460}
1461
1462/**
1463 * @brief Solve linear system of equations using Lapack
1464 *
1465 * @param mat
1466 * @param f
1467 * @return MoFEMErrorCode
1468 */
1470 VectorDouble &f) {
1472 // copy matrix since on output lapack returns factorisation
1473 auto mat_copy = mat;
1474 CHKERR solveLinearSystem(mat_copy, f);
1476}
1477
1478/**
1479 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1480 *
1481 * @param mat input symmetric matrix
1482 * @param eig output eigen values sorted
1483 * @param eigen_vec output matrix of row eigen vectors
1484 * @return MoFEMErrorCode
1485 */
1487 VectorDouble &eig,
1488 MatrixDouble &eigen_vec) {
1490
1491 const size_t M = mat.size1();
1492 const size_t N = mat.size2();
1493
1494 if (M == 0 || M != N)
1495 SETERRQ(
1496 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1497 "The input matrix for eigen value computation is not square %ld != %ld",
1498 M, N);
1499 if (eig.size() != M)
1500 eig.resize(M, false);
1501
1502 eigen_vec = mat;
1503 const int n = M;
1504 const int lda = M;
1505 const int size = (M + 2) * M;
1506 int lwork = size;
1507 double *work = new double[size];
1508
1509 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
1510 &*eig.data().begin(), work, lwork) > 0)
1511 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1512 "The algorithm failed to compute eigenvalues.");
1513
1514 delete[] work;
1516}
1517/**
1518 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1519 *
1520 * @tparam DIM
1521 * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
1522 * @param eig output eigen values sorted
1523 * @return MoFEMErrorCode
1524 */
1525template <int DIM, typename T1, typename T2>
1526inline MoFEMErrorCode
1530
1531 const int n = DIM;
1532 const int lda = DIM;
1533 const int lwork = (DIM + 2) * DIM;
1534 std::array<double, (DIM + 2) * DIM> work;
1535
1536 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
1537 lwork) > 0)
1538 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1539 "The algorithm failed to compute eigenvalues.");
1541}
1542
1543/**
1544 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1545 *
1546 * @tparam DIM
1547 * @param mat input tensor pointer of size DIM x DIM
1548 * @param eig output eigen values sorted
1549 * @param eigen_vec output matrix of row eigen vectors
1550 * @return MoFEMErrorCode
1551 */
1552template <int DIM, typename T1, typename T2, typename T3>
1553inline MoFEMErrorCode
1556 FTensor::Tensor2<T3, DIM, DIM> &eigen_vec) {
1558 for (int ii = 0; ii != DIM; ii++)
1559 for (int jj = 0; jj != DIM; jj++)
1560 eigen_vec(ii, jj) = mat(ii, jj);
1561
1562 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
1563
1565}
1566
1567/**
1568 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
1569 *
1570 * @tparam T
1571 * @param t
1572 * @return double
1573 */
1574template <typename T> static inline auto determinantTensor3by3(T &t) {
1575 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
1576 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
1577 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
1578}
1579
1580/**
1581 * @brief Calculate the determinant of a 2x2 matrix or a tensor of rank 2
1582 *
1583 * @tparam T
1584 * @param t
1585 * @return double
1586 */
1587template <typename T> static inline auto determinantTensor2by2(T &t) {
1588 return t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
1589}
1590
1591template <typename T, int DIM> struct DeterminantTensorImpl;
1592
1593template <typename T> struct DeterminantTensorImpl<T, 3> {
1594 static inline auto get(T &t) { return determinantTensor3by3(t); }
1595};
1596
1597template <typename T> struct DeterminantTensorImpl<T, 2> {
1598 static auto get(T &t) { return determinantTensor2by2(t); }
1599};
1600
1601/**
1602 * @brief Calculate the determinant of a tensor of rank DIM
1603 */
1604template <typename T, int DIM>
1608
1609/**
1610 * @brief Calculate the determinant of a tensor of rank DIM
1611 */
1612template <typename T, int DIM>
1616
1617/**
1618 * \brief Calculate inverse of tensor rank 2 at integration points
1619
1620 */
1621template <int Tensor_Dim, class T, class L, class A>
1622inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
1623 ublas::vector<T, A> &det_data,
1624 ublas::matrix<T, L, A> &inv_jac_data) {
1626 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1627 "Specialization for this template not yet implemented");
1629}
1630
1631template <>
1632inline MoFEMErrorCode
1634 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
1635
1636/**
1637 * \brief Calculate determinant 3 by 3
1638
1639 */
1640template <class T1, class T2>
1646
1647/**
1648 * \brief Calculate determinant 2 by 2
1649
1650 */
1651template <class T1, class T2>
1657
1658/**
1659 * \brief Calculate matrix inverse 3 by 3
1660
1661 */
1662template <class T1, class T2, class T3>
1663inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
1665 const auto inv_det = 1. / det;
1666 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1667 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1668 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1669 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
1670 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1671 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1672 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
1673 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
1674 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1676}
1677
1678/**
1679 * \brief Calculate matrix inverse 2 by 2
1680
1681 */
1682template <class T1, class T2, class T3>
1683inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
1685 const auto inv_det = 1. / det;
1686 inv_t(0, 0) = t(1, 1) * inv_det;
1687 inv_t(0, 1) = -t(0, 1) * inv_det;
1688 inv_t(1, 0) = -t(1, 0) * inv_det;
1689 inv_t(1, 1) = t(0, 0) * inv_det;
1691}
1692
1693#ifdef WITH_ADOL_C
1694
1695/**
1696 * \brief Calculate matrix inverse, specialization for adouble tensor
1697
1698 */
1699template <>
1705 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1706 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1707 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1708 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
1709 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1710 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1711 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
1712 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
1713 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1715}
1716
1717#endif
1718
1719/**
1720 * \brief Calculate matrix inverse, specialization for symmetric tensor
1721
1722 */
1723template <>
1724inline MoFEMErrorCode
1730 const auto inv_det = 1. / det;
1731 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1732 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1733 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1734 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1735 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1736 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1738}
1739
1740#ifdef WITH_ADOL_C
1741
1742/**
1743 * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1744
1745 */
1746template <>
1747inline MoFEMErrorCode
1753 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1754 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1755 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1756 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1757 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1758 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1760}
1761
1762#endif
1763
1764/**
1765 * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1766 tensor
1767
1768 */
1769template <>
1770inline MoFEMErrorCode
1776 const auto inv_det = 1. / det;
1777 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1778 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1779 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1780 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1781 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1782 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1784}
1785
1786template <typename T1, typename T2, typename T3, int DIM>
1788
1789template <typename T1, typename T2, typename T3>
1790struct InvertTensorImpl<T1, T2, T3, 3> {
1791 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
1792 return invertTensor3by3(t, det, inv_t);
1793 }
1794};
1795
1796template <typename T1, typename T2, typename T3>
1797struct InvertTensorImpl<T1, T2, T3, 2> {
1798 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
1799 return invertTensor2by2(t, det, inv_t);
1800 }
1801};
1802
1803template <typename T1, typename T2, typename T3, int DIM>
1804static inline MoFEMErrorCode
1811
1812template <typename T1, typename T2, typename T3, int DIM>
1813static inline MoFEMErrorCode
1820
1821/**
1822 * @brief Extract entity handle form multi-index container
1823 *
1824 */
1826 template <typename Iterator>
1827 static inline EntityHandle extract(const Iterator &it) {
1828 return (*it)->getEnt();
1829 }
1830};
1831
1832/**
1833 * @brief Insert ordered mofem multi-index into range
1834 *
1835 * \note Inserted range has to be ordered.
1836 *
1837 * \code
1838 * auto hi_rit = refEntsPtr->upper_bound(start);
1839 * auto hi_rit = refEntsPtr->upper_bound(end);
1840 * Range to_erase;
1841 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1842 * \endcode
1843 *
1844 * @tparam Iterator
1845 * @param r
1846 * @param begin_iter
1847 * @param end_iter
1848 * @return moab::Range::iterator
1849 */
1850template <typename Extractor, typename Iterator>
1851moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1852 Iterator end_iter) {
1853 moab::Range::iterator hint = r.begin();
1854 while (begin_iter != end_iter) {
1855 size_t j = 0;
1856 auto bi = Extractor::extract(begin_iter);
1857 Iterator pj = begin_iter;
1858 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1859 ++pj;
1860 ++j;
1861 }
1862 hint = r.insert(hint, bi, bi + (j - 1));
1863 begin_iter = pj;
1864 }
1865 return hint;
1866};
1867
1868/**
1869 * @brief Do nothing, used to rebuild database
1870 *
1871 */
1874 template <typename T> inline void operator()(T &e) {}
1875};
1876
1877/**
1878 * @brief Template used to reconstruct multi-index
1879 *
1880 * @tparam MI multi-index
1881 * @tparam Modifier
1882 * @param mi
1883 * @param mo
1884 * @return MoFEMErrorCode
1885 */
1886template <typename MI, typename MO = Modify_change_nothing>
1888 MO &&mo = Modify_change_nothing()) {
1890 for (auto it = mi.begin(); it != mi.end(); ++it) {
1891 if (!const_cast<MI &>(mi).modify(it, mo))
1892 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1893 "Houston we have a problem");
1894 }
1896}
1897
1899 TempMeshset(moab::Interface &moab) : moab(moab) {
1900 rval = moab.create_meshset(MESHSET_SET, meshset);
1902 }
1904 operator EntityHandle() const { return meshset; }
1905 auto get_ptr() { return &meshset; }
1906
1907private:
1909 rval = moab.delete_entities(&meshset, 1);
1911 }
1913 moab::Interface &moab;
1914};
1915
1916/**
1917 * @brief Create smart pointer to temporary meshset
1918 *
1919 */
1920inline auto get_temp_meshset_ptr(moab::Interface &moab) {
1921 return boost::make_shared<TempMeshset>(moab);
1922};
1923
1924inline auto id_from_handle(const EntityHandle h) {
1925 return static_cast<EntityID>(h & MB_ID_MASK);
1926};
1927
1928/**
1929 * @brief get type from entity handle
1930 *
1931 */
1932inline auto type_from_handle(const EntityHandle h) {
1933 return static_cast<EntityType>(h >> MB_ID_WIDTH);
1934};
1935
1936/**
1937 * @brief get entity handle from type and id
1938 *
1939 */
1940inline auto ent_form_type_and_id(const EntityType type, const EntityID id) {
1941 return (static_cast<EntityHandle>(type) << MB_ID_WIDTH) | id;
1942};
1943
1944/**
1945 * @brief get entity dimension form handle
1946 *
1947 */
1949 return moab::CN::Dimension(type_from_handle(h));
1950};
1951
1952/**
1953 * @brief get entity type name from handle
1954 *
1955 */
1957 return moab::CN::EntityTypeName(type_from_handle(h));
1958};
1959
1960/**
1961 * @brief get field bit id from bit number
1962 *
1963 */
1964inline auto field_bit_from_bit_number(const int bit_number) {
1965 return BitFieldId().set(bit_number - 1);
1966};
1967
1968/**
1969 * @brief Insert ranges
1970 *
1971 * @tparam I
1972 * @param f
1973 * @param s
1974 * @param tester
1975 * @param inserter
1976 * @return auto
1977 */
1978template <typename I>
1979auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
1980 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
1982
1983 auto first = f;
1984 while (first != s)
1985 if (tester(first)) {
1986
1987 auto second = first;
1988 ++second;
1989
1990 while (second != s) {
1991 if (tester(second))
1992 ++second;
1993 else
1994 break;
1995 }
1996
1997 CHKERR inserter(first, second);
1998
1999 first = second;
2000 if (first != s)
2001 ++first;
2002
2003 } else {
2004 ++first;
2005 }
2006
2008}
2009
2010/**
2011 * @brief Create Array
2012 *
2013 * See:
2014 * <a
2015 * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
2016 * stack overflow</a>
2017 *
2018 * @tparam Dest
2019 * @tparam Arg
2020 * @param arg
2021 * @return constexpr auto
2022 */
2023template <typename Dest = void, typename... Arg>
2024constexpr auto make_array(Arg &&...arg) {
2025 if constexpr (std::is_same<void, Dest>::value)
2026 return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
2027 {std::forward<Arg>(arg)...}};
2028 else
2029 return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
2030}
2031
2032template <int DIM> using i_FTIndex = FTensor::Index<'i', DIM>;
2033template <int DIM> using j_FTIndex = FTensor::Index<'j', DIM>;
2034template <int DIM> using k_FTIndex = FTensor::Index<'k', DIM>;
2035template <int DIM> using l_FTIndex = FTensor::Index<'l', DIM>;
2036template <int DIM> using m_FTIndex = FTensor::Index<'m', DIM>;
2037template <int DIM> using n_FTIndex = FTensor::Index<'n', DIM>;
2038template <int DIM> using o_FTIndex = FTensor::Index<'o', DIM>;
2039template <int DIM> using p_FTIndex = FTensor::Index<'p', DIM>;
2040template <int DIM> using x_FTIndex = FTensor::Index<'x', DIM>;
2041template <int DIM> using y_FTIndex = FTensor::Index<'x', DIM>;
2042template <int DIM> using I_FTIndex = FTensor::Index<'I', DIM>;
2043template <int DIM> using J_FTIndex = FTensor::Index<'J', DIM>;
2044template <int DIM> using K_FTIndex = FTensor::Index<'K', DIM>;
2045template <int DIM> using L_FTIndex = FTensor::Index<'L', DIM>;
2046template <int DIM> using M_FTIndex = FTensor::Index<'M', DIM>;
2047template <int DIM> using N_FTIndex = FTensor::Index<'N', DIM>;
2048
2049#define FTENSOR_INDEX(DIM, I) I##_FTIndex<DIM> I;
2050
2051} // namespace MoFEM
2052
2053#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.
VecAllocator< double > DoubleAllocator
Definition Types.hpp:62
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition Types.hpp:119
std::vector< T, std::allocator< T > > VecAllocator
Definition Types.hpp:59
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition Types.hpp:42
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
VectorBoundedArray< double, 3 > VectorDouble3
Definition Types.hpp:92
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition Types.hpp:113
UBlasVector< double > VectorDouble
Definition Types.hpp:68
UBlasMatrix< adouble > MatrixADouble
Definition Types.hpp:80
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
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)
boost::shared_ptr< std::vector< T > > ShardVec
Definition Templates.hpp:10
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.
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 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.
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.
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
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
KeyExtractor1::result_type result_type
Definition Templates.hpp:73
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition Templates.hpp:75
KeyExtractor2 key2
Definition Templates.hpp:85
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)