v0.13.1
UserDataOperators.hpp
Go to the documentation of this file.
1/** \file UserDataOperators.hpp
2 * \brief User data Operators
3
4*/
5
6/* This file is part of MoFEM.
7 * MoFEM is free software: you can redistribute it and/or modify it under
8 * the terms of the GNU Lesser General Public License as published by the
9 * Free Software Foundation, either version 3 of the License, or (at your
10 * option) any later version.
11 *
12 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 * License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19
20#ifndef __USER_DATA_OPERATORS_HPP__
21#define __USER_DATA_OPERATORS_HPP__
22
23namespace MoFEM {
24
25/** \name Get values at Gauss pts */
26
27/**@{*/
28
29/** \name Scalar values */
30
31/**@{*/
32
33/** \brief Scalar field values at integration points
34 *
35 */
36template <class T, class A>
39
41 const std::string field_name,
42 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
43 const EntityType zero_type = MBVERTEX)
45 dataPtr(data_ptr), zeroType(zero_type) {
46 if (!dataPtr)
47 THROW_MESSAGE("Pointer is not set");
48 }
49
51 const std::string field_name,
52 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
53 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
55 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
56 if (!dataPtr)
57 THROW_MESSAGE("Pointer is not set");
58 }
59
60 /**
61 * \brief calculate values of scalar field at integration points
62 * @param side side entity number
63 * @param type side entity type
64 * @param data entity data
65 * @return error code
66 */
69
70protected:
71 boost::shared_ptr<ublas::vector<T, A>> dataPtr;
72 const EntityHandle zeroType;
75};
76
77/**
78 * \brief Specialization of member function
79 *
80 */
81template <class T, class A>
85 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented for T = %s",
86 typeid(T).name() // boost::core::demangle(typeid(T).name()).c_str()
87 );
89}
90
91/**
92 * \brief Get value at integration points for scalar field
93 *
94 * \ingroup mofem_forces_and_sources_user_data_operators
95 */
97 : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
98
101
102 /**
103 * \brief calculate values of scalar field at integration points
104 * @param side side entity number
105 * @param type side entity type
106 * @param data entity data
107 * @return error code
108 */
112 VectorDouble &vec = *dataPtr;
113 const size_t nb_gauss_pts = getGaussPts().size2();
114 if (type == zeroType) {
115 vec.resize(nb_gauss_pts, false);
116 vec.clear();
117 }
118
119 const size_t nb_dofs = data.getFieldData().size();
120
121 if (nb_dofs) {
122
123 if (dataVec.use_count()) {
124 dotVector.resize(nb_dofs, false);
125 const double *array;
126 CHKERR VecGetArrayRead(dataVec, &array);
127 const auto &local_indices = data.getLocalIndices();
128 for (int i = 0; i != local_indices.size(); ++i)
129 if (local_indices[i] != -1)
130 dotVector[i] = array[local_indices[i]];
131 else
132 dotVector[i] = 0;
133 CHKERR VecRestoreArrayRead(dataVec, &array);
134 data.getFieldData().swap(dotVector);
135 }
136
137 const size_t nb_base_functions = data.getN().size2();
138 auto base_function = data.getFTensor0N();
139 auto values_at_gauss_pts = getFTensor0FromVec(vec);
140 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
141 auto field_data = data.getFTensor0FieldData();
142 size_t bb = 0;
143 for (; bb != nb_dofs; ++bb) {
144 values_at_gauss_pts += field_data * base_function;
145 ++field_data;
146 ++base_function;
147 }
148 // It is possible to have more base functions than dofs
149 for (; bb != nb_base_functions; ++bb)
150 ++base_function;
151 ++values_at_gauss_pts;
152 }
153
154 if (dataVec.use_count()) {
155 data.getFieldData().swap(dotVector);
156 }
157 }
158
160 }
161};
162
163/**
164 * @brief Get rate of scalar field at integration points
165 *
166 * \ingroup mofem_forces_and_sources_user_data_operators
167 */
168
169template <PetscData::DataContext CTX>
172
174 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
175 const EntityType zero_at_type = MBVERTEX)
178 dataPtr(data_ptr), zeroAtType(zero_at_type) {
179 if (!dataPtr)
180 THROW_MESSAGE("Pointer is not set");
181 }
182
186
187 const size_t nb_gauss_pts = getGaussPts().size2();
188
189 VectorDouble &vec = *dataPtr;
190 if (type == zeroAtType) {
191 vec.resize(nb_gauss_pts, false);
192 vec.clear();
193 }
194
195 auto &local_indices = data.getLocalIndices();
196 const size_t nb_dofs = local_indices.size();
197 if (nb_dofs) {
198
199 const double *array;
200
201 auto get_array = [&](const auto ctx, auto vec) {
203 if ((getFEMethod()->data_ctx & ctx).none())
204 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
205 CHKERR VecGetArrayRead(vec, &array);
207 };
208
209 auto restore_array = [&](auto vec) {
210 return VecRestoreArrayRead(vec, &array);
211 };
212
213 switch (CTX) {
215 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
216 break;
218 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
219 break;
221 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
222 break;
223 default:
224 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
225 "That case is not implemented");
226 }
227
228 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
229 for (int i = 0; i != local_indices.size(); ++i)
230 if (local_indices[i] != -1)
231 dot_dofs_vector[i] = array[local_indices[i]];
232 else
233 dot_dofs_vector[i] = 0;
234
235 switch (CTX) {
237 CHKERR restore_array(getFEMethod()->ts_u);
238 break;
240 CHKERR restore_array(getFEMethod()->ts_u_t);
241 break;
243 CHKERR restore_array(getFEMethod()->ts_u_tt);
244 break;
245 default:
246 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
247 "That case is not implemented");
248 }
249
250 const size_t nb_base_functions = data.getN().size2();
251 auto base_function = data.getFTensor0N();
252 auto values_at_gauss_pts = getFTensor0FromVec(vec);
253
254 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
255 size_t bb = 0;
256 for (; bb != nb_dofs; ++bb) {
257 values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
258 ++base_function;
259 }
260 // Number of dofs can be smaller than number of Tensor_Dim x base
261 // functions
262 for (; bb != nb_base_functions; ++bb)
263 ++base_function;
264 ++values_at_gauss_pts;
265 }
266 }
268 }
269
270private:
271 boost::shared_ptr<VectorDouble> dataPtr;
272 const EntityHandle zeroAtType;
273};
274
279
280/**
281 * \depreacted Name inconstent with other operators
282 *
283 */
285
286/**@}*/
287
288/** \name Vector field values at integration points */
289
290/**@{*/
291
292/** \brief Calculate field values for tenor field rank 1, i.e. vector field
293 *
294 */
295template <int Tensor_Dim, class T, class L, class A>
298
300 const std::string field_name,
301 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
302 const EntityType zero_type = MBVERTEX)
305 dataPtr(data_ptr), zeroType(zero_type) {
306 if (!dataPtr)
307 THROW_MESSAGE("Pointer is not set");
308 }
309
310 /**
311 * \brief calculate values of vector field at integration points
312 * @param side side entity number
313 * @param type side entity type
314 * @param data entity data
315 * @return error code
316 */
319
320protected:
321 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
322 const EntityHandle zeroType;
323};
324
325template <int Tensor_Dim, class T, class L, class A>
328 int side, EntityType type, EntitiesFieldData::EntData &data) {
330 SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
331 "Not implemented for T = %s and dim = %d",
332 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
333 Tensor_Dim);
335}
336
337/** \brief Calculate field values (template specialization) for tensor field
338 * rank 1, i.e. vector field
339 *
340 */
341template <int Tensor_Dim>
345
347 boost::shared_ptr<MatrixDouble> data_ptr,
348 const EntityType zero_type = MBVERTEX)
351 dataPtr(data_ptr), zeroType(zero_type) {
352 if (!dataPtr)
353 THROW_MESSAGE("Pointer is not set");
354 }
355
358
359protected:
360 boost::shared_ptr<MatrixDouble> dataPtr;
361 const EntityHandle zeroType;
362};
363
364/**
365 * \brief Member function specialization calculating values for tenor field rank
366 *
367 */
368template <int Tensor_Dim>
370 Tensor_Dim, double, ublas::row_major,
371 DoubleAllocator>::doWork(int side, EntityType type,
374
375 const size_t nb_gauss_pts = getGaussPts().size2();
376 auto &mat = *dataPtr;
377 if (type == zeroType) {
378 mat.resize(Tensor_Dim, nb_gauss_pts, false);
379 mat.clear();
380 }
381
382 const size_t nb_dofs = data.getFieldData().size();
383 if (nb_dofs) {
384
385 if (nb_gauss_pts) {
386 const size_t nb_base_functions = data.getN().size2();
387 auto base_function = data.getFTensor0N();
388 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
390 const size_t size = nb_dofs / Tensor_Dim;
391 if (nb_dofs % Tensor_Dim) {
392 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
393 "Data inconsistency");
394 }
395 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
396 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
397 size_t bb = 0;
398 for (; bb != size; ++bb) {
399 values_at_gauss_pts(I) += field_data(I) * base_function;
400 ++field_data;
401 ++base_function;
402 }
403 // Number of dofs can be smaller than number of Tensor_Dim x base
404 // functions
405 for (; bb != nb_base_functions; ++bb)
406 ++base_function;
407 ++values_at_gauss_pts;
408 }
409 }
410 }
412}
413
414/** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
415 * field
416 *
417 * \ingroup mofem_forces_and_sources_user_data_operators
418 */
419template <int Tensor_Dim>
422 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
423
425 Tensor_Dim, double, ublas::row_major,
427};
428
429/**@}*/
430
431/** \name Vector field values at integration points */
432
433/**@{*/
434
435/** \brief Calculate field values (template specialization) for tensor field
436 * rank 1, i.e. vector field
437 *
438 */
439template <int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
442
444 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
445 const EntityType zero_type = MBVERTEX)
448 dataPtr(data_ptr), zeroType(zero_type) {
449 if (!dataPtr)
450 THROW_MESSAGE("Pointer is not set");
451 }
452
456
457 // When we move to C++17 add if constexpr()
458 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
459 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
460 "%s coordiante not implemented",
461 CoordinateTypesNames[COORDINATE_SYSTEM]);
462
463 const size_t nb_gauss_pts = getGaussPts().size2();
464 auto &vec = *dataPtr;
465 if (type == zeroType) {
466 vec.resize(nb_gauss_pts, false);
467 vec.clear();
468 }
469
470 const size_t nb_dofs = data.getFieldData().size();
471 if (nb_dofs) {
472
473 if (nb_gauss_pts) {
474 const size_t nb_base_functions = data.getN().size2();
475 auto values_at_gauss_pts = getFTensor0FromVec(vec);
477 const size_t size = nb_dofs / Tensor_Dim;
478#ifndef NDEBUG
479 if (nb_dofs % Tensor_Dim) {
480 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
481 "Number of dofs should multiple of dimensions");
482 }
483#endif
484
485 // When we move to C++17 add if constexpr()
486 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
487 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
488 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
489 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
490 size_t bb = 0;
491 for (; bb != size; ++bb) {
492 values_at_gauss_pts += field_data(I) * diff_base_function(I);
493 ++field_data;
494 ++diff_base_function;
495 }
496 // Number of dofs can be smaller than number of Tensor_Dim x base
497 // functions
498 for (; bb != nb_base_functions; ++bb)
499 ++diff_base_function;
500 ++values_at_gauss_pts;
501 }
502 }
503
504 // When we move to C++17 add if constexpr()
505 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
506 auto t_coords = getFTensor1CoordsAtGaussPts();
507 auto values_at_gauss_pts = getFTensor0FromVec(vec);
508 auto base_function = data.getFTensor0N();
509 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
510 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
511 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
512 size_t bb = 0;
513 for (; bb != size; ++bb) {
514 values_at_gauss_pts += field_data(I) * diff_base_function(I);
515 values_at_gauss_pts +=
516 field_data(0) * base_function / t_coords(0);
517 ++field_data;
518 ++base_function;
519 ++diff_base_function;
520 }
521 // Number of dofs can be smaller than number of Tensor_Dim x base
522 // functions
523 for (; bb != nb_base_functions; ++bb) {
524 ++base_function;
525 ++diff_base_function;
526 }
527 ++values_at_gauss_pts;
528 ++t_coords;
529 }
530 }
531 }
532 }
534 }
535
536protected:
537 boost::shared_ptr<VectorDouble> dataPtr;
538 const EntityHandle zeroType;
539};
540
541/** \brief Approximate field valuse for given petsc vector
542 *
543 * \note Look at PetscData to see what vectors could be extarcted with that user
544 * data operator.
545 *
546 * \ingroup mofem_forces_and_sources_user_data_operators
547 */
548template <int Tensor_Dim, PetscData::DataContext CTX>
551
553 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
554 const EntityType zero_at_type = MBVERTEX)
557 dataPtr(data_ptr), zeroAtType(zero_at_type) {
558 if (!dataPtr)
559 THROW_MESSAGE("Pointer is not set");
560 }
561
565 auto &local_indices = data.getLocalIndices();
566 const size_t nb_dofs = local_indices.size();
567 const size_t nb_gauss_pts = getGaussPts().size2();
568
569 MatrixDouble &mat = *dataPtr;
570 if (type == zeroAtType) {
571 mat.resize(Tensor_Dim, nb_gauss_pts, false);
572 mat.clear();
573 }
574 if (!nb_dofs)
576
577 const double *array;
578
579 auto get_array = [&](const auto ctx, auto vec) {
581 if ((getFEMethod()->data_ctx & ctx).none())
582 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
583 CHKERR VecGetArrayRead(vec, &array);
585 };
586
587 auto restore_array = [&](auto vec) {
588 return VecRestoreArrayRead(vec, &array);
589 };
590
591 switch (CTX) {
593 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
594 break;
596 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
597 break;
599 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
600 break;
601 default:
602 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
603 "That case is not implemented");
604 }
605
606 dotVector.resize(local_indices.size());
607 for (int i = 0; i != local_indices.size(); ++i)
608 if (local_indices[i] != -1)
609 dotVector[i] = array[local_indices[i]];
610 else
611 dotVector[i] = 0;
612
613 switch (CTX) {
615 CHKERR restore_array(getFEMethod()->ts_u);
616 break;
618 CHKERR restore_array(getFEMethod()->ts_u_t);
619 break;
621 CHKERR restore_array(getFEMethod()->ts_u_tt);
622 break;
623 default:
624 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
625 "That case is not implemented");
626 }
627
628 const size_t nb_base_functions = data.getN().size2();
629 auto base_function = data.getFTensor0N();
630 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
631
633 const size_t size = nb_dofs / Tensor_Dim;
634 if (nb_dofs % Tensor_Dim) {
635 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
636 }
637 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
638 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(dotVector);
639 size_t bb = 0;
640 for (; bb != size; ++bb) {
641 values_at_gauss_pts(I) += field_data(I) * base_function;
642 ++field_data;
643 ++base_function;
644 }
645 // Number of dofs can be smaller than number of Tensor_Dim x base
646 // functions
647 for (; bb != nb_base_functions; ++bb)
648 ++base_function;
649 ++values_at_gauss_pts;
650 }
652 }
653
654protected:
655 boost::shared_ptr<MatrixDouble> dataPtr;
656 const EntityHandle zeroAtType;
658};
659
660/** \brief Get time derivatives of values at integration pts for tensor filed
661 * rank 1, i.e. vector field
662 *
663 * \ingroup mofem_forces_and_sources_user_data_operators
664 */
665template <int Tensor_Dim>
669
670/** \brief Get second time derivatives of values at integration pts for tensor
671 * filed rank 1, i.e. vector field
672 *
673 * \ingroup mofem_forces_and_sources_user_data_operators
674 */
675template <int Tensor_Dim>
679
680/**@}*/
681
682/** \name Tensor field values at integration points */
683
684/**@{*/
685
686/** \brief Calculate field values for tenor field rank 2.
687 *
688 */
689template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
692
694 const std::string field_name,
695 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
696 const EntityType zero_type = MBVERTEX)
699 dataPtr(data_ptr), zeroType(zero_type) {
700 if (!dataPtr)
701 THROW_MESSAGE("Pointer is not set");
702 }
703
706
707protected:
708 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
709 const EntityHandle zeroType;
710};
711
712template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
717 SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
718 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
719 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
720 Tensor_Dim0, Tensor_Dim1);
722}
723
724template <int Tensor_Dim0, int Tensor_Dim1>
725struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
728
730 const std::string field_name,
731 boost::shared_ptr<
732 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
733 data_ptr,
734 const EntityType zero_type = MBVERTEX)
737 dataPtr(data_ptr), zeroType(zero_type) {
738 if (!dataPtr)
739 THROW_MESSAGE("Pointer is not set");
740 }
741
743 const std::string field_name,
744 boost::shared_ptr<
745 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
746 data_ptr,
747 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
750 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
751 if (!dataPtr)
752 THROW_MESSAGE("Pointer is not set");
753 }
754
757
758protected:
759 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
761 const EntityHandle zeroType;
764};
765
766template <int Tensor_Dim0, int Tensor_Dim1>
768 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
769 DoubleAllocator>::doWork(int side, EntityType type,
772
773 MatrixDouble &mat = *dataPtr;
774 const size_t nb_gauss_pts = data.getN().size1();
775 if (type == zeroType) {
776 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
777 mat.clear();
778 }
779
780 const size_t nb_dofs = data.getFieldData().size();
781
782 if (dataVec.use_count()) {
783 dotVector.resize(nb_dofs, false);
784 const double *array;
785 CHKERR VecGetArrayRead(dataVec, &array);
786 const auto &local_indices = data.getLocalIndices();
787 for (int i = 0; i != local_indices.size(); ++i)
788 if (local_indices[i] != -1)
789 dotVector[i] = array[local_indices[i]];
790 else
791 dotVector[i] = 0;
792 CHKERR VecRestoreArrayRead(dataVec, &array);
793 data.getFieldData().swap(dotVector);
794 }
795
796 if (nb_dofs) {
797 const size_t nb_base_functions = data.getN().size2();
798 auto base_function = data.getFTensor0N();
799 auto values_at_gauss_pts =
800 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
803 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
804 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
805 auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
806 size_t bb = 0;
807 for (; bb != size; ++bb) {
808 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
809 ++field_data;
810 ++base_function;
811 }
812 for (; bb != nb_base_functions; ++bb)
813 ++base_function;
814 ++values_at_gauss_pts;
815 }
816
817 if (dataVec.use_count()) {
818 data.getFieldData().swap(dotVector);
819 }
820 }
822}
823
824/** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
825 * field
826 *
827 * \ingroup mofem_forces_and_sources_user_data_operators
828 */
829template <int Tensor_Dim0, int Tensor_Dim1>
832 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
833
835 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
837};
838
839/** \brief Get time direvarive values at integration pts for tensor filed rank
840 * 2, i.e. matrix field
841 *
842 * \ingroup mofem_forces_and_sources_user_data_operators
843 */
844template <int Tensor_Dim0, int Tensor_Dim1>
847
849 boost::shared_ptr<MatrixDouble> data_ptr,
850 const EntityType zero_at_type = MBVERTEX)
853 dataPtr(data_ptr), zeroAtType(zero_at_type) {
854 if (!dataPtr)
855 THROW_MESSAGE("Pointer is not set");
856 }
857
861
862 const size_t nb_gauss_pts = getGaussPts().size2();
863 MatrixDouble &mat = *dataPtr;
864 if (type == zeroAtType) {
865 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
866 mat.clear();
867 }
868 const auto &local_indices = data.getLocalIndices();
869 const size_t nb_dofs = local_indices.size();
870 if (nb_dofs) {
871 dotVector.resize(nb_dofs, false);
872 const double *array;
873 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
874 for (size_t i = 0; i != local_indices.size(); ++i)
875 if (local_indices[i] != -1)
876 dotVector[i] = array[local_indices[i]];
877 else
878 dotVector[i] = 0;
879 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
880
881 const size_t nb_base_functions = data.getN().size2();
882
883 auto base_function = data.getFTensor0N();
884 auto values_at_gauss_pts =
885 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
888 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
889 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
890 auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
891 size_t bb = 0;
892 for (; bb != size; ++bb) {
893 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
894 ++field_data;
895 ++base_function;
896 }
897 for (; bb != nb_base_functions; ++bb)
898 ++base_function;
899 ++values_at_gauss_pts;
900 }
901 }
903 }
904
905protected:
906 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
907 EntityType zeroAtType; ///< Zero values at Gauss point at this type
908 VectorDouble dotVector; ///< Keeps temoorary values of time directives
909
910 template <int Dim0, int Dim1> auto getFTensorDotData() {
911 static_assert(Dim0 || !Dim0 || Dim1 || !Dim1, "not implemented");
912 }
913};
914
915template <>
916template <>
919 &dotVector[0], &dotVector[1], &dotVector[2],
920
921 &dotVector[3], &dotVector[4], &dotVector[5],
922
923 &dotVector[6], &dotVector[7], &dotVector[8]);
924}
925
926/**
927 * @brief Calculate symmetric tensor field values at integration pts.
928 *
929 * @tparam Tensor_Dim
930
931 * \ingroup mofem_forces_and_sources_user_data_operators
932 */
933template <int Tensor_Dim>
936
938 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
939 const EntityType zero_type = MBEDGE, const int zero_side = 0)
942 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
943 if (!dataPtr)
944 THROW_MESSAGE("Pointer is not set");
945 }
946
948 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
949 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBEDGE,
950 const int zero_side = 0)
953 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side),
954 dataVec(data_vec) {
955 if (!dataPtr)
956 THROW_MESSAGE("Pointer is not set");
957 }
958
962 MatrixDouble &mat = *dataPtr;
963 const int nb_gauss_pts = getGaussPts().size2();
964 if (type == this->zeroType && side == zeroSide) {
965 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
966 mat.clear();
967 }
968 const int nb_dofs = data.getFieldData().size();
969 if (!nb_dofs)
971
972 if (dataVec.use_count()) {
973 dotVector.resize(nb_dofs, false);
974 const double *array;
975 CHKERR VecGetArrayRead(dataVec, &array);
976 const auto &local_indices = data.getLocalIndices();
977 for (int i = 0; i != local_indices.size(); ++i)
978 if (local_indices[i] != -1)
979 dotVector[i] = array[local_indices[i]];
980 else
981 dotVector[i] = 0;
982 CHKERR VecRestoreArrayRead(dataVec, &array);
983 data.getFieldData().swap(dotVector);
984 }
985
986 const int nb_base_functions = data.getN().size2();
987 auto base_function = data.getFTensor0N();
988 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
991 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
992 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
993 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
994 int bb = 0;
995 for (; bb != size; ++bb) {
996 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
997 ++field_data;
998 ++base_function;
999 }
1000 for (; bb != nb_base_functions; ++bb)
1001 ++base_function;
1002 ++values_at_gauss_pts;
1003 }
1004
1005 if (dataVec.use_count()) {
1006 data.getFieldData().swap(dotVector);
1007 }
1008
1010 }
1011
1012protected:
1013 boost::shared_ptr<MatrixDouble> dataPtr;
1014 const EntityHandle zeroType;
1015 const int zeroSide;
1018};
1019
1020/**
1021 * @brief Calculate symmetric tensor field rates ant integratio pts.
1022 *
1023 * @tparam Tensor_Dim
1024 *
1025 * \ingroup mofem_forces_and_sources_user_data_operators
1026 */
1027template <int Tensor_Dim>
1030
1032 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1033 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1036 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1037 if (!dataPtr)
1038 THROW_MESSAGE("Pointer is not set");
1039 }
1040
1044 const int nb_gauss_pts = getGaussPts().size2();
1045 MatrixDouble &mat = *dataPtr;
1046 if (type == zeroType && side == zeroSide) {
1047 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
1048 mat.clear();
1049 }
1050 auto &local_indices = data.getLocalIndices();
1051 const int nb_dofs = local_indices.size();
1052 if (!nb_dofs)
1054
1055 dotVector.resize(nb_dofs, false);
1056 const double *array;
1057 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1058 for (int i = 0; i != local_indices.size(); ++i)
1059 if (local_indices[i] != -1)
1060 dotVector[i] = array[local_indices[i]];
1061 else
1062 dotVector[i] = 0;
1063 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1064
1065 const int nb_base_functions = data.getN().size2();
1066
1067 auto base_function = data.getFTensor0N();
1068 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1071 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1072 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1073 auto field_data = getFTensorDotData<Tensor_Dim>();
1074 int bb = 0;
1075 for (; bb != size; ++bb) {
1076 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1077 ++field_data;
1078 ++base_function;
1079 }
1080 for (; bb != nb_base_functions; ++bb)
1081 ++base_function;
1082 ++values_at_gauss_pts;
1083 }
1084
1086 }
1087
1088protected:
1089 boost::shared_ptr<MatrixDouble> dataPtr;
1090 const EntityHandle zeroType;
1091 const int zeroSide;
1093
1094 template <int Dim> inline auto getFTensorDotData() {
1095 static_assert(Dim || !Dim, "not implemented");
1096 }
1097};
1098
1099template <>
1100template <>
1101inline auto
1104 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1105 &dotVector[5]);
1106}
1107
1108template <>
1109template <>
1110inline auto
1113 &dotVector[0], &dotVector[1], &dotVector[2]);
1114}
1115
1116/**@}*/
1117
1118/** \name Gradients and Hessian of scalar fields at integration points */
1119
1120/**@{*/
1121
1122/**
1123 * \brief Evaluate field gradient values for scalar field, i.e. gradient is
1124 * tensor rank 1 (vector)
1125 *
1126 */
1127template <int Tensor_Dim, class T, class L, class A>
1129 : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
1130
1133};
1134
1135/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1136 * tensor rank 1 (vector), specialization
1137 *
1138 */
1139template <int Tensor_Dim>
1141 ublas::row_major, DoubleAllocator>
1143 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1144
1146 Tensor_Dim, double, ublas::row_major,
1148
1149 /**
1150 * \brief calculate gradient values of scalar field at integration points
1151 * @param side side entity number
1152 * @param type side entity type
1153 * @param data entity data
1154 * @return error code
1155 */
1158};
1159
1160/**
1161 * \brief Member function specialization calculating scalar field gradients for
1162 * tenor field rank 1
1163 *
1164 */
1165template <int Tensor_Dim>
1167 Tensor_Dim, double, ublas::row_major,
1168 DoubleAllocator>::doWork(int side, EntityType type,
1171
1172 const size_t nb_gauss_pts = this->getGaussPts().size2();
1173 auto &mat = *this->dataPtr;
1174 if (type == this->zeroType) {
1175 mat.resize(Tensor_Dim, nb_gauss_pts, false);
1176 mat.clear();
1177 }
1178
1179 const int nb_dofs = data.getFieldData().size();
1180 if (nb_dofs) {
1181
1182 if (nb_gauss_pts) {
1183 const int nb_base_functions = data.getN().size2();
1184 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
1185 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1186
1187#ifndef NDEBUG
1188 if (nb_dofs > nb_base_functions)
1189 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1190 "Number of base functions inconsistent with number of DOFs "
1191 "(%d > %d)",
1192 nb_dofs, nb_base_functions);
1193
1194 if (data.getDiffN().size2() != nb_base_functions * Tensor_Dim)
1195 SETERRQ2(
1196 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1197 "Number of base functions inconsistent with number of derivatives "
1198 "(%d != %d)",
1199 data.getDiffN().size2(), nb_base_functions);
1200
1201 if (data.getDiffN().size1() != nb_gauss_pts)
1202 SETERRQ2(
1203 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1204 "Number of base functions inconsistent with number of integration "
1205 "pts (%d != %d)",
1206 data.getDiffN().size2(), nb_gauss_pts);
1207
1208#endif
1209
1211 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1212 auto field_data = data.getFTensor0FieldData();
1213 int bb = 0;
1214 for (; bb != nb_dofs; ++bb) {
1215 gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
1216 ++field_data;
1217 ++diff_base_function;
1218 }
1219 // Number of dofs can be smaller than number of base functions
1220 for (; bb < nb_base_functions; ++bb)
1221 ++diff_base_function;
1222 ++gradients_at_gauss_pts;
1223 }
1224 }
1225 }
1226
1228}
1229
1230/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1231 * vector field
1232 *
1233 * \ingroup mofem_forces_and_sources_user_data_operators
1234 */
1235template <int Tensor_Dim>
1238 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1240 Tensor_Dim, double, ublas::row_major,
1241 DoubleAllocator>::OpCalculateScalarFieldGradient_General;
1242};
1243
1244/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1245 * tensor rank 1 (vector), specialization
1246 *
1247 */
1248template <int Tensor_Dim>
1251 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1252
1254 Tensor_Dim, double, ublas::row_major,
1256
1257 /**
1258 * \brief calculate gradient values of scalar field at integration points
1259 * @param side side entity number
1260 * @param type side entity type
1261 * @param data entity data
1262 * @return error code
1263 */
1266};
1267
1268template <int Tensor_Dim>
1270 int side, EntityType type, EntitiesFieldData::EntData &data) {
1272
1273 const size_t nb_gauss_pts = this->getGaussPts().size2();
1274 auto &mat = *this->dataPtr;
1275 if (type == this->zeroType) {
1276 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts, false);
1277 mat.clear();
1278 }
1279
1280 const int nb_dofs = data.getFieldData().size();
1281 if (nb_dofs) {
1282
1283 if (nb_gauss_pts) {
1284 const int nb_base_functions = data.getN().size2();
1285
1286 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1287#ifndef NDEBUG
1288 if (hessian_base.size1() != nb_gauss_pts) {
1289 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1290 "Wrong number of integration pts (%d != %d)",
1291 hessian_base.size1(), nb_gauss_pts);
1292 }
1293 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1294 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1295 "Wrong number of base functions (%d != %d)",
1296 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1297 nb_base_functions);
1298 }
1299 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1300 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1301 "Wrong number of base functions (%d < %d)",
1302 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1303 }
1304#endif
1305
1306 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1307 &*hessian_base.data().begin());
1308
1309 auto hessian_at_gauss_pts =
1310 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1311
1314 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1315 auto field_data = data.getFTensor0FieldData();
1316 int bb = 0;
1317 for (; bb != nb_dofs; ++bb) {
1318 hessian_at_gauss_pts(I, J) +=
1319 field_data * t_diff2_base_function(I, J);
1320 ++field_data;
1321 ++t_diff2_base_function;
1322 }
1323 // Number of dofs can be smaller than number of base functions
1324 for (; bb < nb_base_functions; ++bb) {
1325 ++t_diff2_base_function;
1326 }
1327
1328 ++hessian_at_gauss_pts;
1329 }
1330 }
1331 }
1332
1334}
1335
1336/**}*/
1337
1338/** \name Gradients and hessian of tensor fields at integration points */
1339
1340/**@{*/
1341
1342/**
1343 * \brief Evaluate field gradient values for vector field, i.e. gradient is
1344 * tensor rank 2
1345 *
1346 */
1347template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1349 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1350 L, A> {
1351
1353 Tensor_Dim0, Tensor_Dim1, T, L, A>::OpCalculateTensor2FieldValues_General;
1354};
1355
1356template <int Tensor_Dim0, int Tensor_Dim1>
1357struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1358 ublas::row_major, DoubleAllocator>
1360 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1361
1363 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1365
1366 /**
1367 * \brief calculate values of vector field at integration points
1368 * @param side side entity number
1369 * @param type side entity type
1370 * @param data entity data
1371 * @return error code
1372 */
1375};
1376
1377/**
1378 * \brief Member function specialization calculating vector field gradients for
1379 * tenor field rank 2
1380 *
1381 */
1382template <int Tensor_Dim0, int Tensor_Dim1>
1384 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1385 DoubleAllocator>::doWork(int side, EntityType type,
1388 if (!this->dataPtr)
1389 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1390 "Data pointer not allocated");
1391
1392 const size_t nb_gauss_pts = this->getGaussPts().size2();
1393 auto &mat = *this->dataPtr;
1394 if (type == this->zeroType) {
1395 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1396 mat.clear();
1397 }
1398
1399 if (nb_gauss_pts) {
1400 const size_t nb_dofs = data.getFieldData().size();
1401
1402 if (nb_dofs) {
1403
1404 if (this->dataVec.use_count()) {
1405 this->dotVector.resize(nb_dofs, false);
1406 const double *array;
1407 CHKERR VecGetArrayRead(this->dataVec, &array);
1408 const auto &local_indices = data.getLocalIndices();
1409 for (int i = 0; i != local_indices.size(); ++i)
1410 if (local_indices[i] != -1)
1411 this->dotVector[i] = array[local_indices[i]];
1412 else
1413 this->dotVector[i] = 0;
1414 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1415 data.getFieldData().swap(this->dotVector);
1416 }
1417
1418 const int nb_base_functions = data.getN().size2();
1419 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1420 auto gradients_at_gauss_pts =
1421 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1424 int size = nb_dofs / Tensor_Dim0;
1425 if (nb_dofs % Tensor_Dim0) {
1426 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1427 "Data inconsistency");
1428 }
1429 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1430 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1431 int bb = 0;
1432 for (; bb < size; ++bb) {
1433 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1434 ++field_data;
1435 ++diff_base_function;
1436 }
1437 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1438 // functions
1439 for (; bb != nb_base_functions; ++bb)
1440 ++diff_base_function;
1441 ++gradients_at_gauss_pts;
1442 }
1443
1444 if (this->dataVec.use_count()) {
1445 data.getFieldData().swap(this->dotVector);
1446 }
1447 }
1448 }
1450}
1451
1452/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1453 * vector field
1454 *
1455 * \ingroup mofem_forces_and_sources_user_data_operators
1456 */
1457template <int Tensor_Dim0, int Tensor_Dim1>
1460 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1461
1463 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1464 DoubleAllocator>::OpCalculateVectorFieldGradient_General;
1465};
1466
1467/** \brief Get field gradients time derivative at integration pts for scalar
1468 * filed rank 0, i.e. vector field
1469 *
1470 * \ingroup mofem_forces_and_sources_user_data_operators
1471 */
1472template <int Tensor_Dim0, int Tensor_Dim1>
1475
1477 boost::shared_ptr<MatrixDouble> data_ptr,
1478 const EntityType zero_at_type = MBVERTEX)
1481 dataPtr(data_ptr), zeroAtType(zero_at_type) {
1482 if (!dataPtr)
1483 THROW_MESSAGE("Pointer is not set");
1484 }
1485
1489
1490 const auto &local_indices = data.getLocalIndices();
1491 const int nb_dofs = local_indices.size();
1492 const int nb_gauss_pts = this->getGaussPts().size2();
1493
1494 auto &mat = *this->dataPtr;
1495 if (type == this->zeroAtType) {
1496 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1497 mat.clear();
1498 }
1499 if (!nb_dofs)
1501
1502 dotVector.resize(nb_dofs, false);
1503 const double *array;
1504 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1505 for (int i = 0; i != local_indices.size(); ++i)
1506 if (local_indices[i] != -1)
1507 dotVector[i] = array[local_indices[i]];
1508 else
1509 dotVector[i] = 0;
1510 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1511
1512 const int nb_base_functions = data.getN().size2();
1513 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1514 auto gradients_at_gauss_pts =
1515 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1518 int size = nb_dofs / Tensor_Dim0;
1519 if (nb_dofs % Tensor_Dim0) {
1520 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1521 }
1522
1523 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1524 auto field_data = getFTensorDotData<Tensor_Dim0>();
1525 int bb = 0;
1526 for (; bb < size; ++bb) {
1527 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1528 ++field_data;
1529 ++diff_base_function;
1530 }
1531 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1532 // functions
1533 for (; bb != nb_base_functions; ++bb)
1534 ++diff_base_function;
1535 ++gradients_at_gauss_pts;
1536 }
1538 }
1539
1540private:
1541 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1542 EntityType zeroAtType; ///< Zero values at Gauss point at this type
1543 VectorDouble dotVector; ///< Keeps temoorary values of time directives
1544
1545 template <int Dim> inline auto getFTensorDotData() {
1546 static_assert(Dim || !Dim, "not implemented");
1547 }
1548};
1549
1550template <>
1551template <>
1554 &dotVector[0], &dotVector[1], &dotVector[2]);
1555}
1556
1557template <>
1558template <>
1560 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(&dotVector[0],
1561 &dotVector[1]);
1562}
1563
1564/**
1565 * \brief Evaluate field gradient values for symmetric 2nd order tensor field,
1566 * i.e. gradient is tensor rank 3
1567 *
1568 */
1569template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1571 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1572 L, A> {
1573
1575 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1576 const EntityType zero_type = MBVERTEX)
1577 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
1578 A>(field_name, data_ptr,
1579 zero_type) {}
1580};
1581
1582template <int Tensor_Dim0, int Tensor_Dim1>
1584 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator>
1586 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1587
1589 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1590 const EntityType zero_type = MBVERTEX)
1591 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1592 ublas::row_major,
1594 field_name, data_ptr, zero_type) {}
1595
1596 /**
1597 * \brief calculate values of vector field at integration points
1598 * @param side side entity number
1599 * @param type side entity type
1600 * @param data entity data
1601 * @return error code
1602 */
1605};
1606
1607/**
1608 * \brief Member function specialization calculating tensor field gradients for
1609 * symmetric tensor field rank 2
1610 *
1611 */
1612template <int Tensor_Dim0, int Tensor_Dim1>
1613MoFEMErrorCode OpCalculateTensor2SymmetricFieldGradient_General<
1614 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1615 DoubleAllocator>::doWork(int side, EntityType type,
1618 if (!this->dataPtr)
1619 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1620 "Data pointer not allocated");
1621
1622 const size_t nb_gauss_pts = this->getGaussPts().size2();
1623 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1624 auto &mat = *this->dataPtr;
1625 if (type == this->zeroType) {
1626 mat.resize(msize * Tensor_Dim1, nb_gauss_pts, false);
1627 mat.clear();
1628 }
1629
1630 if (nb_gauss_pts) {
1631 const size_t nb_dofs = data.getFieldData().size();
1632
1633 if (nb_dofs) {
1634
1635 const int nb_base_functions = data.getN().size2();
1636 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1637 auto gradients_at_gauss_pts =
1638 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1642 int size = nb_dofs / msize;
1643 if (nb_dofs % msize) {
1644 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1645 "Data inconsistency");
1646 }
1647 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1648 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim0>();
1649 int bb = 0;
1650 for (; bb < size; ++bb) {
1651 gradients_at_gauss_pts(I, J, K) +=
1652 field_data(I, J) * diff_base_function(K);
1653 ++field_data;
1654 ++diff_base_function;
1655 }
1656 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1657 // functions
1658 for (; bb != nb_base_functions; ++bb)
1659 ++diff_base_function;
1660 ++gradients_at_gauss_pts;
1661 }
1662 }
1663 }
1665}
1666
1667/** \brief Get field gradients at integration pts for symmetric tensorial field
1668 * rank 2
1669 *
1670 * \ingroup mofem_forces_and_sources_user_data_operators
1671 */
1672template <int Tensor_Dim0, int Tensor_Dim1>
1675 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1676
1678 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1679 const EntityType zero_type = MBVERTEX)
1681 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1682 DoubleAllocator>(field_name, data_ptr, zero_type) {}
1683};
1684
1685template <int Tensor_Dim0, int Tensor_Dim1>
1688 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1689
1691 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1693
1694 /**
1695 * \brief calculate values of vector field at integration points
1696 * @param side side entity number
1697 * @param type side entity type
1698 * @param data entity data
1699 * @return error code
1700 */
1703};
1704
1705/**
1706 * \brief Member function specialization calculating vector field gradients for
1707 * tenor field rank 2
1708 *
1709 */
1710template <int Tensor_Dim0, int Tensor_Dim1>
1712 int side, EntityType type, EntitiesFieldData::EntData &data) {
1714 if (!this->dataPtr)
1715 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1716 "Data pointer not allocated");
1717
1718 const size_t nb_gauss_pts = this->getGaussPts().size2();
1719 auto &mat = *this->dataPtr;
1720 if (type == this->zeroType) {
1721 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts, false);
1722 mat.clear();
1723 }
1724
1725 if (nb_gauss_pts) {
1726 const size_t nb_dofs = data.getFieldData().size();
1727
1728 if (nb_dofs) {
1729
1730 if (this->dataVec.use_count()) {
1731 this->dotVector.resize(nb_dofs, false);
1732 const double *array;
1733 CHKERR VecGetArrayRead(this->dataVec, &array);
1734 const auto &local_indices = data.getLocalIndices();
1735 for (int i = 0; i != local_indices.size(); ++i)
1736 if (local_indices[i] != -1)
1737 this->dotVector[i] = array[local_indices[i]];
1738 else
1739 this->dotVector[i] = 0;
1740 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1741 data.getFieldData().swap(this->dotVector);
1742 }
1743
1744 const int nb_base_functions = data.getN().size2();
1745
1746 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1747#ifndef NDEBUG
1748 if (hessian_base.size1() != nb_gauss_pts) {
1749 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1750 "Wrong number of integration pts (%d != %d)",
1751 hessian_base.size1(), nb_gauss_pts);
1752 }
1753 if (hessian_base.size2() !=
1754 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1755 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1756 "Wrong number of base functions (%d != %d)",
1757 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1758 nb_base_functions);
1759 }
1760 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1761 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1762 "Wrong number of base functions (%d < %d)",
1763 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1764 }
1765#endif
1766
1767 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1768 &*hessian_base.data().begin());
1769
1770 auto t_hessian_at_gauss_pts =
1771 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1772
1776
1777 int size = nb_dofs / Tensor_Dim0;
1778#ifndef NDEBUG
1779 if (nb_dofs % Tensor_Dim0) {
1780 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1781 "Data inconsistency");
1782 }
1783#endif
1784
1785 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1786 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1787 int bb = 0;
1788 for (; bb < size; ++bb) {
1789 t_hessian_at_gauss_pts(I, J, K) +=
1790 field_data(I) * t_diff2_base_function(J, K);
1791 ++field_data;
1792 ++t_diff2_base_function;
1793 }
1794 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1795 // functions
1796 for (; bb != nb_base_functions; ++bb)
1797 ++t_diff2_base_function;
1798 ++t_hessian_at_gauss_pts;
1799 }
1800
1801 if (this->dataVec.use_count()) {
1802 data.getFieldData().swap(this->dotVector);
1803 }
1804 }
1805 }
1807}
1808
1809/**@}*/
1810
1811/** \name Transform tensors and vectors */
1812
1813/**@{*/
1814
1815/**
1816 * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1817 * \f$
1818 *
1819 * @tparam DIM
1820 *
1821 * \ingroup mofem_forces_and_sources_user_data_operators
1822 */
1823template <int DIM_01, int DIM_23, int S = 0>
1826
1829
1831 boost::shared_ptr<MatrixDouble> in_mat,
1832 boost::shared_ptr<MatrixDouble> out_mat,
1833 boost::shared_ptr<MatrixDouble> d_mat)
1834 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1835 // Only is run for vertices
1836 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1837 if (!inMat)
1838 THROW_MESSAGE("Pointer for in mat is null");
1839 if (!outMat)
1840 THROW_MESSAGE("Pointer for out mat is null");
1841 if (!dMat)
1842 THROW_MESSAGE("Pointer for tensor mat is null");
1843 }
1844
1847 const size_t nb_gauss_pts = getGaussPts().size2();
1848 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
1849 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
1850 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1851 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
1852 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1853 t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1854 ++t_in;
1855 ++t_out;
1856 }
1858 }
1859
1860private:
1865
1866 boost::shared_ptr<MatrixDouble> inMat;
1867 boost::shared_ptr<MatrixDouble> outMat;
1868 boost::shared_ptr<MatrixDouble> dMat;
1869};
1870
1871template <int DIM>
1873
1876
1878 boost::shared_ptr<MatrixDouble> in_mat,
1879 boost::shared_ptr<MatrixDouble> out_mat)
1880 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1881 // Only is run for vertices
1882 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1883 if (!inMat)
1884 THROW_MESSAGE("Pointer not set for in matrix");
1885 if (!outMat)
1886 THROW_MESSAGE("Pointer not set for in matrix");
1887 }
1888
1891 const size_t nb_gauss_pts = getGaussPts().size2();
1892 auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
1893 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
1894 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
1895 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1896 t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
1897 ++t_in;
1898 ++t_out;
1899 }
1901 }
1902
1903private:
1906 boost::shared_ptr<MatrixDouble> inMat;
1907 boost::shared_ptr<MatrixDouble> outMat;
1908};
1909
1911
1914
1915 OpScaleMatrix(const std::string field_name, const double scale,
1916 boost::shared_ptr<MatrixDouble> in_mat,
1917 boost::shared_ptr<MatrixDouble> out_mat)
1918 : UserOp(field_name, OPROW), scale(scale), inMat(in_mat),
1919 outMat(out_mat) {
1920 // Only is run for vertices
1921 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1922 if (!inMat)
1923 THROW_MESSAGE("Pointer not set for in matrix");
1924 if (!outMat)
1925 THROW_MESSAGE("Pointer not set for in matrix");
1926 }
1927
1930 outMat->resize(inMat->size1(), inMat->size2(), false);
1931 noalias(*outMat) = scale * (*inMat);
1933 }
1934
1935private:
1936 const double scale;
1937 boost::shared_ptr<MatrixDouble> inMat;
1938 boost::shared_ptr<MatrixDouble> outMat;
1939};
1940
1941/**@}*/
1942
1943/** \name H-div/H-curls (Vectorial bases) values at integration points */
1944
1945/**@{*/
1946
1947/** \brief Get vector field for H-div approximation
1948 * \ingroup mofem_forces_and_sources_user_data_operators
1949 */
1950template <int Tensor_Dim0, class T, class L, class A>
1953
1955 const std::string field_name,
1956 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
1957 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1960 dataPtr(data_ptr), zeroType(zero_type), zeroSide(0) {
1961 if (!dataPtr)
1962 THROW_MESSAGE("Pointer is not set");
1963 }
1964
1965 /**
1966 * \brief calculate values of vector field at integration points
1967 * @param side side entity number
1968 * @param type side entity type
1969 * @param data entity data
1970 * @return error code
1971 */
1974
1975private:
1976 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
1977 const EntityHandle zeroType;
1978 const int zeroSide;
1979};
1980
1981template <int Tensor_Dim, class T, class L, class A>
1983 int side, EntityType type, EntitiesFieldData::EntData &data) {
1985 SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1986 "Not implemented for T = %s and dim = %d",
1987 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
1988 Tensor_Dim);
1990}
1991
1992/** \brief Get vector field for H-div approximation
1993 * \ingroup mofem_forces_and_sources_user_data_operators
1994 */
1995template <int Tensor_Dim>
1999
2001 boost::shared_ptr<MatrixDouble> data_ptr,
2002 const EntityType zero_type = MBEDGE,
2003 const int zero_side = 0)
2006 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2007 if (!dataPtr)
2008 THROW_MESSAGE("Pointer is not set");
2009 }
2010
2011 /**
2012 * \brief Calculate values of vector field at integration points
2013 * @param side side entity number
2014 * @param type side entity type
2015 * @param data entity data
2016 * @return error code
2017 */
2020
2021private:
2022 boost::shared_ptr<MatrixDouble> dataPtr;
2023 const EntityHandle zeroType;
2024 const int zeroSide;
2025};
2026
2027template <int Tensor_Dim>
2029 Tensor_Dim, double, ublas::row_major,
2030 DoubleAllocator>::doWork(int side, EntityType type,
2033 const size_t nb_integration_points = this->getGaussPts().size2();
2034 if (type == zeroType && side == zeroSide) {
2035 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2036 dataPtr->clear();
2037 }
2038 const size_t nb_dofs = data.getFieldData().size();
2039 if (!nb_dofs)
2041 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim;
2043 auto t_n_hdiv = data.getFTensor1N<Tensor_Dim>();
2044 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2045 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2046 auto t_dof = data.getFTensor0FieldData();
2047 int bb = 0;
2048 for (; bb != nb_dofs; ++bb) {
2049 t_data(i) += t_n_hdiv(i) * t_dof;
2050 ++t_n_hdiv;
2051 ++t_dof;
2052 }
2053 for (; bb != nb_base_functions; ++bb)
2054 ++t_n_hdiv;
2055 ++t_data;
2056 }
2058}
2059
2060/** \brief Get vector field for H-div approximation
2061 *
2062 * \ingroup mofem_forces_and_sources_user_data_operators
2063 */
2064template <int Tensor_Dim>
2067 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
2069 Tensor_Dim, double, ublas::row_major,
2071};
2072
2073/** \brief Get vector field for H-div approximation
2074 * \ingroup mofem_forces_and_sources_user_data_operators
2075 */
2076template <int Tensor_Dim>
2079
2081 boost::shared_ptr<MatrixDouble> data_ptr,
2082 const EntityType zero_type = MBEDGE,
2083 const int zero_side = 0)
2086 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2087 if (!dataPtr)
2088 THROW_MESSAGE("Pointer is not set");
2089 }
2090
2091 /**
2092 * \brief Calculate values of vector field at integration points
2093 * @param side side entity number
2094 * @param type side entity type
2095 * @param data entity data
2096 * @return error code
2097 */
2100
2101private:
2102 boost::shared_ptr<MatrixDouble> dataPtr;
2103 const EntityHandle zeroType;
2104 const int zeroSide;
2105};
2106
2107template <int Tensor_Dim>
2109 int side, EntityType type, EntitiesFieldData::EntData &data) {
2111 const size_t nb_integration_points = this->getGaussPts().size2();
2112 if (type == zeroType && side == zeroSide) {
2113 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2114 dataPtr->clear();
2115 }
2116
2117 auto &local_indices = data.getIndices();
2118 const size_t nb_dofs = local_indices.size();
2119 if (nb_dofs) {
2120
2121 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2122 const double *array;
2123 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2124 for (size_t i = 0; i != nb_dofs; ++i)
2125 if (local_indices[i] != -1)
2126 dot_dofs_vector[i] = array[local_indices[i]];
2127 else
2128 dot_dofs_vector[i] = 0;
2129 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2130
2131 const size_t nb_base_functions = data.getN().size2() / 3;
2133 auto t_n_hdiv = data.getFTensor1N<Tensor_Dim>();
2134 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2135 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2136 int bb = 0;
2137 for (; bb != nb_dofs; ++bb) {
2138 t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
2139 ++t_n_hdiv;
2140 }
2141 for (; bb != nb_base_functions; ++bb)
2142 ++t_n_hdiv;
2143 ++t_data;
2144 }
2145 }
2147}
2148
2149/**
2150 * @brief Calculate divergence of vector field
2151 * @ingroup mofem_forces_and_sources_user_data_operators
2152 *
2153 * @tparam BASE_DIM
2154 * @tparam SPACE_DIM
2155 */
2156template <int BASE_DIM, int SPACE_DIM>
2159
2161 boost::shared_ptr<VectorDouble> data_ptr,
2162 const EntityType zero_type = MBEDGE,
2163 const int zero_side = 0)
2166 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2167 if (!dataPtr)
2168 THROW_MESSAGE("Pointer is not set");
2169 }
2170
2174 const size_t nb_integration_points = getGaussPts().size2();
2175 if (type == zeroType && side == zeroSide) {
2176 dataPtr->resize(nb_integration_points, false);
2177 dataPtr->clear();
2178 }
2179 const size_t nb_dofs = data.getFieldData().size();
2180 if (!nb_dofs)
2182 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2185 auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2186 auto t_data = getFTensor0FromVec(*dataPtr);
2187 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2188 auto t_dof = data.getFTensor0FieldData();
2189 int bb = 0;
2190 for (; bb != nb_dofs; ++bb) {
2191 t_data += t_dof * t_n_diff_hdiv(j, j);
2192 ++t_n_diff_hdiv;
2193 ++t_dof;
2194 }
2195 for (; bb != nb_base_functions; ++bb)
2196 ++t_n_diff_hdiv;
2197 ++t_data;
2198 }
2200 }
2201
2202private:
2203 boost::shared_ptr<VectorDouble> dataPtr;
2204 const EntityHandle zeroType;
2205 const int zeroSide;
2206};
2207
2208/**
2209 * @brief Calculate gradient of vector field
2210 * @ingroup mofem_forces_and_sources_user_data_operators
2211 *
2212 * @tparam BASE_DIM
2213 * @tparam SPACE_DIM
2214 */
2215template <int BASE_DIM, int SPACE_DIM>
2218
2220 boost::shared_ptr<MatrixDouble> data_ptr,
2221 const EntityType zero_type = MBEDGE,
2222 const int zero_side = 0)
2225 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2226 if (!dataPtr)
2227 THROW_MESSAGE("Pointer is not set");
2228 }
2229
2233 const size_t nb_integration_points = getGaussPts().size2();
2234 if (type == zeroType && side == zeroSide) {
2235 dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
2236 dataPtr->clear();
2237 }
2238 const size_t nb_dofs = data.getFieldData().size();
2239 if (!nb_dofs)
2241 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2244 auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2245 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*dataPtr);
2246 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2247 auto t_dof = data.getFTensor0FieldData();
2248 int bb = 0;
2249 for (; bb != nb_dofs; ++bb) {
2250 t_data(i, j) += t_dof * t_base_diff(i, j);
2251 ++t_base_diff;
2252 ++t_dof;
2253 }
2254 for (; bb != nb_base_functions; ++bb)
2255 ++t_base_diff;
2256 ++t_data;
2257 }
2259 }
2260
2261private:
2262 boost::shared_ptr<MatrixDouble> dataPtr;
2263 const EntityHandle zeroType;
2264 const int zeroSide;
2265};
2266
2267/**
2268 * @brief Calculate gradient of vector field
2269 * @ingroup mofem_forces_and_sources_user_data_operators
2270 *
2271 * @tparam BASE_DIM
2272 * @tparam SPACE_DIM
2273 */
2274template <int BASE_DIM, int SPACE_DIM>
2277
2279 boost::shared_ptr<MatrixDouble> data_ptr,
2280 const EntityType zero_type = MBEDGE,
2281 const int zero_side = 0)
2284 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2285 if (!dataPtr)
2286 THROW_MESSAGE("Pointer is not set");
2287 }
2288
2292 const size_t nb_integration_points = getGaussPts().size2();
2293 if (type == zeroType && side == zeroSide) {
2294 dataPtr->resize(BASE_DIM * SPACE_DIM * SPACE_DIM, nb_integration_points,
2295 false);
2296 dataPtr->clear();
2297 }
2298 const size_t nb_dofs = data.getFieldData().size();
2299 if (!nb_dofs)
2301
2302 const int nb_base_functions = data.getN().size2() / BASE_DIM;
2303 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2304
2305#ifndef NDEBUG
2306 if (hessian_base.size1() != nb_integration_points) {
2307 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2308 "Wrong number of integration pts (%d != %d)",
2309 hessian_base.size1(), nb_integration_points);
2310 }
2311 if (hessian_base.size2() !=
2312 BASE_DIM * nb_base_functions * SPACE_DIM * SPACE_DIM) {
2313 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2314 "Wrong number of base functions (%d != %d)",
2315 hessian_base.size2() / (BASE_DIM * SPACE_DIM * SPACE_DIM),
2316 nb_base_functions);
2317 }
2318 if (hessian_base.size2() < BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM) {
2319 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2320 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2321 BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM);
2322 }
2323#endif
2324
2328
2329 auto t_base_diff2 = data.getFTensor3Diff2N<BASE_DIM, SPACE_DIM, SPACE_DIM>();
2330 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*dataPtr);
2331 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2332 auto t_dof = data.getFTensor0FieldData();
2333 int bb = 0;
2334 for (; bb != nb_dofs; ++bb) {
2335 t_data(i, j, k) += t_dof * t_base_diff2(i, j, k);
2336
2337 ++t_base_diff2;
2338 ++t_dof;
2339 }
2340 for (; bb != nb_base_functions; ++bb)
2341 ++t_base_diff2;
2342 ++t_data;
2343 }
2345 }
2346
2347private:
2348 boost::shared_ptr<MatrixDouble> dataPtr;
2349 const EntityHandle zeroType;
2350 const int zeroSide;
2351};
2352
2353/**
2354 * @brief Calculate divergence of vector field dot
2355 * @ingroup mofem_forces_and_sources_user_data_operators
2356 *
2357 * @tparam Tensor_Dim dimension of space
2358 */
2359template <int Tensor_Dim1, int Tensor_Dim2>
2362
2364 boost::shared_ptr<VectorDouble> data_ptr,
2365 const EntityType zero_type = MBEDGE,
2366 const int zero_side = 0)
2369 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2370 if (!dataPtr)
2371 THROW_MESSAGE("Pointer is not set");
2372 }
2373
2377 const size_t nb_integration_points = getGaussPts().size2();
2378 if (type == zeroType && side == zeroSide) {
2379 dataPtr->resize(nb_integration_points, false);
2380 dataPtr->clear();
2381 }
2382
2383 const auto &local_indices = data.getLocalIndices();
2384 const int nb_dofs = local_indices.size();
2385 if (nb_dofs) {
2386
2387 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2388 const double *array;
2389 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2390 for (size_t i = 0; i != local_indices.size(); ++i)
2391 if (local_indices[i] != -1)
2392 dot_dofs_vector[i] = array[local_indices[i]];
2393 else
2394 dot_dofs_vector[i] = 0;
2395 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2396
2397 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
2399 auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
2400 auto t_data = getFTensor0FromVec(*dataPtr);
2401 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2402 int bb = 0;
2403 for (; bb != nb_dofs; ++bb) {
2404 double div = 0;
2405 for (auto ii = 0; ii != Tensor_Dim2; ++ii)
2406 div += t_n_diff_hdiv(ii, ii);
2407 t_data += dot_dofs_vector[bb] * div;
2408 ++t_n_diff_hdiv;
2409 }
2410 for (; bb != nb_base_functions; ++bb)
2411 ++t_n_diff_hdiv;
2412 ++t_data;
2413 }
2414 }
2416 }
2417
2418private:
2419 boost::shared_ptr<VectorDouble> dataPtr;
2420 const EntityHandle zeroType;
2421 const int zeroSide;
2422};
2423
2424/**
2425 * @brief Calculate curl of vector field
2426 * @ingroup mofem_forces_and_sources_user_data_operators
2427 *
2428 * @tparam Tensor_Dim dimension of space
2429 */
2430template <int Tensor_Dim>
2433
2435 boost::shared_ptr<MatrixDouble> data_ptr,
2436 const EntityType zero_type = MBEDGE,
2437 const int zero_side = 0)
2440 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2441 if (!dataPtr)
2442 THROW_MESSAGE("Pointer is not set");
2443 }
2444
2448 const size_t nb_integration_points = getGaussPts().size2();
2449 if (type == zeroType && side == zeroSide) {
2450 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2451 dataPtr->clear();
2452 }
2453 const size_t nb_dofs = data.getFieldData().size();
2454 if (!nb_dofs)
2456
2457 MatrixDouble curl_mat;
2458
2459 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim;
2463 auto t_n_diff_hcurl = data.getFTensor2DiffN<Tensor_Dim, Tensor_Dim>();
2464 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2465 for (int gg = 0; gg != nb_integration_points; ++gg) {
2466
2467 auto t_dof = data.getFTensor0FieldData();
2468 int bb = 0;
2469 for (; bb != nb_dofs; ++bb) {
2470
2471 t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
2472 ++t_n_diff_hcurl;
2473 ++t_dof;
2474 }
2475
2476 for (; bb != nb_base_functions; ++bb)
2477 ++t_n_diff_hcurl;
2478 ++t_data;
2479 }
2480
2482 }
2483
2484private:
2485 boost::shared_ptr<MatrixDouble> dataPtr;
2486 const EntityHandle zeroType;
2487 const int zeroSide;
2488};
2489
2490/**
2491 * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
2492 * \ingroup mofem_forces_and_sources_user_data_operators
2493 *
2494 * @tparam Tensor_Dim0 rank of the filed
2495 * @tparam Tensor_Dim1 dimension of space
2496 */
2497template <int Tensor_Dim0, int Tensor_Dim1>
2500
2502 boost::shared_ptr<MatrixDouble> data_ptr,
2503 const EntityType zero_type = MBEDGE,
2504 const int zero_side = 0)
2507 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2508 if (!dataPtr)
2509 THROW_MESSAGE("Pointer is not set");
2510 }
2511
2515 const size_t nb_integration_points = getGaussPts().size2();
2516 if (type == zeroType && side == zeroSide) {
2517 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2518 dataPtr->clear();
2519 }
2520 const size_t nb_dofs = data.getFieldData().size();
2521 if (nb_dofs) {
2522 const size_t nb_base_functions = data.getN().size2() / 3;
2525 auto t_n_hvec = data.getFTensor1N<3>();
2526 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2527 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2528 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2529 size_t bb = 0;
2530 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2531 t_data(i, j) += t_dof(i) * t_n_hvec(j);
2532 ++t_n_hvec;
2533 ++t_dof;
2534 }
2535 for (; bb != nb_base_functions; ++bb)
2536 ++t_n_hvec;
2537 ++t_data;
2538 }
2539 }
2541 }
2542
2543private:
2544 boost::shared_ptr<MatrixDouble> dataPtr;
2545 const EntityHandle zeroType;
2546 const int zeroSide;
2547};
2548
2549/**
2550 * @brief Calculate tenor field using tensor base, i.e. Hdiv/Hcurl
2551 * \ingroup mofem_forces_and_sources_user_data_operators
2552 *
2553 * @tparam Tensor_Dim0 rank of the filed
2554 * @tparam Tensor_Dim1 dimension of space
2555 */
2556template <int Tensor_Dim0, int Tensor_Dim1>
2559
2561 boost::shared_ptr<MatrixDouble> data_ptr,
2562 const EntityType zero_type = MBEDGE,
2563 const int zero_side = 0)
2566 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2567 if (!dataPtr)
2568 THROW_MESSAGE("Pointer is not set");
2569 }
2570
2574 const size_t nb_integration_points = getGaussPts().size2();
2575 if (type == zeroType && side == zeroSide) {
2576 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2577 dataPtr->clear();
2578 }
2579 const size_t nb_dofs = data.getFieldData().size();
2580 if (!nb_dofs)
2582 const size_t nb_base_functions =
2583 data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2586 auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2587 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2588 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2589 auto t_dof = data.getFTensor0FieldData();
2590 size_t bb = 0;
2591 for (; bb != nb_dofs; ++bb) {
2592 t_data(i, j) += t_dof * t_n_hten(i, j);
2593 ++t_n_hten;
2594 ++t_dof;
2595 }
2596 for (; bb != nb_base_functions; ++bb)
2597 ++t_n_hten;
2598 ++t_data;
2599 }
2601 }
2602
2603private:
2604 boost::shared_ptr<MatrixDouble> dataPtr;
2605 const EntityHandle zeroType;
2606 const int zeroSide;
2607};
2608
2609/**
2610 * @brief Calculate divergence of tonsorial field using vectorial base
2611 * \ingroup mofem_forces_and_sources_user_data_operators
2612 *
2613 * @tparam Tensor_Dim0 rank of the field
2614 * @tparam Tensor_Dim1 dimension of space
2615 */
2616template <int Tensor_Dim0, int Tensor_Dim1>
2619
2621 boost::shared_ptr<MatrixDouble> data_ptr,
2622 const EntityType zero_type = MBEDGE,
2623 const int zero_side = 0)
2626 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2627 if (!dataPtr)
2628 THROW_MESSAGE("Pointer is not set");
2629 }
2630
2634 const size_t nb_integration_points = getGaussPts().size2();
2635 if (type == zeroType && side == 0) {
2636 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
2637 dataPtr->clear();
2638 }
2639 const size_t nb_dofs = data.getFieldData().size();
2640 if (nb_dofs) {
2641 const size_t nb_base_functions = data.getN().size2() / 3;
2644 auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
2645 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
2646 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2647 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2648 size_t bb = 0;
2649 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2650 double div = t_n_diff_hvec(j, j);
2651 t_data(i) += t_dof(i) * div;
2652 ++t_n_diff_hvec;
2653 ++t_dof;
2654 }
2655 for (; bb < nb_base_functions; ++bb)
2656 ++t_n_diff_hvec;
2657 ++t_data;
2658 }
2659 }
2661 }
2662
2663private:
2664 boost::shared_ptr<MatrixDouble> dataPtr;
2665 const EntityHandle zeroType;
2666 const int zeroSide;
2667};
2668
2669/**
2670 * @brief Calculate trace of vector (Hdiv/Hcurl) space
2671 *
2672 * @tparam Tensor_Dim
2673 * @tparam OpBase
2674 */
2675template <int Tensor_Dim, typename OpBase>
2677
2679 boost::shared_ptr<MatrixDouble> data_ptr,
2680 const EntityType zero_type = MBEDGE,
2681 const int zero_side = 0)
2682 : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
2683 zeroType(zero_type), zeroSide(zero_side) {
2684 if (!dataPtr)
2685 THROW_MESSAGE("Pointer is not set");
2686 }
2687
2691 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2692 if (type == zeroType && side == 0) {
2693 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2694 dataPtr->clear();
2695 }
2696 const size_t nb_dofs = data.getFieldData().size();
2697 if (nb_dofs) {
2698 auto t_normal = OpBase::getFTensor1Normal();
2699 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
2700 const size_t nb_base_functions = data.getN().size2() / 3;
2701 auto t_base = data.getFTensor1N<3>();
2702 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2703 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2704 auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
2705 size_t bb = 0;
2706 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2707 t_data(i) += t_dof(i) * (t_base(j) * t_normal(j));
2708 ++t_base;
2709 ++t_dof;
2710 }
2711 for (; bb < nb_base_functions; ++bb)
2712 ++t_base;
2713 ++t_data;
2714 }
2715 }
2717 }
2718
2719private:
2720 boost::shared_ptr<MatrixDouble> dataPtr;
2721 const EntityHandle zeroType;
2722 const int zeroSide;
2725};
2726
2727template <>
2731
2733 boost::shared_ptr<MatrixDouble> data_ptr,
2734 const EntityType zero_type = MBEDGE,
2735 const int zero_side = 0)
2736 : UserDataOperator(field_name, OPROW), dataPtr(data_ptr),
2737 zeroType(zero_type), zeroSide(zero_side) {
2738 if (!dataPtr)
2739 THROW_MESSAGE("Pointer is not set");
2740 }
2741
2745 const size_t nb_integration_points = getGaussPts().size2();
2746 if (type == zeroType && side == 0) {
2747 dataPtr->resize(3, nb_integration_points, false);
2748 dataPtr->clear();
2749 }
2750 const size_t nb_dofs = data.getFieldData().size();
2751 if (nb_dofs) {
2752 auto t_normal = getFTensor1Normal();
2753 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
2754 const size_t nb_base_functions = data.getN().size2() / 3;
2755 auto t_base = data.getFTensor1N<3>();
2756 auto t_data = getFTensor1FromMat<3>(*dataPtr);
2757 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2758 auto t_dof = data.getFTensor1FieldData<3>();
2759 if (getNormalsAtGaussPts().size1() == nb_integration_points) {
2760 VectorDouble n = getNormalsAtGaussPts(gg);
2761 auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
2762 t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
2763 }
2764 size_t bb = 0;
2765 for (; bb != nb_dofs / 3; ++bb) {
2766 t_data(i) += t_dof(i) * (t_base(j) * t_normal(j));
2767 ++t_base;
2768 ++t_dof;
2769 }
2770 for (; bb < nb_base_functions; ++bb)
2771 ++t_base;
2772 ++t_data;
2773 }
2774 }
2776 }
2777
2778private:
2779 boost::shared_ptr<MatrixDouble> dataPtr;
2780 const EntityHandle zeroType;
2781 const int zeroSide;
2784};
2785
2786/**@}*/
2787
2788/** \name Other operators */
2789
2790/**@{*/
2791
2792/**@}*/
2793
2794/** \name Operators for faces */
2795
2796/**@{*/
2797
2798/** \brief Transform local reference derivatives of shape functions to global
2799derivatives
2800
2801\ingroup mofem_forces_and_sources_tri_element
2802
2803*/
2804template <int DIM, int DERIVATIVE = 1> struct OpSetInvJacSpaceForFaceImpl;
2805
2808
2810 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2811
2812protected:
2813 template <int D1, int D2, int J1, int J2>
2816 size_t nb_functions = diff_n.size2() / D1;
2817 if (nb_functions) {
2818 size_t nb_gauss_pts = diff_n.size1();
2819
2820#ifndef NDEBUG
2821 if (nb_gauss_pts != getGaussPts().size2())
2822 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2823 "Wrong number of Gauss Pts");
2824 if (diff_n.size2() % D1)
2825 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2826 "Number of direvatives of base functions and D1 dimension does "
2827 "not match");
2828#endif
2829
2830 diffNinvJac.resize(diff_n.size1(), diff_n.size2(), false);
2831
2834 auto t_diff_n = getFTensor1FromPtr<D2>(&*diffNinvJac.data().begin());
2835 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2836 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*invJacPtr);
2837 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2838 for (size_t dd = 0; dd != nb_functions; ++dd) {
2839 t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
2840 ++t_diff_n;
2841 ++t_diff_n_ref;
2842 }
2843 }
2844
2845 diff_n.swap(diffNinvJac);
2846 }
2848 }
2849
2850 boost::shared_ptr<MatrixDouble> invJacPtr;
2852};
2853
2854template <>
2857
2859
2860 MoFEMErrorCode doWork(int side, EntityType type,
2862};
2863
2864template <>
2866 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2867
2868 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2869
2870 MoFEMErrorCode doWork(int side, EntityType type,
2872};
2873
2874template <>
2876 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2877
2878 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2879
2880 MoFEMErrorCode doWork(int side, EntityType type,
2882};
2883
2884template <int DERIVARIVE = 1>
2886 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2887 OpSetInvJacH1ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2888 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(H1, inv_jac_ptr) {}
2889};
2890
2891template <int DERIVARIVE = 1>
2893 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2894 OpSetInvJacL2ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2895 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(L2, inv_jac_ptr) {}
2896};
2897
2898template <int DERIVARIVE = 1>
2900 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2902 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2903 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(H1, inv_jac_ptr) {}
2904};
2905
2906template <int DERIVARIVE = 1>
2908 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2910 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2911 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(L2, inv_jac_ptr) {}
2912};
2913
2914/**
2915 * \brief Transform local reference derivatives of shape function to
2916 global derivatives for face
2917
2918 * \ingroup mofem_forces_and_sources_tri_element
2919 */
2920template <int DIM> struct OpSetInvJacHcurlFaceImpl;
2921
2922template <>
2925
2926 OpSetInvJacHcurlFaceImpl(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2928 invJacPtr(inv_jac_ptr) {}
2929
2930 MoFEMErrorCode doWork(int side, EntityType type,
2932
2933protected:
2934 boost::shared_ptr<MatrixDouble> invJacPtr;
2936};
2937
2938template <>
2940 using OpSetInvJacHcurlFaceImpl<2>::OpSetInvJacHcurlFaceImpl;
2941 MoFEMErrorCode doWork(int side, EntityType type,
2943};
2944
2947
2948/**
2949 * @brief Make Hdiv space from Hcurl space in 2d
2950 * @ingroup mofem_forces_and_sources_tri_element
2951 */
2954
2957
2960};
2961
2962/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
2963 *
2964 * \note Hdiv space is generated by Hcurl space in 2d.
2965 *
2966 * Contravariant Piola transformation
2967 * \f[
2968 * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
2969 * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
2970 * =
2971 * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
2972 * \f]
2973 *
2974 * \ingroup mofem_forces_and_sources
2975 *
2976 */
2978
2979template <>
2982
2984 boost::shared_ptr<MatrixDouble> jac_ptr)
2986 jacPtr(jac_ptr) {}
2987
2988 MoFEMErrorCode doWork(int side, EntityType type,
2990
2991protected:
2992 boost::shared_ptr<MatrixDouble> jacPtr;
2995};
2996
2997template <>
3001 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3002
3003 MoFEMErrorCode doWork(int side, EntityType type,
3005};
3006
3011
3012/**@}*/
3013
3014/** \name Operators for edges */
3015
3016/**@{*/
3017
3020
3023
3026
3027private:
3028 std::vector<double> l1;
3029};
3030
3031/**
3032 * @deprecated Name is deprecated and this is added for back compatibility
3033 */
3036
3037/**@}*/
3038
3039/** \name Operator for fat prisms */
3040
3041/**@{*/
3042
3043/**
3044 * @brief Operator for fat prism element updating integration weights in the
3045 * volume.
3046 *
3047 * Jacobian on the distorted element is nonconstant. This operator updates
3048 * integration weight on prism to take into account nonconstat jacobian.
3049 *
3050 * \f[
3051 * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
3052 * \pmb\xi} \right\| \right)
3053 * \f]
3054 * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
3055 * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
3056 * element coordinate.
3057 *
3058 */
3061
3064
3067};
3068
3069/** \brief Calculate inverse of jacobian for face element
3070
3071 It is assumed that face element is XY plane. Applied
3072 only for 2d problems.
3073
3074 FIXME Generalize function for arbitrary face orientation in 3d space
3075 FIXME Calculate to Jacobins for two faces
3076
3077 \ingroup mofem_forces_and_sources_prism_element
3078
3079*/
3082
3083 OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3085 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3086
3089 invJac(inv_jac) {}
3092
3093private:
3094 const boost::shared_ptr<MatrixDouble> invJacPtr;
3096};
3097
3098/** \brief Transform local reference derivatives of shape functions to global
3099derivatives
3100
3101FIXME Generalize to curved shapes
3102FIXME Generalize to case that top and bottom face has different shape
3103
3104\ingroup mofem_forces_and_sources_prism_element
3105
3106*/
3109
3110 OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3112 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3113
3116 invJac(inv_jac) {}
3117
3120
3121private:
3122 const boost::shared_ptr<MatrixDouble> invJacPtr;
3125};
3126
3127// Flat prism
3128
3129/** \brief Calculate inverse of jacobian for face element
3130
3131 It is assumed that face element is XY plane. Applied
3132 only for 2d problems.
3133
3134 FIXME Generalize function for arbitrary face orientation in 3d space
3135 FIXME Calculate to Jacobins for two faces
3136
3137 \ingroup mofem_forces_and_sources_prism_element
3138
3139*/
3142
3145 invJacF3(inv_jac_f3) {}
3148
3149private:
3151};
3152
3153/** \brief Transform local reference derivatives of shape functions to global
3154derivatives
3155
3156FIXME Generalize to curved shapes
3157FIXME Generalize to case that top and bottom face has different shape
3158
3159\ingroup mofem_forces_and_sources_prism_element
3160
3161*/
3164
3167 invJacF3(inv_jac_f3) {}
3168
3171
3172private:
3175};
3176
3177/**@}*/
3178
3179/** \name Operation on matrices at integration points */
3180
3181/**@{*/
3182
3183template <int DIM>
3185
3186 OpInvertMatrix(boost::shared_ptr<MatrixDouble> in_ptr,
3187 boost::shared_ptr<VectorDouble> det_ptr,
3188 boost::shared_ptr<MatrixDouble> out_ptr)
3190 outPtr(out_ptr), detPtr(det_ptr) {}
3191
3194 return doWorkImpl(side, type, data, FTensor::Number<DIM>());
3195 }
3196
3197private:
3198 boost::shared_ptr<MatrixDouble> inPtr;
3199 boost::shared_ptr<MatrixDouble> outPtr;
3200 boost::shared_ptr<VectorDouble> detPtr;
3201
3204 const FTensor::Number<3> &);
3205
3208 const FTensor::Number<2> &);
3209};
3210
3211template <int DIM>
3214 const FTensor::Number<3> &) {
3216
3217 if (!inPtr)
3218 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3219 "Pointer for inPtr matrix not allocated");
3220 if (!detPtr)
3221 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3222 "Pointer for detPtr matrix not allocated");
3223
3224 const auto nb_rows = inPtr->size1();
3225 const auto nb_integration_pts = inPtr->size2();
3226
3227 // Calculate determinant
3228 {
3229 detPtr->resize(nb_integration_pts, false);
3230 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3231 auto t_det = getFTensor0FromVec(*detPtr);
3232 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3233 determinantTensor3by3(t_in, t_det);
3234 ++t_in;
3235 ++t_det;
3236 }
3237 }
3238
3239 // Invert jacobian
3240 if (outPtr) {
3241 outPtr->resize(nb_rows, nb_integration_pts, false);
3242 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3243 auto t_out = getFTensor2FromMat<3, 3>(*outPtr);
3244 auto t_det = getFTensor0FromVec(*detPtr);
3245 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3246 invertTensor3by3(t_in, t_det, t_out);
3247 ++t_in;
3248 ++t_out;
3249 ++t_det;
3250 }
3251 }
3252
3254}
3255
3256template <int DIM>
3259 const FTensor::Number<2> &) {
3261
3262 if (!inPtr)
3263 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3264 "Pointer for inPtr matrix not allocated");
3265 if (!detPtr)
3266 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3267 "Pointer for detPtr matrix not allocated");
3268
3269 const auto nb_rows = inPtr->size1();
3270 const auto nb_integration_pts = inPtr->size2();
3271
3272 // Calculate determinant
3273 {
3274 detPtr->resize(nb_integration_pts, false);
3275 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3276 auto t_det = getFTensor0FromVec(*detPtr);
3277 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3278 determinantTensor2by2(t_in, t_det);
3279 ++t_in;
3280 ++t_det;
3281 }
3282 }
3283
3284 // Invert jacobian
3285 if (outPtr) {
3286 outPtr->resize(nb_rows, nb_integration_pts, false);
3287 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3288 auto t_out = getFTensor2FromMat<2, 2>(*outPtr);
3289 auto t_det = getFTensor0FromVec(*detPtr);
3290 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3291 invertTensor2by2(t_in, t_det, t_out);
3292 ++t_in;
3293 ++t_out;
3294 ++t_det;
3295 }
3296 }
3297
3299}
3300
3301/**@}*/
3302
3303} // namespace MoFEM
3304
3305#endif // __USER_DATA_OPERATORS_HPP__
3306
3307/**
3308 * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
3309 *
3310 * \brief Classes and functions used to evaluate fields at integration pts,
3311 *jacobians, etc..
3312 *
3313 * \ingroup mofem_forces_and_sources
3314 **/
static Index< 'n', 3 > n
static Index< 'J', 3 > J
static Index< 'I', 3 > I
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
constexpr int SPACE_DIM
static const char *const CoordinateTypesNames[]
Coordinate system names.
Definition: definitions.h:121
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FieldSpace
approximation spaces
Definition: definitions.h:95
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ NOSPACE
Definition: definitions.h:96
@ HCURL
field with continuous tangents
Definition: definitions.h:99
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
@ CYLINDRICAL
Definition: definitions.h:130
@ POLAR
Definition: definitions.h:129
@ CARTESIAN
Definition: definitions.h:128
@ SPHERICAL
Definition: definitions.h:131
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#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
constexpr int BASE_DIM
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const double T
EntityType zero_type
@ row_major
Definition: Layout.hpp:13
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VecAllocator< double > DoubleAllocator
Definition: Types.hpp:73
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1248
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1279
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 auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
Definition: Templates.hpp:747
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
double A
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3Diff2N(FieldApproximationBase base)
Get second derivatives of base functions for Hvec space.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(Tensor_Dim *(Tensor_Dim+1))/2 >, Tensor_Dim > getFTensor2SymmetricFieldData()
Return symmetric FTensor rank 2, i.e. matrix from filed data coeffinects.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from filed data coeffinects.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
boost::shared_ptr< VectorDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateDivergenceVectorFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate trace of vector (Hdiv/Hcurl) space.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Index< 'j', Tensor_Dim > j
FTensor::Index< 'i', Tensor_Dim > i
OpCalculateHVecVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
OpCalculateHVecVectorField_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A > > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Get vector field for H-div approximation.
OpCalculateHVecVectorFieldDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate values of vector field at integration points.
boost::shared_ptr< MatrixDouble > dataPtr
Get vector field for H-div approximation.
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecVectorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecVectorHessian(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate curl of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHcurlVectorCurl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate divergence of vector field dot.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergenceDot(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate divergence of vector field.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergence(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFatPrism(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateInvJacForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
const boost::shared_ptr< MatrixDouble > invJacPtr
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFlatPrism(MatrixDouble &inv_jac_f3)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector),...
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate gradient values of scalar field at integration points
Scalar field values at integration points.
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A > > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< ublas::vector< T, A > > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar field at integration points
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A > > data_ptr, const EntityType zero_type=MBVERTEX)
Get rate of scalar field at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateScalarFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
Get value at integration points for scalar field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar field at integration points
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator > > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBVERTEX)
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator > > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate field values for tenor field rank 2.
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A > > data_ptr, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get time direvarive values at integration pts for tensor filed rank 2, i.e. matrix field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateTensor2FieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
EntityType zeroAtType
Zero values at Gauss point at this type.
VectorDouble dotVector
Keeps temoorary values of time directives.
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
OpCalculateTensor2SymmetricFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Evaluate field gradient values for symmetric 2nd order tensor field, i.e. gradient is tensor rank 3.
OpCalculateTensor2SymmetricFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Get field gradients at integration pts for symmetric tensorial field rank 2.
OpCalculateTensor2SymmetricFieldGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate symmetric tensor field rates ant integratio pts.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateTensor2SymmetricFieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate symmetric tensor field values at integration pts.
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBEDGE, const int zero_side=0)
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Evaluate field gradient values for vector field, i.e. gradient is tensor rank 2.
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
OpCalculateVectorFieldGradientDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
EntityType zeroAtType
Zero values at Gauss point at this type.
VectorDouble dotVector
Keeps temoorary values of time directives.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate field values for tenor field rank 1, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A > > data_ptr, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Approximate field valuse for given petsc vector.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateVectorFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
Get values at integration pts for tensor filed rank 1, i.e. vector field.
boost::shared_ptr< VectorDouble > detPtr
OpInvertMatrix(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > det_ptr, boost::shared_ptr< MatrixDouble > out_ptr)
boost::shared_ptr< MatrixDouble > outPtr
boost::shared_ptr< MatrixDouble > inPtr
MoFEMErrorCode doWorkImpl(int side, EntityType type, EntitiesFieldData::EntData &data, const FTensor::Number< 3 > &)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Make Hdiv space from Hcurl space in 2d.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Operator for fat prism element updating integration weights in the volume.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpScaleMatrix(const std::string field_name, const double scale, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< MatrixDouble > outMat
boost::shared_ptr< MatrixDouble > inMat
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetContravariantPiolaTransformOnEdge2D(const FieldSpace space=HCURL)
OpSetContravariantPiolaTransformOnFace2DImpl(boost::shared_ptr< MatrixDouble > jac_ptr)
Apply contravariant (Piola) transfer to Hdiv space on face.
OpSetInvJacH1ForFaceEmbeddedIn3DSpace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
OpSetInvJacH1ForFace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
const boost::shared_ptr< MatrixDouble > invJacPtr
OpSetInvJacH1ForFatPrism(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1ForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1ForFlatPrism(MatrixDouble &inv_jac_f3)
boost::shared_ptr< MatrixDouble > invJacPtr
OpSetInvJacHcurlFaceImpl(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape function to global derivatives for face.
OpSetInvJacL2ForFaceEmbeddedIn3DSpace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
OpSetInvJacL2ForFace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode applyTransform(MatrixDouble &diff_n)
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
boost::shared_ptr< MatrixDouble > invJacPtr
FTensor::Index< 'i', DIM > i
boost::shared_ptr< MatrixDouble > outMat
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSymmetrizeTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
FTensor::Index< 'j', DIM > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'i', DIM_01 > i
OpTensorTimesSymmetricTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'k', DIM_23 > k
boost::shared_ptr< MatrixDouble > dMat
FTensor::Index< 'j', DIM_01 > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'l', DIM_23 > l
boost::shared_ptr< MatrixDouble > outMat
static constexpr Switches CtxSetX
Definition: LoopMethods.hpp:49
static constexpr Switches CtxSetX_TT
Definition: LoopMethods.hpp:51
static constexpr Switches CtxSetX_T
Definition: LoopMethods.hpp:50