v0.14.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>(
260 Tensor_Dim2) + ">: 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,
281 DoubleAllocator>::get(data);
282}
283
284template <int Tensor_Dim1, int Tensor_Dim2>
285inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim1, Tensor_Dim2>
287
288/**
289 * Template specialization for getFTensor2FromMat
290 */
291template <>
296
297template <int Tensor_Dim, int S, class T, class L, class A>
299
300template <int S, class T, class L, class A>
302 static inline auto get(ublas::matrix<T, L, A> &data) {
303#ifndef NDEBUG
304 if (data.size1() != 6)
306 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
307 "of rows should be 6 but is " +
308 boost::lexical_cast<std::string>(data.size1()));
309#endif
311 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
312 &data(5, 0));
313 }
314};
315
316template <int S, class T, class L, class A>
318 static inline auto get(ublas::matrix<T, L, A> &data) {
319#ifndef NDEBUG
320 if (data.size1() != 3)
322 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
323 "of rows should be 3 but is " +
324 boost::lexical_cast<std::string>(data.size1()));
325#endif
327 &data(0, 0), &data(1, 0), &data(2, 0));
328 }
329};
330
331/**
332 * \brief Get symmetric tensor rank 2 (matrix) form data matrix
333 */
334template <int Tensor_Dim, int S, class T, class L, class A>
335static inline auto getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
337}
338
339template <int Tensor_Dim, int S = 1>
340static inline auto getFTensor2SymmetricFromMat(MatrixDouble &data) {
341 return getFTensor2SymmetricFromMat<Tensor_Dim, S, double, ublas::row_major,
342 DoubleAllocator>(data);
343}
344
345template <int Tensor_Dim01, int Tensor_Dim23, int S, class T, class L, class A>
347
348template <int S, class T, class A>
349struct GetFTensor4DdgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
350 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
351#ifndef NDEBUG
352 if (data.size1() != 1)
354 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
355 "of rows should be 1 but is " +
356 boost::lexical_cast<std::string>(data.size1()));
357#endif
358 return FTensor::Ddg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
359 }
360};
361
362template <int S, class T, class A>
363struct GetFTensor4DdgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
364 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
365#ifndef NDEBUG
366 if (data.size1() != 9) {
368 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
369 "of rows should be 9 but is " +
370 boost::lexical_cast<std::string>(data.size1()));
371 }
372#endif
374 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
375 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
376 }
377};
378
379template <int S, class T, class A>
380struct GetFTensor4DdgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
381 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
382#ifndef NDEBUG
383 if (data.size1() != 36) {
384 cerr << data.size1() << endl;
386 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
387 "of rows should be 36 but is " +
388 boost::lexical_cast<std::string>(data.size1()));
389 }
390#endif
392 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
393 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
394 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
395 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
396 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
397 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
398 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
399 &data(35, 0)};
400 }
401};
402
403/**
404 * @brief Get symmetric tensor rank 4 on first two and last indices from
405 * form data matrix
406 *
407 * @tparam Tensor_Dim01 dimension of first two indicies
408 * @tparam Tensor_Dim23 dimension of second two indicies
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_Dim2, int S, class T, class L, class A>
432
433template <int S, class T, class A>
434struct GetFTensor3DgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
435 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
436#ifndef NDEBUG
437 if (data.size1() != 1)
439 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
440 "of rows should be 1 but is " +
441 boost::lexical_cast<std::string>(data.size1()));
442#endif
443 return FTensor::Dg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
444 }
445};
446
447template <int S, class T, class A>
448struct GetFTensor3DgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
449 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
450#ifndef NDEBUG
451 if (data.size1() != 6) {
453 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
454 "of rows should be 6 but is " +
455 boost::lexical_cast<std::string>(data.size1()));
456 }
457#endif
459 &data(0, 0), &data(1, 0), &data(2, 0),
460 &data(3, 0), &data(4, 0), &data(5, 0)};
461 }
462};
463
464template <int S, class T, class A>
465struct GetFTensor3DgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
466 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
467#ifndef NDEBUG
468 if (data.size1() != 18) {
469 cerr << data.size1() << endl;
471 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
472 "of rows should be 18 but is " +
473 boost::lexical_cast<std::string>(data.size1()));
474 }
475#endif
477 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
478 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
479 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
480 &data(15, 0), &data(16, 0), &data(17, 0)};
481 }
482};
483
484/**
485 * @brief Get symmetric tensor rank 3 on the first two indices from
486 * form data matrix
487 *
488 * @tparam Tensor_Dim01 dimension of first two indicies
489 * @tparam Tensor_Dim2 dimension of last index
490 * @tparam T the type of object stored
491 * @tparam L the storage organization
492 * @tparam A the type of Storage array
493 * @param data data container
494 * @return FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
495 */
496template <int Tensor_Dim01, int Tensor_Dim2, int S = 1, class T, class L,
497 class A>
498static inline FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim2>
499getFTensor3DgFromMat(ublas::matrix<T, L, A> &data) {
500 static_assert(!std::is_same<T, T>::value,
501 "Such getFTensor3DgFromMat specialisation is not implemented");
502}
503
504template <int Tensor_Dim01, int Tensor_Dim2, int S = 1>
505static inline auto getFTensor3DgFromMat(MatrixDouble &data) {
506 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, double,
507 ublas::row_major, DoubleAllocator>::get(data);
508}
509
510template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
511 int S, class T, class L, class A>
513
514template <int S, class T, class A>
515struct GetFTensor4FromMatImpl<1, 1, 1, 1, S, T, ublas::row_major, A> {
516 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
517#ifndef NDEBUG
518 if (data.size1() != 1)
520 "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
521 "of rows should be 1 but is " +
522 boost::lexical_cast<std::string>(data.size1()));
523#endif
525 &data(0, 0)};
526 }
527};
528
529template <int S, class T, class A>
530struct GetFTensor4FromMatImpl<2, 2, 2, 2, S, T, ublas::row_major, A> {
531 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
532#ifndef NDEBUG
533 if (data.size1() != 16) {
535 "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
536 "of rows should be 16 but is " +
537 boost::lexical_cast<std::string>(data.size1()));
538 }
539#endif
541 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
542 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
543 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
544 &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
545 }
546};
547
548template <int S, class T, class A>
549struct GetFTensor4FromMatImpl<3, 3, 3, 3, S, T, ublas::row_major, A> {
550 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
551#ifndef NDEBUG
552 if (data.size1() != 81) {
553 cerr << data.size1() << endl;
555 "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
556 "of rows should be 81 but is " +
557 boost::lexical_cast<std::string>(data.size1()));
558 }
559#endif
561 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
562 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
563 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
564 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
565 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
566 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
567 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
568 &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
569 &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
570 &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
571 &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
572 &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
573 &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
574 &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
575 &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
576 &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
577 &data(80, 0)};
578 }
579};
580
581/**
582 * @brief Get tensor rank 4 (non symmetric) form data matrix
583 *
584 * @tparam Tensor_Dim0 dimension of frirst index
585 * @tparam Tensor_Dim1 dimension of second index
586 * @tparam Tensor_Dim2 dimension of third index
587 * @tparam Tensor_Dim3 dimension of fourth index
588 * @tparam T the type of object stored
589 * @tparam L the storage organization
590 * @tparam A the type of Storage array
591 * @param data data container
592 * @return FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
593 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
594 */
595template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
596 int S = 1, class T, class L, class A>
597static inline FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
598 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
599getFTensor4FromMat(ublas::matrix<T, L, A> &data) {
600 static_assert(!std::is_same<T, T>::value,
601 "Such getFTensor4FromMat specialisation is not implemented");
602}
603
604template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
605 int S = 1>
606static inline auto getFTensor4FromMat(MatrixDouble &data) {
607 return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
608 Tensor_Dim3, S, double, ublas::row_major,
609 DoubleAllocator>::get(data);
610}
611
612template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class T,
613 class L, class A>
615
616template <int S, class T, class A>
617struct GetFTensor3FromMatImpl<1, 1, 1, S, T, ublas::row_major, A> {
618 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
619#ifndef NDEBUG
620 if (data.size1() != 1)
622 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
623 "of rows should be 1 but is " +
624 boost::lexical_cast<std::string>(data.size1()));
625#endif
627 &data(0, 0)};
628 }
629};
630
631template <int S, class T, class A>
632struct GetFTensor3FromMatImpl<2, 2, 2, S, T, ublas::row_major, A> {
633 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
634#ifndef NDEBUG
635 if (data.size1() != 8)
637 "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
638 "of rows should be 8 but is " +
639 boost::lexical_cast<std::string>(data.size1()));
640#endif
642 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
643 &data(5, 0), &data(6, 0), &data(7, 0)
644
645 };
646 }
647};
648
649template <int S, class T, class A>
650struct GetFTensor3FromMatImpl<3, 2, 2, S, T, ublas::row_major, A> {
651 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
652#ifndef NDEBUG
653 if (data.size1() != 12)
655 "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
656 "of rows should be 12 but is " +
657 boost::lexical_cast<std::string>(data.size1()));
658#endif
660 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
661 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
662 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
663 }
664};
665
666template <int S, class T, class A>
667struct GetFTensor3FromMatImpl<2, 2, 3, S, T, ublas::row_major, A> {
668 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
669#ifndef NDEBUG
670 if (data.size1() != 12)
672 "getFTensor3FromMat<2, 2, 3>: wrong size of data matrix, number "
673 "of rows should be 12 but is " +
674 boost::lexical_cast<std::string>(data.size1()));
675#endif
677 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
678 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
679 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
680 }
681};
682
683template <int S, class T, class A>
684struct GetFTensor3FromMatImpl<3, 3, 3, S, T, ublas::row_major, A> {
685 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
686#ifndef NDEBUG
687 if (data.size1() != 27)
689 "getFTensor3FromMat<3, 3, 3>: wrong size of data matrix, number "
690 "of rows should be 27 but is " +
691 boost::lexical_cast<std::string>(data.size1()));
692#endif
694 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
695 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
696 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
697 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
698 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
699 &data(25, 0), &data(26, 0)};
700 }
701};
702
703template <int S, class T, class A>
704struct GetFTensor3FromMatImpl<6, 3, 3, S, T, ublas::row_major, A> {
705 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
706#ifndef NDEBUG
707 if (data.size1() != 54)
709 "getFTensor3FromMat<6, 3, 3>: wrong size of data matrix, number "
710 "of rows should be 54 but is " +
711 boost::lexical_cast<std::string>(data.size1()));
712#endif
713 std::array<double *, 54> ptrs;
714 for (auto i = 0; i != 54; ++i)
715 ptrs[i] = &data(i, 0);
717 }
718};
719
720template <int S, class T, class A>
721struct GetFTensor3FromMatImpl<3, 3, 6, S, T, ublas::row_major, A> {
722 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
723#ifndef NDEBUG
724 if (data.size1() != 54)
726 "getFTensor3FromMat<3, 3, 6>: wrong size of data matrix, number "
727 "of rows should be 54 but is " +
728 boost::lexical_cast<std::string>(data.size1()));
729#endif
730 std::array<double *, 54> ptrs;
731 for (auto i = 0; i != 54; ++i)
732 ptrs[i] = &data(i, 0);
734 }
735};
736
737/**
738 * @brief Get tensor rank 3 (non symmetries) form data matrix
739 *
740 * @tparam Tensor_Dim0 dimension of first index
741 * @tparam Tensor_Dim1 dimension of second index
742 * @tparam Tensor_Dim2 dimension of third index
743 * @tparam S shift size
744 * @tparam T the type of object stored
745 * @tparam L the storage organization
746 * @tparam A the type of Storage array
747 * @param data data container
748 * @return FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
749 Tensor_Dim1, Tensor_Dim2>
750 */
751template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1, class T,
752 class L, class A>
753static inline FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
754 Tensor_Dim1, Tensor_Dim2>
755getFTensor3FromMat(ublas::matrix<T, L, A> &data) {
756 static_assert(!std::is_same<T, T>::value,
757 "Such getFTensor3FromMat specialisation is not implemented");
758}
759
760template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1>
761static inline auto getFTensor3FromMat(MatrixDouble &data) {
762 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
763 double, ublas::row_major,
764 DoubleAllocator>::get(data);
765}
766
767template <int DIM, int S, typename T> struct GetFTensor1FromPtrImpl;
768
769template <int S, typename T> struct GetFTensor1FromPtrImpl<2, S, T> {
771 inline static auto get(T *ptr) {
773 &ptr[HVEC1]);
774 }
775};
776
777template <int S, typename T> struct GetFTensor1FromPtrImpl<3, S, T> {
779 inline static auto get(T *ptr) {
781 &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
782 }
783};
784
785template <int S, typename T> struct GetFTensor1FromPtrImpl<6, S, T> {
787 inline static auto get(T *ptr) {
789 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
790 }
791};
792
793/**
794 * @brief Make Tensor1 from pointer
795 *
796 * @tparam DIM
797 * @param ptr
798 * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
799 */
800template <int DIM, int S = DIM>
805
806#ifdef WITH_ADOL_C
807template <int DIM, int S = DIM>
812#endif
813
814template <int DIM, int S = DIM>
816getFTensor1FromPtr(std::complex<double> *ptr) {
818};
819
820template <int DIM1, int DIM2, int S, typename T> struct GetFTensor2FromPtr;
821
822template <int S, typename T> struct GetFTensor2FromPtr<3, 2, S, T> {
824 inline static auto get(T *ptr) {
826 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5]);
827 }
828};
829
830template <int S, typename T> struct GetFTensor2FromPtr<6, 6, S, T> {
832 inline static auto get(T *ptr) {
834 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
835 &ptr[8], &ptr[9], &ptr[10], &ptr[11], &ptr[12], &ptr[13], &ptr[14],
836 &ptr[15], &ptr[16], &ptr[17], &ptr[18], &ptr[19], &ptr[20], &ptr[21],
837 &ptr[22], &ptr[23], &ptr[24], &ptr[25], &ptr[26], &ptr[27], &ptr[28],
838 &ptr[29], &ptr[30], &ptr[31], &ptr[32], &ptr[33], &ptr[34], &ptr[35]);
839 }
840};
841
842template <int S, typename T> struct GetFTensor2FromPtr<3, 3, S, T> {
844 inline static auto get(T *ptr) {
846 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
847 &ptr[8]);
848 }
849};
850
851template <int S, typename T> struct GetFTensor2FromPtr<2, 2, S, T> {
853 inline static auto get(T *ptr) {
854 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 2, 2>(&ptr[0], &ptr[1],
855 &ptr[2], &ptr[3]);
856 }
857};
858
859template <int S, typename T> struct GetFTensor2FromPtr<1, 3, S, T> {
861 inline static auto get(T *ptr) {
862 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 3>(&ptr[0], &ptr[1],
863 &ptr[2]);
864 }
865};
866
867template <int S, typename T> struct GetFTensor2FromPtr<1, 2, S, T> {
869 inline static auto get(T *ptr) {
870 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 2>(&ptr[0], &ptr[1]);
871 }
872};
873
874template <int S, typename T> struct GetFTensor2FromPtr<1, 1, S, T> {
876 inline static auto get(T *ptr) {
877 return FTensor::Tensor2<FTensor::PackPtr<T *, S>, 1, 1>(&ptr[0]);
878 }
879};
880
881/**
882 * @brief Make Tensor2 from pointer
883 *
884 * @tparam DIM
885 * @param ptr
886 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
887 */
888template <int DIM1, int DIM2, int S = DIM1 * DIM2>
889inline auto getFTensor2FromPtr(double *ptr) {
891};
892
893/**
894 * @brief Make Tensor2 from pointer
895 *
896 * @tparam DIM
897 * @param ptr
898 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
899 */
900template <int DIM1, int DIM2, int S = DIM1 * DIM2>
901inline auto getFTensor2FromPtr(std::complex<double> *ptr) {
903};
904
905/**
906 * @brief Make Tensor2 for HVec base from pointer
907 *
908 * @tparam DIM
909 * @param ptr
910 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
911 */
912template <int DIM1, int DIM2>
915 static_assert(DIM1 == DIM1 || DIM2 != DIM2,
916 "Such getFTensor2HVecFromPtr is not implemented");
917};
918
919template <>
921 2> inline getFTensor2HVecFromPtr<3, 2>(double *ptr) {
923 &ptr[HVEC0_0], &ptr[HVEC0_1],
924
925 &ptr[HVEC1_0], &ptr[HVEC1_1],
926
927 &ptr[HVEC2_0], &ptr[HVEC2_1]);
928};
929
930template <>
932 3> inline getFTensor2HVecFromPtr<3, 3>(double *ptr) {
934 &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
935
936 &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
937
938 &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
939};
940
941/*
942 * @brief Make Tensor3 from pointer
943 *
944 * @tparam DIM
945 * @param ptr
946 * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3>, DIM1,
947 * DIM2, DIM3>
948 */
949template <int DIM1, int DIM2, int DIM3>
951 DIM2, DIM3>
952getFTensor3FromPtr(double *ptr) {
953 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
954 "Such getFTensor3FromPtr is not implemented");
955};
956
957template <>
959getFTensor3FromPtr<3, 2, 2>(double *ptr) {
961 &ptr[0], &ptr[1], &ptr[2],
962
963 &ptr[3], &ptr[4], &ptr[5],
964
965 &ptr[6], &ptr[7], &ptr[8],
966
967 &ptr[9], &ptr[10], &ptr[11]
968
969 );
970};
971
972template <>
974getFTensor3FromPtr<3, 3, 3>(double *ptr) {
976 &ptr[0], &ptr[1], &ptr[2], &ptr[3], &ptr[4], &ptr[5], &ptr[6], &ptr[7],
977 &ptr[8], &ptr[9], &ptr[10], &ptr[11], &ptr[12], &ptr[13], &ptr[14],
978 &ptr[15], &ptr[16], &ptr[17], &ptr[18], &ptr[19], &ptr[20], &ptr[21],
979 &ptr[22], &ptr[23], &ptr[24], &ptr[25], &ptr[26]);
980};
981
982/**
983 * @brief Make symmetric Tensor2 from pointer
984 *
985 * @tparam DIM
986 * @param ptr
987 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
988 */
989template <int DIM>
991 FTensor::PackPtr<double *, (DIM * (DIM + 1)) / 2>, DIM>
993 static_assert(DIM, "Such getFTensor2SymmetricFromPtr is not implemented");
994}
995
996template <>
1000 &ptr[0], &ptr[1], &ptr[2],
1001
1002 &ptr[3], &ptr[4],
1003
1004 &ptr[5]);
1005};
1006
1007template <>
1009getFTensor2SymmetricFromPtr<2>(double *ptr) {
1011 &ptr[0], &ptr[1], &ptr[2]);
1012};
1013
1014#ifdef WITH_ADOL_C
1015
1016/**
1017 * @brief Make symmetric Tensor2 from pointer
1018 *
1019 * @tparam DIM
1020 * @param ptr
1021 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
1022 */
1023template <int DIM>
1025 FTensor::PackPtr<adouble *, (DIM * (DIM + 1)) / 2>, DIM>
1027 static_assert(DIM, "Such getFTensor2SymmetricFromPtr is not implemented");
1028}
1029
1030template <>
1034 &ptr[0], &ptr[1], &ptr[2],
1035
1036 &ptr[3], &ptr[4],
1037
1038 &ptr[5]);
1039};
1040
1041template <>
1045 &ptr[0], &ptr[1], &ptr[2]);
1046};
1047
1048#endif
1049
1050/**
1051 * @brief Make symmetric Tensor2 from pointer, taking lower triangle of matrix
1052 *
1053 * @tparam DIM
1054 * @param ptr
1055 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
1056 */
1057template <int DIM>
1060 static_assert(DIM,
1061 "Such getFTensor2SymmetricLowerFromPtr is not implemented");
1062}
1063
1064template <>
1068 &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
1069
1070 &ptr[HVEC1_0], &ptr[HVEC1_1],
1071
1072 &ptr[HVEC2_2]);
1073};
1074
1075template <>
1079 &ptr[0], &ptr[1], &ptr[3]);
1080};
1081
1082template <int DIM, int S> struct GetFTensor1FromArray;
1083
1084template <int S> struct GetFTensor1FromArray<1, S> {
1086 template <typename V> static inline auto get(V &data) {
1087 using T = typename std::remove_reference<decltype(data[0])>::type;
1088 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 1>{&data[0]};
1089 }
1090};
1091
1092template <int S> struct GetFTensor1FromArray<2, S> {
1094 template <typename V> static inline auto get(V &data) {
1095 using T = typename std::remove_reference<decltype(data[0])>::type;
1096 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 2>{&data[0], &data[1]};
1097 }
1098};
1099
1100template <int S> struct GetFTensor1FromArray<3, S> {
1102 template <typename V> static inline auto get(V &data) {
1103 using T = typename std::remove_reference<decltype(data[0])>::type;
1104 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 3>{&data[0], &data[1],
1105 &data[2]};
1106 }
1107};
1108
1109template <int S> struct GetFTensor1FromArray<4, S> {
1111 template <typename V> static inline auto get(V &data) {
1112 using T = typename std::remove_reference<decltype(data[0])>::type;
1113 return FTensor::Tensor1<FTensor::PackPtr<T *, S>, 4>{&data[0], &data[1],
1114 &data[2], &data[3]};
1115 }
1116};
1117
1118template <int S> struct GetFTensor1FromArray<6, S> {
1120 template <typename V> static inline auto get(V &data) {
1121 using T = typename std::remove_reference<decltype(data[0])>::type;
1123 &data[0], &data[1], &data[2], &data[3], &data[4], &data[5]};
1124 }
1125};
1126
1127template <int S> struct GetFTensor1FromArray<9, S> {
1129 template <typename V> static inline auto get(V &data) {
1130 using T = typename std::remove_reference<decltype(data[0])>::type;
1132 &data[0], &data[1], &data[2], &data[3], &data[4],
1133 &data[5], &data[6], &data[7], &data[8]};
1134 }
1135};
1136
1137/**
1138 * @brief Get FTensor1 from array
1139 *
1140 * \todo Generalise for different arrays and data types
1141 *
1142 * @tparam DIM
1143 * @param data
1144 * @return FTensor::Tensor1<FTensor::PackPtr<double *, S>, DIM>
1145 */
1146template <int DIM, int S> inline auto getFTensor1FromArray(VectorDouble &data) {
1148}
1149
1150/** @copydoc getFTensor1FromArray */
1151template <int DIM, int S = 0>
1153
1154template <> inline auto getFTensor1FromArray<3, 0>(VectorDouble3 &data) {
1156}
1157
1158template <int DIM, int S>
1160getFTensor1FromMat(MatrixDouble &data, const size_t rr);
1161
1162template <>
1164getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1165 return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2>{&data(rr + 0, 0),
1166 &data(rr + 1, 0)};
1167}
1168
1169template <>
1171getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1173 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
1174}
1175
1176template <>
1178getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1180 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0), &data(rr + 3, 0)};
1181}
1182
1183template <>
1185getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
1187 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0),
1188 &data(rr + 3, 0), &data(rr + 4, 0), &data(rr + 5, 0),
1189 &data(rr + 6, 0), &data(rr + 7, 0), &data(rr + 8, 0)};
1190}
1191
1192/**
1193 * @brief Get FTensor1 from array
1194 *
1195 * \todo Generalise for diffrent arrays and data types
1196 *
1197 * @tparam DIM
1198 * @param data
1199 * @param rr
1200 * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
1201 */
1202template <int DIM, int S>
1205 static_assert(DIM != DIM, "not implemented");
1207}
1208
1209template <>
1212 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
1213 &data(rr + 1, 1)};
1214}
1215
1216template <>
1220 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
1221}
1222
1223/**
1224 * @brief Get FTensor2 from array
1225 *
1226 * \note Generalise for other data types
1227 *
1228 * @tparam DIM1
1229 * @tparam DIM2
1230 * @tparam S
1231 * @param data
1232 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
1233 */
1234template <int DIM1, int DIM2, int S, class T, class L, class A>
1236
1237template <int DIM1, int DIM2, class T, class L, class A>
1239
1240template <int S, class T, class L, class A>
1241struct GetFTensor2FromArrayImpl<2, 2, S, T, L, A> {
1243 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1244 const size_t cc) {
1246 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1247
1248 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1)};
1249 }
1250};
1251
1252template <int S, class T, class L, class A>
1253struct GetFTensor2FromArrayImpl<3, 3, S, T, L, A> {
1255 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1256 const size_t cc) {
1258 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1259 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1260 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2)};
1261 }
1262};
1263
1264template <class T, class L, class A>
1267 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1268 const size_t cc, const int ss = 0) {
1270 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1),
1271
1272 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), ss);
1273 }
1274};
1275
1276template <class T, class L, class A>
1279 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr,
1280 const size_t cc, const int ss = 0) {
1282 &data(rr + 0, cc + 0), &data(rr + 0, cc + 1), &data(rr + 0, cc + 2),
1283 &data(rr + 1, cc + 0), &data(rr + 1, cc + 1), &data(rr + 1, cc + 2),
1284 &data(rr + 2, cc + 0), &data(rr + 2, cc + 1), &data(rr + 2, cc + 2),
1285 ss);
1286 }
1287};
1288
1289template <int DIM1, int DIM2, int S>
1291getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc = 0) {
1292 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, double, ublas::row_major,
1293 VecAllocator<double>>::get(data, rr, cc);
1294}
1295
1296template <int DIM1, int DIM2>
1298getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc,
1299 const int ss) {
1300 return GetFTensor2FromArrayRawPtrImpl<DIM1, DIM2, double, ublas::row_major,
1301 VecAllocator<double>>::get(data, rr, cc,
1302 ss);
1303}
1304
1305template <int S, typename T, typename L, typename A>
1306inline auto getFTensor2FromArray2by2(ublas::matrix<T, L, A> &data,
1307 const FTensor::Number<S> &,
1308 const size_t rr, const size_t cc = 0) {
1310}
1311
1312template <int S, typename T, typename L, typename A>
1313inline auto getFTensor2FromArray3by3(ublas::matrix<T, L, A> &data,
1314 const FTensor::Number<S> &,
1315 const size_t rr, const size_t cc = 0) {
1317}
1318
1319#ifdef WITH_ADOL_C
1320
1321template <int DIM1, int DIM2, int S>
1322inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
1323 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, adouble, ublas::row_major,
1324 VecAllocator<adouble>>::get(data, rr);
1325}
1326
1327#endif
1328
1329// list of lapack wrappers
1330/**
1331 * @brief compute matrix inverse with lapack dgetri
1332 *
1333 * @param mat input square matrix / output inverse matrix
1334 * @return MoFEMErrorCode
1335 */
1338
1339 const size_t M = mat.size1();
1340 const size_t N = mat.size2();
1341
1342 if (M != N)
1343 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1344 "The input matrix for inverse computation is not square %d != %d",
1345 M, N);
1346
1347 int *ipv = new int[N];
1348 int lwork = N * N;
1349 double *work = new double[lwork];
1350 int info;
1351 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1352 if (info != 0)
1353 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1354 "lapack error info = %d", info);
1355 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1356 if (info != 0)
1357 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1358 "lapack error info = %d", info);
1359
1360 delete[] ipv;
1361 delete[] work;
1362
1364}
1365/**
1366 * @brief solve linear system with lapack dgesv
1367 *
1368 * @param mat input lhs square matrix / output L and U from the factorization
1369 * @param f input rhs vector / output solution vector
1370 * @return MoFEMErrorCode
1371 */
1374
1375 const size_t M = mat.size1();
1376 const size_t N = mat.size2();
1377
1378 if (M == 0 || M != N)
1379 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1380 "The input matrix for inverse computation is not square %d != %d",
1381 M, N);
1382 if (f.size() != M)
1383 f.resize(M, false);
1384
1385 const int nrhs = 1;
1386 int info;
1387 int *ipiv = new int[M];
1388 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1389 &*f.data().begin(), nrhs);
1390
1391 if (info != 0) {
1392 SETERRQ1(PETSC_COMM_SELF, 1, "error lapack solve dgesv info = %d", info);
1393 }
1394
1395 delete[] ipiv;
1397}
1398
1399/**
1400 * @brief Solve linear system of equations using Lapack
1401 *
1402 * @param mat
1403 * @param f
1404 * @return MoFEMErrorCode
1405 */
1407 VectorDouble &f) {
1409 // copy matrix since on output lapack returns factorisation
1410 auto mat_copy = mat;
1411 CHKERR solveLinearSystem(mat_copy, f);
1413}
1414
1415/**
1416 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1417 *
1418 * @param mat input symmetric matrix
1419 * @param eig output eigen values sorted
1420 * @param eigen_vec output matrix of row eigen vectors
1421 * @return MoFEMErrorCode
1422 */
1424 VectorDouble &eig,
1425 MatrixDouble &eigen_vec) {
1427
1428 const size_t M = mat.size1();
1429 const size_t N = mat.size2();
1430
1431 if (M == 0 || M != N)
1432 SETERRQ2(
1433 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1434 "The input matrix for eigen value computation is not square %d != %d",
1435 M, N);
1436 if (eig.size() != M)
1437 eig.resize(M, false);
1438
1439 eigen_vec = mat;
1440 const int n = M;
1441 const int lda = M;
1442 const int size = (M + 2) * M;
1443 int lwork = size;
1444 double *work = new double[size];
1445
1446 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
1447 &*eig.data().begin(), work, lwork) > 0)
1448 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1449 "The algorithm failed to compute eigenvalues.");
1450
1451 delete[] work;
1453}
1454/**
1455 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1456 *
1457 * @tparam DIM
1458 * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
1459 * @param eig output eigen values sorted
1460 * @return MoFEMErrorCode
1461 */
1462template <int DIM, typename T1, typename T2>
1463inline MoFEMErrorCode
1467
1468 const int n = DIM;
1469 const int lda = DIM;
1470 const int lwork = (DIM + 2) * DIM;
1471 std::array<double, (DIM + 2) * DIM> work;
1472
1473 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
1474 lwork) > 0)
1475 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1476 "The algorithm failed to compute eigenvalues.");
1478}
1479
1480/**
1481 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1482 *
1483 * @tparam DIM
1484 * @param mat input tensor pointer of size DIM x DIM
1485 * @param eig output eigen values sorted
1486 * @param eigen_vec output matrix of row eigen vectors
1487 * @return MoFEMErrorCode
1488 */
1489template <int DIM, typename T1, typename T2, typename T3>
1490inline MoFEMErrorCode
1493 FTensor::Tensor2<T3, DIM, DIM> &eigen_vec) {
1495 for (int ii = 0; ii != DIM; ii++)
1496 for (int jj = 0; jj != DIM; jj++)
1497 eigen_vec(ii, jj) = mat(ii, jj);
1498
1499 CHKERR computeEigenValuesSymmetric(eigen_vec, eig);
1500
1502}
1503
1504/**
1505 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
1506 *
1507 * @tparam T
1508 * @param t
1509 * @return double
1510 */
1511template <typename T> static inline auto determinantTensor3by3(T &t) {
1512 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
1513 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
1514 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
1515}
1516
1517/**
1518 * @brief Calculate the determinant of a 2x2 matrix or a tensor of rank 2
1519 *
1520 * @tparam T
1521 * @param t
1522 * @return double
1523 */
1524template <typename T> static inline auto determinantTensor2by2(T &t) {
1525 return t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
1526}
1527
1528template <typename T, int DIM> struct DeterminantTensorImpl;
1529
1530template <typename T> struct DeterminantTensorImpl<T, 3> {
1531 static inline auto get(T &t) { return determinantTensor3by3(t); }
1532};
1533
1534template <typename T> struct DeterminantTensorImpl<T, 2> {
1535 static auto get(T &t) { return determinantTensor2by2(t); }
1536};
1537
1538/**
1539 * @brief Calculate the determinant of a tensor of rank DIM
1540 */
1541template <typename T, int DIM>
1545
1546/**
1547 * @brief Calculate the determinant of a tensor of rank DIM
1548 */
1549template <typename T, int DIM>
1553
1554/**
1555 * \brief Calculate inverse of tensor rank 2 at integration points
1556
1557 */
1558template <int Tensor_Dim, class T, class L, class A>
1559inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
1560 ublas::vector<T, A> &det_data,
1561 ublas::matrix<T, L, A> &inv_jac_data) {
1563 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1564 "Specialization for this template not yet implemented");
1566}
1567
1568template <>
1569inline MoFEMErrorCode
1571 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
1572
1573/**
1574 * \brief Calculate determinant 3 by 3
1575
1576 */
1577template <class T1, class T2>
1583
1584/**
1585 * \brief Calculate determinant 2 by 2
1586
1587 */
1588template <class T1, class T2>
1594
1595/**
1596 * \brief Calculate matrix inverse 3 by 3
1597
1598 */
1599template <class T1, class T2, class T3>
1600inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
1602 const auto inv_det = 1. / det;
1603 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1604 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1605 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1606 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
1607 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1608 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1609 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
1610 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
1611 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1613}
1614
1615/**
1616 * \brief Calculate matrix inverse 2 by 2
1617
1618 */
1619template <class T1, class T2, class T3>
1620inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
1622 const auto inv_det = 1. / det;
1623 inv_t(0, 0) = t(1, 1) * inv_det;
1624 inv_t(0, 1) = -t(0, 1) * inv_det;
1625 inv_t(1, 0) = -t(1, 0) * inv_det;
1626 inv_t(1, 1) = t(0, 0) * inv_det;
1628}
1629
1630#ifdef WITH_ADOL_C
1631
1632/**
1633 * \brief Calculate matrix inverse, specialization for adouble tensor
1634
1635 */
1636template <>
1637inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
1642 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1643 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1644 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1645 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
1646 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1647 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1648 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
1649 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
1650 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1652}
1653
1654#endif
1655
1656/**
1657 * \brief Calculate matrix inverse, specialization for symmetric tensor
1658
1659 */
1660template <>
1661inline MoFEMErrorCode
1662invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1667 const auto inv_det = 1. / det;
1668 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1669 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1670 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1671 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1672 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1673 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1675}
1676
1677#ifdef WITH_ADOL_C
1678
1679/**
1680 * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1681
1682 */
1683template <>
1684inline MoFEMErrorCode
1685invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
1690 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1691 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1692 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1693 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1694 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1695 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1697}
1698
1699#endif
1700
1701/**
1702 * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1703 tensor
1704
1705 */
1706template <>
1707inline MoFEMErrorCode
1708invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1713 const auto inv_det = 1. / det;
1714 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1715 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1716 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1717 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1718 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1719 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1721}
1722
1723template <typename T1, typename T2, typename T3, int DIM>
1725
1726template <typename T1, typename T2, typename T3>
1727struct InvertTensorImpl<T1, T2, T3, 3> {
1728 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
1729 return invertTensor3by3(t, det, inv_t);
1730 }
1731};
1732
1733template <typename T1, typename T2, typename T3>
1734struct InvertTensorImpl<T1, T2, T3, 2> {
1735 inline static MoFEMErrorCode invert(T1 &t, T2 &det, T3 &inv_t) {
1736 return invertTensor2by2(t, det, inv_t);
1737 }
1738};
1739
1740template <typename T1, typename T2, typename T3, int DIM>
1741static inline MoFEMErrorCode
1748
1749template <typename T1, typename T2, typename T3, int DIM>
1750static inline MoFEMErrorCode
1757
1758/**
1759 * @brief Extract entity handle form multi-index container
1760 *
1761 */
1763 template <typename Iterator>
1764 static inline EntityHandle extract(const Iterator &it) {
1765 return (*it)->getEnt();
1766 }
1767};
1768
1769/**
1770 * @brief Insert ordered mofem multi-index into range
1771 *
1772 * \note Inserted range has to be ordered.
1773 *
1774 * \code
1775 * auto hi_rit = refEntsPtr->upper_bound(start);
1776 * auto hi_rit = refEntsPtr->upper_bound(end);
1777 * Range to_erase;
1778 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1779 * \endcode
1780 *
1781 * @tparam Iterator
1782 * @param r
1783 * @param begin_iter
1784 * @param end_iter
1785 * @return moab::Range::iterator
1786 */
1787template <typename Extractor, typename Iterator>
1788moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1789 Iterator end_iter) {
1790 moab::Range::iterator hint = r.begin();
1791 while (begin_iter != end_iter) {
1792 size_t j = 0;
1793 auto bi = Extractor::extract(begin_iter);
1794 Iterator pj = begin_iter;
1795 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1796 ++pj;
1797 ++j;
1798 }
1799 hint = r.insert(hint, bi, bi + (j - 1));
1800 begin_iter = pj;
1801 }
1802 return hint;
1803};
1804
1805/**
1806 * @brief Do nothing, used to rebuild database
1807 *
1808 */
1811 template <typename T> inline void operator()(T &e) {}
1812};
1813
1814/**
1815 * @brief Template used to reconstruct multi-index
1816 *
1817 * @tparam MI multi-index
1818 * @tparam Modifier
1819 * @param mi
1820 * @param mo
1821 * @return MoFEMErrorCode
1822 */
1823template <typename MI, typename MO = Modify_change_nothing>
1825 MO &&mo = Modify_change_nothing()) {
1827 for (auto it = mi.begin(); it != mi.end(); ++it) {
1828 if (!const_cast<MI &>(mi).modify(it, mo))
1829 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1830 "Houston we have a problem");
1831 }
1833}
1834
1836 TempMeshset(moab::Interface &moab) : moab(moab) {
1837 rval = moab.create_meshset(MESHSET_SET, meshset);
1839 }
1841 operator EntityHandle() const { return meshset; }
1842 auto get_ptr() { return &meshset; }
1843
1844private:
1846 rval = moab.delete_entities(&meshset, 1);
1848 }
1850 moab::Interface &moab;
1851};
1852
1853/**
1854 * @brief Create smart pointer to temporary meshset
1855 *
1856 */
1857inline auto get_temp_meshset_ptr(moab::Interface &moab) {
1858 return boost::make_shared<TempMeshset>(moab);
1859};
1860
1861inline auto id_from_handle(const EntityHandle h) {
1862 return static_cast<EntityID>(h & MB_ID_MASK);
1863};
1864
1865/**
1866 * @brief get type from entity handle
1867 *
1868 */
1869inline auto type_from_handle(const EntityHandle h) {
1870 return static_cast<EntityType>(h >> MB_ID_WIDTH);
1871};
1872
1873/**
1874 * @brief get entity handle from type and id
1875 *
1876 */
1877inline auto ent_form_type_and_id(const EntityType type, const EntityID id) {
1878 return (static_cast<EntityHandle>(type) << MB_ID_WIDTH) | id;
1879};
1880
1881/**
1882 * @brief get entity dimension form handle
1883 *
1884 */
1886 return moab::CN::Dimension(type_from_handle(h));
1887};
1888
1889/**
1890 * @brief get entity type name from handle
1891 *
1892 */
1894 return moab::CN::EntityTypeName(type_from_handle(h));
1895};
1896
1897/**
1898 * @brief get field bit id from bit number
1899 *
1900 */
1901inline auto field_bit_from_bit_number(const int bit_number) {
1902 return BitFieldId().set(bit_number - 1);
1903};
1904
1905/**
1906 * @brief Insert ranges
1907 *
1908 * @tparam I
1909 * @param f
1910 * @param s
1911 * @param tester
1912 * @param inserter
1913 * @return auto
1914 */
1915template <typename I>
1916auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
1917 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
1919
1920 auto first = f;
1921 while (first != s)
1922 if (tester(first)) {
1923
1924 auto second = first;
1925 ++second;
1926
1927 while (second != s) {
1928 if (tester(second))
1929 ++second;
1930 else
1931 break;
1932 }
1933
1934 CHKERR inserter(first, second);
1935
1936 first = second;
1937 if (first != s)
1938 ++first;
1939
1940 } else {
1941 ++first;
1942 }
1943
1945}
1946
1947/**
1948 * @brief Create Array
1949 *
1950 * See:
1951 * <a
1952 * href="https://stackoverflow.com/questions/50942556/current-status-of-stdmake-array">See
1953 * stack overflow</a>
1954 *
1955 * @tparam Dest
1956 * @tparam Arg
1957 * @param arg
1958 * @return constexpr auto
1959 */
1960template <typename Dest = void, typename... Arg>
1961constexpr auto make_array(Arg &&...arg) {
1962 if constexpr (std::is_same<void, Dest>::value)
1963 return std::array<std::common_type_t<std::decay_t<Arg>...>, sizeof...(Arg)>{
1964 {std::forward<Arg>(arg)...}};
1965 else
1966 return std::array<Dest, sizeof...(Arg)>{{std::forward<Arg>(arg)...}};
1967}
1968
1969} // namespace MoFEM
1970
1971#endif //__TEMPLATES_HPP__
static Index< 'M', 3 > M
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_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< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
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
const double T
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)
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)
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(DIM *(DIM+1))/2 >, DIM getFTensor2SymmetricFromPtr)(double *ptr)
Make symmetric Tensor2 from pointer.
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.
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::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< 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:59
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)