v0.13.2
Loading...
Searching...
No Matches
UserDataOperators.hpp
Go to the documentation of this file.
1/** \file UserDataOperators.hpp
2 * \brief User data Operators
3
4*/
5
6#ifndef __USER_DATA_OPERATORS_HPP__
7#define __USER_DATA_OPERATORS_HPP__
8
9namespace MoFEM {
10
11/** \name Get values at Gauss pts */
12
13/**@{*/
14
15/** \name Scalar values */
16
17/**@{*/
18
19/** \brief Scalar field values at integration points
20 *
21 */
22template <class T, class A>
25
27 const std::string field_name,
28 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
29 const EntityType zero_type = MBVERTEX)
31 dataPtr(data_ptr), zeroType(zero_type) {
32 if (!dataPtr)
33 THROW_MESSAGE("Pointer is not set");
34 }
35
37 const std::string field_name,
38 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
39 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
41 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
42 if (!dataPtr)
43 THROW_MESSAGE("Pointer is not set");
44 }
45
46 /**
47 * \brief calculate values of scalar field at integration points
48 * @param side side entity number
49 * @param type side entity type
50 * @param data entity data
51 * @return error code
52 */
55
56protected:
57 boost::shared_ptr<ublas::vector<T, A>> dataPtr;
61};
62
63/**
64 * \brief Specialization of member function
65 *
66 */
67template <class T, class A>
69 int side, EntityType type, EntitiesFieldData::EntData &data) {
71 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented for T = %s",
72 typeid(T).name() // boost::core::demangle(typeid(T).name()).c_str()
73 );
75}
76
77/**
78 * \brief Get value at integration points for scalar field
79 *
80 * \ingroup mofem_forces_and_sources_user_data_operators
81 */
83 : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
84
86 double, DoubleAllocator>::OpCalculateScalarFieldValues_General;
87
88 /**
89 * \brief calculate values of scalar field at integration points
90 * @param side side entity number
91 * @param type side entity type
92 * @param data entity data
93 * @return error code
94 */
98 VectorDouble &vec = *dataPtr;
99 const size_t nb_gauss_pts = getGaussPts().size2();
100 if (type == zeroType) {
101 vec.resize(nb_gauss_pts, false);
102 vec.clear();
103 }
104
105 const size_t nb_dofs = data.getFieldData().size();
106
107 if (nb_dofs) {
108
109 if (dataVec.use_count()) {
110 dotVector.resize(nb_dofs, false);
111 const double *array;
112 CHKERR VecGetArrayRead(dataVec, &array);
113 const auto &local_indices = data.getLocalIndices();
114 for (int i = 0; i != local_indices.size(); ++i)
115 if (local_indices[i] != -1)
116 dotVector[i] = array[local_indices[i]];
117 else
118 dotVector[i] = 0;
119 CHKERR VecRestoreArrayRead(dataVec, &array);
120 data.getFieldData().swap(dotVector);
121 }
122
123 const size_t nb_base_functions = data.getN().size2();
124 auto base_function = data.getFTensor0N();
125 auto values_at_gauss_pts = getFTensor0FromVec(vec);
126 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
127 auto field_data = data.getFTensor0FieldData();
128 size_t bb = 0;
129 for (; bb != nb_dofs; ++bb) {
130 values_at_gauss_pts += field_data * base_function;
131 ++field_data;
132 ++base_function;
133 }
134 // It is possible to have more base functions than dofs
135 for (; bb < nb_base_functions; ++bb)
136 ++base_function;
137 ++values_at_gauss_pts;
138 }
139
140 if (dataVec.use_count()) {
141 data.getFieldData().swap(dotVector);
142 }
143 }
144
146 }
147};
148
149/**
150 * @brief Get rate of scalar field at integration points
151 *
152 * \ingroup mofem_forces_and_sources_user_data_operators
153 */
154
155template <PetscData::DataContext CTX>
158
160 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
161 const EntityType zero_at_type = MBVERTEX)
164 dataPtr(data_ptr), zeroAtType(zero_at_type) {
165 if (!dataPtr)
166 THROW_MESSAGE("Pointer is not set");
167 }
168
172
173 const size_t nb_gauss_pts = getGaussPts().size2();
174
175 VectorDouble &vec = *dataPtr;
176 if (type == zeroAtType) {
177 vec.resize(nb_gauss_pts, false);
178 vec.clear();
179 }
180
181 auto &local_indices = data.getLocalIndices();
182 const size_t nb_dofs = local_indices.size();
183 if (nb_dofs) {
184
185 const double *array;
186
187 auto get_array = [&](const auto ctx, auto vec) {
189#ifndef NDEBUG
190 if ((getFEMethod()->data_ctx & ctx).none()) {
191 MOFEM_LOG_CHANNEL("SELF");
192 MOFEM_LOG("SELF", Sev::error)
193 << "In this case filed degrees of freedom are read from vector. "
194 "That usually happens when time solver is used, and acces to "
195 "first or second rates is needed. You probably not set ts_u, "
196 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
197 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
198 "respectively";
199 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
200 }
201#endif
202 CHKERR VecGetArrayRead(vec, &array);
204 };
205
206 auto restore_array = [&](auto vec) {
207 return VecRestoreArrayRead(vec, &array);
208 };
209
210 switch (CTX) {
212 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
213 break;
215 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
216 break;
218 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
219 break;
220 default:
221 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
222 "That case is not implemented");
223 }
224
225 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
226 for (int i = 0; i != local_indices.size(); ++i)
227 if (local_indices[i] != -1)
228 dot_dofs_vector[i] = array[local_indices[i]];
229 else
230 dot_dofs_vector[i] = 0;
231
232 switch (CTX) {
234 CHKERR restore_array(getFEMethod()->ts_u);
235 break;
237 CHKERR restore_array(getFEMethod()->ts_u_t);
238 break;
240 CHKERR restore_array(getFEMethod()->ts_u_tt);
241 break;
242 default:
243 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
244 "That case is not implemented");
245 }
246
247 const size_t nb_base_functions = data.getN().size2();
248 auto base_function = data.getFTensor0N();
249 auto values_at_gauss_pts = getFTensor0FromVec(vec);
250
251 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
252 size_t bb = 0;
253 for (; bb != nb_dofs; ++bb) {
254 values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
255 ++base_function;
256 }
257 // Number of dofs can be smaller than number of Tensor_Dim x base
258 // functions
259 for (; bb < nb_base_functions; ++bb)
260 ++base_function;
261 ++values_at_gauss_pts;
262 }
263 }
265 }
266
267private:
268 boost::shared_ptr<VectorDouble> dataPtr;
270};
271
276
277/**
278 * \deprecated Name inconsistent with other operators
279 *
280 */
282
283/**@}*/
284
285/** \name Vector field values at integration points */
286
287/**@{*/
288
289/** \brief Calculate field values for tenor field rank 1, i.e. vector field
290 *
291 */
292template <int Tensor_Dim, class T, class L, class A>
295
297 const std::string field_name,
298 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
299 const EntityType zero_type = MBVERTEX)
302 dataPtr(data_ptr), zeroType(zero_type) {
303 if (!dataPtr)
304 THROW_MESSAGE("Pointer is not set");
305 }
306
307 /**
308 * \brief calculate values of vector field at integration points
309 * @param side side entity number
310 * @param type side entity type
311 * @param data entity data
312 * @return error code
313 */
314 MoFEMErrorCode doWork(int side, EntityType type,
316
317protected:
318 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
320};
321
322template <int Tensor_Dim, class T, class L, class A>
325 int side, EntityType type, EntitiesFieldData::EntData &data) {
327 SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
328 "Not implemented for T = %s and dim = %d",
329 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
330 Tensor_Dim);
332}
333
334/** \brief Calculate field values (template specialization) for tensor field
335 * rank 1, i.e. vector field
336 *
337 */
338template <int Tensor_Dim>
340 ublas::row_major, DoubleAllocator>
342
344 boost::shared_ptr<MatrixDouble> data_ptr,
345 const EntityType zero_type = MBVERTEX)
348 dataPtr(data_ptr), zeroType(zero_type) {
349 if (!dataPtr)
350 THROW_MESSAGE("Pointer is not set");
351 }
352
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) {
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
965/**
966 * @brief Calculate symmetric tensor field values at integration pts.
967 *
968 * @tparam Tensor_Dim
969
970 * \ingroup mofem_forces_and_sources_user_data_operators
971 */
972template <int Tensor_Dim>
975
977 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
978 const EntityType zero_type = MBEDGE, const int zero_side = 0)
981 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
982 if (!dataPtr)
983 THROW_MESSAGE("Pointer is not set");
984 }
985
987 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
988 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBEDGE,
989 const int zero_side = 0)
992 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side),
993 dataVec(data_vec) {
994 if (!dataPtr)
995 THROW_MESSAGE("Pointer is not set");
996 }
997
1001 MatrixDouble &mat = *dataPtr;
1002 const int nb_gauss_pts = getGaussPts().size2();
1003 if (type == this->zeroType && side == zeroSide) {
1004 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
1005 mat.clear();
1006 }
1007 const int nb_dofs = data.getFieldData().size();
1008 if (!nb_dofs)
1010
1011 if (dataVec.use_count()) {
1012 dotVector.resize(nb_dofs, false);
1013 const double *array;
1014 CHKERR VecGetArrayRead(dataVec, &array);
1015 const auto &local_indices = data.getLocalIndices();
1016 for (int i = 0; i != local_indices.size(); ++i)
1017 if (local_indices[i] != -1)
1018 dotVector[i] = array[local_indices[i]];
1019 else
1020 dotVector[i] = 0;
1021 CHKERR VecRestoreArrayRead(dataVec, &array);
1022 data.getFieldData().swap(dotVector);
1023 }
1024
1025 const int nb_base_functions = data.getN().size2();
1026 auto base_function = data.getFTensor0N();
1027 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1030 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1031 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1032 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
1033 int bb = 0;
1034 for (; bb != size; ++bb) {
1035 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1036 ++field_data;
1037 ++base_function;
1038 }
1039 for (; bb < nb_base_functions; ++bb)
1040 ++base_function;
1041 ++values_at_gauss_pts;
1042 }
1043
1044 if (dataVec.use_count()) {
1045 data.getFieldData().swap(dotVector);
1046 }
1047
1049 }
1050
1051protected:
1052 boost::shared_ptr<MatrixDouble> dataPtr;
1054 const int zeroSide;
1057};
1058
1059/**
1060 * @brief Calculate symmetric tensor field rates ant integratio pts.
1061 *
1062 * @tparam Tensor_Dim
1063 *
1064 * \ingroup mofem_forces_and_sources_user_data_operators
1065 */
1066template <int Tensor_Dim>
1069
1071 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1072 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1075 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1076 if (!dataPtr)
1077 THROW_MESSAGE("Pointer is not set");
1078 }
1079
1083 const int nb_gauss_pts = getGaussPts().size2();
1084 MatrixDouble &mat = *dataPtr;
1085 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1086 if (type == zeroType && side == zeroSide) {
1087 mat.resize(symm_size, nb_gauss_pts, false);
1088 mat.clear();
1089 }
1090 auto &local_indices = data.getLocalIndices();
1091 const int nb_dofs = local_indices.size();
1092 if (!nb_dofs)
1094
1095#ifndef NDEBUG
1096 if ((getFEMethod()->data_ctx & PetscData::CtxSetX_T).none()) {
1097 MOFEM_LOG_CHANNEL("SELF");
1098 MOFEM_LOG("SELF", Sev::error)
1099 << "In this case filed degrees of freedom are read from vector. "
1100 "That usually happens when time solver is used, and acces to "
1101 "first rates is needed. You probably not set "
1102 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1103 "respectively";
1104 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
1105 }
1106#endif
1107
1108 dotVector.resize(nb_dofs, false);
1109 const double *array;
1110 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1111 for (int i = 0; i != local_indices.size(); ++i)
1112 if (local_indices[i] != -1)
1113 dotVector[i] = array[local_indices[i]];
1114 else
1115 dotVector[i] = 0;
1116 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1117
1118 const int nb_base_functions = data.getN().size2();
1119
1120 auto base_function = data.getFTensor0N();
1121 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1124 const int size = nb_dofs / symm_size;
1125 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1126 auto field_data = getFTensorDotData<Tensor_Dim>();
1127 int bb = 0;
1128 for (; bb != size; ++bb) {
1129 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1130 ++field_data;
1131 ++base_function;
1132 }
1133 for (; bb < nb_base_functions; ++bb)
1134 ++base_function;
1135 ++values_at_gauss_pts;
1136 }
1137
1139 }
1140
1141protected:
1142 boost::shared_ptr<MatrixDouble> dataPtr;
1144 const int zeroSide;
1146
1147 template <int Dim> inline auto getFTensorDotData() {
1148 static_assert(Dim || !Dim, "not implemented");
1149 }
1150};
1151
1152template <>
1153template <>
1154inline auto
1157 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1158 &dotVector[5]);
1159}
1160
1161template <>
1162template <>
1163inline auto
1166 &dotVector[0], &dotVector[1], &dotVector[2]);
1167}
1168
1169/**@}*/
1170
1171/** \name Gradients and Hessian of scalar fields at integration points */
1172
1173/**@{*/
1174
1175/**
1176 * \brief Evaluate field gradient values for scalar field, i.e. gradient is
1177 * tensor rank 1 (vector)
1178 *
1179 */
1180template <int Tensor_Dim, class T, class L, class A>
1182 : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
1183
1185 Tensor_Dim, T, L, A>::OpCalculateVectorFieldValues_General;
1186};
1187
1188/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1189 * tensor rank 1 (vector), specialization
1190 *
1191 */
1192template <int Tensor_Dim>
1194 ublas::row_major, DoubleAllocator>
1196 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1197
1199 Tensor_Dim, double, ublas::row_major,
1200 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1201
1202 /**
1203 * \brief calculate gradient values of scalar field at integration points
1204 * @param side side entity number
1205 * @param type side entity type
1206 * @param data entity data
1207 * @return error code
1208 */
1209 MoFEMErrorCode doWork(int side, EntityType type,
1211};
1212
1213/**
1214 * \brief Member function specialization calculating scalar field gradients for
1215 * tenor field rank 1
1216 *
1217 */
1218template <int Tensor_Dim>
1220 Tensor_Dim, double, ublas::row_major,
1221 DoubleAllocator>::doWork(int side, EntityType type,
1224
1225 const size_t nb_gauss_pts = this->getGaussPts().size2();
1226 auto &mat = *this->dataPtr;
1227 if (type == this->zeroType) {
1228 mat.resize(Tensor_Dim, nb_gauss_pts, false);
1229 mat.clear();
1230 }
1231
1232 const int nb_dofs = data.getFieldData().size();
1233 if (nb_dofs) {
1234
1235 if (nb_gauss_pts) {
1236 const int nb_base_functions = data.getN().size2();
1237 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
1238 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1239
1240#ifndef NDEBUG
1241 if (nb_dofs > nb_base_functions)
1242 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1243 "Number of base functions inconsistent with number of DOFs "
1244 "(%d > %d)",
1245 nb_dofs, nb_base_functions);
1246
1247 if (data.getDiffN().size2() != nb_base_functions * Tensor_Dim)
1248 SETERRQ2(
1249 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1250 "Number of base functions inconsistent with number of derivatives "
1251 "(%d != %d)",
1252 data.getDiffN().size2(), nb_base_functions);
1253
1254 if (data.getDiffN().size1() != nb_gauss_pts)
1255 SETERRQ2(
1256 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1257 "Number of base functions inconsistent with number of integration "
1258 "pts (%d != %d)",
1259 data.getDiffN().size2(), nb_gauss_pts);
1260
1261#endif
1262
1264 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1265 auto field_data = data.getFTensor0FieldData();
1266 int bb = 0;
1267 for (; bb != nb_dofs; ++bb) {
1268 gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
1269 ++field_data;
1270 ++diff_base_function;
1271 }
1272 // Number of dofs can be smaller than number of base functions
1273 for (; bb < nb_base_functions; ++bb)
1274 ++diff_base_function;
1275 ++gradients_at_gauss_pts;
1276 }
1277 }
1278 }
1279
1281}
1282
1283/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1284 * vector field
1285 *
1286 * \ingroup mofem_forces_and_sources_user_data_operators
1287 */
1288template <int Tensor_Dim>
1291 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1293 Tensor_Dim, double, ublas::row_major,
1294 DoubleAllocator>::OpCalculateScalarFieldGradient_General;
1295};
1296
1297/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1298 * tensor rank 1 (vector), specialization
1299 *
1300 */
1301template <int Tensor_Dim>
1304 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1305
1307 Tensor_Dim, double, ublas::row_major,
1308 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1309
1310 /**
1311 * \brief calculate gradient values of scalar field at integration points
1312 * @param side side entity number
1313 * @param type side entity type
1314 * @param data entity data
1315 * @return error code
1316 */
1317 MoFEMErrorCode doWork(int side, EntityType type,
1319};
1320
1321template <int Tensor_Dim>
1323 int side, EntityType type, EntitiesFieldData::EntData &data) {
1325
1326 const size_t nb_gauss_pts = this->getGaussPts().size2();
1327 auto &mat = *this->dataPtr;
1328 if (type == this->zeroType) {
1329 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts, false);
1330 mat.clear();
1331 }
1332
1333 const int nb_dofs = data.getFieldData().size();
1334 if (nb_dofs) {
1335
1336 if (nb_gauss_pts) {
1337 const int nb_base_functions = data.getN().size2();
1338
1339 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1340#ifndef NDEBUG
1341 if (hessian_base.size1() != nb_gauss_pts) {
1342 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1343 "Wrong number of integration pts (%d != %d)",
1344 hessian_base.size1(), nb_gauss_pts);
1345 }
1346 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1347 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1348 "Wrong number of base functions (%d != %d)",
1349 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1350 nb_base_functions);
1351 }
1352 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1353 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1354 "Wrong number of base functions (%d < %d)",
1355 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1356 }
1357#endif
1358
1359 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1360 &*hessian_base.data().begin());
1361
1362 auto hessian_at_gauss_pts =
1363 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1364
1367 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1368 auto field_data = data.getFTensor0FieldData();
1369 int bb = 0;
1370 for (; bb != nb_dofs; ++bb) {
1371 hessian_at_gauss_pts(I, J) +=
1372 field_data * t_diff2_base_function(I, J);
1373 ++field_data;
1374 ++t_diff2_base_function;
1375 }
1376 // Number of dofs can be smaller than number of base functions
1377 for (; bb < nb_base_functions; ++bb) {
1378 ++t_diff2_base_function;
1379 }
1380
1381 ++hessian_at_gauss_pts;
1382 }
1383 }
1384 }
1385
1387}
1388
1389/**}*/
1390
1391/** \name Gradients and hessian of tensor fields at integration points */
1392
1393/**@{*/
1394
1395/**
1396 * \brief Evaluate field gradient values for vector field, i.e. gradient is
1397 * tensor rank 2
1398 *
1399 */
1400template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1402 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1403 L, A> {
1404
1406 Tensor_Dim0, Tensor_Dim1, T, L, A>::OpCalculateTensor2FieldValues_General;
1407};
1408
1409template <int Tensor_Dim0, int Tensor_Dim1>
1410struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1411 ublas::row_major, DoubleAllocator>
1413 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1414
1416 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1417 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1418
1419 /**
1420 * \brief calculate values of vector field at integration points
1421 * @param side side entity number
1422 * @param type side entity type
1423 * @param data entity data
1424 * @return error code
1425 */
1426 MoFEMErrorCode doWork(int side, EntityType type,
1428};
1429
1430/**
1431 * \brief Member function specialization calculating vector field gradients for
1432 * tenor field rank 2
1433 *
1434 */
1435template <int Tensor_Dim0, int Tensor_Dim1>
1437 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1438 DoubleAllocator>::doWork(int side, EntityType type,
1441 if (!this->dataPtr)
1442 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1443 "Data pointer not allocated");
1444
1445 const size_t nb_gauss_pts = this->getGaussPts().size2();
1446 auto &mat = *this->dataPtr;
1447 if (type == this->zeroType) {
1448 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1449 mat.clear();
1450 }
1451
1452 if (nb_gauss_pts) {
1453 const size_t nb_dofs = data.getFieldData().size();
1454
1455 if (nb_dofs) {
1456
1457 if (this->dataVec.use_count()) {
1458 this->dotVector.resize(nb_dofs, false);
1459 const double *array;
1460 CHKERR VecGetArrayRead(this->dataVec, &array);
1461 const auto &local_indices = data.getLocalIndices();
1462 for (int i = 0; i != local_indices.size(); ++i)
1463 if (local_indices[i] != -1)
1464 this->dotVector[i] = array[local_indices[i]];
1465 else
1466 this->dotVector[i] = 0;
1467 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1468 data.getFieldData().swap(this->dotVector);
1469 }
1470
1471 const int nb_base_functions = data.getN().size2();
1472 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1473 auto gradients_at_gauss_pts =
1474 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1477 int size = nb_dofs / Tensor_Dim0;
1478 if (nb_dofs % Tensor_Dim0) {
1479 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1480 "Data inconsistency");
1481 }
1482 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1483 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1484 int bb = 0;
1485 for (; bb < size; ++bb) {
1486 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1487 ++field_data;
1488 ++diff_base_function;
1489 }
1490 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1491 // functions
1492 for (; bb != nb_base_functions; ++bb)
1493 ++diff_base_function;
1494 ++gradients_at_gauss_pts;
1495 }
1496
1497 if (this->dataVec.use_count()) {
1498 data.getFieldData().swap(this->dotVector);
1499 }
1500 }
1501 }
1503}
1504
1505/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1506 * vector field
1507 *
1508 * \ingroup mofem_forces_and_sources_user_data_operators
1509 */
1510template <int Tensor_Dim0, int Tensor_Dim1>
1513 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1514
1516 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1517 DoubleAllocator>::OpCalculateVectorFieldGradient_General;
1518};
1519
1520/** \brief Get field gradients time derivative at integration pts for scalar
1521 * filed rank 0, i.e. vector field
1522 *
1523 * \ingroup mofem_forces_and_sources_user_data_operators
1524 */
1525template <int Tensor_Dim0, int Tensor_Dim1>
1528
1530 boost::shared_ptr<MatrixDouble> data_ptr,
1531 const EntityType zero_at_type = MBVERTEX)
1534 dataPtr(data_ptr), zeroAtType(zero_at_type) {
1535 if (!dataPtr)
1536 THROW_MESSAGE("Pointer is not set");
1537 }
1538
1542
1543 const auto &local_indices = data.getLocalIndices();
1544 const int nb_dofs = local_indices.size();
1545 const int nb_gauss_pts = this->getGaussPts().size2();
1546
1547 auto &mat = *this->dataPtr;
1548 if (type == this->zeroAtType) {
1549 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1550 mat.clear();
1551 }
1552 if (!nb_dofs)
1554
1555 dotVector.resize(nb_dofs, false);
1556 const double *array;
1557 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1558 for (int i = 0; i != local_indices.size(); ++i)
1559 if (local_indices[i] != -1)
1560 dotVector[i] = array[local_indices[i]];
1561 else
1562 dotVector[i] = 0;
1563 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1564
1565 const int nb_base_functions = data.getN().size2();
1566 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1567 auto gradients_at_gauss_pts =
1568 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1571 int size = nb_dofs / Tensor_Dim0;
1572 if (nb_dofs % Tensor_Dim0) {
1573 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1574 }
1575
1576 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1577 auto field_data = getFTensorDotData<Tensor_Dim0>();
1578 int bb = 0;
1579 for (; bb < size; ++bb) {
1580 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1581 ++field_data;
1582 ++diff_base_function;
1583 }
1584 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1585 // functions
1586 for (; bb != nb_base_functions; ++bb)
1587 ++diff_base_function;
1588 ++gradients_at_gauss_pts;
1589 }
1591 }
1592
1593private:
1594 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1595 EntityType zeroAtType; ///< Zero values at Gauss point at this type
1596 VectorDouble dotVector; ///< Keeps temoorary values of time directives
1597
1598 template <int Dim> inline auto getFTensorDotData() {
1599 static_assert(Dim || !Dim, "not implemented");
1600 }
1601};
1602
1603template <>
1604template <>
1607 &dotVector[0], &dotVector[1], &dotVector[2]);
1608}
1609
1610template <>
1611template <>
1613 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(&dotVector[0],
1614 &dotVector[1]);
1615}
1616
1617/**
1618 * \brief Evaluate field gradient values for symmetric 2nd order tensor field,
1619 * i.e. gradient is tensor rank 3
1620 *
1621 */
1622template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1624 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1625 L, A> {
1626
1628 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1629 const EntityType zero_type = MBVERTEX)
1630 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
1631 A>(field_name, data_ptr,
1632 zero_type) {}
1633};
1634
1635template <int Tensor_Dim0, int Tensor_Dim1>
1637 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator>
1639 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1640
1642 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1643 const EntityType zero_type = MBVERTEX)
1644 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1645 ublas::row_major,
1647 field_name, data_ptr, zero_type) {}
1648
1649 /**
1650 * \brief calculate values of vector field at integration points
1651 * @param side side entity number
1652 * @param type side entity type
1653 * @param data entity data
1654 * @return error code
1655 */
1656 MoFEMErrorCode doWork(int side, EntityType type,
1658};
1659
1660/**
1661 * \brief Member function specialization calculating tensor field gradients for
1662 * symmetric tensor field rank 2
1663 *
1664 */
1665template <int Tensor_Dim0, int Tensor_Dim1>
1666MoFEMErrorCode OpCalculateTensor2SymmetricFieldGradient_General<
1667 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1668 DoubleAllocator>::doWork(int side, EntityType type,
1671 if (!this->dataPtr)
1672 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1673 "Data pointer not allocated");
1674
1675 const size_t nb_gauss_pts = this->getGaussPts().size2();
1676 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1677 auto &mat = *this->dataPtr;
1678 if (type == this->zeroType) {
1679 mat.resize(msize * Tensor_Dim1, nb_gauss_pts, false);
1680 mat.clear();
1681 }
1682
1683 if (nb_gauss_pts) {
1684 const size_t nb_dofs = data.getFieldData().size();
1685
1686 if (nb_dofs) {
1687
1688 const int nb_base_functions = data.getN().size2();
1689 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1690 auto gradients_at_gauss_pts =
1691 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1695 int size = nb_dofs / msize;
1696 if (nb_dofs % msize) {
1697 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1698 "Data inconsistency");
1699 }
1700 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1701 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim0>();
1702 int bb = 0;
1703 for (; bb < size; ++bb) {
1704 gradients_at_gauss_pts(I, J, K) +=
1705 field_data(I, J) * diff_base_function(K);
1706 ++field_data;
1707 ++diff_base_function;
1708 }
1709 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1710 // functions
1711 for (; bb != nb_base_functions; ++bb)
1712 ++diff_base_function;
1713 ++gradients_at_gauss_pts;
1714 }
1715 }
1716 }
1718}
1719
1720/** \brief Get field gradients at integration pts for symmetric tensorial field
1721 * rank 2
1722 *
1723 * \ingroup mofem_forces_and_sources_user_data_operators
1724 */
1725template <int Tensor_Dim0, int Tensor_Dim1>
1728 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1729
1731 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1732 const EntityType zero_type = MBVERTEX)
1734 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1735 DoubleAllocator>(field_name, data_ptr, zero_type) {}
1736};
1737
1738template <int Tensor_Dim0, int Tensor_Dim1>
1741 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1742
1744 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1745 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1746
1747 /**
1748 * \brief calculate values of vector field at integration points
1749 * @param side side entity number
1750 * @param type side entity type
1751 * @param data entity data
1752 * @return error code
1753 */
1754 MoFEMErrorCode doWork(int side, EntityType type,
1756};
1757
1758/**
1759 * \brief Member function specialization calculating vector field gradients for
1760 * tenor field rank 2
1761 *
1762 */
1763template <int Tensor_Dim0, int Tensor_Dim1>
1765 int side, EntityType type, EntitiesFieldData::EntData &data) {
1767 if (!this->dataPtr)
1768 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1769 "Data pointer not allocated");
1770
1771 const size_t nb_gauss_pts = this->getGaussPts().size2();
1772 auto &mat = *this->dataPtr;
1773 if (type == this->zeroType) {
1774 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts, false);
1775 mat.clear();
1776 }
1777
1778 if (nb_gauss_pts) {
1779 const size_t nb_dofs = data.getFieldData().size();
1780
1781 if (nb_dofs) {
1782
1783 if (this->dataVec.use_count()) {
1784 this->dotVector.resize(nb_dofs, false);
1785 const double *array;
1786 CHKERR VecGetArrayRead(this->dataVec, &array);
1787 const auto &local_indices = data.getLocalIndices();
1788 for (int i = 0; i != local_indices.size(); ++i)
1789 if (local_indices[i] != -1)
1790 this->dotVector[i] = array[local_indices[i]];
1791 else
1792 this->dotVector[i] = 0;
1793 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1794 data.getFieldData().swap(this->dotVector);
1795 }
1796
1797 const int nb_base_functions = data.getN().size2();
1798
1799 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1800#ifndef NDEBUG
1801 if (hessian_base.size1() != nb_gauss_pts) {
1802 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1803 "Wrong number of integration pts (%d != %d)",
1804 hessian_base.size1(), nb_gauss_pts);
1805 }
1806 if (hessian_base.size2() !=
1807 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1808 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1809 "Wrong number of base functions (%d != %d)",
1810 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1811 nb_base_functions);
1812 }
1813 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1814 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1815 "Wrong number of base functions (%d < %d)",
1816 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1817 }
1818#endif
1819
1820 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1821 &*hessian_base.data().begin());
1822
1823 auto t_hessian_at_gauss_pts =
1824 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1825
1829
1830 int size = nb_dofs / Tensor_Dim0;
1831#ifndef NDEBUG
1832 if (nb_dofs % Tensor_Dim0) {
1833 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1834 "Data inconsistency");
1835 }
1836#endif
1837
1838 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1839 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1840 int bb = 0;
1841 for (; bb < size; ++bb) {
1842 t_hessian_at_gauss_pts(I, J, K) +=
1843 field_data(I) * t_diff2_base_function(J, K);
1844 ++field_data;
1845 ++t_diff2_base_function;
1846 }
1847 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1848 // functions
1849 for (; bb != nb_base_functions; ++bb)
1850 ++t_diff2_base_function;
1851 ++t_hessian_at_gauss_pts;
1852 }
1853
1854 if (this->dataVec.use_count()) {
1855 data.getFieldData().swap(this->dotVector);
1856 }
1857 }
1858 }
1860}
1861
1862/**@}*/
1863
1864/** \name Transform tensors and vectors */
1865
1866/**@{*/
1867
1868/**
1869 * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1870 * \f$
1871 *
1872 * @tparam DIM
1873 *
1874 * \ingroup mofem_forces_and_sources_user_data_operators
1875 */
1876template <int DIM_01, int DIM_23, int S = 0>
1879
1882
1884 boost::shared_ptr<MatrixDouble> in_mat,
1885 boost::shared_ptr<MatrixDouble> out_mat,
1886 boost::shared_ptr<MatrixDouble> d_mat)
1887 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1888 // Only is run for vertices
1889 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1890 if (!inMat)
1891 THROW_MESSAGE("Pointer for in mat is null");
1892 if (!outMat)
1893 THROW_MESSAGE("Pointer for out mat is null");
1894 if (!dMat)
1895 THROW_MESSAGE("Pointer for tensor mat is null");
1896 }
1897
1898 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1900 const size_t nb_gauss_pts = getGaussPts().size2();
1901 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
1902 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
1903 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1904 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
1905 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1906 t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1907 ++t_in;
1908 ++t_out;
1909 }
1911 }
1912
1913private:
1918
1919 boost::shared_ptr<MatrixDouble> inMat;
1920 boost::shared_ptr<MatrixDouble> outMat;
1921 boost::shared_ptr<MatrixDouble> dMat;
1922};
1923
1924template <int DIM>
1926
1929
1931 boost::shared_ptr<MatrixDouble> in_mat,
1932 boost::shared_ptr<MatrixDouble> out_mat)
1933 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1934 // Only is run for vertices
1935 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1936 if (!inMat)
1937 THROW_MESSAGE("Pointer not set for in matrix");
1938 if (!outMat)
1939 THROW_MESSAGE("Pointer not set for in matrix");
1940 }
1941
1942 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1944 const size_t nb_gauss_pts = getGaussPts().size2();
1945 auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
1946 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
1947 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
1948 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1949 t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
1950 ++t_in;
1951 ++t_out;
1952 }
1954 }
1955
1956private:
1959 boost::shared_ptr<MatrixDouble> inMat;
1960 boost::shared_ptr<MatrixDouble> outMat;
1961};
1962
1964
1967
1968 OpScaleMatrix(const std::string field_name, const double scale,
1969 boost::shared_ptr<MatrixDouble> in_mat,
1970 boost::shared_ptr<MatrixDouble> out_mat)
1971 : UserOp(field_name, OPROW), scale(scale), inMat(in_mat),
1972 outMat(out_mat) {
1973 // Only is run for vertices
1974 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1975 if (!inMat)
1976 THROW_MESSAGE("Pointer not set for in matrix");
1977 if (!outMat)
1978 THROW_MESSAGE("Pointer not set for in matrix");
1979 }
1980
1981 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1983 outMat->resize(inMat->size1(), inMat->size2(), false);
1984 noalias(*outMat) = scale * (*inMat);
1986 }
1987
1988private:
1989 const double scale;
1990 boost::shared_ptr<MatrixDouble> inMat;
1991 boost::shared_ptr<MatrixDouble> outMat;
1992};
1993
1994/**@}*/
1995
1996/** \name H-div/H-curls (Vectorial bases) values at integration points */
1997
1998/**@{*/
1999
2000/** \brief Get vector field for H-div approximation
2001 * \ingroup mofem_forces_and_sources_user_data_operators
2002 */
2003template <int Base_Dim, int Field_Dim, class T, class L, class A>
2005
2006/** \brief Get vector field for H-div approximation
2007 * \ingroup mofem_forces_and_sources_user_data_operators
2008 */
2009template <int Field_Dim>
2011 ublas::row_major, DoubleAllocator>
2013
2015 boost::shared_ptr<MatrixDouble> data_ptr,
2016 const EntityType zero_type = MBEDGE,
2017 const int zero_side = 0)
2020 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2021 if (!dataPtr)
2022 THROW_MESSAGE("Pointer is not set");
2023 }
2024
2025 /**
2026 * \brief Calculate values of vector field at integration points
2027 * @param side side entity number
2028 * @param type side entity type
2029 * @param data entity data
2030 * @return error code
2031 */
2032 MoFEMErrorCode doWork(int side, EntityType type,
2034
2035private:
2036 boost::shared_ptr<MatrixDouble> dataPtr;
2038 const int zeroSide;
2039};
2040
2041template <int Field_Dim>
2043 3, Field_Dim, double, ublas::row_major,
2044 DoubleAllocator>::doWork(int side, EntityType type,
2047 const size_t nb_integration_points = this->getGaussPts().size2();
2048 if (type == zeroType && side == zeroSide) {
2049 dataPtr->resize(Field_Dim, nb_integration_points, false);
2050 dataPtr->clear();
2051 }
2052 const size_t nb_dofs = data.getFieldData().size();
2053 if (!nb_dofs)
2055 const size_t nb_base_functions = data.getN().size2() / 3;
2057 auto t_n_hdiv = data.getFTensor1N<3>();
2058 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2059 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2060 auto t_dof = data.getFTensor0FieldData();
2061 int bb = 0;
2062 for (; bb != nb_dofs; ++bb) {
2063 t_data(i) += t_n_hdiv(i) * t_dof;
2064 ++t_n_hdiv;
2065 ++t_dof;
2066 }
2067 for (; bb != nb_base_functions; ++bb)
2068 ++t_n_hdiv;
2069 ++t_data;
2070 }
2072}
2073
2074/** \brief Get vector field for H-div approximation
2075 *
2076 * \ingroup mofem_forces_and_sources_user_data_operators
2077 */
2078template <int Base_Dim, int Field_Dim = Base_Dim>
2081 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2083 Base_Dim, Field_Dim, double, ublas::row_major,
2084 DoubleAllocator>::OpCalculateHVecVectorField_General;
2085};
2086
2087
2088
2089/** \brief Get vector field for H-div approximation
2090 * \ingroup mofem_forces_and_sources_user_data_operators
2091 */
2092template <int Base_Dim, int Field_Dim = Base_Dim>
2094
2095template <int Field_Dim>
2098
2100 boost::shared_ptr<MatrixDouble> data_ptr,
2101 const EntityType zero_type = MBEDGE,
2102 const int zero_side = 0)
2105 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2106 if (!dataPtr)
2107 THROW_MESSAGE("Pointer is not set");
2108 }
2109
2110 /**
2111 * \brief Calculate values of vector field at integration points
2112 * @param side side entity number
2113 * @param type side entity type
2114 * @param data entity data
2115 * @return error code
2116 */
2117 MoFEMErrorCode doWork(int side, EntityType type,
2119
2120private:
2121 boost::shared_ptr<MatrixDouble> dataPtr;
2123 const int zeroSide;
2124};
2125
2126template <int Field_Dim>
2128 int side, EntityType type, EntitiesFieldData::EntData &data) {
2130
2131 const size_t nb_integration_points = this->getGaussPts().size2();
2132 if (type == zeroType && side == zeroSide) {
2133 dataPtr->resize(Field_Dim, nb_integration_points, false);
2134 dataPtr->clear();
2135 }
2136
2137 auto &local_indices = data.getIndices();
2138 const size_t nb_dofs = local_indices.size();
2139 if (nb_dofs) {
2140
2141 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2142 const double *array;
2143 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2144 for (size_t i = 0; i != nb_dofs; ++i)
2145 if (local_indices[i] != -1)
2146 dot_dofs_vector[i] = array[local_indices[i]];
2147 else
2148 dot_dofs_vector[i] = 0;
2149 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2150
2151 const size_t nb_base_functions = data.getN().size2() / 3;
2153 auto t_n_hdiv = data.getFTensor1N<3>();
2154 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2155 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2156 int bb = 0;
2157 for (; bb != nb_dofs; ++bb) {
2158 t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
2159 ++t_n_hdiv;
2160 }
2161 for (; bb != nb_base_functions; ++bb)
2162 ++t_n_hdiv;
2163 ++t_data;
2164 }
2165 }
2166
2168}
2169
2170/**
2171 * @brief Calculate divergence of vector field
2172 * @ingroup mofem_forces_and_sources_user_data_operators
2173 *
2174 * @tparam BASE_DIM
2175 * @tparam SPACE_DIM
2176 */
2177template <int BASE_DIM, int SPACE_DIM>
2180
2182 boost::shared_ptr<VectorDouble> data_ptr,
2183 const EntityType zero_type = MBEDGE,
2184 const int zero_side = 0)
2187 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2188 if (!dataPtr)
2189 THROW_MESSAGE("Pointer is not set");
2190 }
2191
2195 const size_t nb_integration_points = getGaussPts().size2();
2196 if (type == zeroType && side == zeroSide) {
2197 dataPtr->resize(nb_integration_points, false);
2198 dataPtr->clear();
2199 }
2200 const size_t nb_dofs = data.getFieldData().size();
2201 if (!nb_dofs)
2203 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2206 auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2207 auto t_data = getFTensor0FromVec(*dataPtr);
2208 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2209 auto t_dof = data.getFTensor0FieldData();
2210 int bb = 0;
2211 for (; bb != nb_dofs; ++bb) {
2212 t_data += t_dof * t_n_diff_hdiv(j, j);
2213 ++t_n_diff_hdiv;
2214 ++t_dof;
2215 }
2216 for (; bb != nb_base_functions; ++bb)
2217 ++t_n_diff_hdiv;
2218 ++t_data;
2219 }
2221 }
2222
2223private:
2224 boost::shared_ptr<VectorDouble> dataPtr;
2226 const int zeroSide;
2227};
2228
2229/**
2230 * @brief Calculate gradient of vector field
2231 * @ingroup mofem_forces_and_sources_user_data_operators
2232 *
2233 * @tparam BASE_DIM
2234 * @tparam SPACE_DIM
2235 */
2236template <int BASE_DIM, int SPACE_DIM>
2239
2241 boost::shared_ptr<MatrixDouble> data_ptr,
2242 const EntityType zero_type = MBEDGE,
2243 const int zero_side = 0)
2246 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2247 if (!dataPtr)
2248 THROW_MESSAGE("Pointer is not set");
2249 }
2250
2254 const size_t nb_integration_points = getGaussPts().size2();
2255 if (type == zeroType && side == zeroSide) {
2256 dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
2257 dataPtr->clear();
2258 }
2259 const size_t nb_dofs = data.getFieldData().size();
2260 if (!nb_dofs)
2262 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2265 auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2266 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*dataPtr);
2267 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2268 auto t_dof = data.getFTensor0FieldData();
2269 int bb = 0;
2270 for (; bb != nb_dofs; ++bb) {
2271 t_data(i, j) += t_dof * t_base_diff(i, j);
2272 ++t_base_diff;
2273 ++t_dof;
2274 }
2275 for (; bb != nb_base_functions; ++bb)
2276 ++t_base_diff;
2277 ++t_data;
2278 }
2280 }
2281
2282private:
2283 boost::shared_ptr<MatrixDouble> dataPtr;
2285 const int zeroSide;
2286};
2287
2288/**
2289 * @brief Calculate gradient of vector field
2290 * @ingroup mofem_forces_and_sources_user_data_operators
2291 *
2292 * @tparam BASE_DIM
2293 * @tparam SPACE_DIM
2294 */
2295template <int BASE_DIM, int SPACE_DIM>
2298
2300 boost::shared_ptr<MatrixDouble> data_ptr,
2301 const EntityType zero_type = MBEDGE,
2302 const int zero_side = 0)
2305 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2306 if (!dataPtr)
2307 THROW_MESSAGE("Pointer is not set");
2308 }
2309
2313 const size_t nb_integration_points = getGaussPts().size2();
2314 if (type == zeroType && side == zeroSide) {
2315 dataPtr->resize(BASE_DIM * SPACE_DIM * SPACE_DIM, nb_integration_points,
2316 false);
2317 dataPtr->clear();
2318 }
2319 const size_t nb_dofs = data.getFieldData().size();
2320 if (!nb_dofs)
2322
2323 const int nb_base_functions = data.getN().size2() / BASE_DIM;
2324 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2325
2326#ifndef NDEBUG
2327 if (hessian_base.size1() != nb_integration_points) {
2328 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2329 "Wrong number of integration pts (%d != %d)",
2330 hessian_base.size1(), nb_integration_points);
2331 }
2332 if (hessian_base.size2() !=
2333 BASE_DIM * nb_base_functions * SPACE_DIM * SPACE_DIM) {
2334 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2335 "Wrong number of base functions (%d != %d)",
2336 hessian_base.size2() / (BASE_DIM * SPACE_DIM * SPACE_DIM),
2337 nb_base_functions);
2338 }
2339 if (hessian_base.size2() < BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM) {
2340 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2341 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2342 BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM);
2343 }
2344#endif
2345
2349
2350 auto t_base_diff2 =
2352 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*dataPtr);
2353 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2354 auto t_dof = data.getFTensor0FieldData();
2355 int bb = 0;
2356 for (; bb != nb_dofs; ++bb) {
2357 t_data(i, j, k) += t_dof * t_base_diff2(i, j, k);
2358
2359 ++t_base_diff2;
2360 ++t_dof;
2361 }
2362 for (; bb != nb_base_functions; ++bb)
2363 ++t_base_diff2;
2364 ++t_data;
2365 }
2367 }
2368
2369private:
2370 boost::shared_ptr<MatrixDouble> dataPtr;
2372 const int zeroSide;
2373};
2374
2375/**
2376 * @brief Calculate divergence of vector field dot
2377 * @ingroup mofem_forces_and_sources_user_data_operators
2378 *
2379 * @tparam Tensor_Dim dimension of space
2380 */
2381template <int Tensor_Dim1, int Tensor_Dim2>
2384
2386 boost::shared_ptr<VectorDouble> data_ptr,
2387 const EntityType zero_type = MBEDGE,
2388 const int zero_side = 0)
2391 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2392 if (!dataPtr)
2393 THROW_MESSAGE("Pointer is not set");
2394 }
2395
2399 const size_t nb_integration_points = getGaussPts().size2();
2400 if (type == zeroType && side == zeroSide) {
2401 dataPtr->resize(nb_integration_points, false);
2402 dataPtr->clear();
2403 }
2404
2405 const auto &local_indices = data.getLocalIndices();
2406 const int nb_dofs = local_indices.size();
2407 if (nb_dofs) {
2408
2409 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2410 const double *array;
2411 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2412 for (size_t i = 0; i != local_indices.size(); ++i)
2413 if (local_indices[i] != -1)
2414 dot_dofs_vector[i] = array[local_indices[i]];
2415 else
2416 dot_dofs_vector[i] = 0;
2417 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2418
2419 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
2421 auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
2422 auto t_data = getFTensor0FromVec(*dataPtr);
2423 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2424 int bb = 0;
2425 for (; bb != nb_dofs; ++bb) {
2426 double div = 0;
2427 for (auto ii = 0; ii != Tensor_Dim2; ++ii)
2428 div += t_n_diff_hdiv(ii, ii);
2429 t_data += dot_dofs_vector[bb] * div;
2430 ++t_n_diff_hdiv;
2431 }
2432 for (; bb != nb_base_functions; ++bb)
2433 ++t_n_diff_hdiv;
2434 ++t_data;
2435 }
2436 }
2438 }
2439
2440private:
2441 boost::shared_ptr<VectorDouble> dataPtr;
2443 const int zeroSide;
2444};
2445
2446/**
2447 * @brief Calculate curl of vector field
2448 * @ingroup mofem_forces_and_sources_user_data_operators
2449 *
2450 * @tparam Base_Dim base function dimension
2451 * @tparam Space_Dim dimension of space
2452 * @tparam Hcurl field dimension
2453 */
2454template <int Base_Dim, int Space_Dim>
2456
2457/**
2458 * @brief Calculate curl of vector field
2459 * @ingroup mofem_forces_and_sources_user_data_operators
2460 *
2461 * @tparam Base_Dim base function dimension
2462 * @tparam Space_Dim dimension of space
2463 * @tparam Hcurl field dimension
2464 */
2465template <>
2468 OpCalculateHcurlVectorCurl(const std::string field_name,
2469 boost::shared_ptr<MatrixDouble> data_ptr,
2470 const EntityType zero_type = MBEDGE,
2471 const int zero_side = 0);
2472 MoFEMErrorCode doWork(int side, EntityType type,
2474
2475private:
2476 boost::shared_ptr<MatrixDouble> dataPtr;
2478 const int zeroSide;
2479};
2480
2481/**
2482 * @brief Calculate curl of vector field
2483 * @ingroup mofem_forces_and_sources_user_data_operators
2484 *
2485 * @tparam Field_Dim dimension of field
2486 * @tparam Space_Dim dimension of space
2487 */
2488template <>
2491
2492 OpCalculateHcurlVectorCurl(const std::string field_name,
2493 boost::shared_ptr<MatrixDouble> data_ptr,
2494 const EntityType zero_type = MBVERTEX,
2495 const int zero_side = 0);
2496
2497 MoFEMErrorCode doWork(int side, EntityType type,
2499
2500private:
2501 boost::shared_ptr<MatrixDouble> dataPtr;
2503 const int zeroSide;
2504};
2505
2506/**
2507 * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
2508 * \ingroup mofem_forces_and_sources_user_data_operators
2509 *
2510 * @tparam Tensor_Dim0 rank of the filed
2511 * @tparam Tensor_Dim1 dimension of space
2512 */
2513template <int Tensor_Dim0, int Tensor_Dim1>
2516
2518 boost::shared_ptr<MatrixDouble> data_ptr,
2519 const EntityType zero_type = MBEDGE,
2520 const int zero_side = 0)
2523 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2524 if (!dataPtr)
2525 THROW_MESSAGE("Pointer is not set");
2526 }
2527
2531 const size_t nb_integration_points = getGaussPts().size2();
2532 if (type == zeroType && side == zeroSide) {
2533 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2534 dataPtr->clear();
2535 }
2536 const size_t nb_dofs = data.getFieldData().size();
2537 if (nb_dofs) {
2538 const size_t nb_base_functions = data.getN().size2() / 3;
2541 auto t_n_hvec = data.getFTensor1N<3>();
2542 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2543 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2544 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2545 size_t bb = 0;
2546 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2547 t_data(i, j) += t_dof(i) * t_n_hvec(j);
2548 ++t_n_hvec;
2549 ++t_dof;
2550 }
2551 for (; bb < nb_base_functions; ++bb)
2552 ++t_n_hvec;
2553 ++t_data;
2554 }
2555 }
2557 }
2558
2559private:
2560 boost::shared_ptr<MatrixDouble> dataPtr;
2562 const int zeroSide;
2563};
2564
2565/**
2566 * @brief Calculate tenor field using tensor base, i.e. Hdiv/Hcurl
2567 * \ingroup mofem_forces_and_sources_user_data_operators
2568 *
2569 * @tparam Tensor_Dim0 rank of the filed
2570 * @tparam Tensor_Dim1 dimension of space
2571 */
2572template <int Tensor_Dim0, int Tensor_Dim1>
2575
2577 boost::shared_ptr<MatrixDouble> data_ptr,
2578 const EntityType zero_type = MBEDGE,
2579 const int zero_side = 0)
2582 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2583 if (!dataPtr)
2584 THROW_MESSAGE("Pointer is not set");
2585 }
2586
2590 const size_t nb_integration_points = getGaussPts().size2();
2591 if (type == zeroType && side == zeroSide) {
2592 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2593 dataPtr->clear();
2594 }
2595 const size_t nb_dofs = data.getFieldData().size();
2596 if (!nb_dofs)
2598 const size_t nb_base_functions =
2599 data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2602 auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2603 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2604 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2605 auto t_dof = data.getFTensor0FieldData();
2606 size_t bb = 0;
2607 for (; bb != nb_dofs; ++bb) {
2608 t_data(i, j) += t_dof * t_n_hten(i, j);
2609 ++t_n_hten;
2610 ++t_dof;
2611 }
2612 for (; bb < nb_base_functions; ++bb)
2613 ++t_n_hten;
2614 ++t_data;
2615 }
2617 }
2618
2619private:
2620 boost::shared_ptr<MatrixDouble> dataPtr;
2622 const int zeroSide;
2623};
2624
2625/**
2626 * @brief Calculate divergence of tonsorial field using vectorial base
2627 * \ingroup mofem_forces_and_sources_user_data_operators
2628 *
2629 * @tparam Tensor_Dim0 rank of the field
2630 * @tparam Tensor_Dim1 dimension of space
2631 */
2632template <int Tensor_Dim0, int Tensor_Dim1>
2635
2637 boost::shared_ptr<MatrixDouble> data_ptr,
2638 const EntityType zero_type = MBEDGE,
2639 const int zero_side = 0)
2642 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2643 if (!dataPtr)
2644 THROW_MESSAGE("Pointer is not set");
2645 }
2646
2650 const size_t nb_integration_points = getGaussPts().size2();
2651 if (type == zeroType && side == 0) {
2652 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
2653 dataPtr->clear();
2654 }
2655 const size_t nb_dofs = data.getFieldData().size();
2656 if (nb_dofs) {
2657 const size_t nb_base_functions = data.getN().size2() / 3;
2660 auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
2661 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
2662 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2663 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2664 size_t bb = 0;
2665 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2666 double div = t_n_diff_hvec(j, j);
2667 t_data(i) += t_dof(i) * div;
2668 ++t_n_diff_hvec;
2669 ++t_dof;
2670 }
2671 for (; bb < nb_base_functions; ++bb)
2672 ++t_n_diff_hvec;
2673 ++t_data;
2674 }
2675 }
2677 }
2678
2679private:
2680 boost::shared_ptr<MatrixDouble> dataPtr;
2682 const int zeroSide;
2683};
2684
2685/**
2686 * @brief Calculate trace of vector (Hdiv/Hcurl) space
2687 *
2688 * @tparam Tensor_Dim
2689 * @tparam OpBase
2690 */
2691template <int Tensor_Dim, typename OpBase>
2693
2695 boost::shared_ptr<MatrixDouble> data_ptr,
2696 const EntityType zero_type = MBEDGE,
2697 const int zero_side = 0)
2698 : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
2699 zeroType(zero_type), zeroSide(zero_side) {
2700 if (!dataPtr)
2701 THROW_MESSAGE("Pointer is not set");
2702 }
2703
2707 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2708 if (type == zeroType && side == 0) {
2709 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2710 dataPtr->clear();
2711 }
2712 const size_t nb_dofs = data.getFieldData().size();
2713 if (nb_dofs) {
2714 auto t_normal = OpBase::getFTensor1Normal();
2715 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
2716 const size_t nb_base_functions = data.getN().size2() / 3;
2717 auto t_base = data.getFTensor1N<3>();
2718 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2719 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2720 auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
2721 size_t bb = 0;
2722 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2723 t_data(i) += t_dof(i) * (t_base(j) * t_normal(j));
2724 ++t_base;
2725 ++t_dof;
2726 }
2727 for (; bb < nb_base_functions; ++bb)
2728 ++t_base;
2729 ++t_data;
2730 }
2731 }
2733 }
2734
2735private:
2736 boost::shared_ptr<MatrixDouble> dataPtr;
2738 const int zeroSide;
2741};
2742
2743template <>
2747
2749 boost::shared_ptr<MatrixDouble> data_ptr,
2750 const EntityType zero_type = MBEDGE,
2751 const int zero_side = 0)
2752 : UserDataOperator(field_name, OPROW), dataPtr(data_ptr),
2753 zeroType(zero_type), zeroSide(zero_side) {
2754 if (!dataPtr)
2755 THROW_MESSAGE("Pointer is not set");
2756 }
2757
2761 const size_t nb_integration_points = getGaussPts().size2();
2762 if (type == zeroType && side == 0) {
2763 dataPtr->resize(3, nb_integration_points, false);
2764 dataPtr->clear();
2765 }
2766 const size_t nb_dofs = data.getFieldData().size();
2767 if (nb_dofs) {
2768 auto t_normal = getFTensor1Normal();
2769 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
2770 const size_t nb_base_functions = data.getN().size2() / 3;
2771 auto t_base = data.getFTensor1N<3>();
2772 auto t_data = getFTensor1FromMat<3>(*dataPtr);
2773 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2774 auto t_dof = data.getFTensor1FieldData<3>();
2775 if (getNormalsAtGaussPts().size1() == nb_integration_points) {
2776 VectorDouble n = getNormalsAtGaussPts(gg);
2777 auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
2778 t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
2779 }
2780 size_t bb = 0;
2781 for (; bb != nb_dofs / 3; ++bb) {
2782 t_data(i) += t_dof(i) * (t_base(j) * t_normal(j));
2783 ++t_base;
2784 ++t_dof;
2785 }
2786 for (; bb < nb_base_functions; ++bb)
2787 ++t_base;
2788 ++t_data;
2789 }
2790 }
2792 }
2793
2794private:
2795 boost::shared_ptr<MatrixDouble> dataPtr;
2797 const int zeroSide;
2800};
2801
2802/**@}*/
2803
2804/** \name Other operators */
2805
2806/**@{*/
2807
2808/**@}*/
2809
2810/** \name Operators for faces */
2811
2812/**@{*/
2813
2814/** \brief Transform local reference derivatives of shape functions to global
2815derivatives
2816
2817\ingroup mofem_forces_and_sources_tri_element
2818
2819*/
2820template <int DIM, int DERIVATIVE = 1> struct OpSetInvJacSpaceForFaceImpl;
2821
2824
2826 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2827
2828protected:
2829 template <int D1, int D2, int J1, int J2>
2832
2833 static_assert(D2 == J2, "Dimension of jacobian and dimension of <out> "
2834 "directive does not match");
2835
2836 size_t nb_functions = diff_n.size2() / D1;
2837 if (nb_functions) {
2838 size_t nb_gauss_pts = diff_n.size1();
2839
2840#ifndef NDEBUG
2841 if (nb_gauss_pts != getGaussPts().size2())
2842 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2843 "Wrong number of Gauss Pts");
2844 if (diff_n.size2() % D1)
2845 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2846 "Number of direvatives of base functions and D1 dimension does "
2847 "not match");
2848#endif
2849
2850 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions, false);
2851
2854 auto t_diff_n = getFTensor1FromPtr<D2>(&*diffNinvJac.data().begin());
2855 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2856 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*invJacPtr);
2857 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2858 for (size_t dd = 0; dd != nb_functions; ++dd) {
2859 t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
2860 ++t_diff_n;
2861 ++t_diff_n_ref;
2862 }
2863 }
2864
2865 diff_n.swap(diffNinvJac);
2866 }
2868 }
2869
2870 boost::shared_ptr<MatrixDouble> invJacPtr;
2872};
2873
2874template <>
2877
2879
2880 MoFEMErrorCode doWork(int side, EntityType type,
2882};
2883
2884template <>
2886 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2887
2888 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2889
2890 MoFEMErrorCode doWork(int side, EntityType type,
2892};
2893
2894template <>
2896 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2897
2898 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2899
2900 MoFEMErrorCode doWork(int side, EntityType type,
2902};
2903
2904template <int DERIVARIVE = 1>
2906 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2907 OpSetInvJacH1ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2908 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(H1, inv_jac_ptr) {}
2909};
2910
2911template <int DERIVARIVE = 1>
2913 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2914 OpSetInvJacL2ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2915 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(L2, inv_jac_ptr) {}
2916};
2917
2918template <int DERIVARIVE = 1>
2920 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2922 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2923 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(H1, inv_jac_ptr) {}
2924};
2925
2926template <int DERIVARIVE = 1>
2928 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2930 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2931 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(L2, inv_jac_ptr) {}
2932};
2933
2934/**
2935 * \brief Transform local reference derivatives of shape function to
2936 global derivatives for face
2937
2938 * \ingroup mofem_forces_and_sources_tri_element
2939 */
2940template <int DIM> struct OpSetInvJacHcurlFaceImpl;
2941
2942template <>
2945
2946 OpSetInvJacHcurlFaceImpl(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2948 invJacPtr(inv_jac_ptr) {}
2949
2950 MoFEMErrorCode doWork(int side, EntityType type,
2952
2953protected:
2954 boost::shared_ptr<MatrixDouble> invJacPtr;
2956};
2957
2958template <>
2960 using OpSetInvJacHcurlFaceImpl<2>::OpSetInvJacHcurlFaceImpl;
2961 MoFEMErrorCode doWork(int side, EntityType type,
2963};
2964
2967
2968/**
2969 * @brief Make Hdiv space from Hcurl space in 2d
2970 * @ingroup mofem_forces_and_sources_tri_element
2971 */
2974
2977
2978 MoFEMErrorCode doWork(int side, EntityType type,
2980};
2981
2982/** \brief Transform Hcurl base fluxes from reference element to physical
2983 * triangle
2984 */
2986
2987/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
2988 *
2989 * Covariant Piola transformation
2990 \f[
2991 \psi_i|_t = J^{-1}_{ij}\hat{\psi}_j\\
2992 \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
2993 = J^{-1}_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
2994 \f]
2995
2996 */
2997template <>
3000
3002 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3003
3004 MoFEMErrorCode doWork(int side, EntityType type,
3006
3007private:
3008 boost::shared_ptr<MatrixDouble> invJacPtr;
3009
3012};
3013
3016
3017/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3018 *
3019 * \note Hdiv space is generated by Hcurl space in 2d.
3020 *
3021 * Contravariant Piola transformation
3022 * \f[
3023 * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
3024 * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3025 * =
3026 * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3027 * \f]
3028 *
3029 * \ingroup mofem_forces_and_sources
3030 *
3031 */
3033
3034template <>
3037
3039 boost::shared_ptr<MatrixDouble> jac_ptr)
3041 jacPtr(jac_ptr) {}
3042
3043 MoFEMErrorCode doWork(int side, EntityType type,
3045
3046protected:
3047 boost::shared_ptr<MatrixDouble> jacPtr;
3050};
3051
3052template <>
3056 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3057
3058 MoFEMErrorCode doWork(int side, EntityType type,
3060};
3061
3066
3067/**@}*/
3068
3069/** \name Operators for edges */
3070
3071/**@{*/
3072
3075
3078
3079 MoFEMErrorCode doWork(int side, EntityType type,
3081
3082};
3083
3084/**
3085 * @deprecated Name is deprecated and this is added for back compatibility
3086 */
3089
3090/**@}*/
3091
3092/** \name Operator for fat prisms */
3093
3094/**@{*/
3095
3096/**
3097 * @brief Operator for fat prism element updating integration weights in the
3098 * volume.
3099 *
3100 * Jacobian on the distorted element is nonconstant. This operator updates
3101 * integration weight on prism to take into account nonconstat jacobian.
3102 *
3103 * \f[
3104 * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
3105 * \pmb\xi} \right\| \right)
3106 * \f]
3107 * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
3108 * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
3109 * element coordinate.
3110 *
3111 */
3114
3117
3118 MoFEMErrorCode doWork(int side, EntityType type,
3120};
3121
3122/** \brief Calculate inverse of jacobian for face element
3123
3124 It is assumed that face element is XY plane. Applied
3125 only for 2d problems.
3126
3127 FIXME Generalize function for arbitrary face orientation in 3d space
3128 FIXME Calculate to Jacobins for two faces
3129
3130 \ingroup mofem_forces_and_sources_prism_element
3131
3132*/
3135
3136 OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3138 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3139
3142 invJac(inv_jac) {}
3143 MoFEMErrorCode doWork(int side, EntityType type,
3145
3146private:
3147 const boost::shared_ptr<MatrixDouble> invJacPtr;
3149};
3150
3151/** \brief Transform local reference derivatives of shape functions to global
3152derivatives
3153
3154FIXME Generalize to curved shapes
3155FIXME Generalize to case that top and bottom face has different shape
3156
3157\ingroup mofem_forces_and_sources_prism_element
3158
3159*/
3162
3163 OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3165 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3166
3169 invJac(inv_jac) {}
3170
3171 MoFEMErrorCode doWork(int side, EntityType type,
3173
3174private:
3175 const boost::shared_ptr<MatrixDouble> invJacPtr;
3178};
3179
3180// Flat prism
3181
3182/** \brief Calculate inverse of jacobian for face element
3183
3184 It is assumed that face element is XY plane. Applied
3185 only for 2d problems.
3186
3187 FIXME Generalize function for arbitrary face orientation in 3d space
3188 FIXME Calculate to Jacobins for two faces
3189
3190 \ingroup mofem_forces_and_sources_prism_element
3191
3192*/
3195
3198 invJacF3(inv_jac_f3) {}
3199 MoFEMErrorCode doWork(int side, EntityType type,
3201
3202private:
3204};
3205
3206/** \brief Transform local reference derivatives of shape functions to global
3207derivatives
3208
3209FIXME Generalize to curved shapes
3210FIXME Generalize to case that top and bottom face has different shape
3211
3212\ingroup mofem_forces_and_sources_prism_element
3213
3214*/
3217
3220 invJacF3(inv_jac_f3) {}
3221
3222 MoFEMErrorCode doWork(int side, EntityType type,
3224
3225private:
3228};
3229
3230/**@}*/
3231
3232/** \name Operation on matrices at integration points */
3233
3234/**@{*/
3235
3236template <int DIM>
3238
3239 OpInvertMatrix(boost::shared_ptr<MatrixDouble> in_ptr,
3240 boost::shared_ptr<VectorDouble> det_ptr,
3241 boost::shared_ptr<MatrixDouble> out_ptr)
3243 outPtr(out_ptr), detPtr(det_ptr) {}
3244
3247 return doWorkImpl(side, type, data, FTensor::Number<DIM>());
3248 }
3249
3250private:
3251 boost::shared_ptr<MatrixDouble> inPtr;
3252 boost::shared_ptr<MatrixDouble> outPtr;
3253 boost::shared_ptr<VectorDouble> detPtr;
3254
3255 MoFEMErrorCode doWorkImpl(int side, EntityType type,
3257 const FTensor::Number<3> &);
3258
3259 MoFEMErrorCode doWorkImpl(int side, EntityType type,
3261 const FTensor::Number<2> &);
3262};
3263
3264template <int DIM>
3267 const FTensor::Number<3> &) {
3269
3270 if (!inPtr)
3271 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3272 "Pointer for inPtr matrix not allocated");
3273 if (!detPtr)
3274 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3275 "Pointer for detPtr matrix not allocated");
3276
3277 const auto nb_rows = inPtr->size1();
3278 const auto nb_integration_pts = inPtr->size2();
3279
3280 // Calculate determinant
3281 {
3282 detPtr->resize(nb_integration_pts, false);
3283 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3284 auto t_det = getFTensor0FromVec(*detPtr);
3285 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3286 determinantTensor3by3(t_in, t_det);
3287 ++t_in;
3288 ++t_det;
3289 }
3290 }
3291
3292 // Invert jacobian
3293 if (outPtr) {
3294 outPtr->resize(nb_rows, nb_integration_pts, false);
3295 auto t_in = getFTensor2FromMat<3, 3>(*inPtr);
3296 auto t_out = getFTensor2FromMat<3, 3>(*outPtr);
3297 auto t_det = getFTensor0FromVec(*detPtr);
3298 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3299 invertTensor3by3(t_in, t_det, t_out);
3300 ++t_in;
3301 ++t_out;
3302 ++t_det;
3303 }
3304 }
3305
3307}
3308
3309template <int DIM>
3312 const FTensor::Number<2> &) {
3314
3315 if (!inPtr)
3316 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3317 "Pointer for inPtr matrix not allocated");
3318 if (!detPtr)
3319 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3320 "Pointer for detPtr matrix not allocated");
3321
3322 const auto nb_rows = inPtr->size1();
3323 const auto nb_integration_pts = inPtr->size2();
3324
3325 // Calculate determinant
3326 {
3327 detPtr->resize(nb_integration_pts, false);
3328 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3329 auto t_det = getFTensor0FromVec(*detPtr);
3330 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3331 determinantTensor2by2(t_in, t_det);
3332 ++t_in;
3333 ++t_det;
3334 }
3335 }
3336
3337 // Invert jacobian
3338 if (outPtr) {
3339 outPtr->resize(nb_rows, nb_integration_pts, false);
3340 auto t_in = getFTensor2FromMat<2, 2>(*inPtr);
3341 auto t_out = getFTensor2FromMat<2, 2>(*outPtr);
3342 auto t_det = getFTensor0FromVec(*detPtr);
3343 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3344 invertTensor2by2(t_in, t_det, t_out);
3345 ++t_in;
3346 ++t_out;
3347 ++t_det;
3348 }
3349 }
3350
3352}
3353
3354/**@}*/
3355
3356} // namespace MoFEM
3357
3358#endif // __USER_DATA_OPERATORS_HPP__
3359
3360/**
3361 * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
3362 *
3363 * \brief Classes and functions used to evaluate fields at integration pts,
3364 *jacobians, etc..
3365 *
3366 * \ingroup mofem_forces_and_sources
3367 **/
static Index< 'J', 3 > J
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
static const char *const CoordinateTypesNames[]
Coordinate system names.
Definition: definitions.h:108
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
@ CYLINDRICAL
Definition: definitions.h:117
@ POLAR
Definition: definitions.h:116
@ CARTESIAN
Definition: definitions.h:115
@ SPHERICAL
Definition: definitions.h:118
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
FTensor::Index< 'n', SPACE_DIM > n
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp: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:1393
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1424
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:1363
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:1352
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
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate trace of vector (Hdiv/Hcurl) space.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Index< 'j', Tensor_Dim > j
FTensor::Index< 'i', Tensor_Dim > i
OpCalculateHVecVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
OpCalculateHVecVectorFieldDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
Get vector field for H-div approximation.
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecVectorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecVectorHessian(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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