v0.13.1
Templates.hpp
Go to the documentation of this file.
1/** \file Templates.hpp
2 * \brief Templates declarations
3 *
4 * MoFEM is free software: you can redistribute it and/or modify it under
5 * the terms of the GNU Lesser General Public License as published by the
6 * Free Software Foundation, either version 3 of the License, or (at your
7 * option) any later version.
8 *
9 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12 * License for more details.
13 *
14 * You should have received a copy of the GNU Lesser General Public
15 *
16 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
17 */
18
19#ifndef __TEMPLATES_HPP__
20#define __TEMPLATES_HPP__
21
22namespace MoFEM {
23
24template <typename T> using ShardVec = boost::shared_ptr<std::vector<T>>;
25
26/**
27 * @brief Get Vector adaptor
28 *
29 * \code
30 *
31 * double *a;
32 * CHKERR VecGetArray(v,&a);
33 *
34 * for(int n = 0; n != nodes; ++n) {
35 *
36 * auto a = getVectorAdaptor(&a[3*n], 3);
37 * double dot = inner_prod(a, a);
38 *
39 * }
40 *
41 * CHKERR VecRetsoreArray(v,&a);
42 * \endcode
43 *
44 */
45template <typename T1> inline auto getVectorAdaptor(T1 ptr, const size_t n) {
46 typedef typename std::remove_pointer<T1>::type T;
48 ublas::shallow_array_adaptor<T>(n, ptr));
49};
50
51/**
52 * @brief Get Matrix adaptor
53 *
54 * \code
55 *
56 * double *a;
57 * CHKERR VecGetArray(v,&a);
58 *
59 * for(int n = 0; n != nodes; ++n) {
60 *
61 * auto F = getMatrixAdaptor(&a[3*3*n], 3, 3);
62 * MatrixDouble C = prod(F, trans(F));
63 *
64 * }
65 *
66 * CHKERR VecRetsoreArray(v,&a);
67 * \endcode
68 *
69 */
70template <typename T1>
71inline auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m) {
72 typedef typename std::remove_pointer<T1>::type T;
74 n, m, ublas::shallow_array_adaptor<T>(n * m, ptr));
75};
76
77/**
78 * This small utility that cascades two key extractors will be
79 * used throughout the boost example
80 * <a
81 * href=http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp>
82 * http://www.boost.org/doc/libs/1_53_0/libs/multi_index/example/complex_structs.cpp
83 * </a>
84 */
85template <class KeyExtractor1, class KeyExtractor2> struct KeyFromKey {
86public:
87 typedef typename KeyExtractor1::result_type result_type;
88
89 KeyFromKey(const KeyExtractor1 &key1_ = KeyExtractor1(),
90 const KeyExtractor2 &key2_ = KeyExtractor2())
91 : key1(key1_), key2(key2_) {}
92
93 template <typename Arg> result_type operator()(Arg &arg) const {
94 return key1(key2(arg));
95 }
96
97private:
98 KeyExtractor1 key1;
99 KeyExtractor2 key2;
100};
101
102template <typename id_type> struct LtBit {
103 inline bool operator()(const id_type &valueA, const id_type &valueB) const {
104 return valueA.to_ulong() < valueB.to_ulong();
105 }
106};
107
108template <typename id_type> struct EqBit {
109 inline bool operator()(const id_type &valueA, const id_type &valueB) const {
110 return valueA.to_ulong() == valueB.to_ulong();
111 }
112};
113
114template <typename id_type> struct HashBit {
115 inline unsigned int operator()(const id_type &value) const {
116 return value.to_ulong();
117 }
118};
119
120template <class X> inline std::string toString(X x) {
121 std::ostringstream buffer;
122 buffer << x;
123 return buffer.str();
124}
125
126template <int S, class T, class A> struct GetFTensor0FromVecImpl {
127 static inline auto get(ublas::vector<T, A> &data) {
128 return FTensor::Tensor0<FTensor::PackPtr<T *, S>>(&*data.data().begin());
129 }
130};
131
132/**
133* \brief Get tensor rank 0 (scalar) form data vector
134
135Example how to use it.
136\code
137VectorDouble vec;
138vec.resize(nb_gauss_pts,false);
139vec.clear();
140auto t0 = getFTensor0FromData(data);
141for(int gg = 0;gg!=nb_gauss_pts;gg++) {
142
143 ++t0;
144}
145\endcode
146
147*/
148template <int S = 1, class T, class A>
149static inline auto getFTensor0FromVec(ublas::vector<T, A> &data) {
151}
152
153template <int Tensor_Dim, int S, class T, class L, class A>
155
156template <int S, class T, class A>
157struct GetFTensor1FromMatImpl<3, S, T, ublas::row_major, A> {
158 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
159#ifndef NDDEBUG
160 if (data.size1() != 3)
162 "getFTensor1FromMat<3>: wrong size of data matrix, number of "
163 "rows should be 3 but is %d" +
164 boost::lexical_cast<std::string>(data.size1()));
165#endif
167 &data(0, 0), &data(1, 0), &data(2, 0));
168 }
169};
170
171template <int S, class T, class A>
172struct GetFTensor1FromMatImpl<2, S, T, ublas::row_major, A> {
173 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
174#ifndef NDDEBUG
175 if (data.size1() != 2)
177 "getFTensor1FromMat<2>: wrong size of data matrix, number of "
178 "rows should be 2 but is %d" +
179 boost::lexical_cast<std::string>(data.size1()));
180#endif
182 &data(1, 0));
183 }
184};
185
186template <int S, class T, class A>
187struct GetFTensor1FromMatImpl<1, S, T, ublas::row_major, A> {
188 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
189#ifndef NDEBUG
190 if (data.size1() != 1)
192 "getFTensor1FromMat<1>: wrong size of data matrix, number of "
193 "rows should be 1 but is %d" +
194 boost::lexical_cast<std::string>(data.size1()));
195#endif
197 }
198};
199
200/**
201 * \brief Get tensor rank 1 (vector) form data matrix
202 */
203template <int Tensor_Dim, int S = 1, class T, class L, class A>
205getFTensor1FromMat(ublas::matrix<T, L, A> &data) {
206 static_assert(!std::is_same<T, T>::value, "not implemented");
207}
208
209/**
210 * \brief Get tensor rank 1 (vector) form data matrix (specialization)
211 */
212template <int Tensor_Dim, int S = 1>
214 return GetFTensor1FromMatImpl<Tensor_Dim, S, double, ublas::row_major,
215 DoubleAllocator>::get(data);
216}
217
218/**
219 * \brief Get tensor rank 2 (matrix) form data matrix
220 */
221template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
222inline FTensor::Tensor2<FTensor::PackPtr<T *, 1>, Tensor_Dim0, Tensor_Dim1>
223getFTensor2FromMat(ublas::matrix<T, L, A> &data) {
224 static_assert(!std::is_same<T, T>::value,
225 "Such getFTensor2FromMat specialisation is not implemented");
226}
227
228/**
229 * Template specialization for getFTensor2FromMat
230 *
231 */
232template <>
235#ifndef NDEBUG
236 if (data.size1() != 9)
237 THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix; numer "
238 "of rows should be 9 but is " +
239 boost::lexical_cast<std::string>(data.size1()));
240#endif
242 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
243 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0));
244}
245
246/**
247 * Template specialization for getFTensor2FromMat
248 */
249template <>
252#ifndef NDEBUG
253 if (data.size1() != 6)
254 THROW_MESSAGE("getFTensor2FromMat<3,3>: wrong size of data matrix, numer "
255 "of rows should be 6 but is " +
256 boost::lexical_cast<std::string>(data.size1()));
257#endif
259 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
260 &data(5, 0));
261}
262
263/**
264 * Template specialization for getFTensor2FromMat
265 */
266template <>
269#ifndef NDEBUG
270 if (data.size1() != 4)
271 THROW_MESSAGE("getFTensor2FromMat<2,2>: wrong size of data matrix, numer "
272 "of rows should be 4 but is " +
273 boost::lexical_cast<std::string>(data.size1()));
274#endif
276 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0));
277}
278
279/**
280 * Template specialization for getFTensor2FromMat
281 */
282template <>
285#ifndef NDEBUG
286 if (data.size1() != 1)
287 THROW_MESSAGE("getFTensor2FromMat<1,1>: wrong size of data matrix, numer "
288 "of rows should be 1 but is " +
289 boost::lexical_cast<std::string>(data.size1()));
290#endif
291 return FTensor::Tensor2<FTensor::PackPtr<double *, 1>, 1, 1>(&data(0, 0));
292}
293
294/**
295 * \brief Get tensor rank 2 (matrix) form data matrix (specialization)
296 */
297template <int Tensor_Dim0, int Tensor_Dim1>
298static inline FTensor::Tensor2<FTensor::PackPtr<double *, 1>, Tensor_Dim0,
299 Tensor_Dim1>
301 return getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
302 DoubleAllocator>(data);
303}
304
305template <int Tensor_Dim, int S, class T, class L, class A>
307
308template <int S, class T, class L, class A>
310 static inline auto get(ublas::matrix<T, L, A> &data) {
311#ifndef NDEBUG
312 if (data.size1() != 6)
314 "getFTensor2SymmetricFromMat<3>: wrong size of data matrix, numer "
315 "of rows should be 6 but is " +
316 boost::lexical_cast<std::string>(data.size1()));
317#endif
319 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
320 &data(5, 0));
321 }
322};
323
324template <int S, class T, class L, class A>
326 static inline auto get(ublas::matrix<T, L, A> &data) {
327#ifndef NDEBUG
328 if (data.size1() != 3)
330 "getFTensor2SymmetricFromMat<2>: wrong size of data matrix, numer "
331 "of rows should be 3 but is " +
332 boost::lexical_cast<std::string>(data.size1()));
333#endif
335 &data(0, 0), &data(1, 0), &data(2, 0));
336 }
337};
338
339/**
340 * \brief Get symmetric tensor rank 2 (matrix) form data matrix
341 */
342template <int Tensor_Dim, int S, class T, class L, class A>
343static inline auto getFTensor2SymmetricFromMat(ublas::matrix<T, L, A> &data) {
345}
346
347template <int Tensor_Dim, int S = 1>
348static inline auto getFTensor2SymmetricFromMat(MatrixDouble &data) {
350 DoubleAllocator>(data);
351}
352
353template <int Tensor_Dim01, int Tensor_Dim23, int S, class T, class L, class A>
355
356template <int S, class T, class A>
357struct GetFTensor4DdgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
358 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
359#ifndef NDEBUG
360 if (data.size1() != 1)
362 "getFTensor4DdgFromMat<1, 1>: wrong size of data matrix, number "
363 "of rows should be 1 but is " +
364 boost::lexical_cast<std::string>(data.size1()));
365#endif
366 return FTensor::Ddg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
367 }
368};
369
370template <int S, class T, class A>
371struct GetFTensor4DdgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
372 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
373#ifndef NDEBUG
374 if (data.size1() != 9) {
376 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
377 "of rows should be 9 but is " +
378 boost::lexical_cast<std::string>(data.size1()));
379 }
380#endif
382 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
383 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0)};
384 }
385};
386
387template <int S, class T, class A>
388struct GetFTensor4DdgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
389 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
390#ifndef NDEBUG
391 if (data.size1() != 36) {
392 cerr << data.size1() << endl;
394 "getFTensor4DdgFromMat<3, 3>: wrong size of data matrix, number "
395 "of rows should be 36 but is " +
396 boost::lexical_cast<std::string>(data.size1()));
397 }
398#endif
400 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
401 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
402 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
403 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
404 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
405 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
406 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
407 &data(35, 0)};
408 }
409};
410
411/**
412 * @brief Get symmetric tensor rank 4 on first two and last indices from
413 * form data matrix
414 *
415 * @tparam Tensor_Dim01 dimension of first two indicies
416 * @tparam Tensor_Dim23 dimension of second two indicies
417 * @tparam T the type of object stored
418 * @tparam L the storage organization
419 * @tparam A the type of Storage array
420 * @param data data container
421 * @return FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
422 */
423template <int Tensor_Dim01, int Tensor_Dim23, int S = 1, class T, class L,
424 class A>
425static inline FTensor::Ddg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim23>
426getFTensor4DdgFromMat(ublas::matrix<T, L, A> &data) {
427 static_assert(!std::is_same<T, T>::value,
428 "Such getFTensor4DdgFromMat specialisation is not implemented");
429}
430
431template <int Tensor_Dim01, int Tensor_Dim23, int S = 1>
432static inline auto getFTensor4DdgFromMat(MatrixDouble &data) {
433 return GetFTensor4DdgFromMatImpl<Tensor_Dim01, Tensor_Dim23, S, double,
435 DoubleAllocator>::get(data);
436}
437
438template <int Tensor_Dim01, int Tensor_Dim2, int S, class T, class L, class A>
440
441template <int S, class T, class A>
442struct GetFTensor3DgFromMatImpl<1, 1, S, T, ublas::row_major, A> {
443 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
444#ifndef NDEBUG
445 if (data.size1() != 1)
447 "getFTensor3DgFromMat<1, 1>: wrong size of data matrix, number "
448 "of rows should be 1 but is " +
449 boost::lexical_cast<std::string>(data.size1()));
450#endif
451 return FTensor::Dg<FTensor::PackPtr<double *, S>, 1, 1>{&data(0, 0)};
452 }
453};
454
455template <int S, class T, class A>
456struct GetFTensor3DgFromMatImpl<2, 2, S, T, ublas::row_major, A> {
457 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
458#ifndef NDEBUG
459 if (data.size1() != 6) {
461 "getFTensor4DdgFromMat<2, 2>: wrong size of data matrix, number "
462 "of rows should be 6 but is " +
463 boost::lexical_cast<std::string>(data.size1()));
464 }
465#endif
467 &data(0, 0), &data(1, 0), &data(2, 0),
468 &data(3, 0), &data(4, 0), &data(5, 0)};
469 }
470};
471
472template <int S, class T, class A>
473struct GetFTensor3DgFromMatImpl<3, 3, S, T, ublas::row_major, A> {
474 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
475#ifndef NDEBUG
476 if (data.size1() != 18) {
477 cerr << data.size1() << endl;
479 "getFTensor3DgFromMat<3, 3>: wrong size of data matrix, number "
480 "of rows should be 18 but is " +
481 boost::lexical_cast<std::string>(data.size1()));
482 }
483#endif
485 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
486 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
487 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
488 &data(15, 0), &data(16, 0), &data(17, 0)};
489 }
490};
491
492/**
493 * @brief Get symmetric tensor rank 3 on the first two indices from
494 * form data matrix
495 *
496 * @tparam Tensor_Dim01 dimension of first two indicies
497 * @tparam Tensor_Dim2 dimension of last index
498 * @tparam T the type of object stored
499 * @tparam L the storage organization
500 * @tparam A the type of Storage array
501 * @param data data container
502 * @return FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, TensorDim23>
503 */
504template <int Tensor_Dim01, int Tensor_Dim2, int S = 1, class T, class L,
505 class A>
506static inline FTensor::Dg<FTensor::PackPtr<T *, 1>, Tensor_Dim01, Tensor_Dim2>
507getFTensor3DgFromMat(ublas::matrix<T, L, A> &data) {
508 static_assert(!std::is_same<T, T>::value,
509 "Such getFTensor3DgFromMat specialisation is not implemented");
510}
511
512template <int Tensor_Dim01, int Tensor_Dim2, int S = 1>
513static inline auto getFTensor3DgFromMat(MatrixDouble &data) {
514 return GetFTensor3DgFromMatImpl<Tensor_Dim01, Tensor_Dim2, S, double,
516}
517
518template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
519 int S, class T, class L, class A>
521
522template <int S, class T, class A>
523struct GetFTensor4FromMatImpl<1, 1, 1, 1, S, T, ublas::row_major, A> {
524 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
525#ifndef NDEBUG
526 if (data.size1() != 1)
528 "getFTensor4FromMat<1, 1, 1, 1>: wrong size of data matrix, number "
529 "of rows should be 1 but is " +
530 boost::lexical_cast<std::string>(data.size1()));
531#endif
533 &data(0, 0)};
534 }
535};
536
537template <int S, class T, class A>
538struct GetFTensor4FromMatImpl<2, 2, 2, 2, S, T, ublas::row_major, A> {
539 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
540#ifndef NDEBUG
541 if (data.size1() != 16) {
543 "getFTensor4FromMat<2, 2, 2, 2>: wrong size of data matrix, number "
544 "of rows should be 16 but is " +
545 boost::lexical_cast<std::string>(data.size1()));
546 }
547#endif
549 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
550 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
551 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0),
552 &data(12, 0), &data(13, 0), &data(14, 0), &data(15, 0)};
553 }
554};
555
556template <int S, class T, class A>
557struct GetFTensor4FromMatImpl<3, 3, 3, 3, S, T, ublas::row_major, A> {
558 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
559#ifndef NDEBUG
560 if (data.size1() != 81) {
561 cerr << data.size1() << endl;
563 "getFTensor4FromMat<3, 3, 3, 3>: wrong size of data matrix, number "
564 "of rows should be 81 but is " +
565 boost::lexical_cast<std::string>(data.size1()));
566 }
567#endif
569 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
570 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
571 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
572 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
573 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
574 &data(25, 0), &data(26, 0), &data(27, 0), &data(28, 0), &data(29, 0),
575 &data(30, 0), &data(31, 0), &data(32, 0), &data(33, 0), &data(34, 0),
576 &data(35, 0), &data(36, 0), &data(37, 0), &data(38, 0), &data(39, 0),
577 &data(40, 0), &data(41, 0), &data(42, 0), &data(43, 0), &data(44, 0),
578 &data(45, 0), &data(46, 0), &data(47, 0), &data(48, 0), &data(49, 0),
579 &data(50, 0), &data(51, 0), &data(52, 0), &data(53, 0), &data(54, 0),
580 &data(55, 0), &data(56, 0), &data(57, 0), &data(58, 0), &data(59, 0),
581 &data(60, 0), &data(61, 0), &data(62, 0), &data(63, 0), &data(64, 0),
582 &data(65, 0), &data(66, 0), &data(67, 0), &data(68, 0), &data(69, 0),
583 &data(70, 0), &data(71, 0), &data(72, 0), &data(73, 0), &data(74, 0),
584 &data(75, 0), &data(76, 0), &data(77, 0), &data(78, 0), &data(79, 0),
585 &data(80, 0)};
586 }
587};
588
589/**
590 * @brief Get tensor rank 4 (non symmetric) form data matrix
591 *
592 * @tparam Tensor_Dim0 dimension of frirst index
593 * @tparam Tensor_Dim1 dimension of second index
594 * @tparam Tensor_Dim2 dimension of third index
595 * @tparam Tensor_Dim3 dimension of fourth index
596 * @tparam T the type of object stored
597 * @tparam L the storage organization
598 * @tparam A the type of Storage array
599 * @param data data container
600 * @return FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
601 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
602 */
603template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
604 int S = 1, class T, class L, class A>
605static inline FTensor::Tensor4<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
606 Tensor_Dim1, Tensor_Dim2, Tensor_Dim3>
607getFTensor4FromMat(ublas::matrix<T, L, A> &data) {
608 static_assert(!std::is_same<T, T>::value,
609 "Such getFTensor4FromMat specialisation is not implemented");
610}
611
612template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int Tensor_Dim3,
613 int S = 1>
614static inline auto getFTensor4FromMat(MatrixDouble &data) {
615 return GetFTensor4FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2,
616 Tensor_Dim3, S, double, ublas::row_major,
617 DoubleAllocator>::get(data);
618}
619
620template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S, class T,
621 class L, class A>
623
624template <int S, class T, class A>
625struct GetFTensor3FromMatImpl<1, 1, 1, S, T, ublas::row_major, A> {
626 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
627#ifndef NDEBUG
628 if (data.size1() != 1)
630 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
631 "of rows should be 1 but is " +
632 boost::lexical_cast<std::string>(data.size1()));
633#endif
635 &data(0, 0)};
636 }
637};
638
639template <int S, class T, class A>
640struct GetFTensor3FromMatImpl<2, 2, 2, S, T, ublas::row_major, A> {
641 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
642#ifndef NDEBUG
643 if (data.size1() != 8)
645 "getFTensor3FromMat<2, 2, 2>: wrong size of data matrix, number "
646 "of rows should be 8 but is " +
647 boost::lexical_cast<std::string>(data.size1()));
648#endif
650 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
651 &data(5, 0), &data(6, 0), &data(7, 0)
652
653 };
654 }
655};
656
657template <int S, class T, class A>
658struct GetFTensor3FromMatImpl<3, 2, 2, S, T, ublas::row_major, A> {
659 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
660#ifndef NDEBUG
661 if (data.size1() != 12)
663 "getFTensor3FromMat<3, 2, 2>: wrong size of data matrix, number "
664 "of rows should be 12 but is " +
665 boost::lexical_cast<std::string>(data.size1()));
666#endif
668 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0),
669 &data(4, 0), &data(5, 0), &data(6, 0), &data(7, 0),
670 &data(8, 0), &data(9, 0), &data(10, 0), &data(11, 0)};
671 }
672};
673
674template <int S, class T, class A>
675struct GetFTensor3FromMatImpl<3, 3, 3, S, T, ublas::row_major, A> {
676 static inline auto get(ublas::matrix<T, ublas::row_major, A> &data) {
677#ifndef NDEBUG
678 if (data.size1() != 26)
680 "getFTensor3FromMat<1, 1, 1>: wrong size of data matrix, number "
681 "of rows should be 8 but is " +
682 boost::lexical_cast<std::string>(data.size1()));
683#endif
685 &data(0, 0), &data(1, 0), &data(2, 0), &data(3, 0), &data(4, 0),
686 &data(5, 0), &data(6, 0), &data(7, 0), &data(8, 0), &data(9, 0),
687 &data(10, 0), &data(11, 0), &data(12, 0), &data(13, 0), &data(14, 0),
688 &data(15, 0), &data(16, 0), &data(17, 0), &data(18, 0), &data(19, 0),
689 &data(20, 0), &data(21, 0), &data(22, 0), &data(23, 0), &data(24, 0),
690 &data(25, 0), &data(26, 0)};
691 }
692};
693
694/**
695 * @brief Get tensor rank 3 (non symmetries) form data matrix
696 *
697 * @tparam Tensor_Dim0 dimension of frirst index
698 * @tparam Tensor_Dim1 dimension of second index
699 * @tparam Tensor_Dim2 dimension of third index
700 * @tparam S shift size
701 * @tparam T the type of object stored
702 * @tparam L the storage organization
703 * @tparam A the type of Storage array
704 * @param data data container
705 * @return FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
706 Tensor_Dim1, Tensor_Dim2>
707 */
708template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1, class T,
709 class L, class A>
710static inline FTensor::Tensor3<FTensor::PackPtr<T *, 1>, Tensor_Dim0,
711 Tensor_Dim1, Tensor_Dim2>
712getFTensor3FromMat(ublas::matrix<T, L, A> &data) {
713 static_assert(!std::is_same<T, T>::value,
714 "Such getFTensor3FromMat specialisation is not implemented");
715}
716
717template <int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2, int S = 1>
718static inline auto getFTensor3FromMat(MatrixDouble &data) {
719 return GetFTensor3FromMatImpl<Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S,
721 DoubleAllocator>::get(data);
722}
723
724/**
725 * @brief Make Tensor1 from pointer
726 *
727 * @tparam DIM
728 * @param ptr
729 * @return FTensor::Tensor2<FTensor::PackPtr<double *, 3 * DIM>, 3, DIM>
730 */
731template <int DIM>
733getFTensor1FromPtr(double *ptr) {
734 static_assert(DIM != 3 && DIM != 2,
735 "Such getFTensor1FromPtr specialization is not implemented");
736};
737
738template <>
740getFTensor1FromPtr<2>(double *ptr) {
742 &ptr[HVEC1]);
743};
744
745template <>
747getFTensor1FromPtr<3>(double *ptr) {
749 &ptr[HVEC0], &ptr[HVEC1], &ptr[HVEC2]);
750};
751
752/**
753 * @brief Make Tensor2 from pointer
754 *
755 * @tparam DIM
756 * @param ptr
757 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
758 */
759template <int DIM1, int DIM2>
761getFTensor2FromPtr(double *ptr) {
762 static_assert(DIM1 == DIM1 || DIM2 != DIM2,
763 "Such getFTensor2FromPtr is not implemented");
764};
765
766template <>
768 2> inline getFTensor2FromPtr<3, 2>(double *ptr) {
770 &ptr[HVEC0_0], &ptr[HVEC0_1],
771
772 &ptr[HVEC1_0], &ptr[HVEC1_1],
773
774 &ptr[HVEC2_0], &ptr[HVEC2_1]);
775};
776
777template <>
779 3> inline getFTensor2FromPtr<3, 3>(double *ptr) {
781 &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
782
783 &ptr[HVEC1_0], &ptr[HVEC1_1], &ptr[HVEC1_2],
784
785 &ptr[HVEC2_0], &ptr[HVEC2_1], &ptr[HVEC2_2]);
786};
787
788template <>
790 2> inline getFTensor2FromPtr<2, 2>(double *ptr) {
792 &ptr[0], &ptr[1], &ptr[2], &ptr[3]);
793};
794
795/*
796 * @brief Make Tensor3 from pointer
797 *
798 * @tparam DIM
799 * @param ptr
800 * @return FTensor::Tensor3<FTensor::PackPtr<double *, DIM1 * DIM2* DIM3>, DIM1,
801 * DIM2, DIM3>
802 */
803template <int DIM1, int DIM2, int DIM3>
805 DIM2, DIM3>
806getFTensor3FromPtr(double *ptr) {
807 static_assert(DIM1 == DIM1 || DIM2 != DIM2 || DIM3 != DIM3,
808 "Such getFTensor2FromPtr is not implemented");
809};
810
811template <>
813getFTensor3FromPtr<3, 2, 2>(double *ptr) {
815 &ptr[0], &ptr[1], &ptr[2],
816
817 &ptr[3], &ptr[4], &ptr[5],
818
819 &ptr[6], &ptr[7], &ptr[8],
820
821 &ptr[9], &ptr[10], &ptr[11]
822
823 );
824};
825
826/**
827 * @brief Make symmetric Tensor2 from pointer, taking lower triangle of matrix
828 *
829 * @tparam DIM
830 * @param ptr
831 * @return FTensor::Tensor2<FTensor::PackPtr<double *, DIM1 * DIM2>, DIM1, DIM2>
832 */
833template <int DIM>
836 static_assert(DIM,
837 "Such getFTensor2SymmetricUpperFromPtr is not implemented");
838}
839
840template <>
844 &ptr[HVEC0_0], &ptr[HVEC0_1], &ptr[HVEC0_2],
845
846 &ptr[HVEC1_0], &ptr[HVEC1_1],
847
848 &ptr[HVEC2_2]);
849};
850
851template <>
855 &ptr[0], &ptr[1], &ptr[3]);
856};
857
858/**
859 * @brief Get FTensor1 from array
860 *
861 * \todo Generalise for diffrent arrays and data types
862 *
863 * @tparam DIM
864 * @param data
865 * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
866 */
867template <int DIM, int S>
870 static_assert(DIM != DIM, "not implemented");
872}
873
874template <>
877 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data[0], &data[1]};
878}
879
880template <>
883 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>{&data[0], &data[1],
884 &data[2]};
885}
886
887template <int DIM, int S>
889getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
890 static_assert(DIM != DIM, "not implemented");
892}
893
894template <>
896getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
897 return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2>{&data(rr + 0, 0),
898 &data(rr + 1, 0)};
899}
900
901template <>
903getFTensor1FromMat(MatrixDouble &data, const size_t rr) {
905 &data(rr + 0, 0), &data(rr + 1, 0), &data(rr + 2, 0)};
906}
907
908/**
909 * @brief Get FTensor1 from array
910 *
911 * \todo Generalise for diffrent arrays and data types
912 *
913 * @tparam DIM
914 * @param data
915 * @param rr
916 * @return FTensor::Tensor1<FTensor::PackPtr<double *, DIM>, DIM>
917 */
918template <int DIM, int S>
920getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
921 static_assert(DIM != DIM, "not implemented");
923}
924
925template <>
927getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
928 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{&data(rr + 0, 0),
929 &data(rr + 1, 1)};
930}
931
932template <>
934getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr) {
936 &data(rr + 0, 0), &data(rr + 1, 1), &data(rr + 2, 2)};
937}
938
939/**
940 * @brief Get FTensor2 from array
941 *
942 * \note Generalise for other data types
943 *
944 * @tparam DIM1
945 * @tparam DIM2
946 * @tparam S
947 * @param data
948 * @return FTensor::Tensor2<FTensor::PackPtr<double *, S>, DIM1, DIM2>
949 */
950template <int DIM1, int DIM2, int S, class T, class L, class A>
952
953template <int S, class T, class L, class A>
954struct GetFTensor2FromArrayImpl<2, 2, S, T, L, A> {
955 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr) {
957 &data(rr + 0, 0), &data(rr + 0, 1), &data(rr + 1, 0), &data(rr + 1, 1)};
958 }
959};
960
961template <int S, class T, class L, class A>
962struct GetFTensor2FromArrayImpl<3, 3, S, T, L, A> {
963 inline static auto get(ublas::matrix<T, L, A> &data, const size_t rr) {
965 &data(rr + 0, 0), &data(rr + 0, 1), &data(rr + 0, 2),
966 &data(rr + 1, 0), &data(rr + 1, 1), &data(rr + 1, 2),
967 &data(rr + 2, 0), &data(rr + 2, 1), &data(rr + 2, 2)};
968 }
969};
970
971template <int DIM1, int DIM2, int S>
973getFTensor2FromArray(MatrixDouble &data, const size_t rr) {
974 return GetFTensor2FromArrayImpl<DIM1, DIM2, S, double, ublas::row_major,
975 VecAllocator<double>>::get(data, rr);
976}
977
978template <int S, typename T, typename L, typename A>
979inline auto getFTensor2FromArray2by2(ublas::matrix<T, L, A> &data,
980 const FTensor::Number<S> &,
981 const size_t rr) {
983}
984
985template <int S, typename T, typename L, typename A>
986inline auto getFTensor2FromArray3by3(ublas::matrix<T, L, A> &data,
987 const FTensor::Number<S> &,
988 const size_t rr) {
990}
991
992#ifdef WITH_ADOL_C
993
994template <int DIM1, int DIM2, int S>
995inline auto getFTensor2FromArray(MatrixADouble &data, const size_t rr) {
997 VecAllocator<adouble>>::get(data, rr);
998}
999
1000#endif
1001
1002// list of lapack wrappers
1003/**
1004 * @brief compute matrix inverse with lapack dgetri
1005 *
1006 * @param mat input square matrix / output inverse matrix
1007 * @return MoFEMErrorCode
1008 */
1011
1012 const size_t M = mat.size1();
1013 const size_t N = mat.size2();
1014
1015 if (M != N)
1016 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1017 "The input matrix for inverse computation is not square %d != %d",
1018 M, N);
1019
1020 int *ipv = new int[N];
1021 int lwork = N * N;
1022 double *work = new double[lwork];
1023 int info;
1024 info = lapack_dgetrf(N, N, &*mat.data().begin(), N, ipv);
1025 if (info != 0)
1026 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1027 "lapack error info = %d", info);
1028 info = lapack_dgetri(N, &*mat.data().begin(), N, ipv, work, lwork);
1029 if (info != 0)
1030 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1031 "lapack error info = %d", info);
1032
1033 delete[] ipv;
1034 delete[] work;
1035
1037}
1038/**
1039 * @brief solve linear system with lapack dgesv
1040 *
1041 * @param mat input lhs square matrix / output L and U from the factorization
1042 * @param f input rhs vector / output solution vector
1043 * @return MoFEMErrorCode
1044 */
1047
1048 const size_t M = mat.size1();
1049 const size_t N = mat.size2();
1050
1051 if (M == 0 || M != N)
1052 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1053 "The input matrix for inverse computation is not square %d != %d",
1054 M, N);
1055 if (f.size() != M)
1056 f.resize(M, false);
1057
1058 const int nrhs = 1;
1059 int info;
1060 int *ipiv = new int[M];
1061 info = lapack_dgesv(M, nrhs, &*mat.data().begin(), M, ipiv,
1062 &*f.data().begin(), nrhs);
1063
1064 if (info != 0) {
1065 SETERRQ1(PETSC_COMM_SELF, 1, "error lapack solve dgesv info = %d", info);
1066 }
1067
1068 delete[] ipiv;
1070}
1071
1072/**
1073 * @brief Solve linear system of equations using Lapack
1074 *
1075 * @param mat
1076 * @param f
1077 * @return MoFEMErrorCode
1078 */
1080 VectorDouble &f) {
1082 // copy matrix since on output lapack returns factorisation
1083 auto mat_copy = mat;
1084 CHKERR solveLinearSystem(mat_copy, f);
1086}
1087
1088/**
1089 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1090 *
1091 * @param mat input symmetric matrix
1092 * @param eig output eigen values sorted
1093 * @param eigen_vec output matrix of row eigen vectors
1094 * @return MoFEMErrorCode
1095 */
1097 VectorDouble &eig,
1098 MatrixDouble &eigen_vec) {
1100
1101 const size_t M = mat.size1();
1102 const size_t N = mat.size2();
1103
1104 if (M == 0 || M != N)
1105 SETERRQ2(
1106 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1107 "The input matrix for eigen value computation is not square %d != %d",
1108 M, N);
1109 if (eig.size() != M)
1110 eig.resize(M, false);
1111
1112 eigen_vec = mat;
1113 const int n = M;
1114 const int lda = M;
1115 const int size = (M + 2) * M;
1116 int lwork = size;
1117 double *work = new double[size];
1118
1119 if (lapack_dsyev('V', 'U', n, &*eigen_vec.data().begin(), lda,
1120 &*eig.data().begin(), work, lwork) > 0)
1121 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1122 "The algorithm failed to compute eigenvalues.");
1123
1124 delete[] work;
1126}
1127/**
1128 * @brief compute eigenvalues of a symmetric matrix using lapack dsyev
1129 *
1130 * @tparam DIM
1131 * @param eigen_vec input / output DIM x DIM matrix of row eigen vectors
1132 * @param eig output eigen values sorted
1133 * @return MoFEMErrorCode
1134 */
1135template <int DIM>
1136inline MoFEMErrorCode
1140
1141 const int n = DIM;
1142 const int lda = DIM;
1143 const int lwork = (DIM + 2) * DIM;
1144 std::array<double, (DIM + 2) * DIM> work;
1145
1146 if (lapack_dsyev('V', 'U', n, &eigen_vec(0, 0), lda, &eig(0), work.data(),
1147 lwork) > 0)
1148 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
1149 "The algorithm failed to compute eigenvalues.");
1151}
1152/**
1153 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1154 *
1155 * @tparam DIM
1156 * @param mat input tensor pointer of size DIM x DIM
1157 * @param eig output eigen values sorted
1158 * @param eigen_vec output matrix of row eigen vectors
1159 * @return MoFEMErrorCode
1160 */
1161template <int DIM>
1167 for (int ii = 0; ii != DIM; ii++)
1168 for (int jj = 0; jj != DIM; jj++)
1169 eigen_vec(ii, jj) = mat(ii, jj);
1170
1171 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1172
1174}
1175
1176/**
1177 * @brief compute eigenvalues of a symmetric tensor using lapack dsyev
1178 *
1179 * @tparam DIM
1180 * @param mat input tensor of size DIM x DIM
1181 * @param eig output eigen values sorted
1182 * @param eigen_vec output matrix of row eigen vectors
1183 * @return MoFEMErrorCode
1184 */
1185template <int DIM>
1186inline MoFEMErrorCode
1191 for (int ii = 0; ii != DIM; ii++)
1192 for (int jj = 0; jj != DIM; jj++)
1193 eigen_vec(ii, jj) = mat(ii, jj);
1194
1195 CHKERR computeEigenValuesSymmetric<DIM>(eigen_vec, eig);
1196
1198}
1199
1200/**
1201 * @brief Calculate the determinant of a 3x3 matrix or a tensor of rank 2
1202 *
1203 * @tparam T
1204 * @param t
1205 * @return double
1206 */
1207template <typename T> static inline auto determinantTensor3by3(T &t) {
1208 return t(0, 0) * t(1, 1) * t(2, 2) + t(1, 0) * t(2, 1) * t(0, 2) +
1209 t(2, 0) * t(0, 1) * t(1, 2) - t(0, 0) * t(2, 1) * t(1, 2) -
1210 t(2, 0) * t(1, 1) * t(0, 2) - t(1, 0) * t(0, 1) * t(2, 2);
1211}
1212
1213/**
1214 * \brief Calculate inverse of tensor rank 2 at integration points
1215
1216 */
1217template <int Tensor_Dim, class T, class L, class A>
1218inline MoFEMErrorCode invertTensor3by3(ublas::matrix<T, L, A> &jac_data,
1219 ublas::vector<T, A> &det_data,
1220 ublas::matrix<T, L, A> &inv_jac_data) {
1222 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1223 "Specialization for this template not yet implemented");
1225}
1226
1227template <>
1228inline MoFEMErrorCode
1229invertTensor3by3<3, double, ublas::row_major, DoubleAllocator>(
1230 MatrixDouble &jac_data, VectorDouble &det_data, MatrixDouble &inv_jac_data);
1231
1232/**
1233 * \brief Calculate determinant 3 by 3
1234
1235 */
1236template <class T1, class T2>
1239 det = determinantTensor3by3(t);
1241}
1242
1243/**
1244 * \brief Calculate determinant 2 by 2
1245
1246 */
1247template <class T1, class T2>
1250 det = t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0);
1252}
1253
1254/**
1255 * \brief Calculate matrix inverse 3 by 3
1256
1257 */
1258template <class T1, class T2, class T3>
1259inline MoFEMErrorCode invertTensor3by3(T1 &t, T2 &det, T3 &inv_t) {
1261 const auto inv_det = 1. / det;
1262 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1263 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1264 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1265 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) * inv_det;
1266 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1267 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1268 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) * inv_det;
1269 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) * inv_det;
1270 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1272}
1273
1274/**
1275 * \brief Calculate matrix inverse 2 by 2
1276
1277 */
1278template <class T1, class T2, class T3>
1279inline MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t) {
1281 const auto inv_det = 1. / det;
1282 inv_t(0, 0) = t(1, 1) * inv_det;
1283 inv_t(0, 1) = -t(0, 1) * inv_det;
1284 inv_t(1, 0) = -t(1, 0) * inv_det;
1285 inv_t(1, 1) = t(0, 0) * inv_det;
1287}
1288
1289#ifdef WITH_ADOL_C
1290
1291/**
1292 * \brief Calculate matrix inverse, specialization for adouble tensor
1293
1294 */
1295template <>
1296inline MoFEMErrorCode invertTensor3by3<FTensor::Tensor2<adouble, 3, 3>, adouble,
1301 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1302 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1303 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1304 inv_t(1, 0) = (t(1, 2) * t(2, 0) - t(1, 0) * t(2, 2)) / det;
1305 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1306 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1307 inv_t(2, 0) = (t(1, 0) * t(2, 1) - t(1, 1) * t(2, 0)) / det;
1308 inv_t(2, 1) = (t(0, 1) * t(2, 0) - t(0, 0) * t(2, 1)) / det;
1309 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1311}
1312
1313#endif
1314
1315/**
1316 * \brief Calculate matrix inverse, specialization for symmetric tensor
1317
1318 */
1319template <>
1320inline MoFEMErrorCode
1321invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1326 const auto inv_det = 1. / det;
1327 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1328 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1329 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1330 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1331 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1332 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1334}
1335
1336#ifdef WITH_ADOL_C
1337
1338/**
1339 * \brief Calculate matrix inverse, specialization for adouble symmetric tensor
1340
1341 */
1342template <>
1343inline MoFEMErrorCode
1344invertTensor3by3<FTensor::Tensor2_symmetric<adouble, 3>, adouble,
1349 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) / det;
1350 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) / det;
1351 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) / det;
1352 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) / det;
1353 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) / det;
1354 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) / det;
1356}
1357
1358#endif
1359
1360/**
1361 * \brief Calculate matrix inverse, specialization for symmetric (pointer)
1362 tensor
1363
1364 */
1365template <>
1366inline MoFEMErrorCode
1367invertTensor3by3<FTensor::Tensor2_symmetric<double, 3>, double,
1372 const auto inv_det = 1. / det;
1373 inv_t(0, 0) = (t(1, 1) * t(2, 2) - t(1, 2) * t(2, 1)) * inv_det;
1374 inv_t(0, 1) = (t(0, 2) * t(2, 1) - t(0, 1) * t(2, 2)) * inv_det;
1375 inv_t(0, 2) = (t(0, 1) * t(1, 2) - t(0, 2) * t(1, 1)) * inv_det;
1376 inv_t(1, 1) = (t(0, 0) * t(2, 2) - t(0, 2) * t(2, 0)) * inv_det;
1377 inv_t(1, 2) = (t(0, 2) * t(1, 0) - t(0, 0) * t(1, 2)) * inv_det;
1378 inv_t(2, 2) = (t(0, 0) * t(1, 1) - t(0, 1) * t(1, 0)) * inv_det;
1380}
1381
1382/**
1383 * @brief Extract entity handle form multi-index container
1384 *
1385 */
1387 template <typename Iterator>
1388 static inline EntityHandle extract(const Iterator &it) {
1389 return (*it)->getEnt();
1390 }
1391};
1392
1393/**
1394 * @brief Insert ordered mofem multi-index into range
1395 *
1396 * \code
1397 * auto hi_rit = refEntsPtr->upper_bound(start);
1398 * auto hi_rit = refEntsPtr->upper_bound(end);
1399 * Range to_erase;
1400 * insertOrdered(to_erase, RefEntExtractor(), rit, hi_rit);
1401 * \endcode
1402 *
1403 * @tparam Iterator
1404 * @param r
1405 * @param begin_iter
1406 * @param end_iter
1407 * @return moab::Range::iterator
1408 */
1409template <typename Extractor, typename Iterator>
1410moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter,
1411 Iterator end_iter) {
1412 moab::Range::iterator hint = r.begin();
1413 while (begin_iter != end_iter) {
1414 size_t j = 0;
1415 auto bi = Extractor::extract(begin_iter);
1416 Iterator pj = begin_iter;
1417 while (pj != end_iter && (bi + j) == Extractor::extract(pj)) {
1418 ++pj;
1419 ++j;
1420 }
1421 hint = r.insert(hint, bi, bi + (j - 1));
1422 begin_iter = pj;
1423 }
1424 return hint;
1425};
1426
1427/**
1428 * @brief Do nothing, used to rebuild database
1429 *
1430 */
1433 template <typename T> inline void operator()(T &e) {}
1434};
1435
1436/**
1437 * @brief Template used to reconstruct multi-index
1438 *
1439 * @tparam MI multi-index
1440 * @tparam Modifier
1441 * @param mi
1442 * @param mo
1443 * @return MoFEMErrorCode
1444 */
1445template <typename MI, typename MO = Modify_change_nothing>
1447 MO &&mo = Modify_change_nothing()) {
1449 for (auto it = mi.begin(); it != mi.end(); ++it) {
1450 if (!const_cast<MI &>(mi).modify(it, mo))
1451 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1452 "Houston we have a problem");
1453 }
1455}
1456
1459 rval = moab.create_meshset(MESHSET_SET, meshset);
1461 }
1463 operator EntityHandle() const { return meshset; }
1464 auto get_ptr() { return &meshset; }
1465
1466private:
1468 rval = moab.delete_entities(&meshset, 1);
1470 }
1471 EntityHandle meshset;
1473};
1474
1475/**
1476 * @brief Create smart pointer to temprary meshset
1477 *
1478 */
1480 return boost::make_shared<TempMeshset>(moab);
1481};
1482
1483/**
1484 * @brief get type from entity handle
1485 *
1486 */
1487inline auto type_from_handle(const EntityHandle h) {
1488 return static_cast<EntityType>(h >> MB_ID_WIDTH);
1489};
1490
1491/**
1492 * @brief get entity dimension form handle
1493 *
1494 */
1495inline auto dimension_from_handle(const EntityHandle h) {
1496 return moab::CN::Dimension(type_from_handle(h));
1497};
1498
1499/**
1500 * @brief get field bit id from bit number
1501 *
1502 */
1503inline auto field_bit_from_bit_number(const int bit_number) {
1504 return BitFieldId().set(bit_number - 1);
1505};
1506
1507/**
1508 * @brief Insert ranges
1509 *
1510 * @tparam I
1511 * @param f
1512 * @param s
1513 * @param tester
1514 * @param inserter
1515 * @return auto
1516 */
1517template <typename I>
1518auto rangeInserter(const I f, const I s, boost::function<bool(I it)> tester,
1519 boost::function<MoFEMErrorCode(I f, I s)> inserter) {
1521
1522 auto first = f;
1523 while (first != s)
1524 if (tester(first)) {
1525
1526 auto second = first;
1527 ++second;
1528
1529 while (second != s) {
1530 if (tester(second))
1531 ++second;
1532 else
1533 break;
1534 }
1535
1536 CHKERR inserter(first, second);
1537
1538 first = second;
1539 if (first != s)
1540 ++first;
1541
1542 } else {
1543 ++first;
1544 }
1545
1547}
1548
1549// template <typename T, typename I>
1550// std::vector<std::pair<T, T>> getPairRange(I s, I e, boost::function<T(I)>){
1551
1552// };
1553
1554} // namespace MoFEM
1555
1556#endif //__TEMPLATES_HPP__
static Index< 'M', 3 > M
static Index< 'n', 3 > n
static Index< 'I', 3 > I
T data[Tensor_Dim]
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MB_ID_WIDTH
Definition: definitions.h:240
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ HVEC0
Definition: definitions.h:199
@ HVEC1
Definition: definitions.h:199
@ HVEC2
Definition: definitions.h:199
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_INVALID_DATA
Definition: definitions.h:49
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
@ HVEC1_1
Definition: definitions.h:209
@ HVEC0_1
Definition: definitions.h:208
@ HVEC1_0
Definition: definitions.h:206
@ HVEC2_1
Definition: definitions.h:210
@ HVEC1_2
Definition: definitions.h:212
@ HVEC2_2
Definition: definitions.h:213
@ HVEC2_0
Definition: definitions.h:207
@ HVEC0_2
Definition: definitions.h:211
@ HVEC0_0
Definition: definitions.h:205
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:188
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:169
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:273
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:197
FTensor::Index< 'j', 3 > j
const double T
@ row_major
Definition: Layout.hpp:13
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::vector< T, std::allocator< T > > VecAllocator
Definition: Types.hpp:70
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
ublas::vector< T, ublas::shallow_array_adaptor< T > > VectorShallowArrayAdaptor
Definition: Types.hpp:125
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasMatrix< adouble > MatrixADouble
Definition: Types.hpp:91
VecAllocator< double > DoubleAllocator
Definition: Types.hpp:73
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition: Types.hpp:132
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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.
Definition: Templates.hpp:712
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:45
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.
Definition: Templates.hpp:507
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1487
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1248
FTensor::Tensor2< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:223
FTensor::Tensor3< FTensor::PackPtr< double *, 12 >, 3, 2, 2 > getFTensor3FromPtr< 3, 2, 2 >(double *ptr)
Definition: Templates.hpp:813
FTensor::Tensor2< FTensor::PackPtr< double *, 4 >, 2, 2 > getFTensor2FromPtr< 2, 2 >(double *ptr)
Definition: Templates.hpp:790
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
Definition: Templates.hpp:1446
moab::Range::iterator insertOrdered(Range &r, Extractor, Iterator begin_iter, Iterator end_iter)
Insert ordered mofem multi-index into range.
Definition: Templates.hpp:1410
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1279
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:986
MoFEMErrorCode computeMatrixInverse(MatrixDouble &mat)
compute matrix inverse with lapack dgetri
Definition: Templates.hpp:1009
auto dimension_from_handle(const EntityHandle h)
get entity dimension form handle
Definition: Templates.hpp:1495
MoFEMErrorCode solveLinearSystem(MatrixDouble &mat, VectorDouble &f)
solve linear system with lapack dgesv
Definition: Templates.hpp:1045
FTensor::Tensor1< FTensor::PackPtr< double *, DIM >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition: Templates.hpp:733
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1218
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.
Definition: Templates.hpp:426
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 4 >, 2 > getFTensor2SymmetricLowerFromPtr< 2 >(double *ptr)
Definition: Templates.hpp:853
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
Definition: Templates.hpp:205
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:71
FTensor::Tensor2< FTensor::PackPtr< double *, DIM1 *DIM2 >, DIM1, DIM2 > getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
Definition: Templates.hpp:761
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
Definition: Templates.hpp:343
auto field_bit_from_bit_number(const int bit_number)
get field bit id from bit number
Definition: Templates.hpp:1503
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
std::string toString(X x)
Definition: Templates.hpp:120
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.
Definition: Templates.hpp:607
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2FromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:779
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArrayDiag(MatrixDouble &data, const size_t rr)
Get FTensor1 from array.
Definition: Templates.hpp:920
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr)
Definition: Templates.hpp:973
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temprary meshset.
Definition: Templates.hpp:1479
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
Definition: Templates.hpp:747
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition: Templates.hpp:1096
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, DIM *DIM >, DIM > getFTensor2SymmetricLowerFromPtr(double *ptr)
Make symmetric Tensor2 from pointer, taking lower triangle of matrix.
Definition: Templates.hpp:835
FTensor::Tensor1< FTensor::PackPtr< double *, 2 >, 2 > getFTensor1FromPtr< 2 >(double *ptr)
Definition: Templates.hpp:740
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
Definition: Templates.hpp:869
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2FromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:768
FTensor::Tensor3< FTensor::PackPtr< double *, DIM1 *DIM2 *DIM3 >, DIM1, DIM2, DIM3 > getFTensor3FromPtr(double *ptr)
Definition: Templates.hpp:806
boost::shared_ptr< std::vector< T > > ShardVec
Definition: Templates.hpp:24
auto getFTensor2FromArray2by2(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:979
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
Definition: Templates.hpp:1518
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *, 9 >, 3 > getFTensor2SymmetricLowerFromPtr< 3 >(double *ptr)
Definition: Templates.hpp:842
const double r
rate factor
double A
double h
convective heat coefficient
constexpr double t
plate stiffness
Definition: plate.cpp:76
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:109
static auto get(ublas::vector< T, A > &data)
Definition: Templates.hpp:127
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:188
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:173
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:158
static auto get(ublas::matrix< T, L, A > &data, const size_t rr)
Definition: Templates.hpp:955
static auto get(ublas::matrix< T, L, A > &data, const size_t rr)
Definition: Templates.hpp:963
Get FTensor2 from array.
Definition: Templates.hpp:951
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:326
static auto get(ublas::matrix< T, L, A > &data)
Definition: Templates.hpp:310
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:443
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:457
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:474
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:626
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:641
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:659
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:676
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:358
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:372
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:389
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:524
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:539
static auto get(ublas::matrix< T, ublas::row_major, A > &data)
Definition: Templates.hpp:558
unsigned int operator()(const id_type &value) const
Definition: Templates.hpp:115
KeyExtractor1 key1
Definition: Templates.hpp:98
result_type operator()(Arg &arg) const
Definition: Templates.hpp:93
KeyFromKey(const KeyExtractor1 &key1_=KeyExtractor1(), const KeyExtractor2 &key2_=KeyExtractor2())
Definition: Templates.hpp:89
KeyExtractor2 key2
Definition: Templates.hpp:99
KeyExtractor1::result_type result_type
Definition: Templates.hpp:87
bool operator()(const id_type &valueA, const id_type &valueB) const
Definition: Templates.hpp:103
Do nothing, used to rebuild database.
Definition: Templates.hpp:1431
Extract entity handle form multi-index container.
Definition: Templates.hpp:1386
static EntityHandle extract(const Iterator &it)
Definition: Templates.hpp:1388
virtual ~TempMeshset()
Definition: Templates.hpp:1462
EntityHandle meshset
Definition: Templates.hpp:1471
moab::Interface & moab
Definition: Templates.hpp:1472
TempMeshset(moab::Interface &moab)
Definition: Templates.hpp:1458