v0.13.1
Loading...
Searching...
No Matches
UserDataOperators.hpp
Go to the documentation of this file.
1/** \file UserDataOperators.hpp
2 * \brief User data Operators
3
4*/
5
6#ifndef __USER_DATA_OPERATORS_HPP__
7#define __USER_DATA_OPERATORS_HPP__
8
9namespace MoFEM {
10
11/** \name Get values at Gauss pts */
12
13/**@{*/
14
15/** \name Scalar values */
16
17/**@{*/
18
19/** \brief Scalar field values at integration points
20 *
21 */
22template <class T, class A>
25
27 const std::string field_name,
28 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
29 const EntityType zero_type = MBVERTEX)
31 dataPtr(data_ptr), zeroType(zero_type) {
32 if (!dataPtr)
33 THROW_MESSAGE("Pointer is not set");
34 }
35
37 const std::string field_name,
38 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
39 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
41 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
42 if (!dataPtr)
43 THROW_MESSAGE("Pointer is not set");
44 }
45
46 /**
47 * \brief calculate values of scalar field at integration points
48 * @param side side entity number
49 * @param type side entity type
50 * @param data entity data
51 * @return error code
52 */
55
56protected:
57 boost::shared_ptr<ublas::vector<T, A>> dataPtr;
61};
62
63/**
64 * \brief Specialization of member function
65 *
66 */
67template <class T, class A>
69 int side, EntityType type, EntitiesFieldData::EntData &data) {
71 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented for T = %s",
72 typeid(T).name() // boost::core::demangle(typeid(T).name()).c_str()
73 );
75}
76
77/**
78 * \brief Get value at integration points for scalar field
79 *
80 * \ingroup mofem_forces_and_sources_user_data_operators
81 */
83 : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
84
86 double, DoubleAllocator>::OpCalculateScalarFieldValues_General;
87
88 /**
89 * \brief calculate values of scalar field at integration points
90 * @param side side entity number
91 * @param type side entity type
92 * @param data entity data
93 * @return error code
94 */
98 VectorDouble &vec = *dataPtr;
99 const size_t nb_gauss_pts = getGaussPts().size2();
100 if (type == zeroType) {
101 vec.resize(nb_gauss_pts, false);
102 vec.clear();
103 }
104
105 const size_t nb_dofs = data.getFieldData().size();
106
107 if (nb_dofs) {
108
109 if (dataVec.use_count()) {
110 dotVector.resize(nb_dofs, false);
111 const double *array;
112 CHKERR VecGetArrayRead(dataVec, &array);
113 const auto &local_indices = data.getLocalIndices();
114 for (int i = 0; i != local_indices.size(); ++i)
115 if (local_indices[i] != -1)
116 dotVector[i] = array[local_indices[i]];
117 else
118 dotVector[i] = 0;
119 CHKERR VecRestoreArrayRead(dataVec, &array);
120 data.getFieldData().swap(dotVector);
121 }
122
123 const size_t nb_base_functions = data.getN().size2();
124 auto base_function = data.getFTensor0N();
125 auto values_at_gauss_pts = getFTensor0FromVec(vec);
126 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
127 auto field_data = data.getFTensor0FieldData();
128 size_t bb = 0;
129 for (; bb != nb_dofs; ++bb) {
130 values_at_gauss_pts += field_data * base_function;
131 ++field_data;
132 ++base_function;
133 }
134 // It is possible to have more base functions than dofs
135 for (; bb != nb_base_functions; ++bb)
136 ++base_function;
137 ++values_at_gauss_pts;
138 }
139
140 if (dataVec.use_count()) {
141 data.getFieldData().swap(dotVector);
142 }
143 }
144
146 }
147};
148
149/**
150 * @brief Get rate of scalar field at integration points
151 *
152 * \ingroup mofem_forces_and_sources_user_data_operators
153 */
154
155template <PetscData::DataContext CTX>
158
160 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
161 const EntityType zero_at_type = MBVERTEX)
164 dataPtr(data_ptr), zeroAtType(zero_at_type) {
165 if (!dataPtr)
166 THROW_MESSAGE("Pointer is not set");
167 }
168
172
173 const size_t nb_gauss_pts = getGaussPts().size2();
174
175 VectorDouble &vec = *dataPtr;
176 if (type == zeroAtType) {
177 vec.resize(nb_gauss_pts, false);
178 vec.clear();
179 }
180
181 auto &local_indices = data.getLocalIndices();
182 const size_t nb_dofs = local_indices.size();
183 if (nb_dofs) {
184
185 const double *array;
186
187 auto get_array = [&](const auto ctx, auto vec) {
189#ifndef NDEBUG
190 if ((getFEMethod()->data_ctx & ctx).none()) {
191 MOFEM_LOG_CHANNEL("SELF");
192 MOFEM_LOG("SELF", Sev::error)
193 << "In this case filed degrees of freedom are read from vector. "
194 "That usually happens when time solver is used, and acces to "
195 "first or second rates is needed. You probably not set ts_u, "
196 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
197 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
198 "respectively";
199 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
200 }
201#endif
202 CHKERR VecGetArrayRead(vec, &array);
204 };
205
206 auto restore_array = [&](auto vec) {
207 return VecRestoreArrayRead(vec, &array);
208 };
209
210 switch (CTX) {
212 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
213 break;
215 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
216 break;
218 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
219 break;
220 default:
221 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
222 "That case is not implemented");
223 }
224
225 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
226 for (int i = 0; i != local_indices.size(); ++i)
227 if (local_indices[i] != -1)
228 dot_dofs_vector[i] = array[local_indices[i]];
229 else
230 dot_dofs_vector[i] = 0;
231
232 switch (CTX) {
234 CHKERR restore_array(getFEMethod()->ts_u);
235 break;
237 CHKERR restore_array(getFEMethod()->ts_u_t);
238 break;
240 CHKERR restore_array(getFEMethod()->ts_u_tt);
241 break;
242 default:
243 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
244 "That case is not implemented");
245 }
246
247 const size_t nb_base_functions = data.getN().size2();
248 auto base_function = data.getFTensor0N();
249 auto values_at_gauss_pts = getFTensor0FromVec(vec);
250
251 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
252 size_t bb = 0;
253 for (; bb != nb_dofs; ++bb) {
254 values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
255 ++base_function;
256 }
257 // Number of dofs can be smaller than number of Tensor_Dim x base
258 // functions
259 for (; bb != nb_base_functions; ++bb)
260 ++base_function;
261 ++values_at_gauss_pts;
262 }
263 }
265 }
266
267private:
268 boost::shared_ptr<VectorDouble> dataPtr;
270};
271
276
277/**
278 * \deprecated Name inconsistent with other operators
279 *
280 */
282
283/**@}*/
284
285/** \name Vector field values at integration points */
286
287/**@{*/
288
289/** \brief Calculate field values for tenor field rank 1, i.e. vector field
290 *
291 */
292template <int Tensor_Dim, class T, class L, class A>
295
297 const std::string field_name,
298 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
299 const EntityType zero_type = MBVERTEX)
302 dataPtr(data_ptr), zeroType(zero_type) {
303 if (!dataPtr)
304 THROW_MESSAGE("Pointer is not set");
305 }
306
307 /**
308 * \brief calculate values of vector field at integration points
309 * @param side side entity number
310 * @param type side entity type
311 * @param data entity data
312 * @return error code
313 */
314 MoFEMErrorCode doWork(int side, EntityType type,
316
317protected:
318 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
320};
321
322template <int Tensor_Dim, class T, class L, class A>
325 int side, EntityType type, EntitiesFieldData::EntData &data) {
327 SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
328 "Not implemented for T = %s and dim = %d",
329 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
330 Tensor_Dim);
332}
333
334/** \brief Calculate field values (template specialization) for tensor field
335 * rank 1, i.e. vector field
336 *
337 */
338template <int Tensor_Dim>
340 ublas::row_major, DoubleAllocator>
342
344 boost::shared_ptr<MatrixDouble> data_ptr,
345 const EntityType zero_type = MBVERTEX)
348 dataPtr(data_ptr), zeroType(zero_type) {
349 if (!dataPtr)
350 THROW_MESSAGE("Pointer is not set");
351 }
352
353 MoFEMErrorCode doWork(int side, EntityType type,
355
356protected:
357 boost::shared_ptr<MatrixDouble> dataPtr;
359};
360
361/**
362 * \brief Member function specialization calculating values for tenor field rank
363 *
364 */
365template <int Tensor_Dim>
367 Tensor_Dim, double, ublas::row_major,
368 DoubleAllocator>::doWork(int side, EntityType type,
371
372 const size_t nb_gauss_pts = getGaussPts().size2();
373 auto &mat = *dataPtr;
374 if (type == zeroType) {
375 mat.resize(Tensor_Dim, nb_gauss_pts, false);
376 mat.clear();
377 }
378
379 const size_t nb_dofs = data.getFieldData().size();
380 if (nb_dofs) {
381
382 if (nb_gauss_pts) {
383 const size_t nb_base_functions = data.getN().size2();
384 auto base_function = data.getFTensor0N();
385 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
387 const size_t size = nb_dofs / Tensor_Dim;
388 if (nb_dofs % Tensor_Dim) {
389 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
390 "Data inconsistency");
391 }
392 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
393 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
394 size_t bb = 0;
395 for (; bb != size; ++bb) {
396 values_at_gauss_pts(I) += field_data(I) * base_function;
397 ++field_data;
398 ++base_function;
399 }
400 // Number of dofs can be smaller than number of Tensor_Dim x base
401 // functions
402 for (; bb != nb_base_functions; ++bb)
403 ++base_function;
404 ++values_at_gauss_pts;
405 }
406 }
407 }
409}
410
411/** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
412 * field
413 *
414 * \ingroup mofem_forces_and_sources_user_data_operators
415 */
416template <int Tensor_Dim>
419 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
420
422 Tensor_Dim, double, ublas::row_major,
423 DoubleAllocator>::OpCalculateVectorFieldValues_General;
424};
425
426/**@}*/
427
428/** \name Vector field values at integration points */
429
430/**@{*/
431
432/** \brief Calculate field values (template specialization) for tensor field
433 * rank 1, i.e. vector field
434 *
435 */
436template <int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
439
441 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
442 const EntityType zero_type = MBVERTEX)
445 dataPtr(data_ptr), zeroType(zero_type) {
446 if (!dataPtr)
447 THROW_MESSAGE("Pointer is not set");
448 }
449
453
454 // When we move to C++17 add if constexpr()
455 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
456 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
457 "%s coordiante not implemented",
458 CoordinateTypesNames[COORDINATE_SYSTEM]);
459
460 const size_t nb_gauss_pts = getGaussPts().size2();
461 auto &vec = *dataPtr;
462 if (type == zeroType) {
463 vec.resize(nb_gauss_pts, false);
464 vec.clear();
465 }
466
467 const size_t nb_dofs = data.getFieldData().size();
468 if (nb_dofs) {
469
470 if (nb_gauss_pts) {
471 const size_t nb_base_functions = data.getN().size2();
472 auto values_at_gauss_pts = getFTensor0FromVec(vec);
474 const size_t size = nb_dofs / Tensor_Dim;
475#ifndef NDEBUG
476 if (nb_dofs % Tensor_Dim) {
477 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
478 "Number of dofs should multiple of dimensions");
479 }
480#endif
481
482 // When we move to C++17 add if constexpr()
483 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
484 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
485 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
486 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
487 size_t bb = 0;
488 for (; bb != size; ++bb) {
489 values_at_gauss_pts += field_data(I) * diff_base_function(I);
490 ++field_data;
491 ++diff_base_function;
492 }
493 // Number of dofs can be smaller than number of Tensor_Dim x base
494 // functions
495 for (; bb != nb_base_functions; ++bb)
496 ++diff_base_function;
497 ++values_at_gauss_pts;
498 }
499 }
500
501 // When we move to C++17 add if constexpr()
502 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
503 auto t_coords = getFTensor1CoordsAtGaussPts();
504 auto values_at_gauss_pts = getFTensor0FromVec(vec);
505 auto base_function = data.getFTensor0N();
506 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
507 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
508 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
509 size_t bb = 0;
510 for (; bb != size; ++bb) {
511 values_at_gauss_pts += field_data(I) * diff_base_function(I);
512 values_at_gauss_pts +=
513 field_data(0) * base_function / t_coords(0);
514 ++field_data;
515 ++base_function;
516 ++diff_base_function;
517 }
518 // Number of dofs can be smaller than number of Tensor_Dim x base
519 // functions
520 for (; bb != nb_base_functions; ++bb) {
521 ++base_function;
522 ++diff_base_function;
523 }
524 ++values_at_gauss_pts;
525 ++t_coords;
526 }
527 }
528 }
529 }
531 }
532
533protected:
534 boost::shared_ptr<VectorDouble> dataPtr;
536};
537
538/** \brief Approximate field valuse for given petsc vector
539 *
540 * \note Look at PetscData to see what vectors could be extarcted with that user
541 * data operator.
542 *
543 * \ingroup mofem_forces_and_sources_user_data_operators
544 */
545template <int Tensor_Dim, PetscData::DataContext CTX>
548
550 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
551 const EntityType zero_at_type = MBVERTEX)
554 dataPtr(data_ptr), zeroAtType(zero_at_type) {
555 if (!dataPtr)
556 THROW_MESSAGE("Pointer is not set");
557 }
558
562 auto &local_indices = data.getLocalIndices();
563 const size_t nb_dofs = local_indices.size();
564 const size_t nb_gauss_pts = getGaussPts().size2();
565
566 MatrixDouble &mat = *dataPtr;
567 if (type == zeroAtType) {
568 mat.resize(Tensor_Dim, nb_gauss_pts, false);
569 mat.clear();
570 }
571 if (!nb_dofs)
573
574 const double *array;
575
576 auto get_array = [&](const auto ctx, auto vec) {
578#ifndef NDEBUG
579 if ((getFEMethod()->data_ctx & ctx).none()) {
580 MOFEM_LOG_CHANNEL("SELF");
581 MOFEM_LOG("SELF", Sev::error)
582 << "In this case filed degrees of freedom are read from vector. "
583 "That usually happens when time solver is used, and acces to "
584 "first or second rates is needed. You probably not set ts_u, "
585 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
586 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
587 "respectively";
588 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
589 }
590#endif
591 CHKERR VecGetArrayRead(vec, &array);
593 };
594
595 auto restore_array = [&](auto vec) {
596 return VecRestoreArrayRead(vec, &array);
597 };
598
599 switch (CTX) {
601 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
602 break;
604 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
605 break;
607 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
608 break;
609 default:
610 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
611 "That case is not implemented");
612 }
613
614 dotVector.resize(local_indices.size());
615 for (int i = 0; i != local_indices.size(); ++i)
616 if (local_indices[i] != -1)
617 dotVector[i] = array[local_indices[i]];
618 else
619 dotVector[i] = 0;
620
621 switch (CTX) {
623 CHKERR restore_array(getFEMethod()->ts_u);
624 break;
626 CHKERR restore_array(getFEMethod()->ts_u_t);
627 break;
629 CHKERR restore_array(getFEMethod()->ts_u_tt);
630 break;
631 default:
632 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
633 "That case is not implemented");
634 }
635
636 const size_t nb_base_functions = data.getN().size2();
637 auto base_function = data.getFTensor0N();
638 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
639
641 const size_t size = nb_dofs / Tensor_Dim;
642 if (nb_dofs % Tensor_Dim) {
643 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
644 }
645 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
646 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(dotVector);
647 size_t bb = 0;
648 for (; bb != size; ++bb) {
649 values_at_gauss_pts(I) += field_data(I) * base_function;
650 ++field_data;
651 ++base_function;
652 }
653 // Number of dofs can be smaller than number of Tensor_Dim x base
654 // functions
655 for (; bb != nb_base_functions; ++bb)
656 ++base_function;
657 ++values_at_gauss_pts;
658 }
660 }
661
662protected:
663 boost::shared_ptr<MatrixDouble> dataPtr;
666};
667
668/** \brief Get time derivatives of values at integration pts for tensor filed
669 * rank 1, i.e. vector field
670 *
671 * \ingroup mofem_forces_and_sources_user_data_operators
672 */
673template <int Tensor_Dim>
677
678/** \brief Get second time derivatives of values at integration pts for tensor
679 * filed rank 1, i.e. vector field
680 *
681 * \ingroup mofem_forces_and_sources_user_data_operators
682 */
683template <int Tensor_Dim>
687
688/**@}*/
689
690/** \name Tensor field values at integration points */
691
692/**@{*/
693
694/** \brief Calculate field values for tenor field rank 2.
695 *
696 */
697template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
700
702 const std::string field_name,
703 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
704 const EntityType zero_type = MBVERTEX)
707 dataPtr(data_ptr), zeroType(zero_type) {
708 if (!dataPtr)
709 THROW_MESSAGE("Pointer is not set");
710 }
711
712 MoFEMErrorCode doWork(int side, EntityType type,
714
715protected:
716 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
718};
719
720template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
723 doWork(int side, EntityType type, EntitiesFieldData::EntData &data) {
725 SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
726 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
727 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
728 Tensor_Dim0, Tensor_Dim1);
730}
731
732template <int Tensor_Dim0, int Tensor_Dim1>
733struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
734 ublas::row_major, DoubleAllocator>
736
738 const std::string field_name,
739 boost::shared_ptr<
740 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
741 data_ptr,
742 const EntityType zero_type = MBVERTEX)
745 dataPtr(data_ptr), zeroType(zero_type) {
746 if (!dataPtr)
747 THROW_MESSAGE("Pointer is not set");
748 }
749
751 const std::string field_name,
752 boost::shared_ptr<
753 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
754 data_ptr,
755 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
758 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
759 if (!dataPtr)
760 THROW_MESSAGE("Pointer is not set");
761 }
762
763 MoFEMErrorCode doWork(int side, EntityType type,
765
766protected:
767 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
772};
773
774template <int Tensor_Dim0, int Tensor_Dim1>
776 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
777 DoubleAllocator>::doWork(int side, EntityType type,
780
781 MatrixDouble &mat = *dataPtr;
782 const size_t nb_gauss_pts = data.getN().size1();
783 if (type == zeroType) {
784 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
785 mat.clear();
786 }
787
788 const size_t nb_dofs = data.getFieldData().size();
789
790 if (dataVec.use_count()) {
791 dotVector.resize(nb_dofs, false);
792 const double *array;
793 CHKERR VecGetArrayRead(dataVec, &array);
794 const auto &local_indices = data.getLocalIndices();
795 for (int i = 0; i != local_indices.size(); ++i)
796 if (local_indices[i] != -1)
797 dotVector[i] = array[local_indices[i]];
798 else
799 dotVector[i] = 0;
800 CHKERR VecRestoreArrayRead(dataVec, &array);
801 data.getFieldData().swap(dotVector);
802 }
803
804 if (nb_dofs) {
805 const size_t nb_base_functions = data.getN().size2();
806 auto base_function = data.getFTensor0N();
807 auto values_at_gauss_pts =
808 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
811 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
812 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
813 auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
814 size_t bb = 0;
815 for (; bb != size; ++bb) {
816 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
817 ++field_data;
818 ++base_function;
819 }
820 for (; bb != nb_base_functions; ++bb)
821 ++base_function;
822 ++values_at_gauss_pts;
823 }
824
825 if (dataVec.use_count()) {
826 data.getFieldData().swap(dotVector);
827 }
828 }
830}
831
832/** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
833 * field
834 *
835 * \ingroup mofem_forces_and_sources_user_data_operators
836 */
837template <int Tensor_Dim0, int Tensor_Dim1>
840 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
841
843 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
844 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
845};
846
847/** \brief Get time direvarive values at integration pts for tensor filed rank
848 * 2, i.e. matrix field
849 *
850 * \ingroup mofem_forces_and_sources_user_data_operators
851 */
852template <int Tensor_Dim0, int Tensor_Dim1>
855
857 boost::shared_ptr<MatrixDouble> data_ptr,
858 const EntityType zero_at_type = MBVERTEX)
861 dataPtr(data_ptr), zeroAtType(zero_at_type) {
862 if (!dataPtr)
863 THROW_MESSAGE("Pointer is not set");
864 }
865
869
870 const size_t nb_gauss_pts = getGaussPts().size2();
871 MatrixDouble &mat = *dataPtr;
872 if (type == zeroAtType) {
873 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
874 mat.clear();
875 }
876 const auto &local_indices = data.getLocalIndices();
877 const size_t nb_dofs = local_indices.size();
878 if (nb_dofs) {
879 dotVector.resize(nb_dofs, false);
880 const double *array;
881 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
882 for (size_t i = 0; i != local_indices.size(); ++i)
883 if (local_indices[i] != -1)
884 dotVector[i] = array[local_indices[i]];
885 else
886 dotVector[i] = 0;
887 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
888
889 const size_t nb_base_functions = data.getN().size2();
890
891 auto base_function = data.getFTensor0N();
892 auto values_at_gauss_pts =
893 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
896 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
897 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
898 auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
899 size_t bb = 0;
900 for (; bb != size; ++bb) {
901 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
902 ++field_data;
903 ++base_function;
904 }
905 for (; bb != nb_base_functions; ++bb)
906 ++base_function;
907 ++values_at_gauss_pts;
908 }
909 }
911 }
912
913protected:
914 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
915 EntityType zeroAtType; ///< Zero values at Gauss point at this type
916 VectorDouble dotVector; ///< Keeps temoorary values of time directives
917
918 template <int Dim0, int Dim1> auto getFTensorDotData() {
919 static_assert(Dim0 || !Dim0 || Dim1 || !Dim1, "not implemented");
920 }
921};
922
923template <>
924template <>
927 &dotVector[0], &dotVector[1], &dotVector[2],
928
929 &dotVector[3], &dotVector[4], &dotVector[5],
930
931 &dotVector[6], &dotVector[7], &dotVector[8]);
932}
933
934/**
935 * @brief Calculate symmetric tensor field values at integration pts.
936 *
937 * @tparam Tensor_Dim
938
939 * \ingroup mofem_forces_and_sources_user_data_operators
940 */
941template <int Tensor_Dim>
944
946 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
947 const EntityType zero_type = MBEDGE, const int zero_side = 0)
950 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
951 if (!dataPtr)
952 THROW_MESSAGE("Pointer is not set");
953 }
954
956 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
957 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBEDGE,
958 const int zero_side = 0)
961 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side),
962 dataVec(data_vec) {
963 if (!dataPtr)
964 THROW_MESSAGE("Pointer is not set");
965 }
966
970 MatrixDouble &mat = *dataPtr;
971 const int nb_gauss_pts = getGaussPts().size2();
972 if (type == this->zeroType && side == zeroSide) {
973 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
974 mat.clear();
975 }
976 const int nb_dofs = data.getFieldData().size();
977 if (!nb_dofs)
979
980 if (dataVec.use_count()) {
981 dotVector.resize(nb_dofs, false);
982 const double *array;
983 CHKERR VecGetArrayRead(dataVec, &array);
984 const auto &local_indices = data.getLocalIndices();
985 for (int i = 0; i != local_indices.size(); ++i)
986 if (local_indices[i] != -1)
987 dotVector[i] = array[local_indices[i]];
988 else
989 dotVector[i] = 0;
990 CHKERR VecRestoreArrayRead(dataVec, &array);
991 data.getFieldData().swap(dotVector);
992 }
993
994 const int nb_base_functions = data.getN().size2();
995 auto base_function = data.getFTensor0N();
996 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
999 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1000 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1001 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
1002 int bb = 0;
1003 for (; bb != size; ++bb) {
1004 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1005 ++field_data;
1006 ++base_function;
1007 }
1008 for (; bb != nb_base_functions; ++bb)
1009 ++base_function;
1010 ++values_at_gauss_pts;
1011 }
1012
1013 if (dataVec.use_count()) {
1014 data.getFieldData().swap(dotVector);
1015 }
1016
1018 }
1019
1020protected:
1021 boost::shared_ptr<MatrixDouble> dataPtr;
1023 const int zeroSide;
1026};
1027
1028/**
1029 * @brief Calculate symmetric tensor field rates ant integratio pts.
1030 *
1031 * @tparam Tensor_Dim
1032 *
1033 * \ingroup mofem_forces_and_sources_user_data_operators
1034 */
1035template <int Tensor_Dim>
1038
1040 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1041 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1044 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1045 if (!dataPtr)
1046 THROW_MESSAGE("Pointer is not set");
1047 }
1048
1052 const int nb_gauss_pts = getGaussPts().size2();
1053 MatrixDouble &mat = *dataPtr;
1054 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1055 if (type == zeroType && side == zeroSide) {
1056 mat.resize(symm_size, nb_gauss_pts, false);
1057 mat.clear();
1058 }
1059 auto &local_indices = data.getLocalIndices();
1060 const int nb_dofs = local_indices.size();
1061 if (!nb_dofs)
1063
1064#ifndef NDEBUG
1065 if ((getFEMethod()->data_ctx & PetscData::CtxSetX_T).none()) {
1066 MOFEM_LOG_CHANNEL("SELF");
1067 MOFEM_LOG("SELF", Sev::error)
1068 << "In this case filed degrees of freedom are read from vector. "
1069 "That usually happens when time solver is used, and acces to "
1070 "first rates is needed. You probably not set "
1071 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1072 "respectively";
1073 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
1074 }
1075#endif
1076
1077 dotVector.resize(nb_dofs, false);
1078 const double *array;
1079 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1080 for (int i = 0; i != local_indices.size(); ++i)
1081 if (local_indices[i] != -1)
1082 dotVector[i] = array[local_indices[i]];
1083 else
1084 dotVector[i] = 0;
1085 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1086
1087 const int nb_base_functions = data.getN().size2();
1088
1089 auto base_function = data.getFTensor0N();
1090 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1093 const int size = nb_dofs / symm_size;
1094 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1095 auto field_data = getFTensorDotData<Tensor_Dim>();
1096 int bb = 0;
1097 for (; bb != size; ++bb) {
1098 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1099 ++field_data;
1100 ++base_function;
1101 }
1102 for (; bb != nb_base_functions; ++bb)
1103 ++base_function;
1104 ++values_at_gauss_pts;
1105 }
1106
1108 }
1109
1110protected:
1111 boost::shared_ptr<MatrixDouble> dataPtr;
1113 const int zeroSide;
1115
1116 template <int Dim> inline auto getFTensorDotData() {
1117 static_assert(Dim || !Dim, "not implemented");
1118 }
1119};
1120
1121template <>
1122template <>
1123inline auto
1126 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1127 &dotVector[5]);
1128}
1129
1130template <>
1131template <>
1132inline auto
1135 &dotVector[0], &dotVector[1], &dotVector[2]);
1136}
1137
1138/**@}*/
1139
1140/** \name Gradients and Hessian of scalar fields at integration points */
1141
1142/**@{*/
1143
1144/**
1145 * \brief Evaluate field gradient values for scalar field, i.e. gradient is
1146 * tensor rank 1 (vector)
1147 *
1148 */
1149template <int Tensor_Dim, class T, class L, class A>
1151 : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
1152
1154 Tensor_Dim, T, L, A>::OpCalculateVectorFieldValues_General;
1155};
1156
1157/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1158 * tensor rank 1 (vector), specialization
1159 *
1160 */
1161template <int Tensor_Dim>
1163 ublas::row_major, DoubleAllocator>
1165 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1166
1168 Tensor_Dim, double, ublas::row_major,
1169 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1170
1171 /**
1172 * \brief calculate gradient values of scalar field at integration points
1173 * @param side side entity number
1174 * @param type side entity type
1175 * @param data entity data
1176 * @return error code
1177 */
1178 MoFEMErrorCode doWork(int side, EntityType type,
1180};
1181
1182/**
1183 * \brief Member function specialization calculating scalar field gradients for
1184 * tenor field rank 1
1185 *
1186 */
1187template <int Tensor_Dim>
1189 Tensor_Dim, double, ublas::row_major,
1190 DoubleAllocator>::doWork(int side, EntityType type,
1193
1194 const size_t nb_gauss_pts = this->getGaussPts().size2();
1195 auto &mat = *this->dataPtr;
1196 if (type == this->zeroType) {
1197 mat.resize(Tensor_Dim, nb_gauss_pts, false);
1198 mat.clear();
1199 }
1200
1201 const int nb_dofs = data.getFieldData().size();
1202 if (nb_dofs) {
1203
1204 if (nb_gauss_pts) {
1205 const int nb_base_functions = data.getN().size2();
1206 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
1207 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1208
1209#ifndef NDEBUG
1210 if (nb_dofs > nb_base_functions)
1211 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1212 "Number of base functions inconsistent with number of DOFs "
1213 "(%d > %d)",
1214 nb_dofs, nb_base_functions);
1215
1216 if (data.getDiffN().size2() != nb_base_functions * Tensor_Dim)
1217 SETERRQ2(
1218 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1219 "Number of base functions inconsistent with number of derivatives "
1220 "(%d != %d)",
1221 data.getDiffN().size2(), nb_base_functions);
1222
1223 if (data.getDiffN().size1() != nb_gauss_pts)
1224 SETERRQ2(
1225 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1226 "Number of base functions inconsistent with number of integration "
1227 "pts (%d != %d)",
1228 data.getDiffN().size2(), nb_gauss_pts);
1229
1230#endif
1231
1233 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1234 auto field_data = data.getFTensor0FieldData();
1235 int bb = 0;
1236 for (; bb != nb_dofs; ++bb) {
1237 gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
1238 ++field_data;
1239 ++diff_base_function;
1240 }
1241 // Number of dofs can be smaller than number of base functions
1242 for (; bb < nb_base_functions; ++bb)
1243 ++diff_base_function;
1244 ++gradients_at_gauss_pts;
1245 }
1246 }
1247 }
1248
1250}
1251
1252/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1253 * vector field
1254 *
1255 * \ingroup mofem_forces_and_sources_user_data_operators
1256 */
1257template <int Tensor_Dim>
1260 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1262 Tensor_Dim, double, ublas::row_major,
1263 DoubleAllocator>::OpCalculateScalarFieldGradient_General;
1264};
1265
1266/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1267 * tensor rank 1 (vector), specialization
1268 *
1269 */
1270template <int Tensor_Dim>
1273 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1274
1276 Tensor_Dim, double, ublas::row_major,
1277 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1278
1279 /**
1280 * \brief calculate gradient values of scalar field at integration points
1281 * @param side side entity number
1282 * @param type side entity type
1283 * @param data entity data
1284 * @return error code
1285 */
1286 MoFEMErrorCode doWork(int side, EntityType type,
1288};
1289
1290template <int Tensor_Dim>
1292 int side, EntityType type, EntitiesFieldData::EntData &data) {
1294
1295 const size_t nb_gauss_pts = this->getGaussPts().size2();
1296 auto &mat = *this->dataPtr;
1297 if (type == this->zeroType) {
1298 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts, false);
1299 mat.clear();
1300 }
1301
1302 const int nb_dofs = data.getFieldData().size();
1303 if (nb_dofs) {
1304
1305 if (nb_gauss_pts) {
1306 const int nb_base_functions = data.getN().size2();
1307
1308 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1309#ifndef NDEBUG
1310 if (hessian_base.size1() != nb_gauss_pts) {
1311 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1312 "Wrong number of integration pts (%d != %d)",
1313 hessian_base.size1(), nb_gauss_pts);
1314 }
1315 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1316 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1317 "Wrong number of base functions (%d != %d)",
1318 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1319 nb_base_functions);
1320 }
1321 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1322 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1323 "Wrong number of base functions (%d < %d)",
1324 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1325 }
1326#endif
1327
1328 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1329 &*hessian_base.data().begin());
1330
1331 auto hessian_at_gauss_pts =
1332 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1333
1336 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1337 auto field_data = data.getFTensor0FieldData();
1338 int bb = 0;
1339 for (; bb != nb_dofs; ++bb) {
1340 hessian_at_gauss_pts(I, J) +=
1341 field_data * t_diff2_base_function(I, J);
1342 ++field_data;
1343 ++t_diff2_base_function;
1344 }
1345 // Number of dofs can be smaller than number of base functions
1346 for (; bb < nb_base_functions; ++bb) {
1347 ++t_diff2_base_function;
1348 }
1349
1350 ++hessian_at_gauss_pts;
1351 }
1352 }
1353 }
1354
1356}
1357
1358/**}*/
1359
1360/** \name Gradients and hessian of tensor fields at integration points */
1361
1362/**@{*/
1363
1364/**
1365 * \brief Evaluate field gradient values for vector field, i.e. gradient is
1366 * tensor rank 2
1367 *
1368 */
1369template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1371 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1372 L, A> {
1373
1375 Tensor_Dim0, Tensor_Dim1, T, L, A>::OpCalculateTensor2FieldValues_General;
1376};
1377
1378template <int Tensor_Dim0, int Tensor_Dim1>
1379struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1380 ublas::row_major, DoubleAllocator>
1382 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1383
1385 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1386 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1387
1388 /**
1389 * \brief calculate values of vector field at integration points
1390 * @param side side entity number
1391 * @param type side entity type
1392 * @param data entity data
1393 * @return error code
1394 */
1395 MoFEMErrorCode doWork(int side, EntityType type,
1397};
1398
1399/**
1400 * \brief Member function specialization calculating vector field gradients for
1401 * tenor field rank 2
1402 *
1403 */
1404template <int Tensor_Dim0, int Tensor_Dim1>
1406 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1407 DoubleAllocator>::doWork(int side, EntityType type,
1410 if (!this->dataPtr)
1411 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1412 "Data pointer not allocated");
1413
1414 const size_t nb_gauss_pts = this->getGaussPts().size2();
1415 auto &mat = *this->dataPtr;
1416 if (type == this->zeroType) {
1417 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1418 mat.clear();
1419 }
1420
1421 if (nb_gauss_pts) {
1422 const size_t nb_dofs = data.getFieldData().size();
1423
1424 if (nb_dofs) {
1425
1426 if (this->dataVec.use_count()) {
1427 this->dotVector.resize(nb_dofs, false);
1428 const double *array;
1429 CHKERR VecGetArrayRead(this->dataVec, &array);
1430 const auto &local_indices = data.getLocalIndices();
1431 for (int i = 0; i != local_indices.size(); ++i)
1432 if (local_indices[i] != -1)
1433 this->dotVector[i] = array[local_indices[i]];
1434 else
1435 this->dotVector[i] = 0;
1436 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1437 data.getFieldData().swap(this->dotVector);
1438 }
1439
1440 const int nb_base_functions = data.getN().size2();
1441 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1442 auto gradients_at_gauss_pts =
1443 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1446 int size = nb_dofs / Tensor_Dim0;
1447 if (nb_dofs % Tensor_Dim0) {
1448 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1449 "Data inconsistency");
1450 }
1451 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1452 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1453 int bb = 0;
1454 for (; bb < size; ++bb) {
1455 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1456 ++field_data;
1457 ++diff_base_function;
1458 }
1459 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1460 // functions
1461 for (; bb != nb_base_functions; ++bb)
1462 ++diff_base_function;
1463 ++gradients_at_gauss_pts;
1464 }
1465
1466 if (this->dataVec.use_count()) {
1467 data.getFieldData().swap(this->dotVector);
1468 }
1469 }
1470 }
1472}
1473
1474/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1475 * vector field
1476 *
1477 * \ingroup mofem_forces_and_sources_user_data_operators
1478 */
1479template <int Tensor_Dim0, int Tensor_Dim1>
1482 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1483
1485 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1486 DoubleAllocator>::OpCalculateVectorFieldGradient_General;
1487};
1488
1489/** \brief Get field gradients time derivative at integration pts for scalar
1490 * filed rank 0, i.e. vector field
1491 *
1492 * \ingroup mofem_forces_and_sources_user_data_operators
1493 */
1494template <int Tensor_Dim0, int Tensor_Dim1>
1497
1499 boost::shared_ptr<MatrixDouble> data_ptr,
1500 const EntityType zero_at_type = MBVERTEX)
1503 dataPtr(data_ptr), zeroAtType(zero_at_type) {
1504 if (!dataPtr)
1505 THROW_MESSAGE("Pointer is not set");
1506 }
1507
1511
1512 const auto &local_indices = data.getLocalIndices();
1513 const int nb_dofs = local_indices.size();
1514 const int nb_gauss_pts = this->getGaussPts().size2();
1515
1516 auto &mat = *this->dataPtr;
1517 if (type == this->zeroAtType) {
1518 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1519 mat.clear();
1520 }
1521 if (!nb_dofs)
1523
1524 dotVector.resize(nb_dofs, false);
1525 const double *array;
1526 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1527 for (int i = 0; i != local_indices.size(); ++i)
1528 if (local_indices[i] != -1)
1529 dotVector[i] = array[local_indices[i]];
1530 else
1531 dotVector[i] = 0;
1532 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1533
1534 const int nb_base_functions = data.getN().size2();
1535 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1536 auto gradients_at_gauss_pts =
1537 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1540 int size = nb_dofs / Tensor_Dim0;
1541 if (nb_dofs % Tensor_Dim0) {
1542 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1543 }
1544
1545 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1546 auto field_data = getFTensorDotData<Tensor_Dim0>();
1547 int bb = 0;
1548 for (; bb < size; ++bb) {
1549 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1550 ++field_data;
1551 ++diff_base_function;
1552 }
1553 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1554 // functions
1555 for (; bb != nb_base_functions; ++bb)
1556 ++diff_base_function;
1557 ++gradients_at_gauss_pts;
1558 }
1560 }
1561
1562private:
1563 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1564 EntityType zeroAtType; ///< Zero values at Gauss point at this type
1565 VectorDouble dotVector; ///< Keeps temoorary values of time directives
1566
1567 template <int Dim> inline auto getFTensorDotData() {
1568 static_assert(Dim || !Dim, "not implemented");
1569 }
1570};
1571
1572template <>
1573template <>
1576 &dotVector[0], &dotVector[1], &dotVector[2]);
1577}
1578
1579template <>
1580template <>
1582 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(&dotVector[0],
1583 &dotVector[1]);
1584}
1585
1586/**
1587 * \brief Evaluate field gradient values for symmetric 2nd order tensor field,
1588 * i.e. gradient is tensor rank 3
1589 *
1590 */
1591template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1593 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1594 L, A> {
1595
1597 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1598 const EntityType zero_type = MBVERTEX)
1599 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
1600 A>(field_name, data_ptr,
1601 zero_type) {}
1602};
1603
1604template <int Tensor_Dim0, int Tensor_Dim1>
1606 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator>
1608 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1609
1611 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1612 const EntityType zero_type = MBVERTEX)
1613 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1614 ublas::row_major,
1616 field_name, data_ptr, zero_type) {}
1617
1618 /**
1619 * \brief calculate values of vector field at integration points
1620 * @param side side entity number
1621 * @param type side entity type
1622 * @param data entity data
1623 * @return error code
1624 */
1625 MoFEMErrorCode doWork(int side, EntityType type,
1627};
1628
1629/**
1630 * \brief Member function specialization calculating tensor field gradients for
1631 * symmetric tensor field rank 2
1632 *
1633 */
1634template <int Tensor_Dim0, int Tensor_Dim1>
1635MoFEMErrorCode OpCalculateTensor2SymmetricFieldGradient_General<
1636 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1637 DoubleAllocator>::doWork(int side, EntityType type,
1640 if (!this->dataPtr)
1641 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1642 "Data pointer not allocated");
1643
1644 const size_t nb_gauss_pts = this->getGaussPts().size2();
1645 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1646 auto &mat = *this->dataPtr;
1647 if (type == this->zeroType) {
1648 mat.resize(msize * Tensor_Dim1, nb_gauss_pts, false);
1649 mat.clear();
1650 }
1651
1652 if (nb_gauss_pts) {
1653 const size_t nb_dofs = data.getFieldData().size();
1654
1655 if (nb_dofs) {
1656
1657 const int nb_base_functions = data.getN().size2();
1658 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1659 auto gradients_at_gauss_pts =
1660 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1664 int size = nb_dofs / msize;
1665 if (nb_dofs % msize) {
1666 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1667 "Data inconsistency");
1668 }
1669 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1670 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim0>();
1671 int bb = 0;
1672 for (; bb < size; ++bb) {
1673 gradients_at_gauss_pts(I, J, K) +=
1674 field_data(I, J) * diff_base_function(K);
1675 ++field_data;
1676 ++diff_base_function;
1677 }
1678 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1679 // functions
1680 for (; bb != nb_base_functions; ++bb)
1681 ++diff_base_function;
1682 ++gradients_at_gauss_pts;
1683 }
1684 }
1685 }
1687}
1688
1689/** \brief Get field gradients at integration pts for symmetric tensorial field
1690 * rank 2
1691 *
1692 * \ingroup mofem_forces_and_sources_user_data_operators
1693 */
1694template <int Tensor_Dim0, int Tensor_Dim1>
1697 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1698
1700 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1701 const EntityType zero_type = MBVERTEX)
1703 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1704 DoubleAllocator>(field_name, data_ptr, zero_type) {}
1705};
1706
1707template <int Tensor_Dim0, int Tensor_Dim1>
1710 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1711
1713 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1714 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1715
1716 /**
1717 * \brief calculate values of vector field at integration points
1718 * @param side side entity number
1719 * @param type side entity type
1720 * @param data entity data
1721 * @return error code
1722 */
1723 MoFEMErrorCode doWork(int side, EntityType type,
1725};
1726
1727/**
1728 * \brief Member function specialization calculating vector field gradients for
1729 * tenor field rank 2
1730 *
1731 */
1732template <int Tensor_Dim0, int Tensor_Dim1>
1734 int side, EntityType type, EntitiesFieldData::EntData &data) {
1736 if (!this->dataPtr)
1737 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1738 "Data pointer not allocated");
1739
1740 const size_t nb_gauss_pts = this->getGaussPts().size2();
1741 auto &mat = *this->dataPtr;
1742 if (type == this->zeroType) {
1743 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts, false);
1744 mat.clear();
1745 }
1746
1747 if (nb_gauss_pts) {
1748 const size_t nb_dofs = data.getFieldData().size();
1749
1750 if (nb_dofs) {
1751
1752 if (this->dataVec.use_count()) {
1753 this->dotVector.resize(nb_dofs, false);
1754 const double *array;
1755 CHKERR VecGetArrayRead(this->dataVec, &array);
1756 const auto &local_indices = data.getLocalIndices();
1757 for (int i = 0; i != local_indices.size(); ++i)
1758 if (local_indices[i] != -1)
1759 this->dotVector[i] = array[local_indices[i]];
1760 else
1761 this->dotVector[i] = 0;
1762 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1763 data.getFieldData().swap(this->dotVector);
1764 }
1765
1766 const int nb_base_functions = data.getN().size2();
1767
1768 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1769#ifndef NDEBUG
1770 if (hessian_base.size1() != nb_gauss_pts) {
1771 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1772 "Wrong number of integration pts (%d != %d)",
1773 hessian_base.size1(), nb_gauss_pts);
1774 }
1775 if (hessian_base.size2() !=
1776 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1777 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1778 "Wrong number of base functions (%d != %d)",
1779 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1780 nb_base_functions);
1781 }
1782 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1783 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1784 "Wrong number of base functions (%d < %d)",
1785 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1786 }
1787#endif
1788
1789 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1790 &*hessian_base.data().begin());
1791
1792 auto t_hessian_at_gauss_pts =
1793 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1794
1798
1799 int size = nb_dofs / Tensor_Dim0;
1800#ifndef NDEBUG
1801 if (nb_dofs % Tensor_Dim0) {
1802 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1803 "Data inconsistency");
1804 }
1805#endif
1806
1807 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1808 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1809 int bb = 0;
1810 for (; bb < size; ++bb) {
1811 t_hessian_at_gauss_pts(I, J, K) +=
1812 field_data(I) * t_diff2_base_function(J, K);
1813 ++field_data;
1814 ++t_diff2_base_function;
1815 }
1816 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1817 // functions
1818 for (; bb != nb_base_functions; ++bb)
1819 ++t_diff2_base_function;
1820 ++t_hessian_at_gauss_pts;
1821 }
1822
1823 if (this->dataVec.use_count()) {
1824 data.getFieldData().swap(this->dotVector);
1825 }
1826 }
1827 }
1829}
1830
1831/**@}*/
1832
1833/** \name Transform tensors and vectors */
1834
1835/**@{*/
1836
1837/**
1838 * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1839 * \f$
1840 *
1841 * @tparam DIM
1842 *
1843 * \ingroup mofem_forces_and_sources_user_data_operators
1844 */
1845template <int DIM_01, int DIM_23, int S = 0>
1848
1851
1853 boost::shared_ptr<MatrixDouble> in_mat,
1854 boost::shared_ptr<MatrixDouble> out_mat,
1855 boost::shared_ptr<MatrixDouble> d_mat)
1856 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1857 // Only is run for vertices
1858 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1859 if (!inMat)
1860 THROW_MESSAGE("Pointer for in mat is null");
1861 if (!outMat)
1862 THROW_MESSAGE("Pointer for out mat is null");
1863 if (!dMat)
1864 THROW_MESSAGE("Pointer for tensor mat is null");
1865 }
1866
1867 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1869 const size_t nb_gauss_pts = getGaussPts().size2();
1870 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
1871 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
1872 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1873 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
1874 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1875 t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1876 ++t_in;
1877 ++t_out;
1878 }
1880 }
1881
1882private:
1887
1888 boost::shared_ptr<MatrixDouble> inMat;
1889 boost::shared_ptr<MatrixDouble> outMat;
1890 boost::shared_ptr<MatrixDouble> dMat;
1891};
1892
1893template <int DIM>
1895
1898
1900 boost::shared_ptr<MatrixDouble> in_mat,
1901 boost::shared_ptr<MatrixDouble> out_mat)
1902 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1903 // Only is run for vertices
1904 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1905 if (!inMat)
1906 THROW_MESSAGE("Pointer not set for in matrix");
1907 if (!outMat)
1908 THROW_MESSAGE("Pointer not set for in matrix");
1909 }
1910
1911 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1913 const size_t nb_gauss_pts = getGaussPts().size2();
1914 auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
1915 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
1916 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
1917 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1918 t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
1919 ++t_in;
1920 ++t_out;
1921 }
1923 }
1924
1925private:
1928 boost::shared_ptr<MatrixDouble> inMat;
1929 boost::shared_ptr<MatrixDouble> outMat;
1930};
1931
1933
1936
1937 OpScaleMatrix(const std::string field_name, const double scale,
1938 boost::shared_ptr<MatrixDouble> in_mat,
1939 boost::shared_ptr<MatrixDouble> out_mat)
1940 : UserOp(field_name, OPROW), scale(scale), inMat(in_mat),
1941 outMat(out_mat) {
1942 // Only is run for vertices
1943 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1944 if (!inMat)
1945 THROW_MESSAGE("Pointer not set for in matrix");
1946 if (!outMat)
1947 THROW_MESSAGE("Pointer not set for in matrix");
1948 }
1949
1950 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1952 outMat->resize(inMat->size1(), inMat->size2(), false);
1953 noalias(*outMat) = scale * (*inMat);
1955 }
1956
1957private:
1958 const double scale;
1959 boost::shared_ptr<MatrixDouble> inMat;
1960 boost::shared_ptr<MatrixDouble> outMat;
1961};
1962
1963/**@}*/
1964
1965/** \name H-div/H-curls (Vectorial bases) values at integration points */
1966
1967/**@{*/
1968
1969/** \brief Get vector field for H-div approximation
1970 * \ingroup mofem_forces_and_sources_user_data_operators
1971 */
1972template <int Base_Dim, int Field_Dim, class T, class L, class A>
1974
1975/** \brief Get vector field for H-div approximation
1976 * \ingroup mofem_forces_and_sources_user_data_operators
1977 */
1978template <int Field_Dim>
1980 ublas::row_major, DoubleAllocator>
1982
1984 boost::shared_ptr<MatrixDouble> data_ptr,
1985 const EntityType zero_type = MBEDGE,
1986 const int zero_side = 0)
1989 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1990 if (!dataPtr)
1991 THROW_MESSAGE("Pointer is not set");
1992 }
1993
1994 /**
1995 * \brief Calculate values of vector field at integration points
1996 * @param side side entity number
1997 * @param type side entity type
1998 * @param data entity data
1999 * @return error code
2000 */
2001 MoFEMErrorCode doWork(int side, EntityType type,
2003
2004private:
2005 boost::shared_ptr<MatrixDouble> dataPtr;
2007 const int zeroSide;
2008};
2009
2010template <int Field_Dim>
2012 3, Field_Dim, double, ublas::row_major,
2013 DoubleAllocator>::doWork(int side, EntityType type,
2016 const size_t nb_integration_points = this->getGaussPts().size2();
2017 if (type == zeroType && side == zeroSide) {
2018 dataPtr->resize(Field_Dim, nb_integration_points, false);
2019 dataPtr->clear();
2020 }
2021 const size_t nb_dofs = data.getFieldData().size();
2022 if (!nb_dofs)
2024 const size_t nb_base_functions = data.getN().size2() / 3;
2026 auto t_n_hdiv = data.getFTensor1N<3>();
2027 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2028 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2029 auto t_dof = data.getFTensor0FieldData();
2030 int bb = 0;
2031 for (; bb != nb_dofs; ++bb) {
2032 t_data(i) += t_n_hdiv(i) * t_dof;
2033 ++t_n_hdiv;
2034 ++t_dof;
2035 }
2036 for (; bb != nb_base_functions; ++bb)
2037 ++t_n_hdiv;
2038 ++t_data;
2039 }
2041}
2042
2043/** \brief Get vector field for H-div approximation
2044 *
2045 * \ingroup mofem_forces_and_sources_user_data_operators
2046 */
2047template <int Base_Dim, int Field_Dim = Base_Dim>
2050 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2052 Base_Dim, Field_Dim, double, ublas::row_major,
2053 DoubleAllocator>::OpCalculateHVecVectorField_General;
2054};
2055
2056
2057
2058/** \brief Get vector field for H-div approximation
2059 * \ingroup mofem_forces_and_sources_user_data_operators
2060 */
2061template <int Base_Dim, int Field_Dim = Base_Dim>
2063
2064template <int Field_Dim>
2067
2069 boost::shared_ptr<MatrixDouble> data_ptr,
2070 const EntityType zero_type = MBEDGE,
2071 const int zero_side = 0)
2074 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2075 if (!dataPtr)
2076 THROW_MESSAGE("Pointer is not set");
2077 }
2078
2079 /**
2080 * \brief Calculate values of vector field at integration points
2081 * @param side side entity number
2082 * @param type side entity type
2083 * @param data entity data
2084 * @return error code
2085 */
2086 MoFEMErrorCode doWork(int side, EntityType type,
2088
2089private:
2090 boost::shared_ptr<MatrixDouble> dataPtr;
2092 const int zeroSide;
2093};
2094
2095template <int Field_Dim>
2097 int side, EntityType type, EntitiesFieldData::EntData &data) {
2099
2100 const size_t nb_integration_points = this->getGaussPts().size2();
2101 if (type == zeroType && side == zeroSide) {
2102 dataPtr->resize(Field_Dim, nb_integration_points, false);
2103 dataPtr->clear();
2104 }
2105
2106 auto &local_indices = data.getIndices();
2107 const size_t nb_dofs = local_indices.size();
2108 if (nb_dofs) {
2109
2110 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2111 const double *array;
2112 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2113 for (size_t i = 0; i != nb_dofs; ++i)
2114 if (local_indices[i] != -1)
2115 dot_dofs_vector[i] = array[local_indices[i]];
2116 else
2117 dot_dofs_vector[i] = 0;
2118 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2119
2120 const size_t nb_base_functions = data.getN().size2() / 3;
2122 auto t_n_hdiv = data.getFTensor1N<3>();
2123 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2124 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2125 int bb = 0;
2126 for (; bb != nb_dofs; ++bb) {
2127 t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
2128 ++t_n_hdiv;
2129 }
2130 for (; bb != nb_base_functions; ++bb)
2131 ++t_n_hdiv;
2132 ++t_data;
2133 }
2134 }
2135
2137}
2138
2139/**
2140 * @brief Calculate divergence of vector field
2141 * @ingroup mofem_forces_and_sources_user_data_operators
2142 *
2143 * @tparam BASE_DIM
2144 * @tparam SPACE_DIM
2145 */
2146template <int BASE_DIM, int SPACE_DIM>
2149
2151 boost::shared_ptr<VectorDouble> data_ptr,
2152 const EntityType zero_type = MBEDGE,
2153 const int zero_side = 0)
2156 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2157 if (!dataPtr)
2158 THROW_MESSAGE("Pointer is not set");
2159 }
2160
2164 const size_t nb_integration_points = getGaussPts().size2();
2165 if (type == zeroType && side == zeroSide) {
2166 dataPtr->resize(nb_integration_points, false);
2167 dataPtr->clear();
2168 }
2169 const size_t nb_dofs = data.getFieldData().size();
2170 if (!nb_dofs)
2172 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2175 auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2176 auto t_data = getFTensor0FromVec(*dataPtr);
2177 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2178 auto t_dof = data.getFTensor0FieldData();
2179 int bb = 0;
2180 for (; bb != nb_dofs; ++bb) {
2181 t_data += t_dof * t_n_diff_hdiv(j, j);
2182 ++t_n_diff_hdiv;
2183 ++t_dof;
2184 }
2185 for (; bb != nb_base_functions; ++bb)
2186 ++t_n_diff_hdiv;
2187 ++t_data;
2188 }
2190 }
2191
2192private:
2193 boost::shared_ptr<VectorDouble> dataPtr;
2195 const int zeroSide;
2196};
2197
2198/**
2199 * @brief Calculate gradient of vector field
2200 * @ingroup mofem_forces_and_sources_user_data_operators
2201 *
2202 * @tparam BASE_DIM
2203 * @tparam SPACE_DIM
2204 */
2205template <int BASE_DIM, int SPACE_DIM>
2208
2210 boost::shared_ptr<MatrixDouble> data_ptr,
2211 const EntityType zero_type = MBEDGE,
2212 const int zero_side = 0)
2215 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2216 if (!dataPtr)
2217 THROW_MESSAGE("Pointer is not set");
2218 }
2219
2223 const size_t nb_integration_points = getGaussPts().size2();
2224 if (type == zeroType && side == zeroSide) {
2225 dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
2226 dataPtr->clear();
2227 }
2228 const size_t nb_dofs = data.getFieldData().size();
2229 if (!nb_dofs)
2231 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2234 auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2235 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*dataPtr);
2236 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2237 auto t_dof = data.getFTensor0FieldData();
2238 int bb = 0;
2239 for (; bb != nb_dofs; ++bb) {
2240 t_data(i, j) += t_dof * t_base_diff(i, j);
2241 ++t_base_diff;
2242 ++t_dof;
2243 }
2244 for (; bb != nb_base_functions; ++bb)
2245 ++t_base_diff;
2246 ++t_data;
2247 }
2249 }
2250
2251private:
2252 boost::shared_ptr<MatrixDouble> dataPtr;
2254 const int zeroSide;
2255};
2256
2257/**
2258 * @brief Calculate gradient of vector field
2259 * @ingroup mofem_forces_and_sources_user_data_operators
2260 *
2261 * @tparam BASE_DIM
2262 * @tparam SPACE_DIM
2263 */
2264template <int BASE_DIM, int SPACE_DIM>
2267
2269 boost::shared_ptr<MatrixDouble> data_ptr,
2270 const EntityType zero_type = MBEDGE,
2271 const int zero_side = 0)
2274 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2275 if (!dataPtr)
2276 THROW_MESSAGE("Pointer is not set");
2277 }
2278
2282 const size_t nb_integration_points = getGaussPts().size2();
2283 if (type == zeroType && side == zeroSide) {
2284 dataPtr->resize(BASE_DIM * SPACE_DIM * SPACE_DIM, nb_integration_points,
2285 false);
2286 dataPtr->clear();
2287 }
2288 const size_t nb_dofs = data.getFieldData().size();
2289 if (!nb_dofs)
2291
2292 const int nb_base_functions = data.getN().size2() / BASE_DIM;
2293 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2294
2295#ifndef NDEBUG
2296 if (hessian_base.size1() != nb_integration_points) {
2297 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2298 "Wrong number of integration pts (%d != %d)",
2299 hessian_base.size1(), nb_integration_points);
2300 }
2301 if (hessian_base.size2() !=
2302 BASE_DIM * nb_base_functions * SPACE_DIM * SPACE_DIM) {
2303 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2304 "Wrong number of base functions (%d != %d)",
2305 hessian_base.size2() / (BASE_DIM * SPACE_DIM * SPACE_DIM),
2306 nb_base_functions);
2307 }
2308 if (hessian_base.size2() < BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM) {
2309 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2310 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2311 BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM);
2312 }
2313#endif
2314
2318
2319 auto t_base_diff2 =
2321 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*dataPtr);
2322 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2323 auto t_dof = data.getFTensor0FieldData();
2324 int bb = 0;
2325 for (; bb != nb_dofs; ++bb) {
2326 t_data(i, j, k) += t_dof * t_base_diff2(i, j, k);
2327
2328 ++t_base_diff2;
2329 ++t_dof;
2330 }
2331 for (; bb != nb_base_functions; ++bb)
2332 ++t_base_diff2;
2333 ++t_data;
2334 }
2336 }
2337
2338private:
2339 boost::shared_ptr<MatrixDouble> dataPtr;
2341 const int zeroSide;
2342};
2343
2344/**
2345 * @brief Calculate divergence of vector field dot
2346 * @ingroup mofem_forces_and_sources_user_data_operators
2347 *
2348 * @tparam Tensor_Dim dimension of space
2349 */
2350template <int Tensor_Dim1, int Tensor_Dim2>
2353
2355 boost::shared_ptr<VectorDouble> data_ptr,
2356 const EntityType zero_type = MBEDGE,
2357 const int zero_side = 0)
2360 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2361 if (!dataPtr)
2362 THROW_MESSAGE("Pointer is not set");
2363 }
2364
2368 const size_t nb_integration_points = getGaussPts().size2();
2369 if (type == zeroType && side == zeroSide) {
2370 dataPtr->resize(nb_integration_points, false);
2371 dataPtr->clear();
2372 }
2373
2374 const auto &local_indices = data.getLocalIndices();
2375 const int nb_dofs = local_indices.size();
2376 if (nb_dofs) {
2377
2378 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2379 const double *array;
2380 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2381 for (size_t i = 0; i != local_indices.size(); ++i)
2382 if (local_indices[i] != -1)
2383 dot_dofs_vector[i] = array[local_indices[i]];
2384 else
2385 dot_dofs_vector[i] = 0;
2386 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2387
2388 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
2390 auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
2391 auto t_data = getFTensor0FromVec(*dataPtr);
2392 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2393 int bb = 0;
2394 for (; bb != nb_dofs; ++bb) {
2395 double div = 0;
2396 for (auto ii = 0; ii != Tensor_Dim2; ++ii)
2397 div += t_n_diff_hdiv(ii, ii);
2398 t_data += dot_dofs_vector[bb] * div;
2399 ++t_n_diff_hdiv;
2400 }
2401 for (; bb != nb_base_functions; ++bb)
2402 ++t_n_diff_hdiv;
2403 ++t_data;
2404 }
2405 }
2407 }
2408
2409private:
2410 boost::shared_ptr<VectorDouble> dataPtr;
2412 const int zeroSide;
2413};
2414
2415/**
2416 * @brief Calculate curl of vector field
2417 * @ingroup mofem_forces_and_sources_user_data_operators
2418 *
2419 * @tparam Tensor_Dim dimension of space
2420 */
2421template <int Tensor_Dim>
2424
2426 boost::shared_ptr<MatrixDouble> data_ptr,
2427 const EntityType zero_type = MBEDGE,
2428 const int zero_side = 0)
2431 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2432 if (!dataPtr)
2433 THROW_MESSAGE("Pointer is not set");
2434 }
2435
2439 const size_t nb_integration_points = getGaussPts().size2();
2440 if (type == zeroType && side == zeroSide) {
2441 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2442 dataPtr->clear();
2443 }
2444 const size_t nb_dofs = data.getFieldData().size();
2445 if (!nb_dofs)
2447
2448 MatrixDouble curl_mat;
2449
2450 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim;
2454 auto t_n_diff_hcurl = data.getFTensor2DiffN<Tensor_Dim, Tensor_Dim>();
2455 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2456 for (int gg = 0; gg != nb_integration_points; ++gg) {
2457
2458 auto t_dof = data.getFTensor0FieldData();
2459 int bb = 0;
2460 for (; bb != nb_dofs; ++bb) {
2461
2462 t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
2463 ++t_n_diff_hcurl;
2464 ++t_dof;
2465 }
2466
2467 for (; bb != nb_base_functions; ++bb)
2468 ++t_n_diff_hcurl;
2469 ++t_data;
2470 }
2471
2473 }
2474
2475private:
2476 boost::shared_ptr<MatrixDouble> dataPtr;
2478 const int zeroSide;
2479};
2480
2481/**
2482 * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
2483 * \ingroup mofem_forces_and_sources_user_data_operators
2484 *
2485 * @tparam Tensor_Dim0 rank of the filed
2486 * @tparam Tensor_Dim1 dimension of space
2487 */
2488template <int Tensor_Dim0, int Tensor_Dim1>
2491
2493 boost::shared_ptr<MatrixDouble> data_ptr,
2494 const EntityType zero_type = MBEDGE,
2495 const int zero_side = 0)
2498 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2499 if (!dataPtr)
2500 THROW_MESSAGE("Pointer is not set");
2501 }
2502
2506 const size_t nb_integration_points = getGaussPts().size2();
2507 if (type == zeroType && side == zeroSide) {
2508 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2509 dataPtr->clear();
2510 }
2511 const size_t nb_dofs = data.getFieldData().size();
2512 if (nb_dofs) {
2513 const size_t nb_base_functions = data.getN().size2() / 3;
2516 auto t_n_hvec = data.getFTensor1N<3>();
2517 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2518 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2519 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2520 size_t bb = 0;
2521 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2522 t_data(i, j) += t_dof(i) * t_n_hvec(j);
2523 ++t_n_hvec;
2524 ++t_dof;
2525 }
2526 for (; bb != nb_base_functions; ++bb)
2527 ++t_n_hvec;
2528 ++t_data;
2529 }
2530 }
2532 }
2533
2534private:
2535 boost::shared_ptr<MatrixDouble> dataPtr;
2537 const int zeroSide;
2538};
2539
2540/**
2541 * @brief Calculate tenor field using tensor base, i.e. Hdiv/Hcurl
2542 * \ingroup mofem_forces_and_sources_user_data_operators
2543 *
2544 * @tparam Tensor_Dim0 rank of the filed
2545 * @tparam Tensor_Dim1 dimension of space
2546 */
2547template <int Tensor_Dim0, int Tensor_Dim1>
2550
2552 boost::shared_ptr<MatrixDouble> data_ptr,
2553 const EntityType zero_type = MBEDGE,
2554 const int zero_side = 0)
2557 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2558 if (!dataPtr)
2559 THROW_MESSAGE("Pointer is not set");
2560 }
2561
2565 const size_t nb_integration_points = getGaussPts().size2();
2566 if (type == zeroType && side == zeroSide) {
2567 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2568 dataPtr->clear();
2569 }
2570 const size_t nb_dofs = data.getFieldData().size();
2571 if (!nb_dofs)
2573 const size_t nb_base_functions =
2574 data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2577 auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2578 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2579 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2580 auto t_dof = data.getFTensor0FieldData();
2581 size_t bb = 0;
2582 for (; bb != nb_dofs; ++bb) {
2583 t_data(i, j) += t_dof * t_n_hten(i, j);
2584 ++t_n_hten;
2585 ++t_dof;
2586 }
2587 for (; bb != nb_base_functions; ++bb)
2588 ++t_n_hten;
2589 ++t_data;
2590 }
2592 }
2593
2594private:
2595 boost::shared_ptr<MatrixDouble> dataPtr;
2597 const int zeroSide;
2598};
2599
2600/**
2601 * @brief Calculate divergence of tonsorial field using vectorial base
2602 * \ingroup mofem_forces_and_sources_user_data_operators
2603 *
2604 * @tparam Tensor_Dim0 rank of the field
2605 * @tparam Tensor_Dim1 dimension of space
2606 */
2607template <int Tensor_Dim0, int Tensor_Dim1>
2610
2612 boost::shared_ptr<MatrixDouble> data_ptr,
2613 const EntityType zero_type = MBEDGE,
2614 const int zero_side = 0)
2617 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2618 if (!dataPtr)
2619 THROW_MESSAGE("Pointer is not set");
2620 }
2621
2625 const size_t nb_integration_points = getGaussPts().size2();
2626 if (type == zeroType && side == 0) {
2627 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
2628 dataPtr->clear();
2629 }
2630 const size_t nb_dofs = data.getFieldData().size();
2631 if (nb_dofs) {
2632 const size_t nb_base_functions = data.getN().size2() / 3;
2635 auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
2636 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
2637 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2638 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2639 size_t bb = 0;
2640 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2641 double div = t_n_diff_hvec(j, j);
2642 t_data(i) += t_dof(i) * div;
2643 ++t_n_diff_hvec;
2644 ++t_dof;
2645 }
2646 for (; bb < nb_base_functions; ++bb)
2647 ++t_n_diff_hvec;
2648 ++t_data;
2649 }
2650 }
2652 }
2653
2654private:
2655 boost::shared_ptr<MatrixDouble> dataPtr;
2657 const int zeroSide;
2658};
2659
2660/**
2661 * @brief Calculate trace of vector (Hdiv/Hcurl) space
2662 *
2663 * @tparam Tensor_Dim
2664 * @tparam OpBase
2665 */
2666template <int Tensor_Dim, typename OpBase>
2668
2670 boost::shared_ptr<MatrixDouble> data_ptr,
2671 const EntityType zero_type = MBEDGE,
2672 const int zero_side = 0)
2673 : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
2674 zeroType(zero_type), zeroSide(zero_side) {
2675 if (!dataPtr)
2676 THROW_MESSAGE("Pointer is not set");
2677 }
2678
2682 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2683 if (type == zeroType && side == 0) {
2684 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2685 dataPtr->clear();
2686 }
2687 const size_t nb_dofs = data.getFieldData().size();
2688 if (nb_dofs) {
2689 auto t_normal = OpBase::getFTensor1Normal();
2690 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
2691 const size_t nb_base_functions = data.getN().size2() / 3;
2692 auto t_base = data.getFTensor1N<3>();
2693 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2694 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2695 auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
2696 size_t bb = 0;
2697 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2698 t_data(i) += t_dof(i) * (t_base(j) * t_normal(j));
2699 ++t_base;
2700 ++t_dof;
2701 }
2702 for (; bb < nb_base_functions; ++bb)
2703 ++t_base;
2704 ++t_data;
2705 }
2706 }
2708 }
2709
2710private:
2711 boost::shared_ptr<MatrixDouble> dataPtr;
2713 const int zeroSide;
2716};
2717
2718template <>
2722
2724 boost::shared_ptr<MatrixDouble> data_ptr,
2725 const EntityType zero_type = MBEDGE,
2726 const int zero_side = 0)
2727 : UserDataOperator(field_name, OPROW), dataPtr(data_ptr),
2728 zeroType(zero_type), zeroSide(zero_side) {
2729 if (!dataPtr)
2730 THROW_MESSAGE("Pointer is not set");
2731 }
2732
2736 const size_t nb_integration_points = getGaussPts().size2();
2737 if (type == zeroType && side == 0) {
2738 dataPtr->resize(3, nb_integration_points, false);
2739 dataPtr->clear();
2740 }
2741 const size_t nb_dofs = data.getFieldData().size();
2742 if (nb_dofs) {
2743 auto t_normal = getFTensor1Normal();
2744 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
2745 const size_t nb_base_functions = data.getN().size2() / 3;
2746 auto t_base = data.getFTensor1N<3>();
2747 auto t_data = getFTensor1FromMat<3>(*dataPtr);
2748 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2749 auto t_dof = data.getFTensor1FieldData<3>();
2750 if (getNormalsAtGaussPts().size1() == nb_integration_points) {
2751 VectorDouble n = getNormalsAtGaussPts(gg);
2752 auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
2753 t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
2754 }
2755 size_t bb = 0;
2756 for (; bb != nb_dofs / 3; ++bb) {
2757 t_data(i) += t_dof(i) * (t_base(j) * t_normal(j));
2758 ++t_base;
2759 ++t_dof;
2760 }
2761 for (; bb < nb_base_functions; ++bb)
2762 ++t_base;
2763 ++t_data;
2764 }
2765 }
2767 }
2768
2769private:
2770 boost::shared_ptr<MatrixDouble> dataPtr;
2772 const int zeroSide;
2775};
2776
2777/**@}*/
2778
2779/** \name Other operators */
2780
2781/**@{*/
2782
2783/**@}*/
2784
2785/** \name Operators for faces */
2786
2787/**@{*/
2788
2789/** \brief Transform local reference derivatives of shape functions to global
2790derivatives
2791
2792\ingroup mofem_forces_and_sources_tri_element
2793
2794*/
2795template <int DIM, int DERIVATIVE = 1> struct OpSetInvJacSpaceForFaceImpl;
2796
2799
2801 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2802
2803protected:
2804 template <int D1, int D2, int J1, int J2>
2807
2808 static_assert(D2 == J2, "Dimension of jacobian and dimension of <out> "
2809 "directive does not match");
2810
2811 size_t nb_functions = diff_n.size2() / D1;
2812 if (nb_functions) {
2813 size_t nb_gauss_pts = diff_n.size1();
2814
2815#ifndef NDEBUG
2816 if (nb_gauss_pts != getGaussPts().size2())
2817 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2818 "Wrong number of Gauss Pts");
2819 if (diff_n.size2() % D1)
2820 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2821 "Number of direvatives of base functions and D1 dimension does "
2822 "not match");
2823#endif
2824
2825 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions, false);
2826
2829 auto t_diff_n = getFTensor1FromPtr<D2>(&*diffNinvJac.data().begin());
2830 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2831 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*invJacPtr);
2832 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2833 for (size_t dd = 0; dd != nb_functions; ++dd) {
2834 t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
2835 ++t_diff_n;
2836 ++t_diff_n_ref;
2837 }
2838 }
2839
2840 diff_n.swap(diffNinvJac);
2841 }
2843 }
2844
2845 boost::shared_ptr<MatrixDouble> invJacPtr;
2847};
2848
2849template <>
2852
2854
2855 MoFEMErrorCode doWork(int side, EntityType type,
2857};
2858
2859template <>
2861 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2862
2863 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2864
2865 MoFEMErrorCode doWork(int side, EntityType type,
2867};
2868
2869template <>
2871 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2872
2873 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2874
2875 MoFEMErrorCode doWork(int side, EntityType type,
2877};
2878
2879template <int DERIVARIVE = 1>
2881 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2882 OpSetInvJacH1ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2883 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(H1, inv_jac_ptr) {}
2884};
2885
2886template <int DERIVARIVE = 1>
2888 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2889 OpSetInvJacL2ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2890 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(L2, inv_jac_ptr) {}
2891};
2892
2893template <int DERIVARIVE = 1>
2895 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2897 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2898 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(H1, inv_jac_ptr) {}
2899};
2900
2901template <int DERIVARIVE = 1>
2903 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2905 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2906 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(L2, inv_jac_ptr) {}
2907};
2908
2909/**
2910 * \brief Transform local reference derivatives of shape function to
2911 global derivatives for face
2912
2913 * \ingroup mofem_forces_and_sources_tri_element
2914 */
2915template <int DIM> struct OpSetInvJacHcurlFaceImpl;
2916
2917template <>
2920
2921 OpSetInvJacHcurlFaceImpl(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2923 invJacPtr(inv_jac_ptr) {}
2924
2925 MoFEMErrorCode doWork(int side, EntityType type,
2927
2928protected:
2929 boost::shared_ptr<MatrixDouble> invJacPtr;
2931};
2932
2933template <>
2935 using OpSetInvJacHcurlFaceImpl<2>::OpSetInvJacHcurlFaceImpl;
2936 MoFEMErrorCode doWork(int side, EntityType type,
2938};
2939
2942
2943/**
2944 * @brief Make Hdiv space from Hcurl space in 2d
2945 * @ingroup mofem_forces_and_sources_tri_element
2946 */
2949
2952
2953 MoFEMErrorCode doWork(int side, EntityType type,
2955};
2956
2957/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
2958 *
2959 * \note Hdiv space is generated by Hcurl space in 2d.
2960 *
2961 * Contravariant Piola transformation
2962 * \f[
2963 * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
2964 * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
2965 * =
2966 * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
2967 * \f]
2968 *
2969 * \ingroup mofem_forces_and_sources
2970 *
2971 */
2973
2974template <>
2977
2979 boost::shared_ptr<MatrixDouble> jac_ptr)
2981 jacPtr(jac_ptr) {}
2982
2983 MoFEMErrorCode doWork(int side, EntityType type,
2985
2986protected:
2987 boost::shared_ptr<MatrixDouble> jacPtr;
2990};
2991
2992template <>
2996 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
2997
2998 MoFEMErrorCode doWork(int side, EntityType type,
3000};
3001
3006
3007/**@}*/
3008
3009/** \name Operators for edges */
3010
3011/**@{*/
3012
3015
3018
3019 MoFEMErrorCode doWork(int side, EntityType type,
3021
3022private:
3023 std::vector<double> l1;
3024};
3025
3026/**
3027 * @deprecated Name is deprecated and this is added for back compatibility
3028 */
3031
3032/**@}*/
3033
3034/** \name Operator for fat prisms */
3035
3036/**@{*/
3037
3038/**
3039 * @brief Operator for fat prism element updating integration weights in the
3040 * volume.
3041 *
3042 * Jacobian on the distorted element is nonconstant. This operator updates
3043 * integration weight on prism to take into account nonconstat jacobian.
3044 *
3045 * \f[
3046 * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
3047 * \pmb\xi} \right\| \right)
3048 * \f]
3049 * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
3050 * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
3051 * element coordinate.
3052 *
3053 */
3056
3059
3060 MoFEMErrorCode doWork(int side, EntityType type,
3062};
3063
3064/** \brief Calculate inverse of jacobian for face element
3065
3066 It is assumed that face element is XY plane. Applied
3067 only for 2d problems.
3068
3069 FIXME Generalize function for arbitrary face orientation in 3d space
3070 FIXME Calculate to Jacobins for two faces
3071
3072 \ingroup mofem_forces_and_sources_prism_element
3073
3074*/
3077
3078 OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3080 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3081
3084 invJac(inv_jac) {}
3085 MoFEMErrorCode doWork(int side, EntityType type,
3087
3088private:
3089 const boost::shared_ptr<MatrixDouble> invJacPtr;
3091};
3092
3093/** \brief Transform local reference derivatives of shape functions to global
3094derivatives
3095
3096FIXME Generalize to curved shapes
3097FIXME Generalize to case that top and bottom face has different shape
3098
3099\ingroup mofem_forces_and_sources_prism_element
3100
3101*/
3104
3105 OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3107 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3108
3111 invJac(inv_jac) {}
3112
3113 MoFEMErrorCode doWork(int side, EntityType type,
3115
3116private:
3117 const boost::shared_ptr<MatrixDouble> invJacPtr;
3120};
3121
3122// Flat prism
3123
3124/** \brief Calculate inverse of jacobian for face element
3125
3126 It is assumed that face element is XY plane. Applied
3127 only for 2d problems.
3128
3129 FIXME Generalize function for arbitrary face orientation in 3d space
3130 FIXME Calculate to Jacobins for two faces
3131
3132 \ingroup mofem_forces_and_sources_prism_element
3133
3134*/
3137
3140 invJacF3(inv_jac_f3) {}
3141 MoFEMErrorCode doWork(int side, EntityType type,
3143
3144private:
3146};
3147
3148/** \brief Transform local reference derivatives of shape functions to global
3149derivatives
3150
3151FIXME Generalize to curved shapes
3152FIXME Generalize to case that top and bottom face has different shape
3153
3154\ingroup mofem_forces_and_sources_prism_element
3155
3156*/
3159
3162 invJacF3(inv_jac_f3) {}
3163
3164 MoFEMErrorCode doWork(int side, EntityType type,
3166
3167private:
3170};
3171
3172/**@}*/
3173
3174/** \name Operation on matrices at integration points */
3175
3176/**@{*/
3177
3178template <int DIM>
3180
3181 OpInvertMatrix(boost::shared_ptr<MatrixDouble> in_ptr,
3182 boost::shared_ptr<VectorDouble> det_ptr,
3183 boost::shared_ptr<MatrixDouble> out_ptr)
3185 outPtr(out_ptr), detPtr(det_ptr) {}
3186
3189 return doWorkImpl(side, type, data, FTensor::Number<DIM>());
3190 }
3191
3192private:
3193 boost::shared_ptr<MatrixDouble> inPtr;
3194 boost::shared_ptr<MatrixDouble> outPtr;
3195 boost::shared_ptr<VectorDouble> detPtr;
3196
3197 MoFEMErrorCode doWorkImpl(int side, EntityType type,
3199 const FTensor::Number<3> &);
3200
3201 MoFEMErrorCode doWorkImpl(int side, EntityType type,
3203 const FTensor::Number<2> &);
3204};
3205
3206template <int DIM>
3209 const FTensor::Number<3> &) {
3211
3212 if (!inPtr)
3213 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3214 "Pointer for inPtr matrix not allocated");
3215 if (!detPtr)
3216 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3217 "Pointer for detPtr matrix not allocated");
3218
3219 const auto nb_rows = inPtr->size1();
3220 const auto nb_integration_pts = inPtr->size2();
3221
3222 // Calculate determinant
3223 {
3224 detPtr->resize(nb_integration_pts, false);
3225 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3226 auto t_det = getFTensor0FromVec(*detPtr);
3227 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3228 determinantTensor3by3(t_in, t_det);
3229 ++t_in;
3230 ++t_det;
3231 }
3232 }
3233
3234 // Invert jacobian
3235 if (outPtr) {
3236 outPtr->resize(nb_rows, nb_integration_pts, false);
3237 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3238 auto t_out = getFTensor2FromMat<3, 3>(*outPtr);
3239 auto t_det = getFTensor0FromVec(*detPtr);
3240 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3241 invertTensor3by3(t_in, t_det, t_out);
3242 ++t_in;
3243 ++t_out;
3244 ++t_det;
3245 }
3246 }
3247
3249}
3250
3251template <int DIM>
3254 const FTensor::Number<2> &) {
3256
3257 if (!inPtr)
3258 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3259 "Pointer for inPtr matrix not allocated");
3260 if (!detPtr)
3261 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3262 "Pointer for detPtr matrix not allocated");
3263
3264 const auto nb_rows = inPtr->size1();
3265 const auto nb_integration_pts = inPtr->size2();
3266
3267 // Calculate determinant
3268 {
3269 detPtr->resize(nb_integration_pts, false);
3270 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3271 auto t_det = getFTensor0FromVec(*detPtr);
3272 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3273 determinantTensor2by2(t_in, t_det);
3274 ++t_in;
3275 ++t_det;
3276 }
3277 }
3278
3279 // Invert jacobian
3280 if (outPtr) {
3281 outPtr->resize(nb_rows, nb_integration_pts, false);
3282 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3283 auto t_out = getFTensor2FromMat<2, 2>(*outPtr);
3284 auto t_det = getFTensor0FromVec(*detPtr);
3285 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3286 invertTensor2by2(t_in, t_det, t_out);
3287 ++t_in;
3288 ++t_out;
3289 ++t_det;
3290 }
3291 }
3292
3294}
3295
3296/**@}*/
3297
3298} // namespace MoFEM
3299
3300#endif // __USER_DATA_OPERATORS_HPP__
3301
3302/**
3303 * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
3304 *
3305 * \brief Classes and functions used to evaluate fields at integration pts,
3306 *jacobians, etc..
3307 *
3308 * \ingroup mofem_forces_and_sources
3309 **/
static Index< 'J', 3 > J
static Index< 'I', 3 > I
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
static const char *const CoordinateTypesNames[]
Coordinate system names.
Definition: definitions.h:108
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
@ CYLINDRICAL
Definition: definitions.h:117
@ POLAR
Definition: definitions.h:116
@ CARTESIAN
Definition: definitions.h:115
@ SPHERICAL
Definition: definitions.h:118
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
FTensor::Index< 'n', SPACE_DIM > n
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:277
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VecAllocator< double > DoubleAllocator
Definition: Types.hpp:62
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
Definition: Templates.hpp:1353
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1384
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:1323
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1312
constexpr AssemblyType A
Definition: plastic.cpp:46
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.
OpCalculateHVecVectorFieldDot(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.
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:39
static constexpr Switches CtxSetX_TT
Definition: LoopMethods.hpp:41
static constexpr Switches CtxSetX_T
Definition: LoopMethods.hpp:40
intrusive_ptr for managing petsc objects