v0.14.0
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 || vec.size() != nb_gauss_pts) {
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 || vec.size() != nb_gauss_pts) {
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
354 boost::shared_ptr<MatrixDouble> data_ptr,
355 SmartPetscObj<Vec> data_vec,
356 const EntityType zero_type = MBVERTEX)
359 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type) {
360 if (!dataPtr)
361 THROW_MESSAGE("Pointer is not set");
362 }
363
364 MoFEMErrorCode doWork(int side, EntityType type,
366
367protected:
368 boost::shared_ptr<MatrixDouble> dataPtr;
372};
373
374/**
375 * \brief Member function specialization calculating values for tenor field rank
376 *
377 */
378template <int Tensor_Dim>
380 Tensor_Dim, double, ublas::row_major,
381 DoubleAllocator>::doWork(int side, EntityType type,
384
385 const size_t nb_gauss_pts = getGaussPts().size2();
386 auto &mat = *dataPtr;
387 if (type == zeroType || mat.size2() != nb_gauss_pts) {
388 mat.resize(Tensor_Dim, nb_gauss_pts, false);
389 mat.clear();
390 }
391
392 const size_t nb_dofs = data.getFieldData().size();
393 if (nb_dofs) {
394
395 if (dataVec.use_count()) {
396 dotVector.resize(nb_dofs, false);
397 const double *array;
398 CHKERR VecGetArrayRead(dataVec, &array);
399 const auto &local_indices = data.getLocalIndices();
400 for (int i = 0; i != local_indices.size(); ++i)
401 if (local_indices[i] != -1)
402 dotVector[i] = array[local_indices[i]];
403 else
404 dotVector[i] = 0;
405 CHKERR VecRestoreArrayRead(dataVec, &array);
406 data.getFieldData().swap(dotVector);
407 }
408
409 if (nb_gauss_pts) {
410 const size_t nb_base_functions = data.getN().size2();
411 auto base_function = data.getFTensor0N();
412 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
414 const size_t size = nb_dofs / Tensor_Dim;
415 if (nb_dofs % Tensor_Dim) {
416 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
417 "Data inconsistency");
418 }
419 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
420 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
421 size_t bb = 0;
422 for (; bb != size; ++bb) {
423 values_at_gauss_pts(I) += field_data(I) * base_function;
424 ++field_data;
425 ++base_function;
426 }
427 // Number of dofs can be smaller than number of Tensor_Dim x base
428 // functions
429 for (; bb < nb_base_functions; ++bb)
430 ++base_function;
431 ++values_at_gauss_pts;
432 }
433 }
434
435 if (dataVec.use_count()) {
436 data.getFieldData().swap(dotVector);
437 }
438 }
440}
441
442/** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
443 * field
444 *
445 * \ingroup mofem_forces_and_sources_user_data_operators
446 */
447template <int Tensor_Dim>
450 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
451
453 Tensor_Dim, double, ublas::row_major,
454 DoubleAllocator>::OpCalculateVectorFieldValues_General;
455};
456
457/**@}*/
458
459/** \name Vector field values at integration points */
460
461/**@{*/
462
463/** \brief Calculate field values (template specialization) for tensor field
464 * rank 1, i.e. vector field
465 *
466 */
467template <int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
470
472 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
473 const EntityType zero_type = MBVERTEX)
476 dataPtr(data_ptr), zeroType(zero_type) {
477 if (!dataPtr)
478 THROW_MESSAGE("Pointer is not set");
479 }
480
484
485 // When we move to C++17 add if constexpr()
486 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
487 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
488 "%s coordiante not implemented",
489 CoordinateTypesNames[COORDINATE_SYSTEM]);
490
491 const size_t nb_gauss_pts = getGaussPts().size2();
492 auto &vec = *dataPtr;
493 if (type == zeroType) {
494 vec.resize(nb_gauss_pts, false);
495 vec.clear();
496 }
497
498 const size_t nb_dofs = data.getFieldData().size();
499 if (nb_dofs) {
500
501 if (nb_gauss_pts) {
502 const size_t nb_base_functions = data.getN().size2();
503 auto values_at_gauss_pts = getFTensor0FromVec(vec);
505 const size_t size = nb_dofs / Tensor_Dim;
506#ifndef NDEBUG
507 if (nb_dofs % Tensor_Dim) {
508 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
509 "Number of dofs should multiple of dimensions");
510 }
511#endif
512
513 // When we move to C++17 add if constexpr()
514 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
515 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
516 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
517 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
518 size_t bb = 0;
519 for (; bb != size; ++bb) {
520 values_at_gauss_pts += field_data(I) * diff_base_function(I);
521 ++field_data;
522 ++diff_base_function;
523 }
524 // Number of dofs can be smaller than number of Tensor_Dim x base
525 // functions
526 for (; bb < nb_base_functions; ++bb)
527 ++diff_base_function;
528 ++values_at_gauss_pts;
529 }
530 }
531
532 // When we move to C++17 add if constexpr()
533 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
534 auto t_coords = getFTensor1CoordsAtGaussPts();
535 auto values_at_gauss_pts = getFTensor0FromVec(vec);
536 auto base_function = data.getFTensor0N();
537 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
538 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
539 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
540 size_t bb = 0;
541 for (; bb != size; ++bb) {
542 values_at_gauss_pts += field_data(I) * diff_base_function(I);
543 values_at_gauss_pts +=
544 field_data(0) * base_function / t_coords(0);
545 ++field_data;
546 ++base_function;
547 ++diff_base_function;
548 }
549 // Number of dofs can be smaller than number of Tensor_Dim x base
550 // functions
551 for (; bb < nb_base_functions; ++bb) {
552 ++base_function;
553 ++diff_base_function;
554 }
555 ++values_at_gauss_pts;
556 ++t_coords;
557 }
558 }
559 }
560 }
562 }
563
564protected:
565 boost::shared_ptr<VectorDouble> dataPtr;
567};
568
569/** \brief Approximate field values for given petsc vector
570 *
571 * \note Look at PetscData to see what vectors could be extracted with that user
572 * data operator.
573 *
574 * \ingroup mofem_forces_and_sources_user_data_operators
575 */
576template <int Tensor_Dim, PetscData::DataContext CTX>
579
581 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
582 const EntityType zero_at_type = MBVERTEX)
585 dataPtr(data_ptr), zeroAtType(zero_at_type) {
586 if (!dataPtr)
587 THROW_MESSAGE("Pointer is not set");
588 }
589
593 auto &local_indices = data.getLocalIndices();
594 const size_t nb_dofs = local_indices.size();
595 const size_t nb_gauss_pts = getGaussPts().size2();
596
597 MatrixDouble &mat = *dataPtr;
598 if (type == zeroAtType) {
599 mat.resize(Tensor_Dim, nb_gauss_pts, false);
600 mat.clear();
601 }
602 if (!nb_dofs)
604
605 const double *array;
606
607 auto get_array = [&](const auto ctx, auto vec) {
609#ifndef NDEBUG
610 if ((getFEMethod()->data_ctx & ctx).none()) {
611 MOFEM_LOG_CHANNEL("SELF");
612 MOFEM_LOG("SELF", Sev::error)
613 << "In this case filed degrees of freedom are read from vector. "
614 "That usually happens when time solver is used, and acces to "
615 "first or second rates is needed. You probably not set ts_u, "
616 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
617 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
618 "respectively";
619 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
620 }
621#endif
622 CHKERR VecGetArrayRead(vec, &array);
624 };
625
626 auto restore_array = [&](auto vec) {
627 return VecRestoreArrayRead(vec, &array);
628 };
629
630 switch (CTX) {
632 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
633 break;
635 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
636 break;
638 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
639 break;
640 default:
641 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
642 "That case is not implemented");
643 }
644
645 dotVector.resize(local_indices.size());
646 for (int i = 0; i != local_indices.size(); ++i)
647 if (local_indices[i] != -1)
648 dotVector[i] = array[local_indices[i]];
649 else
650 dotVector[i] = 0;
651
652 switch (CTX) {
654 CHKERR restore_array(getFEMethod()->ts_u);
655 break;
657 CHKERR restore_array(getFEMethod()->ts_u_t);
658 break;
660 CHKERR restore_array(getFEMethod()->ts_u_tt);
661 break;
662 default:
663 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
664 "That case is not implemented");
665 }
666
667 const size_t nb_base_functions = data.getN().size2();
668 auto base_function = data.getFTensor0N();
669 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
670
672 const size_t size = nb_dofs / Tensor_Dim;
673 if (nb_dofs % Tensor_Dim) {
674 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
675 }
676 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
677 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(dotVector);
678 size_t bb = 0;
679 for (; bb != size; ++bb) {
680 values_at_gauss_pts(I) += field_data(I) * base_function;
681 ++field_data;
682 ++base_function;
683 }
684 // Number of dofs can be smaller than number of Tensor_Dim x base
685 // functions
686 for (; bb < nb_base_functions; ++bb)
687 ++base_function;
688 ++values_at_gauss_pts;
689 }
691 }
692
693protected:
694 boost::shared_ptr<MatrixDouble> dataPtr;
697};
698
699/** \brief Get time derivatives of values at integration pts for tensor filed
700 * rank 1, i.e. vector field
701 *
702 * \ingroup mofem_forces_and_sources_user_data_operators
703 */
704template <int Tensor_Dim>
708
709/** \brief Get second time derivatives of values at integration pts for tensor
710 * filed rank 1, i.e. vector field
711 *
712 * \ingroup mofem_forces_and_sources_user_data_operators
713 */
714template <int Tensor_Dim>
718
719/**@}*/
720
721/** \name Tensor field values at integration points */
722
723/**@{*/
724
725/** \brief Calculate field values for tenor field rank 2.
726 *
727 */
728template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
731
733 const std::string field_name,
734 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
735 const EntityType zero_type = MBVERTEX)
738 dataPtr(data_ptr), zeroType(zero_type) {
739 if (!dataPtr)
740 THROW_MESSAGE("Pointer is not set");
741 }
742
743 MoFEMErrorCode doWork(int side, EntityType type,
745
746protected:
747 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
749};
750
751template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
754 doWork(int side, EntityType type, EntitiesFieldData::EntData &data) {
756 SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
757 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
758 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
759 Tensor_Dim0, Tensor_Dim1);
761}
762
763template <int Tensor_Dim0, int Tensor_Dim1>
764struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
765 ublas::row_major, DoubleAllocator>
767
769 const std::string field_name,
770 boost::shared_ptr<
771 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
772 data_ptr,
773 const EntityType zero_type = MBVERTEX)
776 dataPtr(data_ptr), zeroType(zero_type) {
777 if (!dataPtr)
778 THROW_MESSAGE("Pointer is not set");
779 }
780
782 const std::string field_name,
783 boost::shared_ptr<
784 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
785 data_ptr,
786 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
789 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
790 if (!dataPtr)
791 THROW_MESSAGE("Pointer is not set");
792 }
793
794 MoFEMErrorCode doWork(int side, EntityType type,
796
797protected:
798 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
803};
804
805template <int Tensor_Dim0, int Tensor_Dim1>
807 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
808 DoubleAllocator>::doWork(int side, EntityType type,
811
812 MatrixDouble &mat = *dataPtr;
813 const size_t nb_gauss_pts = data.getN().size1();
814 if (type == zeroType) {
815 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
816 mat.clear();
817 }
818
819 const size_t nb_dofs = data.getFieldData().size();
820
821 if (dataVec.use_count()) {
822 dotVector.resize(nb_dofs, false);
823 const double *array;
824 CHKERR VecGetArrayRead(dataVec, &array);
825 const auto &local_indices = data.getLocalIndices();
826 for (int i = 0; i != local_indices.size(); ++i)
827 if (local_indices[i] != -1)
828 dotVector[i] = array[local_indices[i]];
829 else
830 dotVector[i] = 0;
831 CHKERR VecRestoreArrayRead(dataVec, &array);
832 data.getFieldData().swap(dotVector);
833 }
834
835 if (nb_dofs) {
836 const size_t nb_base_functions = data.getN().size2();
837 auto base_function = data.getFTensor0N();
838 auto values_at_gauss_pts =
839 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
842 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
843 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
844 auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
845 size_t bb = 0;
846 for (; bb != size; ++bb) {
847 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
848 ++field_data;
849 ++base_function;
850 }
851 for (; bb < nb_base_functions; ++bb)
852 ++base_function;
853 ++values_at_gauss_pts;
854 }
855
856 if (dataVec.use_count()) {
857 data.getFieldData().swap(dotVector);
858 }
859 }
861}
862
863/** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
864 * field
865 *
866 * \ingroup mofem_forces_and_sources_user_data_operators
867 */
868template <int Tensor_Dim0, int Tensor_Dim1>
871 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
872
874 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
875 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
876};
877
878/** \brief Get time direvarive values at integration pts for tensor filed rank
879 * 2, i.e. matrix field
880 *
881 * \ingroup mofem_forces_and_sources_user_data_operators
882 */
883template <int Tensor_Dim0, int Tensor_Dim1>
886
888 boost::shared_ptr<MatrixDouble> data_ptr,
889 const EntityType zero_at_type = MBVERTEX)
892 dataPtr(data_ptr), zeroAtType(zero_at_type) {
893 if (!dataPtr)
894 THROW_MESSAGE("Pointer is not set");
895 }
896
900
901 const size_t nb_gauss_pts = getGaussPts().size2();
902 MatrixDouble &mat = *dataPtr;
903 if (type == zeroAtType) {
904 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
905 mat.clear();
906 }
907 const auto &local_indices = data.getLocalIndices();
908 const size_t nb_dofs = local_indices.size();
909 if (nb_dofs) {
910 dotVector.resize(nb_dofs, false);
911 const double *array;
912 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
913 for (size_t i = 0; i != local_indices.size(); ++i)
914 if (local_indices[i] != -1)
915 dotVector[i] = array[local_indices[i]];
916 else
917 dotVector[i] = 0;
918 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
919
920 const size_t nb_base_functions = data.getN().size2();
921
922 auto base_function = data.getFTensor0N();
923 auto values_at_gauss_pts =
924 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
927 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
928 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
929 auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
930 size_t bb = 0;
931 for (; bb != size; ++bb) {
932 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
933 ++field_data;
934 ++base_function;
935 }
936 for (; bb < nb_base_functions; ++bb)
937 ++base_function;
938 ++values_at_gauss_pts;
939 }
940 }
942 }
943
944protected:
945 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
946 EntityType zeroAtType; ///< Zero values at Gauss point at this type
947 VectorDouble dotVector; ///< Keeps temoorary values of time directives
948
949 template <int Dim0, int Dim1> auto getFTensorDotData() {
950 static_assert(Dim0 || !Dim0 || Dim1 || !Dim1, "not implemented");
951 }
952};
953
954template <>
955template <>
958 &dotVector[0], &dotVector[1], &dotVector[2],
959
960 &dotVector[3], &dotVector[4], &dotVector[5],
961
962 &dotVector[6], &dotVector[7], &dotVector[8]);
963}
964
965template <>
966template <>
969 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3]);
970}
971
972/**
973 * @brief Calculate symmetric tensor field values at integration pts.
974 *
975 * @tparam Tensor_Dim
976
977 * \ingroup mofem_forces_and_sources_user_data_operators
978 */
979template <int Tensor_Dim>
982
984 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
985 const EntityType zero_type = MBEDGE, const int zero_side = 0)
988 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
989 if (!dataPtr)
990 THROW_MESSAGE("Pointer is not set");
991 }
992
994 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
995 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBEDGE,
996 const int zero_side = 0)
999 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side),
1000 dataVec(data_vec) {
1001 if (!dataPtr)
1002 THROW_MESSAGE("Pointer is not set");
1003 }
1004
1008 MatrixDouble &mat = *dataPtr;
1009 const int nb_gauss_pts = getGaussPts().size2();
1010 if (type == this->zeroType && side == zeroSide) {
1011 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
1012 mat.clear();
1013 }
1014 const int nb_dofs = data.getFieldData().size();
1015 if (!nb_dofs)
1017
1018 if (dataVec.use_count()) {
1019 dotVector.resize(nb_dofs, false);
1020 const double *array;
1021 CHKERR VecGetArrayRead(dataVec, &array);
1022 const auto &local_indices = data.getLocalIndices();
1023 for (int i = 0; i != local_indices.size(); ++i)
1024 if (local_indices[i] != -1)
1025 dotVector[i] = array[local_indices[i]];
1026 else
1027 dotVector[i] = 0;
1028 CHKERR VecRestoreArrayRead(dataVec, &array);
1029 data.getFieldData().swap(dotVector);
1030 }
1031
1032 const int nb_base_functions = data.getN().size2();
1033 auto base_function = data.getFTensor0N();
1034 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1037 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1038 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1039 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
1040 int bb = 0;
1041 for (; bb != size; ++bb) {
1042 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1043 ++field_data;
1044 ++base_function;
1045 }
1046 for (; bb < nb_base_functions; ++bb)
1047 ++base_function;
1048 ++values_at_gauss_pts;
1049 }
1050
1051 if (dataVec.use_count()) {
1052 data.getFieldData().swap(dotVector);
1053 }
1054
1056 }
1057
1058protected:
1059 boost::shared_ptr<MatrixDouble> dataPtr;
1061 const int zeroSide;
1064};
1065
1066/**
1067 * @brief Calculate symmetric tensor field rates ant integratio pts.
1068 *
1069 * @tparam Tensor_Dim
1070 *
1071 * \ingroup mofem_forces_and_sources_user_data_operators
1072 */
1073template <int Tensor_Dim>
1076
1078 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1079 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1082 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1083 if (!dataPtr)
1084 THROW_MESSAGE("Pointer is not set");
1085 }
1086
1090 const int nb_gauss_pts = getGaussPts().size2();
1091 MatrixDouble &mat = *dataPtr;
1092 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1093 if (type == zeroType && side == zeroSide) {
1094 mat.resize(symm_size, nb_gauss_pts, false);
1095 mat.clear();
1096 }
1097 auto &local_indices = data.getLocalIndices();
1098 const int nb_dofs = local_indices.size();
1099 if (!nb_dofs)
1101
1102#ifndef NDEBUG
1103 if ((getFEMethod()->data_ctx & PetscData::CtxSetX_T).none()) {
1104 MOFEM_LOG_CHANNEL("SELF");
1105 MOFEM_LOG("SELF", Sev::error)
1106 << "In this case filed degrees of freedom are read from vector. "
1107 "That usually happens when time solver is used, and acces to "
1108 "first rates is needed. You probably not set "
1109 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1110 "respectively";
1111 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
1112 }
1113#endif
1114
1115 dotVector.resize(nb_dofs, false);
1116 const double *array;
1117 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1118 for (int i = 0; i != local_indices.size(); ++i)
1119 if (local_indices[i] != -1)
1120 dotVector[i] = array[local_indices[i]];
1121 else
1122 dotVector[i] = 0;
1123 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1124
1125 const int nb_base_functions = data.getN().size2();
1126
1127 auto base_function = data.getFTensor0N();
1128 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1131 const int size = nb_dofs / symm_size;
1132 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1133 auto field_data = getFTensorDotData<Tensor_Dim>();
1134 int bb = 0;
1135 for (; bb != size; ++bb) {
1136 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1137 ++field_data;
1138 ++base_function;
1139 }
1140 for (; bb < nb_base_functions; ++bb)
1141 ++base_function;
1142 ++values_at_gauss_pts;
1143 }
1144
1146 }
1147
1148protected:
1149 boost::shared_ptr<MatrixDouble> dataPtr;
1151 const int zeroSide;
1153
1154 template <int Dim> inline auto getFTensorDotData() {
1155 static_assert(Dim || !Dim, "not implemented");
1156 }
1157};
1158
1159template <>
1160template <>
1161inline auto
1164 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1165 &dotVector[5]);
1166}
1167
1168template <>
1169template <>
1170inline auto
1173 &dotVector[0], &dotVector[1], &dotVector[2]);
1174}
1175
1176/**@}*/
1177
1178/** \name Gradients and Hessian of scalar fields at integration points */
1179
1180/**@{*/
1181
1182/**
1183 * \brief Evaluate field gradient values for scalar field, i.e. gradient is
1184 * tensor rank 1 (vector)
1185 *
1186 */
1187template <int Tensor_Dim, class T, class L, class A>
1189 : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
1190
1192 Tensor_Dim, T, L, A>::OpCalculateVectorFieldValues_General;
1193};
1194
1195/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1196 * tensor rank 1 (vector), specialization
1197 *
1198 */
1199template <int Tensor_Dim>
1201 ublas::row_major, DoubleAllocator>
1203 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1204
1206 Tensor_Dim, double, ublas::row_major,
1207 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1208
1209 /**
1210 * \brief calculate gradient values of scalar field at integration points
1211 * @param side side entity number
1212 * @param type side entity type
1213 * @param data entity data
1214 * @return error code
1215 */
1216 MoFEMErrorCode doWork(int side, EntityType type,
1218};
1219
1220/**
1221 * \brief Member function specialization calculating scalar field gradients for
1222 * tenor field rank 1
1223 *
1224 */
1225template <int Tensor_Dim>
1227 Tensor_Dim, double, ublas::row_major,
1228 DoubleAllocator>::doWork(int side, EntityType type,
1231
1232 const size_t nb_gauss_pts = this->getGaussPts().size2();
1233 auto &mat = *this->dataPtr;
1234 if (type == this->zeroType) {
1235 mat.resize(Tensor_Dim, nb_gauss_pts, false);
1236 mat.clear();
1237 }
1238
1239 const int nb_dofs = data.getFieldData().size();
1240 if (nb_dofs) {
1241
1242 if (nb_gauss_pts) {
1243 const int nb_base_functions = data.getN().size2();
1244 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
1245 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1246
1247#ifndef NDEBUG
1248 if (nb_dofs > nb_base_functions)
1249 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1250 "Number of base functions inconsistent with number of DOFs "
1251 "(%d > %d)",
1252 nb_dofs, nb_base_functions);
1253
1254 if (data.getDiffN().size2() != nb_base_functions * Tensor_Dim)
1255 SETERRQ2(
1256 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1257 "Number of base functions inconsistent with number of derivatives "
1258 "(%d != %d)",
1259 data.getDiffN().size2(), nb_base_functions);
1260
1261 if (data.getDiffN().size1() != nb_gauss_pts)
1262 SETERRQ2(
1263 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1264 "Number of base functions inconsistent with number of integration "
1265 "pts (%d != %d)",
1266 data.getDiffN().size2(), nb_gauss_pts);
1267
1268#endif
1269
1271 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1272 auto field_data = data.getFTensor0FieldData();
1273 int bb = 0;
1274 for (; bb != nb_dofs; ++bb) {
1275 gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
1276 ++field_data;
1277 ++diff_base_function;
1278 }
1279 // Number of dofs can be smaller than number of base functions
1280 for (; bb < nb_base_functions; ++bb)
1281 ++diff_base_function;
1282 ++gradients_at_gauss_pts;
1283 }
1284 }
1285 }
1286
1288}
1289
1290/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1291 * vector field
1292 *
1293 * \ingroup mofem_forces_and_sources_user_data_operators
1294 */
1295template <int Tensor_Dim>
1298 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1300 Tensor_Dim, double, ublas::row_major,
1301 DoubleAllocator>::OpCalculateScalarFieldGradient_General;
1302};
1303
1304/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1305 * tensor rank 1 (vector), specialization
1306 *
1307 */
1308template <int Tensor_Dim>
1311 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1312
1314 Tensor_Dim, double, ublas::row_major,
1315 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1316
1317 /**
1318 * \brief calculate gradient values of scalar field at integration points
1319 * @param side side entity number
1320 * @param type side entity type
1321 * @param data entity data
1322 * @return error code
1323 */
1324 MoFEMErrorCode doWork(int side, EntityType type,
1326};
1327
1328template <int Tensor_Dim>
1330 int side, EntityType type, EntitiesFieldData::EntData &data) {
1332
1333 const size_t nb_gauss_pts = this->getGaussPts().size2();
1334 auto &mat = *this->dataPtr;
1335 if (type == this->zeroType) {
1336 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts, false);
1337 mat.clear();
1338 }
1339
1340 const int nb_dofs = data.getFieldData().size();
1341 if (nb_dofs) {
1342
1343 if (nb_gauss_pts) {
1344 const int nb_base_functions = data.getN().size2();
1345
1346 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1347#ifndef NDEBUG
1348 if (hessian_base.size1() != nb_gauss_pts) {
1349 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1350 "Wrong number of integration pts (%d != %d)",
1351 hessian_base.size1(), nb_gauss_pts);
1352 }
1353 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1354 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1355 "Wrong number of base functions (%d != %d)",
1356 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1357 nb_base_functions);
1358 }
1359 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1360 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1361 "Wrong number of base functions (%d < %d)",
1362 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1363 }
1364#endif
1365
1366 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1367 &*hessian_base.data().begin());
1368
1369 auto hessian_at_gauss_pts =
1370 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1371
1374 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1375 auto field_data = data.getFTensor0FieldData();
1376 int bb = 0;
1377 for (; bb != nb_dofs; ++bb) {
1378 hessian_at_gauss_pts(I, J) +=
1379 field_data * t_diff2_base_function(I, J);
1380 ++field_data;
1381 ++t_diff2_base_function;
1382 }
1383 // Number of dofs can be smaller than number of base functions
1384 for (; bb < nb_base_functions; ++bb) {
1385 ++t_diff2_base_function;
1386 }
1387
1388 ++hessian_at_gauss_pts;
1389 }
1390 }
1391 }
1392
1394}
1395
1396/**}*/
1397
1398/** \name Gradients and hessian of tensor fields at integration points */
1399
1400/**@{*/
1401
1402/**
1403 * \brief Evaluate field gradient values for vector field, i.e. gradient is
1404 * tensor rank 2
1405 *
1406 */
1407template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1409 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1410 L, A> {
1411
1413 Tensor_Dim0, Tensor_Dim1, T, L, A>::OpCalculateTensor2FieldValues_General;
1414};
1415
1416template <int Tensor_Dim0, int Tensor_Dim1>
1417struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1418 ublas::row_major, DoubleAllocator>
1420 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1421
1423 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1424 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1425
1426 /**
1427 * \brief calculate values of vector field at integration points
1428 * @param side side entity number
1429 * @param type side entity type
1430 * @param data entity data
1431 * @return error code
1432 */
1433 MoFEMErrorCode doWork(int side, EntityType type,
1435};
1436
1437/**
1438 * \brief Member function specialization calculating vector field gradients for
1439 * tenor field rank 2
1440 *
1441 */
1442template <int Tensor_Dim0, int Tensor_Dim1>
1444 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1445 DoubleAllocator>::doWork(int side, EntityType type,
1448 if (!this->dataPtr)
1449 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1450 "Data pointer not allocated");
1451
1452 const size_t nb_gauss_pts = this->getGaussPts().size2();
1453 auto &mat = *this->dataPtr;
1454 if (type == this->zeroType) {
1455 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1456 mat.clear();
1457 }
1458
1459 if (nb_gauss_pts) {
1460 const size_t nb_dofs = data.getFieldData().size();
1461
1462 if (nb_dofs) {
1463
1464 if (this->dataVec.use_count()) {
1465 this->dotVector.resize(nb_dofs, false);
1466 const double *array;
1467 CHKERR VecGetArrayRead(this->dataVec, &array);
1468 const auto &local_indices = data.getLocalIndices();
1469 for (int i = 0; i != local_indices.size(); ++i)
1470 if (local_indices[i] != -1)
1471 this->dotVector[i] = array[local_indices[i]];
1472 else
1473 this->dotVector[i] = 0;
1474 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1475 data.getFieldData().swap(this->dotVector);
1476 }
1477
1478 const int nb_base_functions = data.getN().size2();
1479 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1480 auto gradients_at_gauss_pts =
1481 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1484 int size = nb_dofs / Tensor_Dim0;
1485 if (nb_dofs % Tensor_Dim0) {
1486 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1487 "Data inconsistency");
1488 }
1489 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1490 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1491 int bb = 0;
1492 for (; bb < size; ++bb) {
1493 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1494 ++field_data;
1495 ++diff_base_function;
1496 }
1497 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1498 // functions
1499 for (; bb != nb_base_functions; ++bb)
1500 ++diff_base_function;
1501 ++gradients_at_gauss_pts;
1502 }
1503
1504 if (this->dataVec.use_count()) {
1505 data.getFieldData().swap(this->dotVector);
1506 }
1507 }
1508 }
1510}
1511
1512/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1513 * vector field
1514 *
1515 * \ingroup mofem_forces_and_sources_user_data_operators
1516 */
1517template <int Tensor_Dim0, int Tensor_Dim1>
1520 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1521
1523 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1524 DoubleAllocator>::OpCalculateVectorFieldGradient_General;
1525};
1526
1527/** \brief Get field gradients time derivative at integration pts for scalar
1528 * filed rank 0, i.e. vector field
1529 *
1530 * \ingroup mofem_forces_and_sources_user_data_operators
1531 */
1532template <int Tensor_Dim0, int Tensor_Dim1>
1535
1537 boost::shared_ptr<MatrixDouble> data_ptr,
1538 const EntityType zero_at_type = MBVERTEX)
1541 dataPtr(data_ptr), zeroAtType(zero_at_type) {
1542 if (!dataPtr)
1543 THROW_MESSAGE("Pointer is not set");
1544 }
1545
1549
1550 const auto &local_indices = data.getLocalIndices();
1551 const int nb_dofs = local_indices.size();
1552 const int nb_gauss_pts = this->getGaussPts().size2();
1553
1554 auto &mat = *this->dataPtr;
1555 if (type == this->zeroAtType) {
1556 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1557 mat.clear();
1558 }
1559 if (!nb_dofs)
1561
1562 dotVector.resize(nb_dofs, false);
1563 const double *array;
1564 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1565 for (int i = 0; i != local_indices.size(); ++i)
1566 if (local_indices[i] != -1)
1567 dotVector[i] = array[local_indices[i]];
1568 else
1569 dotVector[i] = 0;
1570 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1571
1572 const int nb_base_functions = data.getN().size2();
1573 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1574 auto gradients_at_gauss_pts =
1575 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1578 int size = nb_dofs / Tensor_Dim0;
1579 if (nb_dofs % Tensor_Dim0) {
1580 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1581 }
1582
1583 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1584 auto field_data = getFTensorDotData<Tensor_Dim0>();
1585 int bb = 0;
1586 for (; bb < size; ++bb) {
1587 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1588 ++field_data;
1589 ++diff_base_function;
1590 }
1591 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1592 // functions
1593 for (; bb != nb_base_functions; ++bb)
1594 ++diff_base_function;
1595 ++gradients_at_gauss_pts;
1596 }
1598 }
1599
1600private:
1601 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1602 EntityType zeroAtType; ///< Zero values at Gauss point at this type
1603 VectorDouble dotVector; ///< Keeps temoorary values of time directives
1604
1605 template <int Dim> inline auto getFTensorDotData() {
1606 static_assert(Dim || !Dim, "not implemented");
1607 }
1608};
1609
1610template <>
1611template <>
1614 &dotVector[0], &dotVector[1], &dotVector[2]);
1615}
1616
1617template <>
1618template <>
1620 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(&dotVector[0],
1621 &dotVector[1]);
1622}
1623
1624/**
1625 * \brief Evaluate field gradient values for symmetric 2nd order tensor field,
1626 * i.e. gradient is tensor rank 3
1627 *
1628 */
1629template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1631 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1632 L, A> {
1633
1635 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1636 const EntityType zero_type = MBVERTEX)
1637 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
1638 A>(field_name, data_ptr,
1639 zero_type) {}
1640};
1641
1642template <int Tensor_Dim0, int Tensor_Dim1>
1644 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator>
1646 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1647
1649 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1650 const EntityType zero_type = MBVERTEX)
1651 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1652 ublas::row_major,
1654 field_name, data_ptr, zero_type) {}
1655
1656 /**
1657 * \brief calculate values of vector field at integration points
1658 * @param side side entity number
1659 * @param type side entity type
1660 * @param data entity data
1661 * @return error code
1662 */
1663 MoFEMErrorCode doWork(int side, EntityType type,
1665};
1666
1667/**
1668 * \brief Member function specialization calculating tensor field gradients for
1669 * symmetric tensor field rank 2
1670 *
1671 */
1672template <int Tensor_Dim0, int Tensor_Dim1>
1673MoFEMErrorCode OpCalculateTensor2SymmetricFieldGradient_General<
1674 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1675 DoubleAllocator>::doWork(int side, EntityType type,
1678 if (!this->dataPtr)
1679 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1680 "Data pointer not allocated");
1681
1682 const size_t nb_gauss_pts = this->getGaussPts().size2();
1683 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1684 auto &mat = *this->dataPtr;
1685 if (type == this->zeroType) {
1686 mat.resize(msize * Tensor_Dim1, nb_gauss_pts, false);
1687 mat.clear();
1688 }
1689
1690 if (nb_gauss_pts) {
1691 const size_t nb_dofs = data.getFieldData().size();
1692
1693 if (nb_dofs) {
1694
1695 const int nb_base_functions = data.getN().size2();
1696 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1697 auto gradients_at_gauss_pts =
1698 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1702 int size = nb_dofs / msize;
1703 if (nb_dofs % msize) {
1704 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1705 "Data inconsistency");
1706 }
1707 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1708 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim0>();
1709 int bb = 0;
1710 for (; bb < size; ++bb) {
1711 gradients_at_gauss_pts(I, J, K) +=
1712 field_data(I, J) * diff_base_function(K);
1713 ++field_data;
1714 ++diff_base_function;
1715 }
1716 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1717 // functions
1718 for (; bb != nb_base_functions; ++bb)
1719 ++diff_base_function;
1720 ++gradients_at_gauss_pts;
1721 }
1722 }
1723 }
1725}
1726
1727/** \brief Get field gradients at integration pts for symmetric tensorial field
1728 * rank 2
1729 *
1730 * \ingroup mofem_forces_and_sources_user_data_operators
1731 */
1732template <int Tensor_Dim0, int Tensor_Dim1>
1735 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1736
1738 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1739 const EntityType zero_type = MBVERTEX)
1741 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1742 DoubleAllocator>(field_name, data_ptr, zero_type) {}
1743};
1744
1745template <int Tensor_Dim0, int Tensor_Dim1>
1748 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1749
1751 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1752 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1753
1754 /**
1755 * \brief calculate values of vector field at integration points
1756 * @param side side entity number
1757 * @param type side entity type
1758 * @param data entity data
1759 * @return error code
1760 */
1761 MoFEMErrorCode doWork(int side, EntityType type,
1763};
1764
1765/**
1766 * \brief Member function specialization calculating vector field gradients for
1767 * tenor field rank 2
1768 *
1769 */
1770template <int Tensor_Dim0, int Tensor_Dim1>
1772 int side, EntityType type, EntitiesFieldData::EntData &data) {
1774 if (!this->dataPtr)
1775 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1776 "Data pointer not allocated");
1777
1778 const size_t nb_gauss_pts = this->getGaussPts().size2();
1779 auto &mat = *this->dataPtr;
1780 if (type == this->zeroType) {
1781 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts, false);
1782 mat.clear();
1783 }
1784
1785 if (nb_gauss_pts) {
1786 const size_t nb_dofs = data.getFieldData().size();
1787
1788 if (nb_dofs) {
1789
1790 if (this->dataVec.use_count()) {
1791 this->dotVector.resize(nb_dofs, false);
1792 const double *array;
1793 CHKERR VecGetArrayRead(this->dataVec, &array);
1794 const auto &local_indices = data.getLocalIndices();
1795 for (int i = 0; i != local_indices.size(); ++i)
1796 if (local_indices[i] != -1)
1797 this->dotVector[i] = array[local_indices[i]];
1798 else
1799 this->dotVector[i] = 0;
1800 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1801 data.getFieldData().swap(this->dotVector);
1802 }
1803
1804 const int nb_base_functions = data.getN().size2();
1805
1806 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1807#ifndef NDEBUG
1808 if (hessian_base.size1() != nb_gauss_pts) {
1809 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1810 "Wrong number of integration pts (%d != %d)",
1811 hessian_base.size1(), nb_gauss_pts);
1812 }
1813 if (hessian_base.size2() !=
1814 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1815 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1816 "Wrong number of base functions (%d != %d)",
1817 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1818 nb_base_functions);
1819 }
1820 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1821 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1822 "Wrong number of base functions (%d < %d)",
1823 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1824 }
1825#endif
1826
1827 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1828 &*hessian_base.data().begin());
1829
1830 auto t_hessian_at_gauss_pts =
1831 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1832
1836
1837 int size = nb_dofs / Tensor_Dim0;
1838#ifndef NDEBUG
1839 if (nb_dofs % Tensor_Dim0) {
1840 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1841 "Data inconsistency");
1842 }
1843#endif
1844
1845 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1846 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1847 int bb = 0;
1848 for (; bb < size; ++bb) {
1849 t_hessian_at_gauss_pts(I, J, K) +=
1850 field_data(I) * t_diff2_base_function(J, K);
1851 ++field_data;
1852 ++t_diff2_base_function;
1853 }
1854 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1855 // functions
1856 for (; bb != nb_base_functions; ++bb)
1857 ++t_diff2_base_function;
1858 ++t_hessian_at_gauss_pts;
1859 }
1860
1861 if (this->dataVec.use_count()) {
1862 data.getFieldData().swap(this->dotVector);
1863 }
1864 }
1865 }
1867}
1868
1869/**@}*/
1870
1871/** \name Transform tensors and vectors */
1872
1873/**@{*/
1874
1875/**
1876 * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1877 * \f$
1878 *
1879 * @tparam DIM
1880 *
1881 * \ingroup mofem_forces_and_sources_user_data_operators
1882 */
1883template <int DIM_01, int DIM_23, int S = 0>
1886
1889
1891 boost::shared_ptr<MatrixDouble> in_mat,
1892 boost::shared_ptr<MatrixDouble> out_mat,
1893 boost::shared_ptr<MatrixDouble> d_mat)
1894 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1895 // Only is run for vertices
1896 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1897 if (!inMat)
1898 THROW_MESSAGE("Pointer for in mat is null");
1899 if (!outMat)
1900 THROW_MESSAGE("Pointer for out mat is null");
1901 if (!dMat)
1902 THROW_MESSAGE("Pointer for tensor mat is null");
1903 }
1904
1905 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1907 const size_t nb_gauss_pts = getGaussPts().size2();
1908 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
1909 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
1910 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1911 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
1912 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1913 t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1914 ++t_in;
1915 ++t_out;
1916 }
1918 }
1919
1920private:
1925
1926 boost::shared_ptr<MatrixDouble> inMat;
1927 boost::shared_ptr<MatrixDouble> outMat;
1928 boost::shared_ptr<MatrixDouble> dMat;
1929};
1930
1931template <int DIM>
1933
1936
1938 boost::shared_ptr<MatrixDouble> in_mat,
1939 boost::shared_ptr<MatrixDouble> out_mat)
1940 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1941 // Only is run for vertices
1942 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1943 if (!inMat)
1944 THROW_MESSAGE("Pointer not set for in matrix");
1945 if (!outMat)
1946 THROW_MESSAGE("Pointer not set for in matrix");
1947 }
1948
1949 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1951 const size_t nb_gauss_pts = getGaussPts().size2();
1952 auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
1953 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
1954 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
1955 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1956 t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
1957 ++t_in;
1958 ++t_out;
1959 }
1961 }
1962
1963private:
1966 boost::shared_ptr<MatrixDouble> inMat;
1967 boost::shared_ptr<MatrixDouble> outMat;
1968};
1969
1971
1974
1975 OpScaleMatrix(const std::string field_name, const double scale,
1976 boost::shared_ptr<MatrixDouble> in_mat,
1977 boost::shared_ptr<MatrixDouble> out_mat)
1978 : UserOp(field_name, OPROW), scale(scale), inMat(in_mat),
1979 outMat(out_mat) {
1980 // Only is run for vertices
1981 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1982 if (!inMat)
1983 THROW_MESSAGE("Pointer not set for in matrix");
1984 if (!outMat)
1985 THROW_MESSAGE("Pointer not set for in matrix");
1986 }
1987
1988 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1990 outMat->resize(inMat->size1(), inMat->size2(), false);
1991 noalias(*outMat) = scale * (*inMat);
1993 }
1994
1995private:
1996 const double scale;
1997 boost::shared_ptr<MatrixDouble> inMat;
1998 boost::shared_ptr<MatrixDouble> outMat;
1999};
2000
2001/**@}*/
2002
2003/** \name H-div/H-curls (Vectorial bases) values at integration points */
2004
2005/**@{*/
2006
2007/** \brief Get vector field for H-div approximation
2008 * \ingroup mofem_forces_and_sources_user_data_operators
2009 */
2010template <int Base_Dim, int Field_Dim, class T, class L, class A>
2012
2013/** \brief Get vector field for H-div approximation
2014 * \ingroup mofem_forces_and_sources_user_data_operators
2015 */
2016template <int Field_Dim>
2018 ublas::row_major, DoubleAllocator>
2020
2022 boost::shared_ptr<MatrixDouble> data_ptr,
2023 const EntityType zero_type = MBEDGE,
2024 const int zero_side = 0)
2027 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2028 if (!dataPtr)
2029 THROW_MESSAGE("Pointer is not set");
2030 }
2031
2032 /**
2033 * \brief Calculate values of vector field at integration points
2034 * @param side side entity number
2035 * @param type side entity type
2036 * @param data entity data
2037 * @return error code
2038 */
2039 MoFEMErrorCode doWork(int side, EntityType type,
2041
2042private:
2043 boost::shared_ptr<MatrixDouble> dataPtr;
2045 const int zeroSide;
2046};
2047
2048template <int Field_Dim>
2050 3, Field_Dim, double, ublas::row_major,
2051 DoubleAllocator>::doWork(int side, EntityType type,
2054 const size_t nb_integration_points = this->getGaussPts().size2();
2055 if (type == zeroType && side == zeroSide) {
2056 dataPtr->resize(Field_Dim, nb_integration_points, false);
2057 dataPtr->clear();
2058 }
2059 const size_t nb_dofs = data.getFieldData().size();
2060 if (!nb_dofs)
2062 const size_t nb_base_functions = data.getN().size2() / 3;
2064 auto t_n_hdiv = data.getFTensor1N<3>();
2065 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2066 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2067 auto t_dof = data.getFTensor0FieldData();
2068 int bb = 0;
2069 for (; bb != nb_dofs; ++bb) {
2070 t_data(i) += t_n_hdiv(i) * t_dof;
2071 ++t_n_hdiv;
2072 ++t_dof;
2073 }
2074 for (; bb != nb_base_functions; ++bb)
2075 ++t_n_hdiv;
2076 ++t_data;
2077 }
2079}
2080
2081/** \brief Get vector field for H-div approximation
2082 *
2083 * \ingroup mofem_forces_and_sources_user_data_operators
2084 */
2085template <int Base_Dim, int Field_Dim = Base_Dim>
2088 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2090 Base_Dim, Field_Dim, double, ublas::row_major,
2091 DoubleAllocator>::OpCalculateHVecVectorField_General;
2092};
2093
2094
2095
2096/** \brief Get vector field for H-div approximation
2097 * \ingroup mofem_forces_and_sources_user_data_operators
2098 */
2099template <int Base_Dim, int Field_Dim = Base_Dim>
2101
2102template <int Field_Dim>
2105
2107 boost::shared_ptr<MatrixDouble> data_ptr,
2108 const EntityType zero_type = MBEDGE,
2109 const int zero_side = 0)
2112 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2113 if (!dataPtr)
2114 THROW_MESSAGE("Pointer is not set");
2115 }
2116
2117 /**
2118 * \brief Calculate values of vector field at integration points
2119 * @param side side entity number
2120 * @param type side entity type
2121 * @param data entity data
2122 * @return error code
2123 */
2124 MoFEMErrorCode doWork(int side, EntityType type,
2126
2127private:
2128 boost::shared_ptr<MatrixDouble> dataPtr;
2130 const int zeroSide;
2131};
2132
2133template <int Field_Dim>
2135 int side, EntityType type, EntitiesFieldData::EntData &data) {
2137
2138 const size_t nb_integration_points = this->getGaussPts().size2();
2139 if (type == zeroType && side == zeroSide) {
2140 dataPtr->resize(Field_Dim, nb_integration_points, false);
2141 dataPtr->clear();
2142 }
2143
2144 auto &local_indices = data.getIndices();
2145 const size_t nb_dofs = local_indices.size();
2146 if (nb_dofs) {
2147
2148 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2149 const double *array;
2150 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2151 for (size_t i = 0; i != nb_dofs; ++i)
2152 if (local_indices[i] != -1)
2153 dot_dofs_vector[i] = array[local_indices[i]];
2154 else
2155 dot_dofs_vector[i] = 0;
2156 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2157
2158 const size_t nb_base_functions = data.getN().size2() / 3;
2160 auto t_n_hdiv = data.getFTensor1N<3>();
2161 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2162 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2163 int bb = 0;
2164 for (; bb != nb_dofs; ++bb) {
2165 t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
2166 ++t_n_hdiv;
2167 }
2168 for (; bb != nb_base_functions; ++bb)
2169 ++t_n_hdiv;
2170 ++t_data;
2171 }
2172 }
2173
2175}
2176
2177/**
2178 * @brief Calculate divergence of vector field
2179 * @ingroup mofem_forces_and_sources_user_data_operators
2180 *
2181 * @tparam BASE_DIM
2182 * @tparam SPACE_DIM
2183 */
2184template <int BASE_DIM, int SPACE_DIM>
2187
2189 boost::shared_ptr<VectorDouble> data_ptr,
2190 const EntityType zero_type = MBEDGE,
2191 const int zero_side = 0)
2194 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2195 if (!dataPtr)
2196 THROW_MESSAGE("Pointer is not set");
2197 }
2198
2202 const size_t nb_integration_points = getGaussPts().size2();
2203 if (type == zeroType && side == zeroSide) {
2204 dataPtr->resize(nb_integration_points, false);
2205 dataPtr->clear();
2206 }
2207 const size_t nb_dofs = data.getFieldData().size();
2208 if (!nb_dofs)
2210 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2213 auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2214 auto t_data = getFTensor0FromVec(*dataPtr);
2215 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2216 auto t_dof = data.getFTensor0FieldData();
2217 int bb = 0;
2218 for (; bb != nb_dofs; ++bb) {
2219 t_data += t_dof * t_n_diff_hdiv(j, j);
2220 ++t_n_diff_hdiv;
2221 ++t_dof;
2222 }
2223 for (; bb != nb_base_functions; ++bb)
2224 ++t_n_diff_hdiv;
2225 ++t_data;
2226 }
2228 }
2229
2230private:
2231 boost::shared_ptr<VectorDouble> dataPtr;
2233 const int zeroSide;
2234};
2235
2236/**
2237 * @brief Calculate gradient of vector field
2238 * @ingroup mofem_forces_and_sources_user_data_operators
2239 *
2240 * @tparam BASE_DIM
2241 * @tparam SPACE_DIM
2242 */
2243template <int BASE_DIM, int SPACE_DIM>
2246
2248 boost::shared_ptr<MatrixDouble> data_ptr,
2249 const EntityType zero_type = MBEDGE,
2250 const int zero_side = 0)
2253 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2254 if (!dataPtr)
2255 THROW_MESSAGE("Pointer is not set");
2256 }
2257
2261 const size_t nb_integration_points = getGaussPts().size2();
2262 if (type == zeroType && side == zeroSide) {
2263 dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
2264 dataPtr->clear();
2265 }
2266 const size_t nb_dofs = data.getFieldData().size();
2267 if (!nb_dofs)
2269 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2272 auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2273 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*dataPtr);
2274 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2275 auto t_dof = data.getFTensor0FieldData();
2276 int bb = 0;
2277 for (; bb != nb_dofs; ++bb) {
2278 t_data(i, j) += t_dof * t_base_diff(i, j);
2279 ++t_base_diff;
2280 ++t_dof;
2281 }
2282 for (; bb != nb_base_functions; ++bb)
2283 ++t_base_diff;
2284 ++t_data;
2285 }
2287 }
2288
2289private:
2290 boost::shared_ptr<MatrixDouble> dataPtr;
2292 const int zeroSide;
2293};
2294
2295/**
2296 * @brief Calculate gradient of vector field
2297 * @ingroup mofem_forces_and_sources_user_data_operators
2298 *
2299 * @tparam BASE_DIM
2300 * @tparam SPACE_DIM
2301 */
2302template <int BASE_DIM, int SPACE_DIM>
2305
2307 boost::shared_ptr<MatrixDouble> data_ptr,
2308 const EntityType zero_type = MBEDGE,
2309 const int zero_side = 0)
2312 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2313 if (!dataPtr)
2314 THROW_MESSAGE("Pointer is not set");
2315 }
2316
2320 const size_t nb_integration_points = getGaussPts().size2();
2321 if (type == zeroType && side == zeroSide) {
2322 dataPtr->resize(BASE_DIM * SPACE_DIM * SPACE_DIM, nb_integration_points,
2323 false);
2324 dataPtr->clear();
2325 }
2326 const size_t nb_dofs = data.getFieldData().size();
2327 if (!nb_dofs)
2329
2330 const int nb_base_functions = data.getN().size2() / BASE_DIM;
2331
2332#ifndef NDEBUG
2333 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2334 if (hessian_base.size1() != nb_integration_points) {
2335 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2336 "Wrong number of integration pts (%d != %d)",
2337 hessian_base.size1(), nb_integration_points);
2338 }
2339 if (hessian_base.size2() !=
2340 BASE_DIM * nb_base_functions * SPACE_DIM * SPACE_DIM) {
2341 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2342 "Wrong number of base functions (%d != %d)",
2343 hessian_base.size2() / (BASE_DIM * SPACE_DIM * SPACE_DIM),
2344 nb_base_functions);
2345 }
2346 if (hessian_base.size2() < BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM) {
2347 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2348 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2349 BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM);
2350 }
2351#endif
2352
2356
2357 auto t_base_diff2 =
2359 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*dataPtr);
2360 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2361 auto t_dof = data.getFTensor0FieldData();
2362 int bb = 0;
2363 for (; bb != nb_dofs; ++bb) {
2364 t_data(i, j, k) += t_dof * t_base_diff2(i, j, k);
2365
2366 ++t_base_diff2;
2367 ++t_dof;
2368 }
2369 for (; bb != nb_base_functions; ++bb)
2370 ++t_base_diff2;
2371 ++t_data;
2372 }
2374 }
2375
2376private:
2377 boost::shared_ptr<MatrixDouble> dataPtr;
2379 const int zeroSide;
2380};
2381
2382/**
2383 * @brief Calculate divergence of vector field dot
2384 * @ingroup mofem_forces_and_sources_user_data_operators
2385 *
2386 * @tparam Tensor_Dim dimension of space
2387 */
2388template <int Tensor_Dim1, int Tensor_Dim2>
2391
2393 boost::shared_ptr<VectorDouble> data_ptr,
2394 const EntityType zero_type = MBEDGE,
2395 const int zero_side = 0)
2398 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2399 if (!dataPtr)
2400 THROW_MESSAGE("Pointer is not set");
2401 }
2402
2406 const size_t nb_integration_points = getGaussPts().size2();
2407 if (type == zeroType && side == zeroSide) {
2408 dataPtr->resize(nb_integration_points, false);
2409 dataPtr->clear();
2410 }
2411
2412 const auto &local_indices = data.getLocalIndices();
2413 const int nb_dofs = local_indices.size();
2414 if (nb_dofs) {
2415
2416 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2417 const double *array;
2418 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2419 for (size_t i = 0; i != local_indices.size(); ++i)
2420 if (local_indices[i] != -1)
2421 dot_dofs_vector[i] = array[local_indices[i]];
2422 else
2423 dot_dofs_vector[i] = 0;
2424 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2425
2426 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
2428 auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
2429 auto t_data = getFTensor0FromVec(*dataPtr);
2430 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2431 int bb = 0;
2432 for (; bb != nb_dofs; ++bb) {
2433 double div = 0;
2434 for (auto ii = 0; ii != Tensor_Dim2; ++ii)
2435 div += t_n_diff_hdiv(ii, ii);
2436 t_data += dot_dofs_vector[bb] * div;
2437 ++t_n_diff_hdiv;
2438 }
2439 for (; bb != nb_base_functions; ++bb)
2440 ++t_n_diff_hdiv;
2441 ++t_data;
2442 }
2443 }
2445 }
2446
2447private:
2448 boost::shared_ptr<VectorDouble> dataPtr;
2450 const int zeroSide;
2451};
2452
2453/**
2454 * @brief Calculate curl of vector field
2455 * @ingroup mofem_forces_and_sources_user_data_operators
2456 *
2457 * @tparam Base_Dim base function dimension
2458 * @tparam Space_Dim dimension of space
2459 * @tparam Hcurl field dimension
2460 */
2461template <int Base_Dim, int Space_Dim>
2463
2464/**
2465 * @brief Calculate curl of vector field
2466 * @ingroup mofem_forces_and_sources_user_data_operators
2467 *
2468 * @tparam Base_Dim base function dimension
2469 * @tparam Space_Dim dimension of space
2470 * @tparam Hcurl field dimension
2471 */
2472template <>
2475 OpCalculateHcurlVectorCurl(const std::string field_name,
2476 boost::shared_ptr<MatrixDouble> data_ptr,
2477 const EntityType zero_type = MBEDGE,
2478 const int zero_side = 0);
2479 MoFEMErrorCode doWork(int side, EntityType type,
2481
2482private:
2483 boost::shared_ptr<MatrixDouble> dataPtr;
2485 const int zeroSide;
2486};
2487
2488/**
2489 * @brief Calculate curl of vector field
2490 * @ingroup mofem_forces_and_sources_user_data_operators
2491 *
2492 * @tparam Field_Dim dimension of field
2493 * @tparam Space_Dim dimension of space
2494 */
2495template <>
2498
2499 OpCalculateHcurlVectorCurl(const std::string field_name,
2500 boost::shared_ptr<MatrixDouble> data_ptr,
2501 const EntityType zero_type = MBVERTEX,
2502 const int zero_side = 0);
2503
2504 MoFEMErrorCode doWork(int side, EntityType type,
2506
2507private:
2508 boost::shared_ptr<MatrixDouble> dataPtr;
2510 const int zeroSide;
2511};
2512
2513/**
2514 * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
2515 * \ingroup mofem_forces_and_sources_user_data_operators
2516 *
2517 * @tparam Tensor_Dim0 rank of the filed
2518 * @tparam Tensor_Dim1 dimension of space
2519 */
2520template <int Tensor_Dim0, int Tensor_Dim1>
2523
2525 boost::shared_ptr<MatrixDouble> data_ptr,
2526 const EntityType zero_type = MBEDGE,
2527 const int zero_side = 0)
2530 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2531 if (!dataPtr)
2532 THROW_MESSAGE("Pointer is not set");
2533 }
2534
2538 const size_t nb_integration_points = getGaussPts().size2();
2539 if (type == zeroType && side == zeroSide) {
2540 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2541 dataPtr->clear();
2542 }
2543 const size_t nb_dofs = data.getFieldData().size();
2544 if (nb_dofs) {
2545 const size_t nb_base_functions = data.getN().size2() / 3;
2548 auto t_n_hvec = data.getFTensor1N<3>();
2549 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2550 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2551 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2552 size_t bb = 0;
2553 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2554 t_data(i, j) += t_dof(i) * t_n_hvec(j);
2555 ++t_n_hvec;
2556 ++t_dof;
2557 }
2558 for (; bb < nb_base_functions; ++bb)
2559 ++t_n_hvec;
2560 ++t_data;
2561 }
2562 }
2564 }
2565
2566private:
2567 boost::shared_ptr<MatrixDouble> dataPtr;
2569 const int zeroSide;
2570};
2571
2572/**
2573 * @brief Calculate tenor field using tensor base, i.e. Hdiv/Hcurl
2574 * \ingroup mofem_forces_and_sources_user_data_operators
2575 *
2576 * @tparam Tensor_Dim0 rank of the filed
2577 * @tparam Tensor_Dim1 dimension of space
2578 */
2579template <int Tensor_Dim0, int Tensor_Dim1>
2582
2584 boost::shared_ptr<MatrixDouble> data_ptr,
2585 const EntityType zero_type = MBEDGE,
2586 const int zero_side = 0)
2589 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2590 if (!dataPtr)
2591 THROW_MESSAGE("Pointer is not set");
2592 }
2593
2597 const size_t nb_integration_points = getGaussPts().size2();
2598 if (type == zeroType && side == zeroSide) {
2599 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2600 dataPtr->clear();
2601 }
2602 const size_t nb_dofs = data.getFieldData().size();
2603 if (!nb_dofs)
2605 const size_t nb_base_functions =
2606 data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2609 auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2610 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2611 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2612 auto t_dof = data.getFTensor0FieldData();
2613 size_t bb = 0;
2614 for (; bb != nb_dofs; ++bb) {
2615 t_data(i, j) += t_dof * t_n_hten(i, j);
2616 ++t_n_hten;
2617 ++t_dof;
2618 }
2619 for (; bb < nb_base_functions; ++bb)
2620 ++t_n_hten;
2621 ++t_data;
2622 }
2624 }
2625
2626private:
2627 boost::shared_ptr<MatrixDouble> dataPtr;
2629 const int zeroSide;
2630};
2631
2632/**
2633 * @brief Calculate divergence of tonsorial field using vectorial base
2634 * \ingroup mofem_forces_and_sources_user_data_operators
2635 *
2636 * @tparam Tensor_Dim0 rank of the field
2637 * @tparam Tensor_Dim1 dimension of space
2638 */
2639template <int Tensor_Dim0, int Tensor_Dim1,
2640 CoordinateTypes CoordSys = CARTESIAN>
2643
2645 boost::shared_ptr<MatrixDouble> data_ptr,
2646 const EntityType zero_type = MBEDGE,
2647 const int zero_side = 0)
2650 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2651 if (!dataPtr)
2652 THROW_MESSAGE("Pointer is not set");
2653 }
2654
2658 const size_t nb_integration_points = getGaussPts().size2();
2659 if (type == zeroType && side == 0) {
2660 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
2661 dataPtr->clear();
2662 }
2663 const size_t nb_dofs = data.getFieldData().size();
2664 if (nb_dofs) {
2665 const size_t nb_base_functions = data.getN().size2() / 3;
2668 auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
2669 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
2670 auto t_base = data.getFTensor1N<3>();
2671 auto t_coords = getFTensor1CoordsAtGaussPts();
2672 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2673 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2674 size_t bb = 0;
2675 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2676 double div = t_n_diff_hvec(j, j);
2677 t_data(i) += t_dof(i) * div;
2678 if constexpr (CoordSys == CYLINDRICAL) {
2679 t_data(i) += t_dof(i) * (t_base(0) / t_coords(0));
2680 }
2681 ++t_n_diff_hvec;
2682 ++t_dof;
2683 ++t_base;
2684 }
2685 for (; bb < nb_base_functions; ++bb) {
2686 ++t_base;
2687 ++t_n_diff_hvec;
2688 }
2689 ++t_data;
2690 ++t_coords;
2691 }
2692 }
2694 }
2695
2696private:
2697 boost::shared_ptr<MatrixDouble> dataPtr;
2699 const int zeroSide;
2700};
2701
2702/**
2703 * @brief Calculate trace of vector (Hdiv/Hcurl) space
2704 *
2705 * @tparam Tensor_Dim
2706 * @tparam OpBase
2707 */
2708template <int Tensor_Dim, typename OpBase>
2710
2712 boost::shared_ptr<MatrixDouble> data_ptr,
2713 const EntityType zero_type = MBEDGE,
2714 const int zero_side = 0)
2715 : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
2716 zeroType(zero_type), zeroSide(zero_side) {
2717 if (!dataPtr)
2718 THROW_MESSAGE("Pointer is not set");
2719 }
2720
2724 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2725 if (type == zeroType && side == 0) {
2726 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2727 dataPtr->clear();
2728 }
2729 const size_t nb_dofs = data.getFieldData().size();
2730 if (nb_dofs) {
2731 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
2732 const size_t nb_base_functions = data.getN().size2() / 3;
2733 auto t_base = data.getFTensor1N<3>();
2734 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2735 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2736 FTensor::Tensor1<double, Tensor_Dim> t_normalized_normal;
2737 t_normalized_normal(j) = t_normal(j);
2738 t_normalized_normal.normalize();
2739 auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
2740 size_t bb = 0;
2741 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2742 t_data(i) += t_dof(i) * (t_base(j) * t_normalized_normal(j));
2743 ++t_base;
2744 ++t_dof;
2745 }
2746 for (; bb < nb_base_functions; ++bb) {
2747 ++t_base;
2748 }
2749 ++t_data;
2750 ++t_normal;
2751 }
2752 }
2754 }
2755
2756private:
2757 boost::shared_ptr<MatrixDouble> dataPtr;
2759 const int zeroSide;
2762};
2763
2764/**@}*/
2765
2766/** \name Other operators */
2767
2768/**@{*/
2769
2770/**@}*/
2771
2772/** \name Operators for faces */
2773
2774/**@{*/
2775
2776/** \brief Transform local reference derivatives of shape functions to global
2777derivatives
2778
2779\ingroup mofem_forces_and_sources_tri_element
2780
2781*/
2782template <int DIM, int DERIVATIVE = 1> struct OpSetInvJacSpaceForFaceImpl;
2783
2786
2788 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2789
2790protected:
2791 template <int D1, int D2, int J1, int J2>
2794
2795 static_assert(D2 == J2, "Dimension of jacobian and dimension of <out> "
2796 "directive does not match");
2797
2798 size_t nb_functions = diff_n.size2() / D1;
2799 if (nb_functions) {
2800 size_t nb_gauss_pts = diff_n.size1();
2801
2802#ifndef NDEBUG
2803 if (nb_gauss_pts != getGaussPts().size2())
2804 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2805 "Wrong number of Gauss Pts");
2806 if (diff_n.size2() % D1)
2807 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2808 "Number of direvatives of base functions and D1 dimension does "
2809 "not match");
2810#endif
2811
2812 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions, false);
2813
2816 auto t_diff_n = getFTensor1FromPtr<D2>(&*diffNinvJac.data().begin());
2817 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2818 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*invJacPtr);
2819 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2820 for (size_t dd = 0; dd != nb_functions; ++dd) {
2821 t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
2822 ++t_diff_n;
2823 ++t_diff_n_ref;
2824 }
2825 }
2826
2827 diff_n.swap(diffNinvJac);
2828 }
2830 }
2831
2832 boost::shared_ptr<MatrixDouble> invJacPtr;
2834};
2835
2836template <>
2839
2841
2842 MoFEMErrorCode doWork(int side, EntityType type,
2844};
2845
2846template <>
2848 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2849
2850 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2851
2852 MoFEMErrorCode doWork(int side, EntityType type,
2854};
2855
2856template <>
2858 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2859
2860 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2861
2862 MoFEMErrorCode doWork(int side, EntityType type,
2864};
2865
2866template <int DERIVARIVE = 1>
2868 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2869 OpSetInvJacH1ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2870 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(H1, inv_jac_ptr) {}
2871};
2872
2873template <int DERIVARIVE = 1>
2875 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2876 OpSetInvJacL2ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2877 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(L2, inv_jac_ptr) {}
2878};
2879
2880template <int DERIVARIVE = 1>
2882 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2884 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2885 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(H1, inv_jac_ptr) {}
2886};
2887
2888template <int DERIVARIVE = 1>
2890 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2892 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2893 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(L2, inv_jac_ptr) {}
2894};
2895
2896/**
2897 * \brief Transform local reference derivatives of shape function to
2898 global derivatives for face
2899
2900 * \ingroup mofem_forces_and_sources_tri_element
2901 */
2902template <int DIM> struct OpSetInvJacHcurlFaceImpl;
2903
2904template <>
2907
2908 OpSetInvJacHcurlFaceImpl(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2910 invJacPtr(inv_jac_ptr) {}
2911
2912 MoFEMErrorCode doWork(int side, EntityType type,
2914
2915protected:
2916 boost::shared_ptr<MatrixDouble> invJacPtr;
2918};
2919
2920template <>
2922 using OpSetInvJacHcurlFaceImpl<2>::OpSetInvJacHcurlFaceImpl;
2923 MoFEMErrorCode doWork(int side, EntityType type,
2925};
2926
2929
2930/**
2931 * @brief Make Hdiv space from Hcurl space in 2d
2932 * @ingroup mofem_forces_and_sources_tri_element
2933 */
2936
2939
2940 MoFEMErrorCode doWork(int side, EntityType type,
2942};
2943
2944/** \brief Transform Hcurl base fluxes from reference element to physical
2945 * triangle
2946 */
2948
2949/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
2950 *
2951 * Covariant Piola transformation
2952 \f[
2953 \psi_i|_t = J^{-1}_{ij}\hat{\psi}_j\\
2954 \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
2955 = J^{-1}_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
2956 \f]
2957
2958 */
2959template <>
2962
2964 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2965
2966 MoFEMErrorCode doWork(int side, EntityType type,
2968
2969private:
2970 boost::shared_ptr<MatrixDouble> invJacPtr;
2971
2974};
2975
2978
2979/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
2980 *
2981 * \note Hdiv space is generated by Hcurl space in 2d.
2982 *
2983 * Contravariant Piola transformation
2984 * \f[
2985 * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
2986 * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
2987 * =
2988 * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
2989 * \f]
2990 *
2991 * \ingroup mofem_forces_and_sources
2992 *
2993 */
2995
2996template <>
2999
3001 boost::shared_ptr<MatrixDouble> jac_ptr)
3003 jacPtr(jac_ptr) {}
3004
3005 MoFEMErrorCode doWork(int side, EntityType type,
3007
3008protected:
3009 boost::shared_ptr<MatrixDouble> jacPtr;
3012};
3013
3014template <>
3018 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3019
3020 MoFEMErrorCode doWork(int side, EntityType type,
3022};
3023
3028
3029/**@}*/
3030
3031/** \name Operators for edges */
3032
3033/**@{*/
3034
3037
3040
3041 MoFEMErrorCode doWork(int side, EntityType type,
3043
3044};
3045
3046/**
3047 * @deprecated Name is deprecated and this is added for back compatibility
3048 */
3051
3052/**@}*/
3053
3054/** \name Operator for fat prisms */
3055
3056/**@{*/
3057
3058/**
3059 * @brief Operator for fat prism element updating integration weights in the
3060 * volume.
3061 *
3062 * Jacobian on the distorted element is nonconstant. This operator updates
3063 * integration weight on prism to take into account nonconstat jacobian.
3064 *
3065 * \f[
3066 * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
3067 * \pmb\xi} \right\| \right)
3068 * \f]
3069 * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
3070 * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
3071 * element coordinate.
3072 *
3073 */
3076
3079
3080 MoFEMErrorCode doWork(int side, EntityType type,
3082};
3083
3084/** \brief Calculate inverse of jacobian for face element
3085
3086 It is assumed that face element is XY plane. Applied
3087 only for 2d problems.
3088
3089 FIXME Generalize function for arbitrary face orientation in 3d space
3090 FIXME Calculate to Jacobins for two faces
3091
3092 \ingroup mofem_forces_and_sources_prism_element
3093
3094*/
3097
3098 OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3100 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3101
3104 invJac(inv_jac) {}
3105 MoFEMErrorCode doWork(int side, EntityType type,
3107
3108private:
3109 const boost::shared_ptr<MatrixDouble> invJacPtr;
3111};
3112
3113/** \brief Transform local reference derivatives of shape functions to global
3114derivatives
3115
3116FIXME Generalize to curved shapes
3117FIXME Generalize to case that top and bottom face has different shape
3118
3119\ingroup mofem_forces_and_sources_prism_element
3120
3121*/
3124
3125 OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3127 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3128
3131 invJac(inv_jac) {}
3132
3133 MoFEMErrorCode doWork(int side, EntityType type,
3135
3136private:
3137 const boost::shared_ptr<MatrixDouble> invJacPtr;
3140};
3141
3142// Flat prism
3143
3144/** \brief Calculate inverse of jacobian for face element
3145
3146 It is assumed that face element is XY plane. Applied
3147 only for 2d problems.
3148
3149 FIXME Generalize function for arbitrary face orientation in 3d space
3150 FIXME Calculate to Jacobins for two faces
3151
3152 \ingroup mofem_forces_and_sources_prism_element
3153
3154*/
3157
3160 invJacF3(inv_jac_f3) {}
3161 MoFEMErrorCode doWork(int side, EntityType type,
3163
3164private:
3166};
3167
3168/** \brief Transform local reference derivatives of shape functions to global
3169derivatives
3170
3171FIXME Generalize to curved shapes
3172FIXME Generalize to case that top and bottom face has different shape
3173
3174\ingroup mofem_forces_and_sources_prism_element
3175
3176*/
3179
3182 invJacF3(inv_jac_f3) {}
3183
3184 MoFEMErrorCode doWork(int side, EntityType type,
3186
3187private:
3190};
3191
3192/**@}*/
3193
3194/** \name Operation on matrices at integration points */
3195
3196/**@{*/
3197
3198template <int DIM>
3200
3201 OpInvertMatrix(boost::shared_ptr<MatrixDouble> in_ptr,
3202 boost::shared_ptr<VectorDouble> det_ptr,
3203 boost::shared_ptr<MatrixDouble> out_ptr)
3205 outPtr(out_ptr), detPtr(det_ptr) {}
3206
3209 return doWorkImpl(side, type, data, FTensor::Number<DIM>());
3210 }
3211
3212private:
3213 boost::shared_ptr<MatrixDouble> inPtr;
3214 boost::shared_ptr<MatrixDouble> outPtr;
3215 boost::shared_ptr<VectorDouble> detPtr;
3216
3217 MoFEMErrorCode doWorkImpl(int side, EntityType type,
3219 const FTensor::Number<3> &);
3220
3221 MoFEMErrorCode doWorkImpl(int side, EntityType type,
3223 const FTensor::Number<2> &);
3224};
3225
3226template <int DIM>
3229 const FTensor::Number<3> &) {
3231
3232 if (!inPtr)
3233 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3234 "Pointer for inPtr matrix not allocated");
3235 if (!detPtr)
3236 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3237 "Pointer for detPtr matrix not allocated");
3238
3239 const auto nb_rows = inPtr->size1();
3240 const auto nb_integration_pts = inPtr->size2();
3241
3242 // Calculate determinant
3243 {
3244 detPtr->resize(nb_integration_pts, false);
3245 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3246 auto t_det = getFTensor0FromVec(*detPtr);
3247 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3248 determinantTensor3by3(t_in, t_det);
3249 ++t_in;
3250 ++t_det;
3251 }
3252 }
3253
3254 // Invert jacobian
3255 if (outPtr) {
3256 outPtr->resize(nb_rows, nb_integration_pts, false);
3257 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3258 auto t_out = getFTensor2FromMat<3, 3>(*outPtr);
3259 auto t_det = getFTensor0FromVec(*detPtr);
3260 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3261 invertTensor3by3(t_in, t_det, t_out);
3262 ++t_in;
3263 ++t_out;
3264 ++t_det;
3265 }
3266 }
3267
3269}
3270
3271template <int DIM>
3274 const FTensor::Number<2> &) {
3276
3277 if (!inPtr)
3278 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3279 "Pointer for inPtr matrix not allocated");
3280 if (!detPtr)
3281 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3282 "Pointer for detPtr matrix not allocated");
3283
3284 const auto nb_rows = inPtr->size1();
3285 const auto nb_integration_pts = inPtr->size2();
3286
3287 // Calculate determinant
3288 {
3289 detPtr->resize(nb_integration_pts, false);
3290 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3291 auto t_det = getFTensor0FromVec(*detPtr);
3292 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3293 determinantTensor2by2(t_in, t_det);
3294 ++t_in;
3295 ++t_det;
3296 }
3297 }
3298
3299 // Invert jacobian
3300 if (outPtr) {
3301 outPtr->resize(nb_rows, nb_integration_pts, false);
3302 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3303 auto t_out = getFTensor2FromMat<2, 2>(*outPtr);
3304 auto t_det = getFTensor0FromVec(*detPtr);
3305 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3306 invertTensor2by2(t_in, t_det, t_out);
3307 ++t_in;
3308 ++t_out;
3309 ++t_det;
3310 }
3311 }
3312
3314}
3315
3316/**@}*/
3317
3318} // namespace MoFEM
3319
3320#endif // __USER_DATA_OPERATORS_HPP__
3321
3322/**
3323 * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
3324 *
3325 * \brief Classes and functions used to evaluate fields at integration pts,
3326 *jacobians, etc..
3327 *
3328 * \ingroup mofem_forces_and_sources
3329 **/
static Index< 'J', 3 > J
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
Tensor1< T, Tensor_Dim > normalize()
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
CoordinateTypes
Coodinate system.
Definition: definitions.h:114
@ 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
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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:1516
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1547
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:1486
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:1475
constexpr IntegrationType I
constexpr AssemblyType A
constexpr auto field_name
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3Diff2N(FieldApproximationBase base)
Get second derivatives of base functions for Hvec space.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(Tensor_Dim *(Tensor_Dim+1))/2 >, Tensor_Dim > getFTensor2SymmetricFieldData()
Return symmetric FTensor rank 2, i.e. matrix from filed data coeffinects.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from filed data coeffinects.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
boost::shared_ptr< VectorDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateDivergenceVectorFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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.
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)
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Calculate curl of vector field.
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)
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, 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 values 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.
Apply contravariant (Piola) transfer to Hdiv space on face.
Transform Hcurl base fluxes from reference element to physical triangle.
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