v0.14.0
Loading...
Searching...
No Matches
UserDataOperators.hpp
Go to the documentation of this file.
1/** \file UserDataOperators.hpp
2 * \brief User data Operators
3
4*/
5
6#ifndef __USER_DATA_OPERATORS_HPP__
7#define __USER_DATA_OPERATORS_HPP__
8
9namespace MoFEM {
10
11/** \name Get values at Gauss pts */
12
13/**@{*/
14
15/** \name Scalar values */
16
17/**@{*/
18
19/** \brief Scalar field values at integration points
20 *
21 */
22template <class T, class A>
25
27 const std::string field_name,
28 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
29 const EntityType zero_type = MBVERTEX)
31 dataPtr(data_ptr), zeroType(zero_type) {
32 if (!dataPtr)
33 THROW_MESSAGE("Pointer is not set");
34 }
35
37 const std::string field_name,
38 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
39 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
41 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
42 if (!dataPtr)
43 THROW_MESSAGE("Pointer is not set");
44 }
45
46 /**
47 * \brief calculate values of scalar field at integration points
48 * @param side side entity number
49 * @param type side entity type
50 * @param data entity data
51 * @return error code
52 */
55
56protected:
57 boost::shared_ptr<ublas::vector<T, A>> dataPtr;
61};
62
63/**
64 * \brief Specialization of member function
65 *
66 */
67template <class T, class A>
69 int side, EntityType type, EntitiesFieldData::EntData &data) {
71 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented for T = %s",
72 typeid(T).name() // boost::core::demangle(typeid(T).name()).c_str()
73 );
75}
76
77/**
78 * \brief Get value at integration points for scalar field
79 *
80 * \ingroup mofem_forces_and_sources_user_data_operators
81 */
83 : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
84
86 double, DoubleAllocator>::OpCalculateScalarFieldValues_General;
87
88 /**
89 * \brief calculate values of scalar field at integration points
90 * @param side side entity number
91 * @param type side entity type
92 * @param data entity data
93 * @return error code
94 */
98 VectorDouble &vec = *dataPtr;
99 const size_t nb_gauss_pts = getGaussPts().size2();
100 if (type == zeroType || vec.size() != nb_gauss_pts) {
101 vec.resize(nb_gauss_pts, false);
102 vec.clear();
103 }
104
105 const size_t nb_dofs = data.getFieldData().size();
106
107 if (nb_dofs) {
108
109 if (dataVec.use_count()) {
110 dotVector.resize(nb_dofs, false);
111 const double *array;
112 CHKERR VecGetArrayRead(dataVec, &array);
113 const auto &local_indices = data.getLocalIndices();
114 for (int i = 0; i != local_indices.size(); ++i)
115 if (local_indices[i] != -1)
116 dotVector[i] = array[local_indices[i]];
117 else
118 dotVector[i] = 0;
119 CHKERR VecRestoreArrayRead(dataVec, &array);
120 data.getFieldData().swap(dotVector);
121 }
122
123 const size_t nb_base_functions = data.getN().size2();
124 auto base_function = data.getFTensor0N();
125 auto values_at_gauss_pts = getFTensor0FromVec(vec);
126 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
127 auto field_data = data.getFTensor0FieldData();
128 size_t bb = 0;
129 for (; bb != nb_dofs; ++bb) {
130 values_at_gauss_pts += field_data * base_function;
131 ++field_data;
132 ++base_function;
133 }
134 // It is possible to have more base functions than dofs
135 for (; bb < nb_base_functions; ++bb)
136 ++base_function;
137 ++values_at_gauss_pts;
138 }
139
140 if (dataVec.use_count()) {
141 data.getFieldData().swap(dotVector);
142 }
143 }
144
146 }
147};
148
149/**
150 * @brief Get rate of scalar field at integration points
151 *
152 * \ingroup mofem_forces_and_sources_user_data_operators
153 */
154
155template <PetscData::DataContext CTX>
158
160 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
161 const EntityType zero_at_type = MBVERTEX)
164 dataPtr(data_ptr), zeroAtType(zero_at_type) {
165 if (!dataPtr)
166 THROW_MESSAGE("Pointer is not set");
167 }
168
172
173 const size_t nb_gauss_pts = getGaussPts().size2();
174
175 VectorDouble &vec = *dataPtr;
176 if (type == zeroAtType || vec.size() != nb_gauss_pts) {
177 vec.resize(nb_gauss_pts, false);
178 vec.clear();
179 }
180
181 auto &local_indices = data.getLocalIndices();
182 const size_t nb_dofs = local_indices.size();
183 if (nb_dofs) {
184
185 const double *array;
186
187 auto get_array = [&](const auto ctx, auto vec) {
189#ifndef NDEBUG
190 if ((getFEMethod()->data_ctx & ctx).none()) {
191 MOFEM_LOG_CHANNEL("SELF");
192 MOFEM_LOG("SELF", Sev::error)
193 << "In this case filed degrees of freedom are read from vector. "
194 "That usually happens when time solver is used, and acces to "
195 "first or second rates is needed. You probably not set ts_u, "
196 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
197 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
198 "respectively";
199 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
200 }
201#endif
202 CHKERR VecGetArrayRead(vec, &array);
204 };
205
206 auto restore_array = [&](auto vec) {
207 return VecRestoreArrayRead(vec, &array);
208 };
209
210 switch (CTX) {
212 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
213 break;
215 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
216 break;
218 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
219 break;
220 default:
221 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
222 "That case is not implemented");
223 }
224
225 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
226 for (int i = 0; i != local_indices.size(); ++i)
227 if (local_indices[i] != -1)
228 dot_dofs_vector[i] = array[local_indices[i]];
229 else
230 dot_dofs_vector[i] = 0;
231
232 switch (CTX) {
234 CHKERR restore_array(getFEMethod()->ts_u);
235 break;
237 CHKERR restore_array(getFEMethod()->ts_u_t);
238 break;
240 CHKERR restore_array(getFEMethod()->ts_u_tt);
241 break;
242 default:
243 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
244 "That case is not implemented");
245 }
246
247 const size_t nb_base_functions = data.getN().size2();
248 auto base_function = data.getFTensor0N();
249 auto values_at_gauss_pts = getFTensor0FromVec(vec);
250
251 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
252 size_t bb = 0;
253 for (; bb != nb_dofs; ++bb) {
254 values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
255 ++base_function;
256 }
257 // Number of dofs can be smaller than number of Tensor_Dim x base
258 // functions
259 for (; bb < nb_base_functions; ++bb)
260 ++base_function;
261 ++values_at_gauss_pts;
262 }
263 }
265 }
266
267private:
268 boost::shared_ptr<VectorDouble> dataPtr;
270};
271
276
277/**
278 * \deprecated Name inconsistent with other operators
279 *
280 */
282
283/**@}*/
284
285/** \name Vector field values at integration points */
286
287/**@{*/
288
289/** \brief Calculate field values for tenor field rank 1, i.e. vector field
290 *
291 */
292template <int Tensor_Dim, class T, class L, class A>
295
297 const std::string field_name,
298 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
299 const EntityType zero_type = MBVERTEX)
302 dataPtr(data_ptr), zeroType(zero_type) {
303 if (!dataPtr)
304 THROW_MESSAGE("Pointer is not set");
305 }
306
307 /**
308 * \brief calculate values of vector field at integration points
309 * @param side side entity number
310 * @param type side entity type
311 * @param data entity data
312 * @return error code
313 */
314 MoFEMErrorCode doWork(int side, EntityType type,
316
317protected:
318 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
320};
321
322template <int Tensor_Dim, class T, class L, class A>
325 int side, EntityType type, EntitiesFieldData::EntData &data) {
327 SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
328 "Not implemented for T = %s and dim = %d",
329 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
330 Tensor_Dim);
332}
333
334/** \brief Calculate field values (template specialization) for tensor field
335 * rank 1, i.e. vector field
336 *
337 */
338template <int Tensor_Dim>
340 ublas::row_major, DoubleAllocator>
342
344 boost::shared_ptr<MatrixDouble> data_ptr,
345 const EntityType zero_type = MBVERTEX)
348 dataPtr(data_ptr), zeroType(zero_type) {
349 if (!dataPtr)
350 THROW_MESSAGE("Pointer is not set");
351 }
352
354 boost::shared_ptr<MatrixDouble> data_ptr,
355 SmartPetscObj<Vec> data_vec,
356 const EntityType zero_type = MBVERTEX)
359 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type) {
360 if (!dataPtr)
361 THROW_MESSAGE("Pointer is not set");
362 }
363
364 MoFEMErrorCode doWork(int side, EntityType type,
366
367protected:
368 boost::shared_ptr<MatrixDouble> dataPtr;
372};
373
374/**
375 * \brief Member function specialization calculating values for tenor field rank
376 *
377 */
378template <int Tensor_Dim>
380 Tensor_Dim, double, ublas::row_major,
381 DoubleAllocator>::doWork(int side, EntityType type,
384
385 const size_t nb_gauss_pts = getGaussPts().size2();
386 auto &mat = *dataPtr;
387 if (type == zeroType || mat.size2() != nb_gauss_pts) {
388 mat.resize(Tensor_Dim, nb_gauss_pts, false);
389 mat.clear();
390 }
391
392 const size_t nb_dofs = data.getFieldData().size();
393 if (nb_dofs) {
394
395 if (dataVec.use_count()) {
396 dotVector.resize(nb_dofs, false);
397 const double *array;
398 CHKERR VecGetArrayRead(dataVec, &array);
399 const auto &local_indices = data.getLocalIndices();
400 for (int i = 0; i != local_indices.size(); ++i)
401 if (local_indices[i] != -1)
402 dotVector[i] = array[local_indices[i]];
403 else
404 dotVector[i] = 0;
405 CHKERR VecRestoreArrayRead(dataVec, &array);
406 data.getFieldData().swap(dotVector);
407 }
408
409 if (nb_gauss_pts) {
410 const size_t nb_base_functions = data.getN().size2();
411 auto base_function = data.getFTensor0N();
412 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
414 const size_t size = nb_dofs / Tensor_Dim;
415 if (nb_dofs % Tensor_Dim) {
416 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
417 "Data inconsistency");
418 }
419 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
420 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
421
422#ifndef NDEBUG
423 if (field_data.l2() != field_data.l2()) {
424 MOFEM_LOG("SELF", Sev::error) << "field data: " << field_data;
425 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
426 "Wrong number in coefficients");
427 }
428#endif
429
430 size_t bb = 0;
431 for (; bb != size; ++bb) {
432
433#ifndef NDEBUG
434 if (base_function != base_function) {
435 MOFEM_LOG("SELF", Sev::error) << "base finction: " << base_function;
436 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
437 "Wrong number number in base functions");
438 }
439#endif
440
441 values_at_gauss_pts(I) += field_data(I) * base_function;
442 ++field_data;
443 ++base_function;
444 }
445 // Number of dofs can be smaller than number of Tensor_Dim x base
446 // functions
447 for (; bb < nb_base_functions; ++bb)
448 ++base_function;
449 ++values_at_gauss_pts;
450 }
451 }
452
453 if (dataVec.use_count()) {
454 data.getFieldData().swap(dotVector);
455 }
456 }
458}
459
460/** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
461 * field
462 *
463 * \ingroup mofem_forces_and_sources_user_data_operators
464 */
465template <int Tensor_Dim>
468 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
469
471 Tensor_Dim, double, ublas::row_major,
472 DoubleAllocator>::OpCalculateVectorFieldValues_General;
473};
474
475/**@}*/
476
477/** \name Vector field values at integration points */
478
479/**@{*/
480
481/** \brief Calculate field values (template specialization) for tensor field
482 * rank 1, i.e. vector field
483 *
484 */
485template <int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
488
490 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
491 const EntityType zero_type = MBVERTEX)
494 dataPtr(data_ptr), zeroType(zero_type) {
495 if (!dataPtr)
496 THROW_MESSAGE("Pointer is not set");
497 }
498
502
503 // When we move to C++17 add if constexpr()
504 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
505 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
506 "%s coordiante not implemented",
507 CoordinateTypesNames[COORDINATE_SYSTEM]);
508
509 const size_t nb_gauss_pts = getGaussPts().size2();
510 auto &vec = *dataPtr;
511 if (type == zeroType) {
512 vec.resize(nb_gauss_pts, false);
513 vec.clear();
514 }
515
516 const size_t nb_dofs = data.getFieldData().size();
517 if (nb_dofs) {
518
519 if (nb_gauss_pts) {
520 const size_t nb_base_functions = data.getN().size2();
521 auto values_at_gauss_pts = getFTensor0FromVec(vec);
523 const size_t size = nb_dofs / Tensor_Dim;
524#ifndef NDEBUG
525 if (nb_dofs % Tensor_Dim) {
526 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
527 "Number of dofs should multiple of dimensions");
528 }
529#endif
530
531 // When we move to C++17 add if constexpr()
532 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
533 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
534 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
535 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
536 size_t bb = 0;
537 for (; bb != size; ++bb) {
538 values_at_gauss_pts += field_data(I) * diff_base_function(I);
539 ++field_data;
540 ++diff_base_function;
541 }
542 // Number of dofs can be smaller than number of Tensor_Dim x base
543 // functions
544 for (; bb < nb_base_functions; ++bb)
545 ++diff_base_function;
546 ++values_at_gauss_pts;
547 }
548 }
549
550 // When we move to C++17 add if constexpr()
551 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
552 auto t_coords = getFTensor1CoordsAtGaussPts();
553 auto values_at_gauss_pts = getFTensor0FromVec(vec);
554 auto base_function = data.getFTensor0N();
555 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
556 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
557 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
558 size_t bb = 0;
559 for (; bb != size; ++bb) {
560 values_at_gauss_pts += field_data(I) * diff_base_function(I);
561 values_at_gauss_pts +=
562 field_data(0) * base_function / t_coords(0);
563 ++field_data;
564 ++base_function;
565 ++diff_base_function;
566 }
567 // Number of dofs can be smaller than number of Tensor_Dim x base
568 // functions
569 for (; bb < nb_base_functions; ++bb) {
570 ++base_function;
571 ++diff_base_function;
572 }
573 ++values_at_gauss_pts;
574 ++t_coords;
575 }
576 }
577 }
578 }
580 }
581
582protected:
583 boost::shared_ptr<VectorDouble> dataPtr;
585};
586
587/** \brief Approximate field values for given petsc vector
588 *
589 * \note Look at PetscData to see what vectors could be extracted with that user
590 * data operator.
591 *
592 * \ingroup mofem_forces_and_sources_user_data_operators
593 */
594template <int Tensor_Dim, PetscData::DataContext CTX>
597
599 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
600 const EntityType zero_at_type = MBVERTEX)
603 dataPtr(data_ptr), zeroAtType(zero_at_type) {
604 if (!dataPtr)
605 THROW_MESSAGE("Pointer is not set");
606 }
607
611 auto &local_indices = data.getLocalIndices();
612 const size_t nb_dofs = local_indices.size();
613 const size_t nb_gauss_pts = getGaussPts().size2();
614
615 MatrixDouble &mat = *dataPtr;
616 if (type == zeroAtType) {
617 mat.resize(Tensor_Dim, nb_gauss_pts, false);
618 mat.clear();
619 }
620 if (!nb_dofs)
622
623 const double *array;
624
625 auto get_array = [&](const auto ctx, auto vec) {
627#ifndef NDEBUG
628 if ((getFEMethod()->data_ctx & ctx).none()) {
629 MOFEM_LOG_CHANNEL("SELF");
630 MOFEM_LOG("SELF", Sev::error)
631 << "In this case filed degrees of freedom are read from vector. "
632 "That usually happens when time solver is used, and acces to "
633 "first or second rates is needed. You probably not set ts_u, "
634 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
635 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
636 "respectively";
637 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
638 }
639#endif
640 CHKERR VecGetArrayRead(vec, &array);
642 };
643
644 auto restore_array = [&](auto vec) {
645 return VecRestoreArrayRead(vec, &array);
646 };
647
648 switch (CTX) {
650 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
651 break;
653 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
654 break;
656 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
657 break;
658 default:
659 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
660 "That case is not implemented");
661 }
662
663 dotVector.resize(local_indices.size());
664 for (int i = 0; i != local_indices.size(); ++i)
665 if (local_indices[i] != -1)
666 dotVector[i] = array[local_indices[i]];
667 else
668 dotVector[i] = 0;
669
670 switch (CTX) {
672 CHKERR restore_array(getFEMethod()->ts_u);
673 break;
675 CHKERR restore_array(getFEMethod()->ts_u_t);
676 break;
678 CHKERR restore_array(getFEMethod()->ts_u_tt);
679 break;
680 default:
681 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
682 "That case is not implemented");
683 }
684
685 const size_t nb_base_functions = data.getN().size2();
686 auto base_function = data.getFTensor0N();
687 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
688
690 const size_t size = nb_dofs / Tensor_Dim;
691 if (nb_dofs % Tensor_Dim) {
692 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
693 }
694 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
695 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(dotVector);
696 size_t bb = 0;
697 for (; bb != size; ++bb) {
698 values_at_gauss_pts(I) += field_data(I) * base_function;
699 ++field_data;
700 ++base_function;
701 }
702 // Number of dofs can be smaller than number of Tensor_Dim x base
703 // functions
704 for (; bb < nb_base_functions; ++bb)
705 ++base_function;
706 ++values_at_gauss_pts;
707 }
709 }
710
711protected:
712 boost::shared_ptr<MatrixDouble> dataPtr;
715};
716
717/** \brief Get time derivatives of values at integration pts for tensor filed
718 * rank 1, i.e. vector field
719 *
720 * \ingroup mofem_forces_and_sources_user_data_operators
721 */
722template <int Tensor_Dim>
726
727/** \brief Get second time derivatives of values at integration pts for tensor
728 * filed rank 1, i.e. vector field
729 *
730 * \ingroup mofem_forces_and_sources_user_data_operators
731 */
732template <int Tensor_Dim>
736
737/**@}*/
738
739/** \name Tensor field values at integration points */
740
741/**@{*/
742
743/** \brief Calculate field values for tenor field rank 2.
744 *
745 */
746template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
749
751 const std::string field_name,
752 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
753 const EntityType zero_type = MBVERTEX)
756 dataPtr(data_ptr), zeroType(zero_type) {
757 if (!dataPtr)
758 THROW_MESSAGE("Pointer is not set");
759 }
760
761 MoFEMErrorCode doWork(int side, EntityType type,
763
764protected:
765 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
767};
768
769template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
772 doWork(int side, EntityType type, EntitiesFieldData::EntData &data) {
774 SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
775 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
776 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
777 Tensor_Dim0, Tensor_Dim1);
779}
780
781template <int Tensor_Dim0, int Tensor_Dim1>
782struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
783 ublas::row_major, DoubleAllocator>
785
787 const std::string field_name,
788 boost::shared_ptr<
789 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
790 data_ptr,
791 const EntityType zero_type = MBVERTEX)
794 dataPtr(data_ptr), zeroType(zero_type) {
795 if (!dataPtr)
796 THROW_MESSAGE("Pointer is not set");
797 }
798
800 const std::string field_name,
801 boost::shared_ptr<
802 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
803 data_ptr,
804 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
807 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
808 if (!dataPtr)
809 THROW_MESSAGE("Pointer is not set");
810 }
811
812 MoFEMErrorCode doWork(int side, EntityType type,
814
815protected:
816 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
821};
822
823template <int Tensor_Dim0, int Tensor_Dim1>
825 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
826 DoubleAllocator>::doWork(int side, EntityType type,
829
830 MatrixDouble &mat = *dataPtr;
831 const size_t nb_gauss_pts = data.getN().size1();
832 if (type == zeroType) {
833 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
834 mat.clear();
835 }
836
837 const size_t nb_dofs = data.getFieldData().size();
838
839 if (dataVec.use_count()) {
840 dotVector.resize(nb_dofs, false);
841 const double *array;
842 CHKERR VecGetArrayRead(dataVec, &array);
843 const auto &local_indices = data.getLocalIndices();
844 for (int i = 0; i != local_indices.size(); ++i)
845 if (local_indices[i] != -1)
846 dotVector[i] = array[local_indices[i]];
847 else
848 dotVector[i] = 0;
849 CHKERR VecRestoreArrayRead(dataVec, &array);
850 data.getFieldData().swap(dotVector);
851 }
852
853 if (nb_dofs) {
854 const size_t nb_base_functions = data.getN().size2();
855 auto base_function = data.getFTensor0N();
856 auto values_at_gauss_pts =
857 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
860 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
861 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
862 auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
863 size_t bb = 0;
864 for (; bb != size; ++bb) {
865 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
866 ++field_data;
867 ++base_function;
868 }
869 for (; bb < nb_base_functions; ++bb)
870 ++base_function;
871 ++values_at_gauss_pts;
872 }
873
874 if (dataVec.use_count()) {
875 data.getFieldData().swap(dotVector);
876 }
877 }
879}
880
881/** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
882 * field
883 *
884 * \ingroup mofem_forces_and_sources_user_data_operators
885 */
886template <int Tensor_Dim0, int Tensor_Dim1>
889 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
890
892 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
893 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
894};
895
896/** \brief Get time direvarive values at integration pts for tensor filed rank
897 * 2, i.e. matrix field
898 *
899 * \ingroup mofem_forces_and_sources_user_data_operators
900 */
901template <int Tensor_Dim0, int Tensor_Dim1>
904
906 boost::shared_ptr<MatrixDouble> data_ptr,
907 const EntityType zero_at_type = MBVERTEX)
910 dataPtr(data_ptr), zeroAtType(zero_at_type) {
911 if (!dataPtr)
912 THROW_MESSAGE("Pointer is not set");
913 }
914
918
919 const size_t nb_gauss_pts = getGaussPts().size2();
920 MatrixDouble &mat = *dataPtr;
921 if (type == zeroAtType) {
922 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
923 mat.clear();
924 }
925 const auto &local_indices = data.getLocalIndices();
926 const size_t nb_dofs = local_indices.size();
927 if (nb_dofs) {
928 dotVector.resize(nb_dofs, false);
929 const double *array;
930 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
931 for (size_t i = 0; i != local_indices.size(); ++i)
932 if (local_indices[i] != -1)
933 dotVector[i] = array[local_indices[i]];
934 else
935 dotVector[i] = 0;
936 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
937
938 const size_t nb_base_functions = data.getN().size2();
939
940 auto base_function = data.getFTensor0N();
941 auto values_at_gauss_pts =
942 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
945 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
946 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
947 auto field_data =
948 getFTensor2FromPtr<Tensor_Dim0, Tensor_Dim1>(&*dotVector.begin());
949 size_t bb = 0;
950 for (; bb != size; ++bb) {
951 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
952 ++field_data;
953 ++base_function;
954 }
955 for (; bb < nb_base_functions; ++bb)
956 ++base_function;
957 ++values_at_gauss_pts;
958 }
959 }
961 }
962
963protected:
964 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
965 EntityType zeroAtType; ///< Zero values at Gauss point at this type
966 VectorDouble dotVector; ///< Keeps temoorary values of time directives
967
968};
969
970/**
971 * @brief Calculate symmetric tensor field values at integration pts.
972 *
973 * @tparam Tensor_Dim
974
975 * \ingroup mofem_forces_and_sources_user_data_operators
976 */
977template <int Tensor_Dim>
980
982 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
983 const EntityType zero_type = MBEDGE, const int zero_side = 0)
986 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
987 if (!dataPtr)
988 THROW_MESSAGE("Pointer is not set");
989 }
990
992 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
993 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBEDGE,
994 const int zero_side = 0)
997 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side),
998 dataVec(data_vec) {
999 if (!dataPtr)
1000 THROW_MESSAGE("Pointer is not set");
1001 }
1002
1006 MatrixDouble &mat = *dataPtr;
1007 const int nb_gauss_pts = getGaussPts().size2();
1008 if (type == this->zeroType && side == zeroSide) {
1009 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
1010 mat.clear();
1011 }
1012 const int nb_dofs = data.getFieldData().size();
1013 if (!nb_dofs)
1015
1016 if (dataVec.use_count()) {
1017 dotVector.resize(nb_dofs, false);
1018 const double *array;
1019 CHKERR VecGetArrayRead(dataVec, &array);
1020 const auto &local_indices = data.getLocalIndices();
1021 for (int i = 0; i != local_indices.size(); ++i)
1022 if (local_indices[i] != -1)
1023 dotVector[i] = array[local_indices[i]];
1024 else
1025 dotVector[i] = 0;
1026 CHKERR VecRestoreArrayRead(dataVec, &array);
1027 data.getFieldData().swap(dotVector);
1028 }
1029
1030 const int nb_base_functions = data.getN().size2();
1031 auto base_function = data.getFTensor0N();
1032 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1035 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1036 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1037 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
1038 int bb = 0;
1039 for (; bb != size; ++bb) {
1040 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1041 ++field_data;
1042 ++base_function;
1043 }
1044 for (; bb < nb_base_functions; ++bb)
1045 ++base_function;
1046 ++values_at_gauss_pts;
1047 }
1048
1049 if (dataVec.use_count()) {
1050 data.getFieldData().swap(dotVector);
1051 }
1052
1054 }
1055
1056protected:
1057 boost::shared_ptr<MatrixDouble> dataPtr;
1059 const int zeroSide;
1062};
1063
1064/**
1065 * @brief Calculate symmetric tensor field rates ant integratio pts.
1066 *
1067 * @tparam Tensor_Dim
1068 *
1069 * \ingroup mofem_forces_and_sources_user_data_operators
1070 */
1071template <int Tensor_Dim>
1074
1076 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1077 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1080 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1081 if (!dataPtr)
1082 THROW_MESSAGE("Pointer is not set");
1083 }
1084
1088 const int nb_gauss_pts = getGaussPts().size2();
1089 MatrixDouble &mat = *dataPtr;
1090 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1091 if (type == zeroType && side == zeroSide) {
1092 mat.resize(symm_size, nb_gauss_pts, false);
1093 mat.clear();
1094 }
1095 auto &local_indices = data.getLocalIndices();
1096 const int nb_dofs = local_indices.size();
1097 if (!nb_dofs)
1099
1100#ifndef NDEBUG
1101 if ((getFEMethod()->data_ctx & PetscData::CtxSetX_T).none()) {
1102 MOFEM_LOG_CHANNEL("SELF");
1103 MOFEM_LOG("SELF", Sev::error)
1104 << "In this case filed degrees of freedom are read from vector. "
1105 "That usually happens when time solver is used, and acces to "
1106 "first rates is needed. You probably not set "
1107 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1108 "respectively";
1109 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
1110 }
1111#endif
1112
1113 dotVector.resize(nb_dofs, false);
1114 const double *array;
1115 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1116 for (int i = 0; i != local_indices.size(); ++i)
1117 if (local_indices[i] != -1)
1118 dotVector[i] = array[local_indices[i]];
1119 else
1120 dotVector[i] = 0;
1121 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1122
1123 const int nb_base_functions = data.getN().size2();
1124
1125 auto base_function = data.getFTensor0N();
1126 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1129 const int size = nb_dofs / symm_size;
1130 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1131 auto field_data = getFTensorDotData<Tensor_Dim>();
1132 int bb = 0;
1133 for (; bb != size; ++bb) {
1134 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1135 ++field_data;
1136 ++base_function;
1137 }
1138 for (; bb < nb_base_functions; ++bb)
1139 ++base_function;
1140 ++values_at_gauss_pts;
1141 }
1142
1144 }
1145
1146protected:
1147 boost::shared_ptr<MatrixDouble> dataPtr;
1149 const int zeroSide;
1151
1152 template <int Dim> inline auto getFTensorDotData() {
1153 static_assert(Dim || !Dim, "not implemented");
1154 }
1155};
1156
1157template <>
1158template <>
1159inline auto
1162 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1163 &dotVector[5]);
1164}
1165
1166template <>
1167template <>
1168inline auto
1171 &dotVector[0], &dotVector[1], &dotVector[2]);
1172}
1173
1174/**@}*/
1175
1176/** \name Gradients and Hessian of scalar fields at integration points */
1177
1178/**@{*/
1179
1180/**
1181 * \brief Evaluate field gradient values for scalar field, i.e. gradient is
1182 * tensor rank 1 (vector)
1183 *
1184 */
1185template <int Tensor_Dim, class T, class L, class A>
1187 : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
1188
1190 Tensor_Dim, T, L, A>::OpCalculateVectorFieldValues_General;
1191};
1192
1193/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1194 * tensor rank 1 (vector), specialization
1195 *
1196 */
1197template <int Tensor_Dim>
1199 ublas::row_major, DoubleAllocator>
1201 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1202
1204 Tensor_Dim, double, ublas::row_major,
1205 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1206
1207 /**
1208 * \brief calculate gradient values of scalar field at integration points
1209 * @param side side entity number
1210 * @param type side entity type
1211 * @param data entity data
1212 * @return error code
1213 */
1214 MoFEMErrorCode doWork(int side, EntityType type,
1216};
1217
1218/**
1219 * \brief Member function specialization calculating scalar field gradients for
1220 * tenor field rank 1
1221 *
1222 */
1223template <int Tensor_Dim>
1225 Tensor_Dim, double, ublas::row_major,
1226 DoubleAllocator>::doWork(int side, EntityType type,
1229
1230 const size_t nb_gauss_pts = this->getGaussPts().size2();
1231 auto &mat = *this->dataPtr;
1232 if (type == this->zeroType) {
1233 mat.resize(Tensor_Dim, nb_gauss_pts, false);
1234 mat.clear();
1235 }
1236
1237 const int nb_dofs = data.getFieldData().size();
1238 if (nb_dofs) {
1239
1240 if (nb_gauss_pts) {
1241 const int nb_base_functions = data.getN().size2();
1242 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
1243 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1244
1245#ifndef NDEBUG
1246 if (nb_dofs > nb_base_functions)
1247 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1248 "Number of base functions inconsistent with number of DOFs "
1249 "(%d > %d)",
1250 nb_dofs, nb_base_functions);
1251
1252 if (data.getDiffN().size2() != nb_base_functions * Tensor_Dim)
1253 SETERRQ2(
1254 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1255 "Number of base functions inconsistent with number of derivatives "
1256 "(%d != %d)",
1257 data.getDiffN().size2(), nb_base_functions);
1258
1259 if (data.getDiffN().size1() != nb_gauss_pts)
1260 SETERRQ2(
1261 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1262 "Number of base functions inconsistent with number of integration "
1263 "pts (%d != %d)",
1264 data.getDiffN().size2(), nb_gauss_pts);
1265
1266#endif
1267
1269 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1270 auto field_data = data.getFTensor0FieldData();
1271 int bb = 0;
1272 for (; bb != nb_dofs; ++bb) {
1273 gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
1274 ++field_data;
1275 ++diff_base_function;
1276 }
1277 // Number of dofs can be smaller than number of base functions
1278 for (; bb < nb_base_functions; ++bb)
1279 ++diff_base_function;
1280 ++gradients_at_gauss_pts;
1281 }
1282 }
1283 }
1284
1286}
1287
1288/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1289 * vector field
1290 *
1291 * \ingroup mofem_forces_and_sources_user_data_operators
1292 */
1293template <int Tensor_Dim>
1296 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1298 Tensor_Dim, double, ublas::row_major,
1299 DoubleAllocator>::OpCalculateScalarFieldGradient_General;
1300};
1301
1302/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1303 * tensor rank 1 (vector), specialization
1304 *
1305 */
1306template <int Tensor_Dim>
1309 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1310
1312 Tensor_Dim, double, ublas::row_major,
1313 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1314
1315 /**
1316 * \brief calculate gradient values of scalar field at integration points
1317 * @param side side entity number
1318 * @param type side entity type
1319 * @param data entity data
1320 * @return error code
1321 */
1322 MoFEMErrorCode doWork(int side, EntityType type,
1324};
1325
1326template <int Tensor_Dim>
1328 int side, EntityType type, EntitiesFieldData::EntData &data) {
1330
1331 const size_t nb_gauss_pts = this->getGaussPts().size2();
1332 auto &mat = *this->dataPtr;
1333 if (type == this->zeroType) {
1334 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts, false);
1335 mat.clear();
1336 }
1337
1338 const int nb_dofs = data.getFieldData().size();
1339 if (nb_dofs) {
1340
1341 if (nb_gauss_pts) {
1342 const int nb_base_functions = data.getN().size2();
1343
1344 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1345#ifndef NDEBUG
1346 if (hessian_base.size1() != nb_gauss_pts) {
1347 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1348 "Wrong number of integration pts (%d != %d)",
1349 hessian_base.size1(), nb_gauss_pts);
1350 }
1351 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1352 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1353 "Wrong number of base functions (%d != %d)",
1354 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1355 nb_base_functions);
1356 }
1357 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1358 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1359 "Wrong number of base functions (%d < %d)",
1360 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1361 }
1362#endif
1363
1364 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1365 &*hessian_base.data().begin());
1366
1367 auto hessian_at_gauss_pts =
1368 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1369
1372 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1373 auto field_data = data.getFTensor0FieldData();
1374 int bb = 0;
1375 for (; bb != nb_dofs; ++bb) {
1376 hessian_at_gauss_pts(I, J) +=
1377 field_data * t_diff2_base_function(I, J);
1378 ++field_data;
1379 ++t_diff2_base_function;
1380 }
1381 // Number of dofs can be smaller than number of base functions
1382 for (; bb < nb_base_functions; ++bb) {
1383 ++t_diff2_base_function;
1384 }
1385
1386 ++hessian_at_gauss_pts;
1387 }
1388 }
1389 }
1390
1392}
1393
1394/**}*/
1395
1396/** \name Gradients and hessian of tensor fields at integration points */
1397
1398/**@{*/
1399
1400/**
1401 * \brief Evaluate field gradient values for vector field, i.e. gradient is
1402 * tensor rank 2
1403 *
1404 */
1405template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1407 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1408 L, A> {
1409
1411 Tensor_Dim0, Tensor_Dim1, T, L, A>::OpCalculateTensor2FieldValues_General;
1412};
1413
1414template <int Tensor_Dim0, int Tensor_Dim1>
1415struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1416 ublas::row_major, DoubleAllocator>
1418 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1419
1421 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1422 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1423
1424 /**
1425 * \brief calculate values of vector field at integration points
1426 * @param side side entity number
1427 * @param type side entity type
1428 * @param data entity data
1429 * @return error code
1430 */
1431 MoFEMErrorCode doWork(int side, EntityType type,
1433};
1434
1435/**
1436 * \brief Member function specialization calculating vector field gradients for
1437 * tenor field rank 2
1438 *
1439 */
1440template <int Tensor_Dim0, int Tensor_Dim1>
1442 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1443 DoubleAllocator>::doWork(int side, EntityType type,
1446 if (!this->dataPtr)
1447 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1448 "Data pointer not allocated");
1449
1450 const size_t nb_gauss_pts = this->getGaussPts().size2();
1451 auto &mat = *this->dataPtr;
1452 if (type == this->zeroType) {
1453 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1454 mat.clear();
1455 }
1456
1457 if (nb_gauss_pts) {
1458 const size_t nb_dofs = data.getFieldData().size();
1459
1460 if (nb_dofs) {
1461
1462 if (this->dataVec.use_count()) {
1463 this->dotVector.resize(nb_dofs, false);
1464 const double *array;
1465 CHKERR VecGetArrayRead(this->dataVec, &array);
1466 const auto &local_indices = data.getLocalIndices();
1467 for (int i = 0; i != local_indices.size(); ++i)
1468 if (local_indices[i] != -1)
1469 this->dotVector[i] = array[local_indices[i]];
1470 else
1471 this->dotVector[i] = 0;
1472 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1473 data.getFieldData().swap(this->dotVector);
1474 }
1475
1476 const int nb_base_functions = data.getN().size2();
1477 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1478 auto gradients_at_gauss_pts =
1479 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1482 int size = nb_dofs / Tensor_Dim0;
1483 if (nb_dofs % Tensor_Dim0) {
1484 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1485 "Data inconsistency");
1486 }
1487 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1488 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1489
1490#ifndef NDEBUG
1491 if (field_data.l2() != field_data.l2()) {
1492 MOFEM_LOG("SELF", Sev::error) << "field data " << field_data;
1493 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1494 "Wrong number in coefficients");
1495 }
1496#endif
1497
1498 int bb = 0;
1499 for (; bb < size; ++bb) {
1500#ifndef SINGULARITY
1501#ifndef NDEBUG
1502 if (diff_base_function.l2() != diff_base_function.l2()) {
1503 MOFEM_LOG("SELF", Sev::error)
1504 << "diff_base_function: " << diff_base_function;
1505 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1506 "Wrong number number in base functions");
1507 }
1508#endif
1509#endif
1510
1511 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1512 ++field_data;
1513 ++diff_base_function;
1514 }
1515 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1516 // functions
1517 for (; bb != nb_base_functions; ++bb)
1518 ++diff_base_function;
1519 ++gradients_at_gauss_pts;
1520 }
1521
1522 if (this->dataVec.use_count()) {
1523 data.getFieldData().swap(this->dotVector);
1524 }
1525 }
1526 }
1528}
1529
1530/** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1531 * vector field
1532 *
1533 * \ingroup mofem_forces_and_sources_user_data_operators
1534 */
1535template <int Tensor_Dim0, int Tensor_Dim1>
1538 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1539
1541 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1542 DoubleAllocator>::OpCalculateVectorFieldGradient_General;
1543};
1544
1545/** \brief Get field gradients time derivative at integration pts for scalar
1546 * filed rank 0, i.e. vector field
1547 *
1548 * \ingroup mofem_forces_and_sources_user_data_operators
1549 */
1550template <int Tensor_Dim0, int Tensor_Dim1>
1553
1555 boost::shared_ptr<MatrixDouble> data_ptr,
1556 const EntityType zero_at_type = MBVERTEX)
1559 dataPtr(data_ptr), zeroAtType(zero_at_type) {
1560 if (!dataPtr)
1561 THROW_MESSAGE("Pointer is not set");
1562 }
1563
1567
1568 const auto &local_indices = data.getLocalIndices();
1569 const int nb_dofs = local_indices.size();
1570 const int nb_gauss_pts = this->getGaussPts().size2();
1571
1572 auto &mat = *this->dataPtr;
1573 if (type == this->zeroAtType) {
1574 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1575 mat.clear();
1576 }
1577 if (!nb_dofs)
1579
1580 dotVector.resize(nb_dofs, false);
1581 const double *array;
1582 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1583 for (int i = 0; i != local_indices.size(); ++i)
1584 if (local_indices[i] != -1)
1585 dotVector[i] = array[local_indices[i]];
1586 else
1587 dotVector[i] = 0;
1588 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1589
1590 const int nb_base_functions = data.getN().size2();
1591 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1592 auto gradients_at_gauss_pts =
1593 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1596 int size = nb_dofs / Tensor_Dim0;
1597 if (nb_dofs % Tensor_Dim0) {
1598 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1599 }
1600
1601 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1602 auto field_data = getFTensor1FromPtr<Tensor_Dim0>(&*dotVector.begin());
1603 int bb = 0;
1604 for (; bb < size; ++bb) {
1605 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1606 ++field_data;
1607 ++diff_base_function;
1608 }
1609 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1610 // functions
1611 for (; bb != nb_base_functions; ++bb)
1612 ++diff_base_function;
1613 ++gradients_at_gauss_pts;
1614 }
1616 }
1617
1618private:
1619 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1620 EntityType zeroAtType; ///< Zero values at Gauss point at this type
1621 VectorDouble dotVector; ///< Keeps temoorary values of time directives
1622
1623};
1624
1625/**
1626 * \brief Evaluate field gradient values for symmetric 2nd order tensor field,
1627 * i.e. gradient is tensor rank 3
1628 *
1629 */
1630template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1632 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1633 L, A> {
1634
1636 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1637 const EntityType zero_type = MBVERTEX)
1638 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
1639 A>(field_name, data_ptr,
1640 zero_type) {}
1641};
1642
1643template <int Tensor_Dim0, int Tensor_Dim1>
1645 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator>
1647 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1648
1650 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1651 const EntityType zero_type = MBVERTEX)
1652 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1653 ublas::row_major,
1655 field_name, data_ptr, zero_type) {}
1656
1657 /**
1658 * \brief calculate values of vector field at integration points
1659 * @param side side entity number
1660 * @param type side entity type
1661 * @param data entity data
1662 * @return error code
1663 */
1664 MoFEMErrorCode doWork(int side, EntityType type,
1666};
1667
1668/**
1669 * \brief Member function specialization calculating tensor field gradients for
1670 * symmetric tensor field rank 2
1671 *
1672 */
1673template <int Tensor_Dim0, int Tensor_Dim1>
1674MoFEMErrorCode OpCalculateTensor2SymmetricFieldGradient_General<
1675 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1676 DoubleAllocator>::doWork(int side, EntityType type,
1679 if (!this->dataPtr)
1680 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1681 "Data pointer not allocated");
1682
1683 const size_t nb_gauss_pts = this->getGaussPts().size2();
1684 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1685 auto &mat = *this->dataPtr;
1686 if (type == this->zeroType) {
1687 mat.resize(msize * Tensor_Dim1, nb_gauss_pts, false);
1688 mat.clear();
1689 }
1690
1691 if (nb_gauss_pts) {
1692 const size_t nb_dofs = data.getFieldData().size();
1693
1694 if (nb_dofs) {
1695
1696 const int nb_base_functions = data.getN().size2();
1697 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1698 auto gradients_at_gauss_pts =
1699 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1703 int size = nb_dofs / msize;
1704 if (nb_dofs % msize) {
1705 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1706 "Data inconsistency");
1707 }
1708 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1709 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim0>();
1710 int bb = 0;
1711 for (; bb < size; ++bb) {
1712 gradients_at_gauss_pts(I, J, K) +=
1713 field_data(I, J) * diff_base_function(K);
1714 ++field_data;
1715 ++diff_base_function;
1716 }
1717 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1718 // functions
1719 for (; bb != nb_base_functions; ++bb)
1720 ++diff_base_function;
1721 ++gradients_at_gauss_pts;
1722 }
1723 }
1724 }
1726}
1727
1728/** \brief Get field gradients at integration pts for symmetric tensorial field
1729 * rank 2
1730 *
1731 * \ingroup mofem_forces_and_sources_user_data_operators
1732 */
1733template <int Tensor_Dim0, int Tensor_Dim1>
1736 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1737
1739 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1740 const EntityType zero_type = MBVERTEX)
1742 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1743 DoubleAllocator>(field_name, data_ptr, zero_type) {}
1744};
1745
1746template <int Tensor_Dim0, int Tensor_Dim1>
1749 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1750
1752 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1753 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1754
1755 /**
1756 * \brief calculate values of vector field at integration points
1757 * @param side side entity number
1758 * @param type side entity type
1759 * @param data entity data
1760 * @return error code
1761 */
1762 MoFEMErrorCode doWork(int side, EntityType type,
1764};
1765
1766/**
1767 * \brief Member function specialization calculating vector field gradients for
1768 * tenor field rank 2
1769 *
1770 */
1771template <int Tensor_Dim0, int Tensor_Dim1>
1773 int side, EntityType type, EntitiesFieldData::EntData &data) {
1775 if (!this->dataPtr)
1776 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1777 "Data pointer not allocated");
1778
1779 const size_t nb_gauss_pts = this->getGaussPts().size2();
1780 auto &mat = *this->dataPtr;
1781 if (type == this->zeroType) {
1782 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts, false);
1783 mat.clear();
1784 }
1785
1786 if (nb_gauss_pts) {
1787 const size_t nb_dofs = data.getFieldData().size();
1788
1789 if (nb_dofs) {
1790
1791 if (this->dataVec.use_count()) {
1792 this->dotVector.resize(nb_dofs, false);
1793 const double *array;
1794 CHKERR VecGetArrayRead(this->dataVec, &array);
1795 const auto &local_indices = data.getLocalIndices();
1796 for (int i = 0; i != local_indices.size(); ++i)
1797 if (local_indices[i] != -1)
1798 this->dotVector[i] = array[local_indices[i]];
1799 else
1800 this->dotVector[i] = 0;
1801 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1802 data.getFieldData().swap(this->dotVector);
1803 }
1804
1805 const int nb_base_functions = data.getN().size2();
1806
1807 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1808#ifndef NDEBUG
1809 if (hessian_base.size1() != nb_gauss_pts) {
1810 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1811 "Wrong number of integration pts (%d != %d)",
1812 hessian_base.size1(), nb_gauss_pts);
1813 }
1814 if (hessian_base.size2() !=
1815 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1816 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1817 "Wrong number of base functions (%d != %d)",
1818 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1819 nb_base_functions);
1820 }
1821 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1822 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1823 "Wrong number of base functions (%d < %d)",
1824 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1825 }
1826#endif
1827
1828 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1829 &*hessian_base.data().begin());
1830
1831 auto t_hessian_at_gauss_pts =
1832 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1833
1837
1838 int size = nb_dofs / Tensor_Dim0;
1839#ifndef NDEBUG
1840 if (nb_dofs % Tensor_Dim0) {
1841 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1842 "Data inconsistency");
1843 }
1844#endif
1845
1846 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1847 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1848 int bb = 0;
1849 for (; bb < size; ++bb) {
1850 t_hessian_at_gauss_pts(I, J, K) +=
1851 field_data(I) * t_diff2_base_function(J, K);
1852 ++field_data;
1853 ++t_diff2_base_function;
1854 }
1855 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1856 // functions
1857 for (; bb != nb_base_functions; ++bb)
1858 ++t_diff2_base_function;
1859 ++t_hessian_at_gauss_pts;
1860 }
1861
1862 if (this->dataVec.use_count()) {
1863 data.getFieldData().swap(this->dotVector);
1864 }
1865 }
1866 }
1868}
1869
1870/**@}*/
1871
1872/** \name Transform tensors and vectors */
1873
1874/**@{*/
1875
1876/**
1877 * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1878 * \f$
1879 *
1880 * @tparam DIM
1881 *
1882 * \ingroup mofem_forces_and_sources_user_data_operators
1883 */
1884template <int DIM_01, int DIM_23, int S = 0>
1887
1890
1891 /**
1892 * @deprecated Do not use this constriuctor
1893 */
1895 boost::shared_ptr<MatrixDouble> in_mat,
1896 boost::shared_ptr<MatrixDouble> out_mat,
1897 boost::shared_ptr<MatrixDouble> d_mat)
1898 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1899 // Only is run for vertices
1900 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1901 if (!inMat)
1902 THROW_MESSAGE("Pointer for in mat is null");
1903 if (!outMat)
1904 THROW_MESSAGE("Pointer for out mat is null");
1905 if (!dMat)
1906 THROW_MESSAGE("Pointer for tensor mat is null");
1907 }
1908
1909 OpTensorTimesSymmetricTensor(boost::shared_ptr<MatrixDouble> in_mat,
1910 boost::shared_ptr<MatrixDouble> out_mat,
1911 boost::shared_ptr<MatrixDouble> d_mat)
1912 : UserOp(NOSPACE, OPSPACE), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1913 // Only is run for vertices
1914 if (!inMat)
1915 THROW_MESSAGE("Pointer for in mat is null");
1916 if (!outMat)
1917 THROW_MESSAGE("Pointer for out mat is null");
1918 if (!dMat)
1919 THROW_MESSAGE("Pointer for tensor mat is null");
1920 }
1921
1922 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1924 const size_t nb_gauss_pts = getGaussPts().size2();
1925 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
1926 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
1927 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1928 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
1929 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1930 t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1931 ++t_in;
1932 ++t_out;
1933 }
1935 }
1936
1937private:
1942
1943 boost::shared_ptr<MatrixDouble> inMat;
1944 boost::shared_ptr<MatrixDouble> outMat;
1945 boost::shared_ptr<MatrixDouble> dMat;
1946};
1947
1948template <int DIM>
1950
1953
1954 /**
1955 * @deprecated Do not use this constructor
1956 */
1958 boost::shared_ptr<MatrixDouble> in_mat,
1959 boost::shared_ptr<MatrixDouble> out_mat)
1960 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1961 // Only is run for vertices
1962 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1963 if (!inMat)
1964 THROW_MESSAGE("Pointer not set for in matrix");
1965 if (!outMat)
1966 THROW_MESSAGE("Pointer not set for in matrix");
1967 }
1968
1969 OpSymmetrizeTensor(boost::shared_ptr<MatrixDouble> in_mat,
1970 boost::shared_ptr<MatrixDouble> out_mat)
1971 : UserOp(NOSPACE, OPSPACE), inMat(in_mat), outMat(out_mat) {
1972 // Only is run for vertices
1973 if (!inMat)
1974 THROW_MESSAGE("Pointer not set for in matrix");
1975 if (!outMat)
1976 THROW_MESSAGE("Pointer not set for in matrix");
1977 }
1978
1979 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1981 const size_t nb_gauss_pts = getGaussPts().size2();
1982 auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
1983 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
1984 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
1985 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1986 t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
1987 ++t_in;
1988 ++t_out;
1989 }
1991 }
1992
1993private:
1996 boost::shared_ptr<MatrixDouble> inMat;
1997 boost::shared_ptr<MatrixDouble> outMat;
1998};
1999
2001
2004
2005 OpScaleMatrix(const std::string field_name, const double scale,
2006 boost::shared_ptr<MatrixDouble> in_mat,
2007 boost::shared_ptr<MatrixDouble> out_mat)
2008 : UserOp(field_name, OPROW), scale(scale), inMat(in_mat),
2009 outMat(out_mat) {
2010 // Only is run for vertices
2011 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
2012 if (!inMat)
2013 THROW_MESSAGE("Pointer not set for in matrix");
2014 if (!outMat)
2015 THROW_MESSAGE("Pointer not set for in matrix");
2016 }
2017
2018 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
2020 outMat->resize(inMat->size1(), inMat->size2(), false);
2021 noalias(*outMat) = scale * (*inMat);
2023 }
2024
2025private:
2026 const double scale;
2027 boost::shared_ptr<MatrixDouble> inMat;
2028 boost::shared_ptr<MatrixDouble> outMat;
2029};
2030
2031/**@}*/
2032
2033/** \name H-div/H-curls (Vectorial bases) values at integration points */
2034
2035/**@{*/
2036
2037/** \brief Get vector field for H-div approximation
2038 * \ingroup mofem_forces_and_sources_user_data_operators
2039 */
2040template <int Base_Dim, int Field_Dim, class T, class L, class A>
2042
2043/** \brief Get vector field for H-div approximation
2044 * \ingroup mofem_forces_and_sources_user_data_operators
2045 */
2046template <int Field_Dim>
2048 ublas::row_major, DoubleAllocator>
2050
2052 boost::shared_ptr<MatrixDouble> data_ptr,
2053 const EntityType zero_type = MBEDGE,
2054 const int zero_side = 0)
2057 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2058 if (!dataPtr)
2059 THROW_MESSAGE("Pointer is not set");
2060 }
2061
2062 /**
2063 * \brief Calculate values of vector field at integration points
2064 * @param side side entity number
2065 * @param type side entity type
2066 * @param data entity data
2067 * @return error code
2068 */
2069 MoFEMErrorCode doWork(int side, EntityType type,
2071
2072private:
2073 boost::shared_ptr<MatrixDouble> dataPtr;
2075 const int zeroSide;
2076};
2077
2078template <int Field_Dim>
2080 3, Field_Dim, double, ublas::row_major,
2081 DoubleAllocator>::doWork(int side, EntityType type,
2084 const size_t nb_integration_points = this->getGaussPts().size2();
2085 if (type == zeroType && side == zeroSide) {
2086 dataPtr->resize(Field_Dim, nb_integration_points, false);
2087 dataPtr->clear();
2088 }
2089 const size_t nb_dofs = data.getFieldData().size();
2090 if (!nb_dofs)
2092 const size_t nb_base_functions = data.getN().size2() / 3;
2094 auto t_n_hdiv = data.getFTensor1N<3>();
2095 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2096 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2097 auto t_dof = data.getFTensor0FieldData();
2098 int bb = 0;
2099 for (; bb != nb_dofs; ++bb) {
2100 t_data(i) += t_n_hdiv(i) * t_dof;
2101 ++t_n_hdiv;
2102 ++t_dof;
2103 }
2104 for (; bb != nb_base_functions; ++bb)
2105 ++t_n_hdiv;
2106 ++t_data;
2107 }
2109}
2110
2111/** \brief Get vector field for H-div approximation
2112 *
2113 * \ingroup mofem_forces_and_sources_user_data_operators
2114 */
2115template <int Base_Dim, int Field_Dim = Base_Dim>
2118 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2120 Base_Dim, Field_Dim, double, ublas::row_major,
2121 DoubleAllocator>::OpCalculateHVecVectorField_General;
2122};
2123
2124
2125
2126/** \brief Get vector field for H-div approximation
2127 * \ingroup mofem_forces_and_sources_user_data_operators
2128 */
2129template <int Base_Dim, int Field_Dim = Base_Dim>
2131
2132template <int Field_Dim>
2135
2137 boost::shared_ptr<MatrixDouble> data_ptr,
2138 const EntityType zero_type = MBEDGE,
2139 const int zero_side = 0)
2142 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2143 if (!dataPtr)
2144 THROW_MESSAGE("Pointer is not set");
2145 }
2146
2147 /**
2148 * \brief Calculate values of vector field at integration points
2149 * @param side side entity number
2150 * @param type side entity type
2151 * @param data entity data
2152 * @return error code
2153 */
2154 MoFEMErrorCode doWork(int side, EntityType type,
2156
2157private:
2158 boost::shared_ptr<MatrixDouble> dataPtr;
2160 const int zeroSide;
2161};
2162
2163template <int Field_Dim>
2165 int side, EntityType type, EntitiesFieldData::EntData &data) {
2167
2168 const size_t nb_integration_points = this->getGaussPts().size2();
2169 if (type == zeroType && side == zeroSide) {
2170 dataPtr->resize(Field_Dim, nb_integration_points, false);
2171 dataPtr->clear();
2172 }
2173
2174 auto &local_indices = data.getIndices();
2175 const size_t nb_dofs = local_indices.size();
2176 if (nb_dofs) {
2177
2178 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2179 const double *array;
2180 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2181 for (size_t i = 0; i != nb_dofs; ++i)
2182 if (local_indices[i] != -1)
2183 dot_dofs_vector[i] = array[local_indices[i]];
2184 else
2185 dot_dofs_vector[i] = 0;
2186 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2187
2188 const size_t nb_base_functions = data.getN().size2() / 3;
2190 auto t_n_hdiv = data.getFTensor1N<3>();
2191 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2192 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2193 int bb = 0;
2194 for (; bb != nb_dofs; ++bb) {
2195 t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
2196 ++t_n_hdiv;
2197 }
2198 for (; bb != nb_base_functions; ++bb)
2199 ++t_n_hdiv;
2200 ++t_data;
2201 }
2202 }
2203
2205}
2206
2207/**
2208 * @brief Calculate divergence of vector field
2209 * @ingroup mofem_forces_and_sources_user_data_operators
2210 *
2211 * @tparam BASE_DIM
2212 * @tparam SPACE_DIM
2213 */
2214template <int BASE_DIM, int SPACE_DIM>
2217
2219 boost::shared_ptr<VectorDouble> data_ptr,
2220 const EntityType zero_type = MBEDGE,
2221 const int zero_side = 0)
2224 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2225 if (!dataPtr)
2226 THROW_MESSAGE("Pointer is not set");
2227 }
2228
2232 const size_t nb_integration_points = getGaussPts().size2();
2233 if (type == zeroType && side == zeroSide) {
2234 dataPtr->resize(nb_integration_points, false);
2235 dataPtr->clear();
2236 }
2237 const size_t nb_dofs = data.getFieldData().size();
2238 if (!nb_dofs)
2240 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2243 auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2244 auto t_data = getFTensor0FromVec(*dataPtr);
2245 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2246 auto t_dof = data.getFTensor0FieldData();
2247 int bb = 0;
2248 for (; bb != nb_dofs; ++bb) {
2249 t_data += t_dof * t_n_diff_hdiv(j, j);
2250 ++t_n_diff_hdiv;
2251 ++t_dof;
2252 }
2253 for (; bb != nb_base_functions; ++bb)
2254 ++t_n_diff_hdiv;
2255 ++t_data;
2256 }
2258 }
2259
2260private:
2261 boost::shared_ptr<VectorDouble> dataPtr;
2263 const int zeroSide;
2264};
2265
2266/**
2267 * @brief Calculate gradient of vector field
2268 * @ingroup mofem_forces_and_sources_user_data_operators
2269 *
2270 * @tparam BASE_DIM
2271 * @tparam SPACE_DIM
2272 */
2273template <int BASE_DIM, int SPACE_DIM>
2276
2278 boost::shared_ptr<MatrixDouble> data_ptr,
2279 const EntityType zero_type = MBEDGE,
2280 const int zero_side = 0)
2283 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2284 if (!dataPtr)
2285 THROW_MESSAGE("Pointer is not set");
2286 }
2287
2291 const size_t nb_integration_points = getGaussPts().size2();
2292 if (type == zeroType && side == zeroSide) {
2293 dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
2294 dataPtr->clear();
2295 }
2296 const size_t nb_dofs = data.getFieldData().size();
2297 if (!nb_dofs)
2299 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2302 auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2303 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*dataPtr);
2304 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2305 auto t_dof = data.getFTensor0FieldData();
2306 int bb = 0;
2307 for (; bb != nb_dofs; ++bb) {
2308 t_data(i, j) += t_dof * t_base_diff(i, j);
2309 ++t_base_diff;
2310 ++t_dof;
2311 }
2312 for (; bb != nb_base_functions; ++bb)
2313 ++t_base_diff;
2314 ++t_data;
2315 }
2317 }
2318
2319private:
2320 boost::shared_ptr<MatrixDouble> dataPtr;
2322 const int zeroSide;
2323};
2324
2325/**
2326 * @brief Calculate gradient of vector field
2327 * @ingroup mofem_forces_and_sources_user_data_operators
2328 *
2329 * @tparam BASE_DIM
2330 * @tparam SPACE_DIM
2331 */
2332template <int BASE_DIM, int SPACE_DIM>
2335
2337 boost::shared_ptr<MatrixDouble> data_ptr,
2338 const EntityType zero_type = MBEDGE,
2339 const int zero_side = 0)
2342 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2343 if (!dataPtr)
2344 THROW_MESSAGE("Pointer is not set");
2345 }
2346
2350 const size_t nb_integration_points = getGaussPts().size2();
2351 if (type == zeroType && side == zeroSide) {
2352 dataPtr->resize(BASE_DIM * SPACE_DIM * SPACE_DIM, nb_integration_points,
2353 false);
2354 dataPtr->clear();
2355 }
2356 const size_t nb_dofs = data.getFieldData().size();
2357 if (!nb_dofs)
2359
2360 const int nb_base_functions = data.getN().size2() / BASE_DIM;
2361
2362#ifndef NDEBUG
2363 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2364 if (hessian_base.size1() != nb_integration_points) {
2365 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2366 "Wrong number of integration pts (%d != %d)",
2367 hessian_base.size1(), nb_integration_points);
2368 }
2369 if (hessian_base.size2() !=
2370 BASE_DIM * nb_base_functions * SPACE_DIM * SPACE_DIM) {
2371 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2372 "Wrong number of base functions (%d != %d)",
2373 hessian_base.size2() / (BASE_DIM * SPACE_DIM * SPACE_DIM),
2374 nb_base_functions);
2375 }
2376 if (hessian_base.size2() < BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM) {
2377 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2378 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2379 BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM);
2380 }
2381#endif
2382
2386
2387 auto t_base_diff2 =
2389 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*dataPtr);
2390 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2391 auto t_dof = data.getFTensor0FieldData();
2392 int bb = 0;
2393 for (; bb != nb_dofs; ++bb) {
2394 t_data(i, j, k) += t_dof * t_base_diff2(i, j, k);
2395
2396 ++t_base_diff2;
2397 ++t_dof;
2398 }
2399 for (; bb != nb_base_functions; ++bb)
2400 ++t_base_diff2;
2401 ++t_data;
2402 }
2404 }
2405
2406private:
2407 boost::shared_ptr<MatrixDouble> dataPtr;
2409 const int zeroSide;
2410};
2411
2412/**
2413 * @brief Calculate divergence of vector field dot
2414 * @ingroup mofem_forces_and_sources_user_data_operators
2415 *
2416 * @tparam Tensor_Dim dimension of space
2417 */
2418template <int Tensor_Dim1, int Tensor_Dim2>
2421
2423 boost::shared_ptr<VectorDouble> data_ptr,
2424 const EntityType zero_type = MBEDGE,
2425 const int zero_side = 0)
2428 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2429 if (!dataPtr)
2430 THROW_MESSAGE("Pointer is not set");
2431 }
2432
2436 const size_t nb_integration_points = getGaussPts().size2();
2437 if (type == zeroType && side == zeroSide) {
2438 dataPtr->resize(nb_integration_points, false);
2439 dataPtr->clear();
2440 }
2441
2442 const auto &local_indices = data.getLocalIndices();
2443 const int nb_dofs = local_indices.size();
2444 if (nb_dofs) {
2445
2446 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2447 const double *array;
2448 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2449 for (size_t i = 0; i != local_indices.size(); ++i)
2450 if (local_indices[i] != -1)
2451 dot_dofs_vector[i] = array[local_indices[i]];
2452 else
2453 dot_dofs_vector[i] = 0;
2454 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2455
2456 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
2458 auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
2459 auto t_data = getFTensor0FromVec(*dataPtr);
2460 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2461 int bb = 0;
2462 for (; bb != nb_dofs; ++bb) {
2463 double div = 0;
2464 for (auto ii = 0; ii != Tensor_Dim2; ++ii)
2465 div += t_n_diff_hdiv(ii, ii);
2466 t_data += dot_dofs_vector[bb] * div;
2467 ++t_n_diff_hdiv;
2468 }
2469 for (; bb != nb_base_functions; ++bb)
2470 ++t_n_diff_hdiv;
2471 ++t_data;
2472 }
2473 }
2475 }
2476
2477private:
2478 boost::shared_ptr<VectorDouble> dataPtr;
2480 const int zeroSide;
2481};
2482
2483/**
2484 * @brief Calculate curl of vector field
2485 * @ingroup mofem_forces_and_sources_user_data_operators
2486 *
2487 * @tparam Base_Dim base function dimension
2488 * @tparam Space_Dim dimension of space
2489 * @tparam Hcurl field dimension
2490 */
2491template <int Base_Dim, int Space_Dim>
2493
2494/**
2495 * @brief Calculate curl of vector field
2496 * @ingroup mofem_forces_and_sources_user_data_operators
2497 *
2498 * @tparam Base_Dim base function dimension
2499 * @tparam Space_Dim dimension of space
2500 * @tparam Hcurl field dimension
2501 */
2502template <>
2505 OpCalculateHcurlVectorCurl(const std::string field_name,
2506 boost::shared_ptr<MatrixDouble> data_ptr,
2507 const EntityType zero_type = MBEDGE,
2508 const int zero_side = 0);
2509 MoFEMErrorCode doWork(int side, EntityType type,
2511
2512private:
2513 boost::shared_ptr<MatrixDouble> dataPtr;
2515 const int zeroSide;
2516};
2517
2518/**
2519 * @brief Calculate curl of vector field
2520 * @ingroup mofem_forces_and_sources_user_data_operators
2521 *
2522 * @tparam Field_Dim dimension of field
2523 * @tparam Space_Dim dimension of space
2524 */
2525template <>
2528
2529 OpCalculateHcurlVectorCurl(const std::string field_name,
2530 boost::shared_ptr<MatrixDouble> data_ptr,
2531 const EntityType zero_type = MBVERTEX,
2532 const int zero_side = 0);
2533
2534 MoFEMErrorCode doWork(int side, EntityType type,
2536
2537private:
2538 boost::shared_ptr<MatrixDouble> dataPtr;
2540 const int zeroSide;
2541};
2542
2543/**
2544 * @brief Calculate curl of vector field
2545 * @ingroup mofem_forces_and_sources_user_data_operators
2546 *
2547 * @tparam Field_Dim dimension of field
2548 * @tparam Space_Dim dimension of space
2549 */
2550template <>
2553
2554 OpCalculateHcurlVectorCurl(const std::string field_name,
2555 boost::shared_ptr<MatrixDouble> data_ptr,
2556 const EntityType zero_type = MBVERTEX,
2557 const int zero_side = 0);
2558
2559 MoFEMErrorCode doWork(int side, EntityType type,
2561
2562private:
2563 boost::shared_ptr<MatrixDouble> dataPtr;
2565 const int zeroSide;
2566};
2567
2568/**
2569 * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
2570 * \ingroup mofem_forces_and_sources_user_data_operators
2571 *
2572 * @tparam Tensor_Dim0 rank of the filed
2573 * @tparam Tensor_Dim1 dimension of space
2574 */
2575template <int Tensor_Dim0, int Tensor_Dim1>
2578
2580 boost::shared_ptr<MatrixDouble> data_ptr,
2581 const EntityType zero_type = MBEDGE,
2582 const int zero_side = 0)
2585 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2586 if (!dataPtr)
2587 THROW_MESSAGE("Pointer is not set");
2588 }
2589
2593 const size_t nb_integration_points = getGaussPts().size2();
2594 if (type == zeroType && side == zeroSide) {
2595 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2596 dataPtr->clear();
2597 }
2598 const size_t nb_dofs = data.getFieldData().size();
2599 if (nb_dofs) {
2600 const size_t nb_base_functions = data.getN().size2() / 3;
2603 auto t_n_hvec = data.getFTensor1N<3>();
2604 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2605 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2606 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2607 size_t bb = 0;
2608 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2609 t_data(i, j) += t_dof(i) * t_n_hvec(j);
2610 ++t_n_hvec;
2611 ++t_dof;
2612 }
2613 for (; bb < nb_base_functions; ++bb)
2614 ++t_n_hvec;
2615 ++t_data;
2616 }
2617 }
2619 }
2620
2621private:
2622 boost::shared_ptr<MatrixDouble> dataPtr;
2624 const int zeroSide;
2625};
2626
2627/**
2628 * @brief Calculate tenor field using tensor base, i.e. Hdiv/Hcurl
2629 * \ingroup mofem_forces_and_sources_user_data_operators
2630 *
2631 * @tparam Tensor_Dim0 rank of the filed
2632 * @tparam Tensor_Dim1 dimension of space
2633 */
2634template <int Tensor_Dim0, int Tensor_Dim1>
2637
2639 boost::shared_ptr<MatrixDouble> data_ptr,
2640 const EntityType zero_type = MBEDGE,
2641 const int zero_side = 0)
2644 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2645 if (!dataPtr)
2646 THROW_MESSAGE("Pointer is not set");
2647 }
2648
2652 const size_t nb_integration_points = getGaussPts().size2();
2653 if (type == zeroType && side == zeroSide) {
2654 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2655 dataPtr->clear();
2656 }
2657 const size_t nb_dofs = data.getFieldData().size();
2658 if (!nb_dofs)
2660 const size_t nb_base_functions =
2661 data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2664 auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2665 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2666 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2667 auto t_dof = data.getFTensor0FieldData();
2668 size_t bb = 0;
2669 for (; bb != nb_dofs; ++bb) {
2670 t_data(i, j) += t_dof * t_n_hten(i, j);
2671 ++t_n_hten;
2672 ++t_dof;
2673 }
2674 for (; bb < nb_base_functions; ++bb)
2675 ++t_n_hten;
2676 ++t_data;
2677 }
2679 }
2680
2681private:
2682 boost::shared_ptr<MatrixDouble> dataPtr;
2684 const int zeroSide;
2685};
2686
2687/**
2688 * @brief Calculate divergence of tonsorial field using vectorial base
2689 * \ingroup mofem_forces_and_sources_user_data_operators
2690 *
2691 * @tparam Tensor_Dim0 rank of the field
2692 * @tparam Tensor_Dim1 dimension of space
2693 */
2694template <int Tensor_Dim0, int Tensor_Dim1,
2695 CoordinateTypes CoordSys = CARTESIAN>
2698
2700 boost::shared_ptr<MatrixDouble> data_ptr,
2701 const EntityType zero_type = MBEDGE,
2702 const int zero_side = 0)
2705 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2706 if (!dataPtr)
2707 THROW_MESSAGE("Pointer is not set");
2708 }
2709
2713 const size_t nb_integration_points = getGaussPts().size2();
2714 if (type == zeroType && side == 0) {
2715 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
2716 dataPtr->clear();
2717 }
2718 const size_t nb_dofs = data.getFieldData().size();
2719 if (nb_dofs) {
2720 const size_t nb_base_functions = data.getN().size2() / 3;
2723 auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
2724 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
2725 auto t_base = data.getFTensor1N<3>();
2726 auto t_coords = getFTensor1CoordsAtGaussPts();
2727 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2728 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2729 size_t bb = 0;
2730 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2731 double div = t_n_diff_hvec(j, j);
2732 t_data(i) += t_dof(i) * div;
2733 if constexpr (CoordSys == CYLINDRICAL) {
2734 t_data(i) += t_dof(i) * (t_base(0) / t_coords(0));
2735 }
2736 ++t_n_diff_hvec;
2737 ++t_dof;
2738 ++t_base;
2739 }
2740 for (; bb < nb_base_functions; ++bb) {
2741 ++t_base;
2742 ++t_n_diff_hvec;
2743 }
2744 ++t_data;
2745 ++t_coords;
2746 }
2747 }
2749 }
2750
2751private:
2752 boost::shared_ptr<MatrixDouble> dataPtr;
2754 const int zeroSide;
2755};
2756
2757/**
2758 * @brief Calculate trace of vector (Hdiv/Hcurl) space
2759 *
2760 * @tparam Tensor_Dim
2761 * @tparam OpBase
2762 */
2763template <int Tensor_Dim, typename OpBase>
2765
2767 boost::shared_ptr<MatrixDouble> data_ptr,
2768 const EntityType zero_type = MBEDGE,
2769 const int zero_side = 0)
2770 : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
2771 zeroType(zero_type), zeroSide(zero_side) {
2772 if (!dataPtr)
2773 THROW_MESSAGE("Pointer is not set");
2774 }
2775
2779 const size_t nb_integration_points = OpBase::getGaussPts().size2();
2780 if (type == zeroType && side == 0) {
2781 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
2782 dataPtr->clear();
2783 }
2784 const size_t nb_dofs = data.getFieldData().size();
2785 if (nb_dofs) {
2786 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
2787 const size_t nb_base_functions = data.getN().size2() / 3;
2788 auto t_base = data.getFTensor1N<3>();
2789 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2790 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2791 FTensor::Tensor1<double, Tensor_Dim> t_normalized_normal;
2792 t_normalized_normal(j) = t_normal(j);
2793 t_normalized_normal.normalize();
2794 auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
2795 size_t bb = 0;
2796 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2797 t_data(i) += t_dof(i) * (t_base(j) * t_normalized_normal(j));
2798 ++t_base;
2799 ++t_dof;
2800 }
2801 for (; bb < nb_base_functions; ++bb) {
2802 ++t_base;
2803 }
2804 ++t_data;
2805 ++t_normal;
2806 }
2807 }
2809 }
2810
2811private:
2812 boost::shared_ptr<MatrixDouble> dataPtr;
2814 const int zeroSide;
2817};
2818
2819/**@}*/
2820
2821/** \name Other operators */
2822
2823/**@{*/
2824
2825/**@}*/
2826
2827/** \name Operators for faces */
2828
2829/**@{*/
2830
2831/** \brief Transform local reference derivatives of shape functions to global
2832derivatives
2833
2834\ingroup mofem_forces_and_sources_tri_element
2835
2836*/
2837template <int DIM, int DERIVATIVE = 1> struct OpSetInvJacSpaceForFaceImpl;
2838
2841
2843 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
2844
2845protected:
2846 template <int D1, int D2, int J1, int J2>
2849
2850 static_assert(D2 == J2, "Dimension of jacobian and dimension of <out> "
2851 "directive does not match");
2852
2853 size_t nb_functions = diff_n.size2() / D1;
2854 if (nb_functions) {
2855 size_t nb_gauss_pts = diff_n.size1();
2856
2857#ifndef NDEBUG
2858 if (nb_gauss_pts != getGaussPts().size2())
2859 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2860 "Wrong number of Gauss Pts");
2861 if (diff_n.size2() % D1)
2862 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2863 "Number of direvatives of base functions and D1 dimension does "
2864 "not match");
2865#endif
2866
2867 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions, false);
2868
2871 auto t_diff_n = getFTensor1FromPtr<D2>(&*diffNinvJac.data().begin());
2872 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
2873 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*invJacPtr);
2874 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
2875 for (size_t dd = 0; dd != nb_functions; ++dd) {
2876 t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
2877 ++t_diff_n;
2878 ++t_diff_n_ref;
2879 }
2880 }
2881
2882 diff_n.swap(diffNinvJac);
2883 }
2885 }
2886
2887 boost::shared_ptr<MatrixDouble> invJacPtr;
2889};
2890
2891template <>
2894
2896
2897 MoFEMErrorCode doWork(int side, EntityType type,
2899};
2900
2901template <>
2903 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2904
2905 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2906
2907 MoFEMErrorCode doWork(int side, EntityType type,
2909};
2910
2911template <>
2913 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
2914
2915 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
2916
2917 MoFEMErrorCode doWork(int side, EntityType type,
2919};
2920
2921template <int DERIVARIVE = 1>
2923 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2924 OpSetInvJacH1ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2925 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(H1, inv_jac_ptr) {}
2926};
2927
2928template <int DERIVARIVE = 1>
2930 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
2931 OpSetInvJacL2ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2932 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(L2, inv_jac_ptr) {}
2933};
2934
2935template <int DERIVARIVE = 1>
2937 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2939 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2940 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(H1, inv_jac_ptr) {}
2941};
2942
2943template <int DERIVARIVE = 1>
2945 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
2947 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2948 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(L2, inv_jac_ptr) {}
2949};
2950
2951/**
2952 * \brief Transform local reference derivatives of shape function to
2953 global derivatives for face
2954
2955 * \ingroup mofem_forces_and_sources_tri_element
2956 */
2957template <int DIM> struct OpSetInvJacHcurlFaceImpl;
2958
2959template <>
2962
2963 OpSetInvJacHcurlFaceImpl(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2965 invJacPtr(inv_jac_ptr) {}
2966
2967 MoFEMErrorCode doWork(int side, EntityType type,
2969
2970protected:
2971 boost::shared_ptr<MatrixDouble> invJacPtr;
2973};
2974
2975template <>
2977 using OpSetInvJacHcurlFaceImpl<2>::OpSetInvJacHcurlFaceImpl;
2978 MoFEMErrorCode doWork(int side, EntityType type,
2980};
2981
2984
2985/**
2986 * @brief Make Hdiv space from Hcurl space in 2d
2987 * @ingroup mofem_forces_and_sources_tri_element
2988 */
2991
2994
2995 MoFEMErrorCode doWork(int side, EntityType type,
2997};
2998
2999/** \brief Transform Hcurl base fluxes from reference element to physical
3000 * triangle
3001 */
3003
3004/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3005 *
3006 * Covariant Piola transformation
3007 \f[
3008 \psi_i|_t = J^{-1}_{ij}\hat{\psi}_j\\
3009 \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3010 = J^{-1}_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3011 \f]
3012
3013 */
3014template <>
3017
3019 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3020
3021 MoFEMErrorCode doWork(int side, EntityType type,
3023
3024private:
3025 boost::shared_ptr<MatrixDouble> invJacPtr;
3026
3029};
3030
3033
3034/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3035 *
3036 * \note Hdiv space is generated by Hcurl space in 2d.
3037 *
3038 * Contravariant Piola transformation
3039 * \f[
3040 * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
3041 * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3042 * =
3043 * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3044 * \f]
3045 *
3046 * \ingroup mofem_forces_and_sources
3047 *
3048 */
3050
3051template <>
3054
3056 boost::shared_ptr<MatrixDouble> jac_ptr)
3058 jacPtr(jac_ptr) {}
3059
3060 MoFEMErrorCode doWork(int side, EntityType type,
3062
3063protected:
3064 boost::shared_ptr<MatrixDouble> jacPtr;
3067};
3068
3069template <>
3073 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3074
3075 MoFEMErrorCode doWork(int side, EntityType type,
3077};
3078
3083
3084/**@}*/
3085
3086/** \name Operators for edges */
3087
3088/**@{*/
3089
3092
3095
3096 MoFEMErrorCode doWork(int side, EntityType type,
3098
3099};
3100
3101/**
3102 * @deprecated Name is deprecated and this is added for back compatibility
3103 */
3106
3107/**@}*/
3108
3109/** \name Operator for fat prisms */
3110
3111/**@{*/
3112
3113/**
3114 * @brief Operator for fat prism element updating integration weights in the
3115 * volume.
3116 *
3117 * Jacobian on the distorted element is nonconstant. This operator updates
3118 * integration weight on prism to take into account nonconstat jacobian.
3119 *
3120 * \f[
3121 * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
3122 * \pmb\xi} \right\| \right)
3123 * \f]
3124 * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
3125 * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
3126 * element coordinate.
3127 *
3128 */
3131
3134
3135 MoFEMErrorCode doWork(int side, EntityType type,
3137};
3138
3139/** \brief Calculate inverse of jacobian for face element
3140
3141 It is assumed that face element is XY plane. Applied
3142 only for 2d problems.
3143
3144 FIXME Generalize function for arbitrary face orientation in 3d space
3145 FIXME Calculate to Jacobins for two faces
3146
3147 \ingroup mofem_forces_and_sources_prism_element
3148
3149*/
3152
3153 OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3155 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3156
3159 invJac(inv_jac) {}
3160 MoFEMErrorCode doWork(int side, EntityType type,
3162
3163private:
3164 const boost::shared_ptr<MatrixDouble> invJacPtr;
3166};
3167
3168/** \brief Transform local reference derivatives of shape functions to global
3169derivatives
3170
3171FIXME Generalize to curved shapes
3172FIXME Generalize to case that top and bottom face has different shape
3173
3174\ingroup mofem_forces_and_sources_prism_element
3175
3176*/
3179
3180 OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3182 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3183
3186 invJac(inv_jac) {}
3187
3188 MoFEMErrorCode doWork(int side, EntityType type,
3190
3191private:
3192 const boost::shared_ptr<MatrixDouble> invJacPtr;
3195};
3196
3197// Flat prism
3198
3199/** \brief Calculate inverse of jacobian for face element
3200
3201 It is assumed that face element is XY plane. Applied
3202 only for 2d problems.
3203
3204 FIXME Generalize function for arbitrary face orientation in 3d space
3205 FIXME Calculate to Jacobins for two faces
3206
3207 \ingroup mofem_forces_and_sources_prism_element
3208
3209*/
3212
3215 invJacF3(inv_jac_f3) {}
3216 MoFEMErrorCode doWork(int side, EntityType type,
3218
3219private:
3221};
3222
3223/** \brief Transform local reference derivatives of shape functions to global
3224derivatives
3225
3226FIXME Generalize to curved shapes
3227FIXME Generalize to case that top and bottom face has different shape
3228
3229\ingroup mofem_forces_and_sources_prism_element
3230
3231*/
3234
3237 invJacF3(inv_jac_f3) {}
3238
3239 MoFEMErrorCode doWork(int side, EntityType type,
3241
3242private:
3245};
3246
3247/**@}*/
3248
3249/** \name Operation on matrices at integration points */
3250
3251/**@{*/
3252
3253template <int DIM>
3255
3256 OpInvertMatrix(boost::shared_ptr<MatrixDouble> in_ptr,
3257 boost::shared_ptr<VectorDouble> det_ptr,
3258 boost::shared_ptr<MatrixDouble> out_ptr)
3260 outPtr(out_ptr), detPtr(det_ptr) {}
3261
3262 MoFEMErrorCode doWork(int side, EntityType type,
3264
3265private:
3266 boost::shared_ptr<MatrixDouble> inPtr;
3267 boost::shared_ptr<MatrixDouble> outPtr;
3268 boost::shared_ptr<VectorDouble> detPtr;
3269
3270};
3271
3272template <int DIM>
3276
3277 if (!inPtr)
3278 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3279 "Pointer for inPtr matrix not allocated");
3280 if (!detPtr)
3281 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3282 "Pointer for detPtr matrix not allocated");
3283
3284 const auto nb_rows = inPtr->size1();
3285 const auto nb_integration_pts = inPtr->size2();
3286
3287 // Calculate determinant
3288 {
3289 detPtr->resize(nb_integration_pts, false);
3290 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3291 auto t_det = getFTensor0FromVec(*detPtr);
3292 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3293 t_det = determinantTensor(t_in);
3294 ++t_in;
3295 ++t_det;
3296 }
3297 }
3298
3299 // Invert jacobian
3300 if (outPtr) {
3301 outPtr->resize(nb_rows, nb_integration_pts, false);
3302 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3303 auto t_out = getFTensor2FromMat<DIM, DIM>(*outPtr);
3304 auto t_det = getFTensor0FromVec(*detPtr);
3305 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3306 CHKERR invertTensor(t_in, t_det, t_out);
3307 ++t_in;
3308 ++t_out;
3309 ++t_det;
3310 }
3311 }
3312
3314}
3315
3316/**@}*/
3317
3318/** \brief Calculates the trace of an input matrix
3319
3320\ingroup mofem_forces_and_sources
3321
3322*/
3323
3324template <int DIM>
3326
3327 OpCalculateTraceFromMat(boost::shared_ptr<MatrixDouble> in_ptr,
3328 boost::shared_ptr<VectorDouble> out_ptr)
3330 outPtr(out_ptr) {}
3331
3332 MoFEMErrorCode doWork(int side, EntityType type,
3333 DataForcesAndSourcesCore::EntData &data);
3334
3335private:
3337 boost::shared_ptr<MatrixDouble> inPtr;
3338 boost::shared_ptr<VectorDouble> outPtr;
3339
3340};
3341
3342template <int DIM>
3345 DataForcesAndSourcesCore::EntData &data) {
3347
3348 if (!inPtr)
3349 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3350 "Pointer for inPtr matrix not allocated");
3351
3352 const auto nb_integration_pts = inPtr->size2();
3353 // Invert jacobian
3354 if (outPtr) {
3355 outPtr->resize(nb_integration_pts, false);
3356 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3357 auto t_out = getFTensor0FromVec(*outPtr);
3358
3359 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3360 t_out = t_in(i, i);
3361 ++t_in;
3362 ++t_out;
3363 }
3364 }
3365
3367}
3368
3369
3370/**@}*/
3371
3372/** \brief Calculates the trace of an input matrix
3373
3374\ingroup mofem_forces_and_sources
3375
3376*/
3377
3378template <int DIM>
3380
3381 OpCalculateTraceFromSymmMat(boost::shared_ptr<MatrixDouble> in_ptr,
3382 boost::shared_ptr<VectorDouble> out_ptr)
3384 outPtr(out_ptr) {}
3385
3386 MoFEMErrorCode doWork(int side, EntityType type,
3387 DataForcesAndSourcesCore::EntData &data);
3388
3389private:
3391 boost::shared_ptr<MatrixDouble> inPtr;
3392 boost::shared_ptr<VectorDouble> outPtr;
3393
3394};
3395
3396template <int DIM>
3399 DataForcesAndSourcesCore::EntData &data) {
3401
3402 if (!inPtr)
3403 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3404 "Pointer for inPtr matrix not allocated");
3405
3406 const auto nb_integration_pts = inPtr->size2();
3407 // Invert jacobian
3408 if (outPtr) {
3409 outPtr->resize(nb_integration_pts, false);
3410 auto t_in = getFTensor2SymmetricFromMat<DIM>(*inPtr);
3411 auto t_out = getFTensor0FromVec(*outPtr);
3412
3413 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3414 t_out = t_in(i, i);
3415 ++t_in;
3416 ++t_out;
3417 }
3418 }
3419
3421}
3422
3423
3424} // namespace MoFEM
3425
3426#endif // __USER_DATA_OPERATORS_HPP__
3427
3428/**
3429 * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
3430 *
3431 * \brief Classes and functions used to evaluate fields at integration pts,
3432 *jacobians, etc..
3433 *
3434 * \ingroup mofem_forces_and_sources
3435 **/
static Index< 'J', 3 > J
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
Tensor1< T, Tensor_Dim > normalize()
static const char *const CoordinateTypesNames[]
Coordinate system names.
Definition: definitions.h:108
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
CoordinateTypes
Coodinate system.
Definition: definitions.h:114
@ CYLINDRICAL
Definition: definitions.h:117
@ POLAR
Definition: definitions.h:116
@ CARTESIAN
Definition: definitions.h:115
@ SPHERICAL
Definition: definitions.h:118
#define DEPRECATED
Definition: definitions.h:17
#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
constexpr int BASE_DIM
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
Definition: Templates.hpp:1742
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
Definition: Templates.hpp:1542
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
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
@ OPSPACE
operator do Work is execute on space data
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
boost::shared_ptr< VectorDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateDivergenceVectorFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate trace of vector (Hdiv/Hcurl) space.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Index< 'j', Tensor_Dim > j
FTensor::Index< 'i', Tensor_Dim > i
OpCalculateHVecVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
OpCalculateHVecVectorFieldDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
Get vector field for H-div approximation.
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecVectorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecVectorHessian(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
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.
Calculates the trace of an input matrix.
boost::shared_ptr< MatrixDouble > inPtr
FTensor::Index< 'i', DIM > i
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< VectorDouble > outPtr
OpCalculateTraceFromMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Calculates the trace of an input matrix.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< VectorDouble > outPtr
boost::shared_ptr< MatrixDouble > inPtr
OpCalculateTraceFromSymmMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
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 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
OpSymmetrizeTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'j', DIM > j
DEPRECATED OpSymmetrizeTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'i', DIM_01 > i
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
OpTensorTimesSymmetricTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
FTensor::Index< 'j', DIM_01 > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'l', DIM_23 > l
DEPRECATED OpTensorTimesSymmetricTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
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