v0.15.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 */
53 MoFEMErrorCode doWork(int side, EntityType type,
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 SETERRQ(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
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 */
95 MoFEMErrorCode doWork(int side, EntityType type,
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
169 MoFEMErrorCode doWork(int side, EntityType type,
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 field 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 SETERRQ(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);
413 FTensor::Index<'I', Tensor_Dim> I;
414 const size_t size = nb_dofs / Tensor_Dim;
415 if (nb_dofs % Tensor_Dim) {
416 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
417 "Nb. of DOFs is inconsistent with Tensor_Dim");
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 SINGULARITY
434 #ifndef NDEBUG
435 if (base_function != base_function) {
436 MOFEM_LOG("SELF", Sev::error) << "base function: " << base_function;
437 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
438 "Wrong number number in base functions");
439 }
440 #endif
441 #endif
442
443 values_at_gauss_pts(I) += field_data(I) * base_function;
444 ++field_data;
445 ++base_function;
446 }
447 // Number of dofs can be smaller than number of Tensor_Dim x base
448 // functions
449 for (; bb < nb_base_functions; ++bb)
450 ++base_function;
451 ++values_at_gauss_pts;
452 }
453 }
454
455 if (dataVec.use_count()) {
456 data.getFieldData().swap(dotVector);
457 }
458 }
460}
461
462/** \brief Get values at integration pts for tensor field rank 1, i.e. vector
463 * field
464 *
465 * \ingroup mofem_forces_and_sources_user_data_operators
466 */
467template <int Tensor_Dim>
470 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
471
473 Tensor_Dim, double, ublas::row_major,
475};
476
477/**@}*/
478
479/** \name Vector field values at integration points */
480
481/**@{*/
482
483/** \brief Calculate field values (template specialization) for tensor field
484 * rank 1, i.e. vector field
485 *
486 */
487template <int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
490
492 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
493 const EntityType zero_type = MBVERTEX)
496 dataPtr(data_ptr), zeroType(zero_type) {
497 if (!dataPtr)
498 THROW_MESSAGE("Pointer is not set");
499 }
500
501 MoFEMErrorCode doWork(int side, EntityType type,
504
505 // When we move to C++17 add if constexpr()
506 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
507 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
508 "%s coordiante not implemented",
509 CoordinateTypesNames[COORDINATE_SYSTEM]);
510
511 const size_t nb_gauss_pts = getGaussPts().size2();
512 auto &vec = *dataPtr;
513 if (type == zeroType) {
514 vec.resize(nb_gauss_pts, false);
515 vec.clear();
516 }
517
518 const size_t nb_dofs = data.getFieldData().size();
519 if (nb_dofs) {
520
521 if (nb_gauss_pts) {
522 const size_t nb_base_functions = data.getN().size2();
523 auto values_at_gauss_pts = getFTensor0FromVec(vec);
524 FTensor::Index<'I', Tensor_Dim> I;
525 const size_t size = nb_dofs / Tensor_Dim;
526 #ifndef NDEBUG
527 if (nb_dofs % Tensor_Dim) {
528 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
529 "Number of dofs should multiple of dimensions");
530 }
531 #endif
532
533 // When we move to C++17 add if constexpr()
534 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
535 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
536 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
537 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
538 size_t bb = 0;
539 for (; bb != size; ++bb) {
540 values_at_gauss_pts += field_data(I) * diff_base_function(I);
541 ++field_data;
542 ++diff_base_function;
543 }
544 // Number of dofs can be smaller than number of Tensor_Dim x base
545 // functions
546 for (; bb < nb_base_functions; ++bb)
547 ++diff_base_function;
548 ++values_at_gauss_pts;
549 }
550 }
551
552 // When we move to C++17 add if constexpr()
553 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
554 auto t_coords = getFTensor1CoordsAtGaussPts();
555 auto values_at_gauss_pts = getFTensor0FromVec(vec);
556 auto base_function = data.getFTensor0N();
557 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
558 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
559 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
560 size_t bb = 0;
561 for (; bb != size; ++bb) {
562 values_at_gauss_pts += field_data(I) * diff_base_function(I);
563 values_at_gauss_pts +=
564 base_function * (field_data(0) / t_coords(0));
565 ++field_data;
566 ++base_function;
567 ++diff_base_function;
568 }
569 // Number of dofs can be smaller than number of Tensor_Dim x base
570 // functions
571 for (; bb < nb_base_functions; ++bb) {
572 ++base_function;
573 ++diff_base_function;
574 }
575 ++values_at_gauss_pts;
576 ++t_coords;
577 }
578 }
579 }
580 }
582 }
583
584protected:
585 boost::shared_ptr<VectorDouble> dataPtr;
587};
588
589/** \brief Approximate field values for given petsc vector
590 *
591 * \note Look at PetscData to see what vectors could be extracted with that user
592 * data operator.
593 *
594 * \ingroup mofem_forces_and_sources_user_data_operators
595 */
596template <int Tensor_Dim, PetscData::DataContext CTX>
599
601 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
602 const EntityType zero_at_type = MBVERTEX, bool throw_error = true)
605 dataPtr(data_ptr), zeroAtType(zero_at_type), throwError(throw_error) {
606 if (!dataPtr)
607 THROW_MESSAGE("Pointer is not set");
608 }
609
610 MoFEMErrorCode doWork(int side, EntityType type,
613
614 auto &local_indices = data.getLocalIndices();
615 const size_t nb_dofs = local_indices.size();
616 const size_t nb_gauss_pts = getGaussPts().size2();
617
618 MatrixDouble &mat = *dataPtr;
619 if (type == zeroAtType) {
620 mat.resize(Tensor_Dim, nb_gauss_pts, false);
621 mat.clear();
622 }
623 if (!nb_dofs)
625
626 if (!throwError) {
627 if ((getFEMethod()->data_ctx & PetscData::Switches(CTX)).none()) {
629 }
630 }
631
632 const double *array;
633
634 auto get_array = [&](const auto ctx, auto vec) {
636 #ifndef NDEBUG
637 if ((getFEMethod()->data_ctx & ctx).none()) {
638 MOFEM_LOG_CHANNEL("SELF");
639 MOFEM_LOG("SELF", Sev::error)
640 << "In this case field degrees of freedom are read from vector. "
641 "That usually happens when time solver is used, and access to "
642 "first or second rates is needed. You probably not set ts_u, "
643 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
644 "data_ctx to CTX_SET_X, CTX_SET_DX, CTX_SET_X_T, or "
645 "CTX_SET_X_TT respectively";
646 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
647 }
648 #endif
649 CHKERR VecGetArrayRead(vec, &array);
651 };
652
653 auto restore_array = [&](auto vec) {
654 return VecRestoreArrayRead(vec, &array);
655 };
656
657 switch (CTX) {
659 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
660 break;
662 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->dx);
663 break;
665 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
666 break;
668 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
669 break;
670 default:
671 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
672 "That case is not implemented");
673 }
674
675 dotVector.resize(local_indices.size());
676 for (int i = 0; i != local_indices.size(); ++i)
677 if (local_indices[i] != -1)
678 dotVector[i] = array[local_indices[i]];
679 else
680 dotVector[i] = 0;
681
682 switch (CTX) {
684 CHKERR restore_array(getFEMethod()->ts_u);
685 break;
687 CHKERR restore_array(getFEMethod()->dx);
688 break;
690 CHKERR restore_array(getFEMethod()->ts_u_t);
691 break;
693 CHKERR restore_array(getFEMethod()->ts_u_tt);
694 break;
695 default:
696 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
697 "That case is not implemented");
698 }
699
700 const size_t nb_base_functions = data.getN().size2();
701 auto base_function = data.getFTensor0N();
702 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
703
704 FTensor::Index<'I', Tensor_Dim> I;
705 const size_t size = nb_dofs / Tensor_Dim;
706 if (nb_dofs % Tensor_Dim) {
707 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
708 }
709 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
711 size_t bb = 0;
712 for (; bb != size; ++bb) {
713 values_at_gauss_pts(I) += field_data(I) * base_function;
714 ++field_data;
715 ++base_function;
716 }
717 // Number of dofs can be smaller than number of Tensor_Dim x base
718 // functions
719 for (; bb < nb_base_functions; ++bb)
720 ++base_function;
721 ++values_at_gauss_pts;
722 }
724 }
725
726protected:
727 boost::shared_ptr<MatrixDouble> dataPtr;
731};
732
733/** \brief Get rate of values at integration pts for tensor field
734 * rank 1, i.e. vector field
735 *
736 * \ingroup mofem_forces_and_sources_user_data_operators
737 */
738template <int Tensor_Dim>
742
743/** \brief Get second rate of values at integration pts for tensor
744 * field rank 1, i.e. vector field
745 *
746 * \ingroup mofem_forces_and_sources_user_data_operators
747 */
748template <int Tensor_Dim>
752
753/** \brief Get second time second update vector at integration pts for tensor
754 * field rank 1, i.e. vector field
755 *
756 * \ingroup mofem_forces_and_sources_user_data_operators
757 */
758template <int Tensor_Dim>
762
763/**@}*/
764
765/** \name Tensor field values at integration points */
766
767/**@{*/
768
769/** \brief Calculate field values for tenor field rank 2.
770 *
771 */
772template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
775
777 const std::string field_name,
778 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
779 const EntityType zero_type = MBVERTEX)
782 dataPtr(data_ptr), zeroType(zero_type) {
783 if (!dataPtr)
784 THROW_MESSAGE("Pointer is not set");
785 }
786
787 MoFEMErrorCode doWork(int side, EntityType type,
789
790protected:
791 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
793};
794
795template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
798 doWork(int side, EntityType type, EntitiesFieldData::EntData &data) {
800 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
801 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
802 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
803 Tensor_Dim0, Tensor_Dim1);
805}
806
807template <int Tensor_Dim0, int Tensor_Dim1>
808struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
809 ublas::row_major, DoubleAllocator>
811
813 const std::string field_name,
814 boost::shared_ptr<
815 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
816 data_ptr,
817 const EntityType zero_type = MBVERTEX)
820 dataPtr(data_ptr), zeroType(zero_type) {
821 if (!dataPtr)
822 THROW_MESSAGE("Pointer is not set");
823 }
824
826 const std::string field_name,
827 boost::shared_ptr<
828 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
829 data_ptr,
830 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
833 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
834 if (!dataPtr)
835 THROW_MESSAGE("Pointer is not set");
836 }
837
838 MoFEMErrorCode doWork(int side, EntityType type,
840
841protected:
842 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
847};
848
849template <int Tensor_Dim0, int Tensor_Dim1>
851 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
852 DoubleAllocator>::doWork(int side, EntityType type,
855
856 MatrixDouble &mat = *dataPtr;
857 const size_t nb_gauss_pts = data.getN().size1();
858 if (type == zeroType) {
859 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
860 mat.clear();
861 }
862
863 const size_t nb_dofs = data.getFieldData().size();
864
865 if (dataVec.use_count()) {
866 dotVector.resize(nb_dofs, false);
867 const double *array;
868 CHKERR VecGetArrayRead(dataVec, &array);
869 const auto &local_indices = data.getLocalIndices();
870 for (int i = 0; i != local_indices.size(); ++i)
871 if (local_indices[i] != -1)
872 dotVector[i] = array[local_indices[i]];
873 else
874 dotVector[i] = 0;
875 CHKERR VecRestoreArrayRead(dataVec, &array);
876 data.getFieldData().swap(dotVector);
877 }
878
879 if (nb_dofs) {
880 const size_t nb_base_functions = data.getN().size2();
881 auto base_function = data.getFTensor0N();
882 auto values_at_gauss_pts =
884 FTensor::Index<'i', Tensor_Dim0> i;
885 FTensor::Index<'j', Tensor_Dim1> j;
886 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
887 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
888 auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
889 size_t bb = 0;
890 for (; bb != size; ++bb) {
891 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
892 ++field_data;
893 ++base_function;
894 }
895 for (; bb < nb_base_functions; ++bb)
896 ++base_function;
897 ++values_at_gauss_pts;
898 }
899
900 if (dataVec.use_count()) {
901 data.getFieldData().swap(dotVector);
902 }
903 }
905}
906
907/** \brief Get values at integration pts for tensor field rank 2, i.e. matrix
908 * field
909 *
910 * \ingroup mofem_forces_and_sources_user_data_operators
911 */
912template <int Tensor_Dim0, int Tensor_Dim1>
915 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
916
918 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
920};
921
922/** \brief Get time direvarive values at integration pts for tensor field rank
923 * 2, i.e. matrix field
924 *
925 * \ingroup mofem_forces_and_sources_user_data_operators
926 */
927template <int Tensor_Dim0, int Tensor_Dim1>
930
932 boost::shared_ptr<MatrixDouble> data_ptr,
933 const EntityType zero_at_type = MBVERTEX)
936 dataPtr(data_ptr), zeroAtType(zero_at_type) {
937 if (!dataPtr)
938 THROW_MESSAGE("Pointer is not set");
939 }
940
941 MoFEMErrorCode doWork(int side, EntityType type,
944
945 const size_t nb_gauss_pts = getGaussPts().size2();
946 MatrixDouble &mat = *dataPtr;
947 if (type == zeroAtType) {
948 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
949 mat.clear();
950 }
951 const auto &local_indices = data.getLocalIndices();
952 const size_t nb_dofs = local_indices.size();
953 if (nb_dofs) {
954 dotVector.resize(nb_dofs, false);
955 const double *array;
956 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
957 for (size_t i = 0; i != local_indices.size(); ++i)
958 if (local_indices[i] != -1)
959 dotVector[i] = array[local_indices[i]];
960 else
961 dotVector[i] = 0;
962 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
963
964 const size_t nb_base_functions = data.getN().size2();
965
966 auto base_function = data.getFTensor0N();
967 auto values_at_gauss_pts =
969 FTensor::Index<'i', Tensor_Dim0> i;
970 FTensor::Index<'j', Tensor_Dim1> j;
971 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
972 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
973 auto field_data =
975 size_t bb = 0;
976 for (; bb != size; ++bb) {
977 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
978 ++field_data;
979 ++base_function;
980 }
981 for (; bb < nb_base_functions; ++bb)
982 ++base_function;
983 ++values_at_gauss_pts;
984 }
985 }
987 }
988
989protected:
990 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
991 EntityType zeroAtType; ///< Zero values at Gauss point at this type
992 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
993};
994
995/**
996 * @brief Calculate symmetric tensor field values at integration pts.
997 *
998 * @tparam Tensor_Dim
999
1000 * \ingroup mofem_forces_and_sources_user_data_operators
1001 */
1002template <int Tensor_Dim>
1005
1007 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1008 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1011 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1012 if (!dataPtr)
1013 THROW_MESSAGE("Pointer is not set");
1014 }
1015
1017 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1018 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBEDGE,
1019 const int zero_side = 0)
1022 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side),
1023 dataVec(data_vec) {
1024 if (!dataPtr)
1025 THROW_MESSAGE("Pointer is not set");
1026 }
1027
1028 MoFEMErrorCode doWork(int side, EntityType type,
1031 MatrixDouble &mat = *dataPtr;
1032 const int nb_gauss_pts = getGaussPts().size2();
1033 if (type == this->zeroType && side == zeroSide) {
1034 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
1035 mat.clear();
1036 }
1037 const int nb_dofs = data.getFieldData().size();
1038 if (!nb_dofs)
1040
1041 if (dataVec.use_count()) {
1042 dotVector.resize(nb_dofs, false);
1043 const double *array;
1044 CHKERR VecGetArrayRead(dataVec, &array);
1045 const auto &local_indices = data.getLocalIndices();
1046 for (int i = 0; i != local_indices.size(); ++i)
1047 if (local_indices[i] != -1)
1048 dotVector[i] = array[local_indices[i]];
1049 else
1050 dotVector[i] = 0;
1051 CHKERR VecRestoreArrayRead(dataVec, &array);
1052 data.getFieldData().swap(dotVector);
1053 }
1054
1055 const int nb_base_functions = data.getN().size2();
1056 auto base_function = data.getFTensor0N();
1057 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1058 FTensor::Index<'i', Tensor_Dim> i;
1059 FTensor::Index<'j', Tensor_Dim> j;
1060 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1061 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1062 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
1063 int bb = 0;
1064 for (; bb != size; ++bb) {
1065 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1066 ++field_data;
1067 ++base_function;
1068 }
1069 for (; bb < nb_base_functions; ++bb)
1070 ++base_function;
1071 ++values_at_gauss_pts;
1072 }
1073
1074 if (dataVec.use_count()) {
1075 data.getFieldData().swap(dotVector);
1076 }
1077
1079 }
1080
1081protected:
1082 boost::shared_ptr<MatrixDouble> dataPtr;
1084 const int zeroSide;
1087};
1088
1089/**
1090 * @brief Calculate symmetric tensor field rates ant integratio pts.
1091 *
1092 * @tparam Tensor_Dim
1093 *
1094 * \ingroup mofem_forces_and_sources_user_data_operators
1095 */
1096template <int Tensor_Dim>
1099
1101 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1102 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1105 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1106 if (!dataPtr)
1107 THROW_MESSAGE("Pointer is not set");
1108 }
1109
1110 MoFEMErrorCode doWork(int side, EntityType type,
1113 const int nb_gauss_pts = getGaussPts().size2();
1114 MatrixDouble &mat = *dataPtr;
1115 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1116 if (type == zeroType && side == zeroSide) {
1117 mat.resize(symm_size, nb_gauss_pts, false);
1118 mat.clear();
1119 }
1120 auto &local_indices = data.getLocalIndices();
1121 const int nb_dofs = local_indices.size();
1122 if (!nb_dofs)
1124
1125 #ifndef NDEBUG
1126 if ((getFEMethod()->data_ctx & PetscData::CtxSetX_T).none()) {
1127 MOFEM_LOG_CHANNEL("SELF");
1128 MOFEM_LOG("SELF", Sev::error)
1129 << "In this case field degrees of freedom are read from vector. "
1130 "That usually happens when time solver is used, and acces to "
1131 "first rates is needed. You probably not set "
1132 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1133 "respectively";
1134 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
1135 }
1136 #endif
1137
1138 dotVector.resize(nb_dofs, false);
1139 const double *array;
1140 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1141 for (int i = 0; i != local_indices.size(); ++i)
1142 if (local_indices[i] != -1)
1143 dotVector[i] = array[local_indices[i]];
1144 else
1145 dotVector[i] = 0;
1146 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1147
1148 const int nb_base_functions = data.getN().size2();
1149
1150 auto base_function = data.getFTensor0N();
1151 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1152 FTensor::Index<'i', Tensor_Dim> i;
1153 FTensor::Index<'j', Tensor_Dim> j;
1154 const int size = nb_dofs / symm_size;
1155 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1156 auto field_data = getFTensorDotData<Tensor_Dim>();
1157 int bb = 0;
1158 for (; bb != size; ++bb) {
1159 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1160 ++field_data;
1161 ++base_function;
1162 }
1163 for (; bb < nb_base_functions; ++bb)
1164 ++base_function;
1165 ++values_at_gauss_pts;
1166 }
1167
1169 }
1170
1171protected:
1172 boost::shared_ptr<MatrixDouble> dataPtr;
1174 const int zeroSide;
1176
1177 template <int Dim> inline auto getFTensorDotData() {
1178 static_assert(Dim || !Dim, "not implemented");
1179 }
1180};
1181
1182template <>
1183template <>
1184inline auto
1187 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1188 &dotVector[5]);
1189}
1190
1191template <>
1192template <>
1193inline auto
1196 &dotVector[0], &dotVector[1], &dotVector[2]);
1197}
1198
1199/**@}*/
1200
1201/** \name Gradients and Hessian of scalar fields at integration points */
1202
1203/**@{*/
1204
1205/**
1206 * \brief Evaluate field gradient values for scalar field, i.e. gradient is
1207 * tensor rank 1 (vector)
1208 *
1209 */
1210template <int Tensor_Dim, class T, class L, class A>
1217
1218/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1219 * tensor rank 1 (vector), specialization
1220 *
1221 */
1222template <int Tensor_Dim>
1224 ublas::row_major, DoubleAllocator>
1226 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1227
1229 Tensor_Dim, double, ublas::row_major,
1231
1232 /**
1233 * \brief calculate gradient values of scalar field at integration points
1234 * @param side side entity number
1235 * @param type side entity type
1236 * @param data entity data
1237 * @return error code
1238 */
1239 MoFEMErrorCode doWork(int side, EntityType type,
1241};
1242
1243/**
1244 * \brief Member function specialization calculating scalar field gradients for
1245 * tenor field rank 1
1246 *
1247 */
1248template <int Tensor_Dim>
1250 Tensor_Dim, double, ublas::row_major,
1251 DoubleAllocator>::doWork(int side, EntityType type,
1254
1255 const size_t nb_gauss_pts = this->getGaussPts().size2();
1256 auto &mat = *this->dataPtr;
1257 if (type == this->zeroType) {
1258 mat.resize(Tensor_Dim, nb_gauss_pts, false);
1259 mat.clear();
1260 }
1261
1262 const int nb_dofs = data.getFieldData().size();
1263 if (nb_dofs) {
1264
1265 if (nb_gauss_pts) {
1266 const int nb_base_functions = data.getN().size2();
1267 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
1268 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1269
1270 #ifndef NDEBUG
1271 if (nb_dofs > nb_base_functions)
1272 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1273 "Number of base functions inconsistent with number of DOFs "
1274 "(%d > %d)",
1275 nb_dofs, nb_base_functions);
1276
1277 if (data.getDiffN().size2() != nb_base_functions * Tensor_Dim)
1278 SETERRQ(
1279 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1280 "Number of base functions inconsistent with number of derivatives "
1281 "(%lu != %d)",
1282 data.getDiffN().size2(), nb_base_functions);
1283
1284 if (data.getDiffN().size1() != nb_gauss_pts)
1285 SETERRQ(
1286 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1287 "Number of base functions inconsistent with number of integration "
1288 "pts (%lu != %lu)",
1289 data.getDiffN().size2(), nb_gauss_pts);
1290
1291 #endif
1292
1293 FTensor::Index<'I', Tensor_Dim> I;
1294 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1295 auto field_data = data.getFTensor0FieldData();
1296 int bb = 0;
1297 for (; bb != nb_dofs; ++bb) {
1298 gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
1299 ++field_data;
1300 ++diff_base_function;
1301 }
1302 // Number of dofs can be smaller than number of base functions
1303 for (; bb < nb_base_functions; ++bb)
1304 ++diff_base_function;
1305 ++gradients_at_gauss_pts;
1306 }
1307 }
1308 }
1309
1311}
1312
1313/** \brief Get field gradients at integration pts for scalar field rank 0, i.e.
1314 * vector field
1315 *
1316 * \ingroup mofem_forces_and_sources_user_data_operators
1317 */
1318template <int Tensor_Dim>
1321 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1323 Tensor_Dim, double, ublas::row_major,
1324 DoubleAllocator>::OpCalculateScalarFieldGradient_General;
1325};
1326
1327/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1328 * tensor rank 1 (vector), specialization
1329 *
1330 */
1331template <int Tensor_Dim>
1334 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1335
1337 Tensor_Dim, double, ublas::row_major,
1339
1340 /**
1341 * \brief calculate gradient values of scalar field at integration points
1342 * @param side side entity number
1343 * @param type side entity type
1344 * @param data entity data
1345 * @return error code
1346 */
1347 MoFEMErrorCode doWork(int side, EntityType type,
1349};
1350
1351template <int Tensor_Dim>
1353 int side, EntityType type, EntitiesFieldData::EntData &data) {
1355
1356 const size_t nb_gauss_pts = this->getGaussPts().size2();
1357 auto &mat = *this->dataPtr;
1358 if (type == this->zeroType) {
1359 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts, false);
1360 mat.clear();
1361 }
1362
1363 const int nb_dofs = data.getFieldData().size();
1364 if (nb_dofs) {
1365
1366 if (nb_gauss_pts) {
1367 const int nb_base_functions = data.getN().size2();
1368
1369 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1370 #ifndef NDEBUG
1371 if (hessian_base.size1() != nb_gauss_pts) {
1372 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1373 "Wrong number of integration pts (%ld != %d)",
1374 hessian_base.size1(), nb_gauss_pts);
1375 }
1376 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1377 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1378 "Wrong number of base functions (%ld != %d)",
1379 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1380 nb_base_functions);
1381 }
1382 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1383 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1384 "Wrong number of base functions (%ld < %d)",
1385 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1386 }
1387 #endif
1388
1389 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1390 &*hessian_base.data().begin());
1391
1392 auto hessian_at_gauss_pts =
1394
1395 FTensor::Index<'I', Tensor_Dim> I;
1396 FTensor::Index<'J', Tensor_Dim> J;
1397 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1398 auto field_data = data.getFTensor0FieldData();
1399 int bb = 0;
1400 for (; bb != nb_dofs; ++bb) {
1401 hessian_at_gauss_pts(I, J) +=
1402 field_data * t_diff2_base_function(I, J);
1403 ++field_data;
1404 ++t_diff2_base_function;
1405 }
1406 // Number of dofs can be smaller than number of base functions
1407 for (; bb < nb_base_functions; ++bb) {
1408 ++t_diff2_base_function;
1409 }
1410
1411 ++hessian_at_gauss_pts;
1412 }
1413 }
1414 }
1415
1417}
1418
1419/**}*/
1420
1421/** \name Gradients and hessian of tensor fields at integration points */
1422
1423/**@{*/
1424
1425/**
1426 * \brief Evaluate field gradient values for vector field, i.e. gradient is
1427 * tensor rank 2
1428 *
1429 * \tparam Tensor_Dim0 Dimension of the vector field
1430 * \tparam Tensor_Dim1 Dimension of the gradient (usually spatial dimension)
1431 * \tparam S Storage order for the field data
1432 * \tparam T Data type
1433 * \tparam L Layout type
1434 * \tparam A Allocator type
1435 *
1436 */
1437template <int Tensor_Dim0, int Tensor_Dim1, int S, class T, class L, class A>
1439 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1440 L, A> {
1441
1443 Tensor_Dim0, Tensor_Dim1, T, L, A>::OpCalculateTensor2FieldValues_General;
1444};
1445
1446template <int Tensor_Dim0, int Tensor_Dim1, int S>
1448 Tensor_Dim0, Tensor_Dim1, S, double, ublas::row_major, DoubleAllocator>
1450 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1451
1453 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1455
1456 /**
1457 * \brief calculate values of vector field at integration points
1458 * @param side side entity number
1459 * @param type side entity type
1460 * @param data entity data
1461 * @return error code
1462 */
1463 MoFEMErrorCode doWork(int side, EntityType type,
1465};
1466
1467/**
1468 * \brief Member function specialization calculating vector field gradients for
1469 * tenor field rank 2
1470 *
1471 */
1472template <int Tensor_Dim0, int Tensor_Dim1, int S>
1474 Tensor_Dim0, Tensor_Dim1, S, double, ublas::row_major,
1475 DoubleAllocator>::doWork(int side, EntityType type,
1478 if (!this->dataPtr)
1479 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1480 "Data pointer not allocated");
1481
1482 const size_t nb_gauss_pts = this->getGaussPts().size2();
1483 auto &mat = *this->dataPtr;
1484 if (type == this->zeroType) {
1485 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1486 mat.clear();
1487 }
1488
1489 if (nb_gauss_pts) {
1490 const size_t nb_dofs = data.getFieldData().size();
1491
1492 if (nb_dofs) {
1493
1494 if (this->dataVec.use_count()) {
1495 this->dotVector.resize(nb_dofs, false);
1496 const double *array;
1497 CHKERR VecGetArrayRead(this->dataVec, &array);
1498 const auto &local_indices = data.getLocalIndices();
1499 for (int i = 0; i != local_indices.size(); ++i)
1500 if (local_indices[i] != -1)
1501 this->dotVector[i] = array[local_indices[i]];
1502 else
1503 this->dotVector[i] = 0;
1504 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1505 data.getFieldData().swap(this->dotVector);
1506 }
1507
1508 const int nb_base_functions = data.getN().size2();
1509 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1510 auto gradients_at_gauss_pts =
1512 FTensor::Index<'I', Tensor_Dim0> I;
1513 FTensor::Index<'J', Tensor_Dim1> J;
1514 int size = nb_dofs / Tensor_Dim0;
1515 if (nb_dofs % Tensor_Dim0) {
1516 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1517 "Data inconsistency");
1518 }
1519 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1520 auto field_data = getFTensor1FromPtr<Tensor_Dim0, S>(
1521 data.getFieldData().data().data());
1522
1523 #ifndef NDEBUG
1524 if (field_data.l2() != field_data.l2()) {
1525 MOFEM_LOG("SELF", Sev::error) << "field data " << field_data;
1526 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1527 "Wrong number in coefficients");
1528 }
1529 #endif
1530
1531 int bb = 0;
1532 for (; bb < size; ++bb) {
1533 #ifndef SINGULARITY
1534 #ifndef NDEBUG
1535 if (diff_base_function.l2() != diff_base_function.l2()) {
1536 MOFEM_LOG("SELF", Sev::error)
1537 << "diff_base_function: " << diff_base_function;
1538 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1539 "Wrong number number in base functions");
1540 }
1541 #endif
1542 #endif
1543
1544 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1545 ++field_data;
1546 ++diff_base_function;
1547 }
1548 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1549 // functions
1550 for (; bb != nb_base_functions; ++bb)
1551 ++diff_base_function;
1552 ++gradients_at_gauss_pts;
1553 }
1554
1555 if (this->dataVec.use_count()) {
1556 data.getFieldData().swap(this->dotVector);
1557 }
1558 }
1559 }
1561}
1562
1563/** \brief Get field gradients at integration pts for scalar field rank 0, i.e.
1564 * vector field
1565 *
1566 * \tparam Tensor_Dim0 Dimension of the vector field
1567 * \tparam Tensor_Dim1 Dimension of the gradient (usually spatial dimension)
1568 * \tparam S Stride in the storage of the vector field, default is Tensor_Dim0
1569 *
1570 * \ingroup mofem_forces_and_sources_user_data_operators
1571 */
1572template <int Tensor_Dim0, int Tensor_Dim1, int S = Tensor_Dim0>
1574 : public OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, S,
1575 double, ublas::row_major,
1576 DoubleAllocator> {
1577
1579 Tensor_Dim0, Tensor_Dim1, S, double, ublas::row_major,
1580 DoubleAllocator>::OpCalculateVectorFieldGradient_General;
1581};
1582
1583/** \brief Get field gradients time derivative at integration pts for scalar
1584 * field rank 0, i.e. vector field
1585 *
1586 * \ingroup mofem_forces_and_sources_user_data_operators
1587 */
1588template <int Tensor_Dim0, int Tensor_Dim1>
1591
1593 boost::shared_ptr<MatrixDouble> data_ptr,
1594 const EntityType zero_at_type = MBVERTEX)
1597 dataPtr(data_ptr), zeroAtType(zero_at_type) {
1598 if (!dataPtr)
1599 THROW_MESSAGE("Pointer is not set");
1600 }
1601
1602 MoFEMErrorCode doWork(int side, EntityType type,
1605
1606 const auto &local_indices = data.getLocalIndices();
1607 const int nb_dofs = local_indices.size();
1608 const int nb_gauss_pts = this->getGaussPts().size2();
1609
1610 auto &mat = *this->dataPtr;
1611 if (type == this->zeroAtType) {
1612 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1613 mat.clear();
1614 }
1615 if (!nb_dofs)
1617
1618 dotVector.resize(nb_dofs, false);
1619 const double *array;
1620 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1621 for (int i = 0; i != local_indices.size(); ++i)
1622 if (local_indices[i] != -1)
1623 dotVector[i] = array[local_indices[i]];
1624 else
1625 dotVector[i] = 0;
1626 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1627
1628 const int nb_base_functions = data.getN().size2();
1629 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1630 auto gradients_at_gauss_pts =
1632 FTensor::Index<'I', Tensor_Dim0> I;
1633 FTensor::Index<'J', Tensor_Dim1> J;
1634 int size = nb_dofs / Tensor_Dim0;
1635 if (nb_dofs % Tensor_Dim0) {
1636 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1637 }
1638
1639 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1640 auto field_data = getFTensor1FromPtr<Tensor_Dim0>(&*dotVector.begin());
1641 int bb = 0;
1642 for (; bb < size; ++bb) {
1643 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1644 ++field_data;
1645 ++diff_base_function;
1646 }
1647 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1648 // functions
1649 for (; bb != nb_base_functions; ++bb)
1650 ++diff_base_function;
1651 ++gradients_at_gauss_pts;
1652 }
1654 }
1655
1656private:
1657 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1658 EntityType zeroAtType; ///< Zero values at Gauss point at this type
1659 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
1660};
1661
1662/**
1663 * \brief Evaluate field gradient values for symmetric 2nd order tensor field,
1664 * i.e. gradient is tensor rank 3
1665 *
1666 */
1667template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1669 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1670 L, A> {
1671
1673 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1674 const EntityType zero_type = MBVERTEX)
1675 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
1676 A>(field_name, data_ptr,
1677 zero_type) {}
1678};
1679
1680template <int Tensor_Dim0, int Tensor_Dim1>
1682 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator>
1684 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1685
1687 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1688 const EntityType zero_type = MBVERTEX)
1689 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1690 ublas::row_major,
1692 field_name, data_ptr, zero_type) {}
1693
1694 /**
1695 * \brief calculate values of vector field at integration points
1696 * @param side side entity number
1697 * @param type side entity type
1698 * @param data entity data
1699 * @return error code
1700 */
1701 MoFEMErrorCode doWork(int side, EntityType type,
1703};
1704
1705/**
1706 * \brief Member function specialization calculating tensor field gradients for
1707 * symmetric tensor field rank 2
1708 *
1709 */
1710template <int Tensor_Dim0, int Tensor_Dim1>
1712 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1713 DoubleAllocator>::doWork(int side, EntityType type,
1716 if (!this->dataPtr)
1717 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1718 "Data pointer not allocated");
1719
1720 const size_t nb_gauss_pts = this->getGaussPts().size2();
1721 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1722 auto &mat = *this->dataPtr;
1723 if (type == this->zeroType) {
1724 mat.resize(msize * Tensor_Dim1, nb_gauss_pts, false);
1725 mat.clear();
1726 }
1727
1728 if (nb_gauss_pts) {
1729 const size_t nb_dofs = data.getFieldData().size();
1730
1731 if (nb_dofs) {
1732
1733 const int nb_base_functions = data.getN().size2();
1734 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1735 auto gradients_at_gauss_pts =
1737 FTensor::Index<'I', Tensor_Dim0> I;
1738 FTensor::Index<'J', Tensor_Dim0> J;
1739 FTensor::Index<'K', Tensor_Dim1> K;
1740 int size = nb_dofs / msize;
1741 if (nb_dofs % msize) {
1742 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1743 "Data inconsistency");
1744 }
1745 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1746 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim0>();
1747 int bb = 0;
1748 for (; bb < size; ++bb) {
1749 gradients_at_gauss_pts(I, J, K) +=
1750 field_data(I, J) * diff_base_function(K);
1751 ++field_data;
1752 ++diff_base_function;
1753 }
1754 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1755 // functions
1756 for (; bb != nb_base_functions; ++bb)
1757 ++diff_base_function;
1758 ++gradients_at_gauss_pts;
1759 }
1760 }
1761 }
1763}
1764
1765/** \brief Get field gradients at integration pts for symmetric tensorial field
1766 * rank 2
1767 *
1768 * \ingroup mofem_forces_and_sources_user_data_operators
1769 */
1770template <int Tensor_Dim0, int Tensor_Dim1>
1773 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1774
1776 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1777 const EntityType zero_type = MBVERTEX)
1779 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1780 DoubleAllocator>(field_name, data_ptr, zero_type) {}
1781};
1782
1783template <int Tensor_Dim0, int Tensor_Dim1>
1786 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1787
1789 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1791
1792 /**
1793 * \brief calculate values of vector field at integration points
1794 * @param side side entity number
1795 * @param type side entity type
1796 * @param data entity data
1797 * @return error code
1798 */
1799 MoFEMErrorCode doWork(int side, EntityType type,
1801};
1802
1803/**
1804 * \brief Member function specialization calculating vector field gradients for
1805 * tenor field rank 2
1806 *
1807 */
1808template <int Tensor_Dim0, int Tensor_Dim1>
1810 int side, EntityType type, EntitiesFieldData::EntData &data) {
1812 if (!this->dataPtr)
1813 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1814 "Data pointer not allocated");
1815
1816 const size_t nb_gauss_pts = this->getGaussPts().size2();
1817 auto &mat = *this->dataPtr;
1818 if (type == this->zeroType) {
1819 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts, false);
1820 mat.clear();
1821 }
1822
1823 if (nb_gauss_pts) {
1824 const size_t nb_dofs = data.getFieldData().size();
1825
1826 if (nb_dofs) {
1827
1828 if (this->dataVec.use_count()) {
1829 this->dotVector.resize(nb_dofs, false);
1830 const double *array;
1831 CHKERR VecGetArrayRead(this->dataVec, &array);
1832 const auto &local_indices = data.getLocalIndices();
1833 for (int i = 0; i != local_indices.size(); ++i)
1834 if (local_indices[i] != -1)
1835 this->dotVector[i] = array[local_indices[i]];
1836 else
1837 this->dotVector[i] = 0;
1838 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1839 data.getFieldData().swap(this->dotVector);
1840 }
1841
1842 const int nb_base_functions = data.getN().size2();
1843
1844 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1845 #ifndef NDEBUG
1846 if (hessian_base.size1() != nb_gauss_pts) {
1847 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1848 "Wrong number of integration pts (%d != %d)",
1849 hessian_base.size1(), nb_gauss_pts);
1850 }
1851 if (hessian_base.size2() !=
1852 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1853 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1854 "Wrong number of base functions (%d != %d)",
1855 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1856 nb_base_functions);
1857 }
1858 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1859 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1860 "Wrong number of base functions (%d < %d)",
1861 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1862 }
1863 #endif
1864
1865 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1866 &*hessian_base.data().begin());
1867
1868 auto t_hessian_at_gauss_pts =
1870
1871 FTensor::Index<'I', Tensor_Dim0> I;
1872 FTensor::Index<'J', Tensor_Dim1> J;
1873 FTensor::Index<'K', Tensor_Dim1> K;
1874
1875 int size = nb_dofs / Tensor_Dim0;
1876 #ifndef NDEBUG
1877 if (nb_dofs % Tensor_Dim0) {
1878 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1879 "Data inconsistency");
1880 }
1881 #endif
1882
1883 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1884 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1885 int bb = 0;
1886 for (; bb < size; ++bb) {
1887 t_hessian_at_gauss_pts(I, J, K) +=
1888 field_data(I) * t_diff2_base_function(J, K);
1889 ++field_data;
1890 ++t_diff2_base_function;
1891 }
1892 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1893 // functions
1894 for (; bb != nb_base_functions; ++bb)
1895 ++t_diff2_base_function;
1896 ++t_hessian_at_gauss_pts;
1897 }
1898
1899 if (this->dataVec.use_count()) {
1900 data.getFieldData().swap(this->dotVector);
1901 }
1902 }
1903 }
1905}
1906
1907/**@}*/
1908
1909/** \name Transform tensors and vectors */
1910
1911/**@{*/
1912
1913/**
1914 * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1915 * \f$
1916 *
1917 * @tparam DIM
1918 *
1919 * \ingroup mofem_forces_and_sources_user_data_operators
1920 */
1921template <int DIM_01, int DIM_23, int S = 0>
1924
1927
1928 /**
1929 * @deprecated Do not use this constructor
1930 */
1933 boost::shared_ptr<MatrixDouble> in_mat,
1934 boost::shared_ptr<MatrixDouble> out_mat,
1935 boost::shared_ptr<MatrixDouble> d_mat)
1936 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1937 // Only is run for vertices
1938 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1939 if (!inMat)
1940 THROW_MESSAGE("Pointer for in mat is null");
1941 if (!outMat)
1942 THROW_MESSAGE("Pointer for out mat is null");
1943 if (!dMat)
1944 THROW_MESSAGE("Pointer for tensor mat is null");
1945 }
1946
1947 OpTensorTimesSymmetricTensor(boost::shared_ptr<MatrixDouble> in_mat,
1948 boost::shared_ptr<MatrixDouble> out_mat,
1949 boost::shared_ptr<MatrixDouble> d_mat)
1950 : UserOp(NOSPACE, OPSPACE), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1951 // Only is run for vertices
1952 if (!inMat)
1953 THROW_MESSAGE("Pointer for in mat is null");
1954 if (!outMat)
1955 THROW_MESSAGE("Pointer for out mat is null");
1956 if (!dMat)
1957 THROW_MESSAGE("Pointer for tensor mat is null");
1958 }
1959
1960 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1962 const size_t nb_gauss_pts = getGaussPts().size2();
1965 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1967 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1968 t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1969 ++t_in;
1970 ++t_out;
1971 }
1973 }
1974
1975private:
1976 FTensor::Index<'i', DIM_01> i;
1977 FTensor::Index<'j', DIM_01> j;
1978 FTensor::Index<'k', DIM_23> k;
1979 FTensor::Index<'l', DIM_23> l;
1980
1981 boost::shared_ptr<MatrixDouble> inMat;
1982 boost::shared_ptr<MatrixDouble> outMat;
1983 boost::shared_ptr<MatrixDouble> dMat;
1984};
1985
1986template <int DIM>
1988
1991
1992 /**
1993 * @deprecated Do not use this constructor
1994 */
1996 boost::shared_ptr<MatrixDouble> in_mat,
1997 boost::shared_ptr<MatrixDouble> out_mat)
1998 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1999 // Only is run for vertices
2000 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
2001 if (!inMat)
2002 THROW_MESSAGE("Pointer not set for in matrix");
2003 if (!outMat)
2004 THROW_MESSAGE("Pointer not set for in matrix");
2005 }
2006
2007 OpSymmetrizeTensor(boost::shared_ptr<MatrixDouble> in_mat,
2008 boost::shared_ptr<MatrixDouble> out_mat)
2009 : UserOp(NOSPACE, OPSPACE), inMat(in_mat), outMat(out_mat) {
2010 // Only is run for vertices
2011 if (!inMat)
2012 THROW_MESSAGE("Pointer not set for in matrix");
2013 if (!outMat)
2014 THROW_MESSAGE("Pointer not set for in matrix");
2015 }
2016
2017 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
2019 const size_t nb_gauss_pts = getGaussPts().size2();
2020 auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
2021 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
2023 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
2024 t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
2025 ++t_in;
2026 ++t_out;
2027 }
2029 }
2030
2031private:
2034 boost::shared_ptr<MatrixDouble> inMat;
2035 boost::shared_ptr<MatrixDouble> outMat;
2036};
2037
2039
2042
2043 /**
2044 * @deprecated Do not use this constructor
2045 */
2046 DEPRECATED OpScaleMatrix(const std::string field_name, const double scale,
2047 boost::shared_ptr<MatrixDouble> in_mat,
2048 boost::shared_ptr<MatrixDouble> out_mat)
2049 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
2050 scalePtr = boost::make_shared<double>(scale);
2051 // Only is run for vertices
2052 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
2053 if (!inMat)
2054 THROW_MESSAGE("Pointer not set for in matrix");
2055 if (!outMat)
2056 THROW_MESSAGE("Pointer not set for in matrix");
2057 }
2058
2059 OpScaleMatrix(boost::shared_ptr<double> scale_ptr,
2060 boost::shared_ptr<MatrixDouble> in_mat,
2061 boost::shared_ptr<MatrixDouble> out_mat)
2062 : UserOp(NOSPACE, OPSPACE), scalePtr(scale_ptr), inMat(in_mat),
2063 outMat(out_mat) {
2064 if (!scalePtr)
2065 THROW_MESSAGE("Pointer not set for scale");
2066 if (!inMat)
2067 THROW_MESSAGE("Pointer not set for in matrix");
2068 if (!outMat)
2069 THROW_MESSAGE("Pointer not set for in matrix");
2070 }
2071
2072 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
2074 outMat->resize(inMat->size1(), inMat->size2(), false);
2075 noalias(*outMat) = (*scalePtr) * (*inMat);
2077 }
2078
2079private:
2080 boost::shared_ptr<double> scalePtr;
2081 boost::shared_ptr<MatrixDouble> inMat;
2082 boost::shared_ptr<MatrixDouble> outMat;
2083};
2084
2085/**@}*/
2086
2087/** \name H-div/H-curls (Vectorial bases) values at integration points */
2088
2089/**@{*/
2090
2091/** \brief Get vector field for H-div approximation
2092 * \ingroup mofem_forces_and_sources_user_data_operators
2093 */
2094template <int Base_Dim, int Field_Dim, class T, class L, class A>
2096
2097/** \brief Get vector field for H-div approximation
2098 * \ingroup mofem_forces_and_sources_user_data_operators
2099 */
2100template <int Field_Dim>
2102 ublas::row_major, DoubleAllocator>
2104
2106 boost::shared_ptr<MatrixDouble> data_ptr,
2107 SmartPetscObj<Vec> data_vec,
2108 const EntityType zero_type = MBEDGE,
2109 const int zero_side = 0)
2112 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type),
2113 zeroSide(zero_side) {
2114 if (!dataPtr)
2115 THROW_MESSAGE("Pointer is not set");
2116 }
2117
2119 boost::shared_ptr<MatrixDouble> data_ptr,
2120 const EntityType zero_type = MBEDGE,
2121 const int zero_side = 0)
2123 field_name, data_ptr, SmartPetscObj<Vec>(), zero_type, zero_side) {}
2124
2125 /**
2126 * \brief Calculate values of vector field at integration points
2127 * @param side side entity number
2128 * @param type side entity type
2129 * @param data entity data
2130 * @return error code
2131 */
2132 MoFEMErrorCode doWork(int side, EntityType type,
2134
2135private:
2136 boost::shared_ptr<MatrixDouble> dataPtr;
2139 const int zeroSide;
2141};
2142
2143template <int Field_Dim>
2145 3, Field_Dim, double, ublas::row_major,
2146 DoubleAllocator>::doWork(int side, EntityType type,
2149 const size_t nb_integration_points = this->getGaussPts().size2();
2150 if (type == zeroType && side == zeroSide) {
2151 dataPtr->resize(Field_Dim, nb_integration_points, false);
2152 dataPtr->clear();
2153 }
2154 const size_t nb_dofs = data.getFieldData().size();
2155 if (!nb_dofs)
2157
2158 if (dataVec.use_count()) {
2159 dotVector.resize(nb_dofs, false);
2160 const double *array;
2161 CHKERR VecGetArrayRead(dataVec, &array);
2162 const auto &local_indices = data.getLocalIndices();
2163 for (int i = 0; i != local_indices.size(); ++i)
2164 if (local_indices[i] != -1)
2165 dotVector[i] = array[local_indices[i]];
2166 else
2167 dotVector[i] = 0;
2168 CHKERR VecRestoreArrayRead(dataVec, &array);
2169 data.getFieldData().swap(dotVector);
2170 }
2171
2172 const size_t nb_base_functions = data.getN().size2() / 3;
2173 FTensor::Index<'i', Field_Dim> i;
2174 auto t_n_hdiv = data.getFTensor1N<3>();
2175 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2176 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2177 auto t_dof = data.getFTensor0FieldData();
2178 int bb = 0;
2179 for (; bb != nb_dofs; ++bb) {
2180 t_data(i) += t_n_hdiv(i) * t_dof;
2181 ++t_n_hdiv;
2182 ++t_dof;
2183 }
2184 for (; bb != nb_base_functions; ++bb)
2185 ++t_n_hdiv;
2186 ++t_data;
2187 }
2188
2189 if (dataVec.use_count()) {
2190 data.getFieldData().swap(dotVector);
2191 }
2193}
2194
2195/** \brief Get vector field for H-div approximation
2196 *
2197 * \ingroup mofem_forces_and_sources_user_data_operators
2198 */
2199template <int Base_Dim, int Field_Dim = Base_Dim>
2202 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2204 Base_Dim, Field_Dim, double, ublas::row_major,
2205 DoubleAllocator>::OpCalculateHVecVectorField_General;
2206};
2207
2208/** \brief Get vector field for H-div approximation
2209 * \ingroup mofem_forces_and_sources_user_data_operators
2210 */
2211template <int Base_Dim, int Field_Dim = Base_Dim>
2213
2214template <int Field_Dim>
2217
2219 boost::shared_ptr<MatrixDouble> 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
2229 /**
2230 * \brief Calculate values of vector field at integration points
2231 * @param side side entity number
2232 * @param type side entity type
2233 * @param data entity data
2234 * @return error code
2235 */
2236 MoFEMErrorCode doWork(int side, EntityType type,
2238
2239private:
2240 boost::shared_ptr<MatrixDouble> dataPtr;
2242 const int zeroSide;
2243};
2244
2245template <int Field_Dim>
2247 int side, EntityType type, EntitiesFieldData::EntData &data) {
2249
2250 const size_t nb_integration_points = this->getGaussPts().size2();
2251 if (type == zeroType && side == zeroSide) {
2252 dataPtr->resize(Field_Dim, nb_integration_points, false);
2253 dataPtr->clear();
2254 }
2255
2256 auto &local_indices = data.getIndices();
2257 const size_t nb_dofs = local_indices.size();
2258 if (nb_dofs) {
2259
2260 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2261 const double *array;
2262 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2263 for (size_t i = 0; i != nb_dofs; ++i)
2264 if (local_indices[i] != -1)
2265 dot_dofs_vector[i] = array[local_indices[i]];
2266 else
2267 dot_dofs_vector[i] = 0;
2268 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2269
2270 const size_t nb_base_functions = data.getN().size2() / 3;
2271 FTensor::Index<'i', Field_Dim> i;
2272 auto t_n_hdiv = data.getFTensor1N<3>();
2273 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2274 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2275 int bb = 0;
2276 for (; bb != nb_dofs; ++bb) {
2277 t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
2278 ++t_n_hdiv;
2279 }
2280 for (; bb != nb_base_functions; ++bb)
2281 ++t_n_hdiv;
2282 ++t_data;
2283 }
2284 }
2285
2287}
2288
2289/**
2290 * @brief Calculate divergence of vector field
2291 * @ingroup mofem_forces_and_sources_user_data_operators
2292 *
2293 * @tparam BASE_DIM
2294 * @tparam SPACE_DIM
2295 */
2296template <int BASE_DIM, int SPACE_DIM>
2299
2301 boost::shared_ptr<VectorDouble> data_ptr,
2302 const EntityType zero_type = MBEDGE,
2303 const int zero_side = 0)
2306 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2307 if (!dataPtr)
2308 THROW_MESSAGE("Pointer is not set");
2309 }
2310
2311 MoFEMErrorCode doWork(int side, EntityType type,
2314 const size_t nb_integration_points = getGaussPts().size2();
2315 if (type == zeroType && side == zeroSide) {
2316 dataPtr->resize(nb_integration_points, false);
2317 dataPtr->clear();
2318 }
2319 const size_t nb_dofs = data.getFieldData().size();
2320 if (!nb_dofs)
2322 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2325 auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2326 auto t_data = getFTensor0FromVec(*dataPtr);
2327 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2328 auto t_dof = data.getFTensor0FieldData();
2329 int bb = 0;
2330 for (; bb != nb_dofs; ++bb) {
2331 t_data += t_dof * t_n_diff_hdiv(j, j);
2332 ++t_n_diff_hdiv;
2333 ++t_dof;
2334 }
2335 for (; bb != nb_base_functions; ++bb)
2336 ++t_n_diff_hdiv;
2337 ++t_data;
2338 }
2340 }
2341
2342private:
2343 boost::shared_ptr<VectorDouble> dataPtr;
2345 const int zeroSide;
2346};
2347
2348/**
2349 * @brief Calculate gradient of vector field
2350 * @ingroup mofem_forces_and_sources_user_data_operators
2351 *
2352 * @tparam BASE_DIM
2353 * @tparam SPACE_DIM
2354 */
2355template <int BASE_DIM, int SPACE_DIM>
2358
2360 boost::shared_ptr<MatrixDouble> data_ptr,
2361 const EntityType zero_type = MBEDGE,
2362 const int zero_side = 0)
2365 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2366 if (!dataPtr)
2367 THROW_MESSAGE("Pointer is not set");
2368 }
2369
2370 MoFEMErrorCode doWork(int side, EntityType type,
2373 const size_t nb_integration_points = getGaussPts().size2();
2374 if (type == zeroType && side == zeroSide) {
2375 dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
2376 dataPtr->clear();
2377 }
2378 const size_t nb_dofs = data.getFieldData().size();
2379 if (!nb_dofs)
2381 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2384 auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2386 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2387 auto t_dof = data.getFTensor0FieldData();
2388 int bb = 0;
2389 for (; bb != nb_dofs; ++bb) {
2390 t_data(i, j) += t_dof * t_base_diff(i, j);
2391 ++t_base_diff;
2392 ++t_dof;
2393 }
2394 for (; bb != nb_base_functions; ++bb)
2395 ++t_base_diff;
2396 ++t_data;
2397 }
2399 }
2400
2401private:
2402 boost::shared_ptr<MatrixDouble> dataPtr;
2404 const int zeroSide;
2405};
2406
2407/**
2408 * @brief Calculate gradient of vector field
2409 * @ingroup mofem_forces_and_sources_user_data_operators
2410 *
2411 * @tparam BASE_DIM
2412 * @tparam SPACE_DIM
2413 */
2414template <int BASE_DIM, int SPACE_DIM>
2417
2419 boost::shared_ptr<MatrixDouble> data_ptr,
2420 const EntityType zero_type = MBEDGE,
2421 const int zero_side = 0)
2424 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2425 if (!dataPtr)
2426 THROW_MESSAGE("Pointer is not set");
2427 }
2428
2429 MoFEMErrorCode doWork(int side, EntityType type,
2432 const size_t nb_integration_points = getGaussPts().size2();
2433 if (type == zeroType && side == zeroSide) {
2434 dataPtr->resize(BASE_DIM * SPACE_DIM * SPACE_DIM, nb_integration_points,
2435 false);
2436 dataPtr->clear();
2437 }
2438 const size_t nb_dofs = data.getFieldData().size();
2439 if (!nb_dofs)
2441
2442 const int nb_base_functions = data.getN().size2() / BASE_DIM;
2443
2444 #ifndef NDEBUG
2445 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2446 if (hessian_base.size1() != nb_integration_points) {
2447 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2448 "Wrong number of integration pts (%d != %d)",
2449 hessian_base.size1(), nb_integration_points);
2450 }
2451 if (hessian_base.size2() !=
2452 BASE_DIM * nb_base_functions * SPACE_DIM * SPACE_DIM) {
2453 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2454 "Wrong number of base functions (%d != %d)",
2455 hessian_base.size2() / (BASE_DIM * SPACE_DIM * SPACE_DIM),
2456 nb_base_functions);
2457 }
2458 if (hessian_base.size2() < BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM) {
2459 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2460 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2461 BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM);
2462 }
2463 #endif
2464
2468
2469 auto t_base_diff2 =
2472 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2473 auto t_dof = data.getFTensor0FieldData();
2474 int bb = 0;
2475 for (; bb != nb_dofs; ++bb) {
2476 t_data(i, j, k) += t_dof * t_base_diff2(i, j, k);
2477
2478 ++t_base_diff2;
2479 ++t_dof;
2480 }
2481 for (; bb != nb_base_functions; ++bb)
2482 ++t_base_diff2;
2483 ++t_data;
2484 }
2486 }
2487
2488private:
2489 boost::shared_ptr<MatrixDouble> dataPtr;
2491 const int zeroSide;
2492};
2493
2494/**
2495 * @brief Calculate divergence of vector field dot
2496 * @ingroup mofem_forces_and_sources_user_data_operators
2497 *
2498 * @tparam Tensor_Dim dimension of space
2499 */
2500template <int Tensor_Dim1, int Tensor_Dim2>
2503
2505 boost::shared_ptr<VectorDouble> data_ptr,
2506 const EntityType zero_type = MBEDGE,
2507 const int zero_side = 0)
2510 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2511 if (!dataPtr)
2512 THROW_MESSAGE("Pointer is not set");
2513 }
2514
2515 MoFEMErrorCode doWork(int side, EntityType type,
2518 const size_t nb_integration_points = getGaussPts().size2();
2519 if (type == zeroType && side == zeroSide) {
2520 dataPtr->resize(nb_integration_points, false);
2521 dataPtr->clear();
2522 }
2523
2524 const auto &local_indices = data.getLocalIndices();
2525 const int nb_dofs = local_indices.size();
2526 if (nb_dofs) {
2527
2528 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2529 const double *array;
2530 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2531 for (size_t i = 0; i != local_indices.size(); ++i)
2532 if (local_indices[i] != -1)
2533 dot_dofs_vector[i] = array[local_indices[i]];
2534 else
2535 dot_dofs_vector[i] = 0;
2536 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2537
2538 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
2539 FTensor::Index<'i', Tensor_Dim1> i;
2540 auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
2541 auto t_data = getFTensor0FromVec(*dataPtr);
2542 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2543 int bb = 0;
2544 for (; bb != nb_dofs; ++bb) {
2545 double div = 0;
2546 for (auto ii = 0; ii != Tensor_Dim2; ++ii)
2547 div += t_n_diff_hdiv(ii, ii);
2548 t_data += dot_dofs_vector[bb] * div;
2549 ++t_n_diff_hdiv;
2550 }
2551 for (; bb != nb_base_functions; ++bb)
2552 ++t_n_diff_hdiv;
2553 ++t_data;
2554 }
2555 }
2557 }
2558
2559private:
2560 boost::shared_ptr<VectorDouble> dataPtr;
2562 const int zeroSide;
2563};
2564
2565/**
2566 * @brief Calculate curl of vector field
2567 * @ingroup mofem_forces_and_sources_user_data_operators
2568 *
2569 * @tparam Base_Dim base function dimension
2570 * @tparam Space_Dim dimension of space
2571 * @tparam Hcurl field dimension
2572 */
2573template <int Base_Dim, int Space_Dim> struct OpCalculateHcurlVectorCurl;
2574
2575/**
2576 * @brief Calculate curl of vector field
2577 * @ingroup mofem_forces_and_sources_user_data_operators
2578 *
2579 * @tparam Base_Dim base function dimension
2580 * @tparam Space_Dim dimension of space
2581 * @tparam Hcurl field dimension
2582 */
2583template <>
2586 OpCalculateHcurlVectorCurl(const std::string field_name,
2587 boost::shared_ptr<MatrixDouble> data_ptr,
2588 const EntityType zero_type = MBEDGE,
2589 const int zero_side = 0);
2590 MoFEMErrorCode doWork(int side, EntityType type,
2592
2593private:
2594 boost::shared_ptr<MatrixDouble> dataPtr;
2596 const int zeroSide;
2597};
2598
2599/**
2600 * @brief Calculate curl of vector field
2601 * @ingroup mofem_forces_and_sources_user_data_operators
2602 *
2603 * @tparam Field_Dim dimension of field
2604 * @tparam Space_Dim dimension of space
2605 */
2606template <>
2609
2610 OpCalculateHcurlVectorCurl(const std::string field_name,
2611 boost::shared_ptr<MatrixDouble> data_ptr,
2612 const EntityType zero_type = MBVERTEX,
2613 const int zero_side = 0);
2614
2615 MoFEMErrorCode doWork(int side, EntityType type,
2617
2618private:
2619 boost::shared_ptr<MatrixDouble> dataPtr;
2621 const int zeroSide;
2622};
2623
2624/**
2625 * @brief Calculate curl of vector field
2626 * @ingroup mofem_forces_and_sources_user_data_operators
2627 *
2628 * @tparam Field_Dim dimension of field
2629 * @tparam Space_Dim dimension of space
2630 */
2631template <>
2634
2635 OpCalculateHcurlVectorCurl(const std::string field_name,
2636 boost::shared_ptr<MatrixDouble> data_ptr,
2637 const EntityType zero_type = MBVERTEX,
2638 const int zero_side = 0);
2639
2640 MoFEMErrorCode doWork(int side, EntityType type,
2642
2643private:
2644 boost::shared_ptr<MatrixDouble> dataPtr;
2646 const int zeroSide;
2647};
2648
2649/**
2650 * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
2651 * \ingroup mofem_forces_and_sources_user_data_operators
2652 *
2653 * @tparam Tensor_Dim0 rank of the field
2654 * @tparam Tensor_Dim1 dimension of space
2655 */
2656template <int Tensor_Dim0, int Tensor_Dim1>
2659
2661 boost::shared_ptr<MatrixDouble> data_ptr,
2662 boost::shared_ptr<double> scale_ptr,
2664 const EntityType zero_type = MBEDGE,
2665 const int zero_side = 0)
2668 dataPtr(data_ptr), scalePtr(scale_ptr), dataVec(data_vec),
2669 zeroType(zero_type), zeroSide(zero_side) {
2670 if (!dataPtr)
2671 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
2672 }
2673
2675 boost::shared_ptr<MatrixDouble> data_ptr,
2676 const EntityType zero_type = MBEDGE,
2677 const int zero_side = 0)
2678 : OpCalculateHVecTensorField(field_name, data_ptr, nullptr,
2679 SmartPetscObj<Vec>(), zero_type, zero_side) {
2680 }
2681
2682 MoFEMErrorCode doWork(int side, EntityType type,
2685 const size_t nb_integration_points = getGaussPts().size2();
2686 if (type == zeroType && side == zeroSide) {
2687 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2688 dataPtr->clear();
2689 }
2690 const size_t nb_dofs = data.getFieldData().size();
2691 if (nb_dofs) {
2692
2693 if (dataVec.use_count()) {
2694 dotVector.resize(nb_dofs, false);
2695 const double *array;
2696 CHKERR VecGetArrayRead(dataVec, &array);
2697 const auto &local_indices = data.getLocalIndices();
2698 for (int i = 0; i != local_indices.size(); ++i)
2699 if (local_indices[i] != -1)
2700 dotVector[i] = array[local_indices[i]];
2701 else
2702 dotVector[i] = 0;
2703 CHKERR VecRestoreArrayRead(dataVec, &array);
2704 data.getFieldData().swap(dotVector);
2705 }
2706
2707 double scale = (scalePtr) ? *scalePtr : 1.0;
2708 const size_t nb_base_functions = data.getN().size2() / 3;
2709 FTensor::Index<'i', Tensor_Dim0> i;
2710 FTensor::Index<'j', Tensor_Dim1> j;
2711 auto t_n_hvec = data.getFTensor1N<3>();
2713 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2714 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2715 size_t bb = 0;
2716 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2717 t_data(i, j) += (scale * t_dof(i)) * t_n_hvec(j);
2718 ++t_n_hvec;
2719 ++t_dof;
2720 }
2721 for (; bb < nb_base_functions; ++bb)
2722 ++t_n_hvec;
2723 ++t_data;
2724 }
2725
2726 if (dataVec.use_count()) {
2727 data.getFieldData().swap(dotVector);
2728 }
2729 }
2731 }
2732
2733private:
2734 boost::shared_ptr<MatrixDouble> dataPtr;
2735 boost::shared_ptr<double> scalePtr;
2738 const int zeroSide;
2739 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
2740};
2741
2742/** \brief Get tensor field for H-div approximation
2743 * \ingroup mofem_forces_and_sources_user_data_operators
2744 *
2745 * \warning This operator is not tested
2746 */
2747template <int Tensor_Dim0, int Tensor_Dim1>
2750
2752
2754 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
2755 SmartPetscObj<Vec> data_vec, EntityType broken_type,
2756 boost::shared_ptr<Range> broken_range_ptr = nullptr,
2757 boost::shared_ptr<double> scale_ptr = nullptr,
2758 const EntityType zero_type = MBEDGE, const int zero_side = 0)
2760 dataPtr(data_ptr), dataVec(data_vec), brokenType(broken_type),
2761 brokenRangePtr(broken_range_ptr), zeroType(zero_type) {
2762 if (!dataPtr)
2763 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
2764 }
2765
2766 /**
2767 * \brief Calculate values of vector field at integration points
2768 * @param side side entity number
2769 * @param type side entity type
2770 * @param data entity data
2771 * @return error code
2772 */
2773 MoFEMErrorCode doWork(int side, EntityType type,
2775
2776private:
2777 boost::shared_ptr<MatrixDouble> dataPtr;
2779 EntityType brokenType;
2780 boost::shared_ptr<Range> brokenRangePtr;
2781 boost::shared_ptr<double> scalePtr;
2784};
2785
2786template <int Tensor_Dim0, int Tensor_Dim1>
2789 int side, EntityType type, EntitiesFieldData::EntData &data) {
2791 const size_t nb_integration_points = OP::getGaussPts().size2();
2792 if (type == zeroType) {
2793 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2794 dataPtr->clear();
2795 }
2796 const size_t nb_dofs = data.getFieldData().size();
2797 if (!nb_dofs)
2799
2800 if (dataVec.use_count()) {
2801 dotVector.resize(nb_dofs, false);
2802 const double *array;
2803 CHKERR VecGetArrayRead(dataVec, &array);
2804 const auto &local_indices = data.getLocalIndices();
2805 for (int i = 0; i != local_indices.size(); ++i)
2806 if (local_indices[i] != -1)
2807 dotVector[i] = array[local_indices[i]];
2808 else
2809 dotVector[i] = 0;
2810 CHKERR VecRestoreArrayRead(dataVec, &array);
2811 data.getFieldData().swap(dotVector);
2812 }
2813
2814 /**
2815 * @brief Get side face dofs
2816 *
2817 * Find which base functions on borken space have adjacent given entity type
2818 * and are in the range ptr if given.
2819 *
2820 */
2821 auto get_get_side_face_dofs = [&]() {
2822 auto fe_type = OP::getFEType();
2823
2824 BaseFunction::DofsSideMap &side_dof_map =
2825 data.getFieldEntities()[0]->getDofSideMap().at(fe_type);
2826 std::vector<int> side_face_dofs;
2827 side_face_dofs.reserve(data.getIndices().size() / Tensor_Dim0);
2828
2829 for (
2830
2831 auto it = side_dof_map.get<1>().begin();
2832 it != side_dof_map.get<1>().end(); ++it
2833
2834 ) {
2835 if ((Tensor_Dim0 * it->dof) >= data.getIndices().size()) {
2836 break;
2837 }
2838 if (it->type == brokenType) {
2839 if (brokenRangePtr) {
2840 auto ent = OP::getSideEntity(it->side, brokenType);
2841 if (brokenRangePtr->find(ent) != brokenRangePtr->end()) {
2842 side_face_dofs.push_back(it->dof);
2843 }
2844 } else {
2845 side_face_dofs.push_back(it->dof);
2846 }
2847 }
2848 }
2849
2850 return side_face_dofs;
2851 };
2852
2853 auto side_face_dofs = get_get_side_face_dofs();
2854
2855 FTensor::Index<'i', Tensor_Dim0> i;
2856 FTensor::Index<'j', Tensor_Dim1> j;
2857 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2858 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2859 for (auto b : side_face_dofs) {
2860 auto t_row_base = data.getFTensor1N<3>(gg, b);
2861 auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(data.getFieldData().data() +
2862 b * Tensor_Dim0);
2863 t_data(i, j) += t_dof(i) * t_row_base(j);
2864 }
2865 ++t_data;
2866 }
2867 *dataPtr *= (scalePtr) ? *scalePtr : 1.0;
2868
2869 if (dataVec.use_count()) {
2870 data.getFieldData().swap(dotVector);
2871 }
2872
2874}
2875
2876/**
2877 * @brief Calculate tenor field using tensor base, i.e. Hdiv/Hcurl
2878 * \ingroup mofem_forces_and_sources_user_data_operators
2879 *
2880 * @tparam Tensor_Dim0 rank of the field
2881 * @tparam Tensor_Dim1 dimension of space
2882 */
2883template <int Tensor_Dim0, int Tensor_Dim1>
2886
2888 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
2889 boost::shared_ptr<double> scale_ptr,
2891 const EntityType zero_type = MBEDGE, const int zero_side = 0)
2894 dataPtr(data_ptr), scalePtr(scale_ptr), dataVec(data_vec),
2895 zeroType(zero_type), zeroSide(zero_side) {
2896 if (!dataPtr)
2897 THROW_MESSAGE("Pointer is not set");
2898 }
2899
2901 boost::shared_ptr<MatrixDouble> data_ptr,
2902 const EntityType zero_type = MBEDGE,
2903 const int zero_side = 0)
2904 : OpCalculateHTensorTensorField(field_name, data_ptr, nullptr,
2905 SmartPetscObj<Vec>(), zero_type,
2906 zero_side) {}
2907
2908 MoFEMErrorCode doWork(int side, EntityType type,
2911 const size_t nb_integration_points = getGaussPts().size2();
2912 if (type == zeroType && side == zeroSide) {
2913 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2914 dataPtr->clear();
2915 }
2916 const size_t nb_dofs = data.getFieldData().size();
2917 if (!nb_dofs)
2919
2920 if (dataVec.use_count()) {
2921 dotVector.resize(nb_dofs, false);
2922 const double *array;
2923 CHKERR VecGetArrayRead(dataVec, &array);
2924 const auto &local_indices = data.getLocalIndices();
2925 for (int i = 0; i != local_indices.size(); ++i)
2926 if (local_indices[i] != -1)
2927 dotVector[i] = array[local_indices[i]];
2928 else
2929 dotVector[i] = 0;
2930 CHKERR VecRestoreArrayRead(dataVec, &array);
2931 data.getFieldData().swap(dotVector);
2932 }
2933
2934 double scale = (scalePtr) ? *scalePtr : 1.0;
2935 const size_t nb_base_functions =
2936 data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2937 FTensor::Index<'i', Tensor_Dim0> i;
2938 FTensor::Index<'j', Tensor_Dim1> j;
2939 auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2941 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2942 auto t_dof = data.getFTensor0FieldData();
2943 size_t bb = 0;
2944 for (; bb != nb_dofs; ++bb) {
2945 t_data(i, j) += (scale * t_dof) * t_n_hten(i, j);
2946 ++t_n_hten;
2947 ++t_dof;
2948 }
2949 for (; bb < nb_base_functions; ++bb)
2950 ++t_n_hten;
2951 ++t_data;
2952 }
2953
2954 if (dataVec.use_count()) {
2955 data.getFieldData().swap(dotVector);
2956 }
2957
2959 }
2960
2961private:
2962 boost::shared_ptr<MatrixDouble> dataPtr;
2963 boost::shared_ptr<double> scalePtr;
2966 const int zeroSide;
2967 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
2968};
2969
2970/**
2971 * @brief Calculate divergence of tonsorial field using vectorial base
2972 * \ingroup mofem_forces_and_sources_user_data_operators
2973 *
2974 * @tparam Tensor_Dim0 rank of the field
2975 * @tparam Tensor_Dim1 dimension of space
2976 */
2977template <int Tensor_Dim0, int Tensor_Dim1,
2978 CoordinateTypes CoordSys = CARTESIAN>
2981
2983 boost::shared_ptr<MatrixDouble> data_ptr,
2984 SmartPetscObj<Vec> data_vec,
2985 const EntityType zero_type = MBEDGE,
2986 const int zero_side = 0)
2989 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type),
2990 zeroSide(zero_side) {
2991 if (!dataPtr)
2992 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
2993 }
2994
2996 boost::shared_ptr<MatrixDouble> data_ptr,
2997 const EntityType zero_type = MBEDGE,
2998 const int zero_side = 0)
3000 field_name, data_ptr, SmartPetscObj<Vec>(), zero_type, zero_side) {}
3001
3002 MoFEMErrorCode doWork(int side, EntityType type,
3005 const size_t nb_integration_points = getGaussPts().size2();
3006 if (type == zeroType && side == zeroSide) {
3007 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
3008 dataPtr->clear();
3009 }
3010 const size_t nb_dofs = data.getFieldData().size();
3011 if (nb_dofs) {
3012
3013 if (dataVec.use_count()) {
3014 dotVector.resize(nb_dofs, false);
3015 const double *array;
3016 CHKERR VecGetArrayRead(dataVec, &array);
3017 const auto &local_indices = data.getLocalIndices();
3018 for (int i = 0; i != local_indices.size(); ++i)
3019 if (local_indices[i] != -1)
3020 dotVector[i] = array[local_indices[i]];
3021 else
3022 dotVector[i] = 0;
3023 CHKERR VecRestoreArrayRead(dataVec, &array);
3024 data.getFieldData().swap(dotVector);
3025 }
3026
3027 const size_t nb_base_functions = data.getN().size2() / 3;
3028 FTensor::Index<'i', Tensor_Dim0> i;
3029 FTensor::Index<'j', Tensor_Dim1> j;
3030 auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
3032 auto t_base = data.getFTensor1N<3>();
3033 auto t_coords = getFTensor1CoordsAtGaussPts();
3034 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3035 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
3036 size_t bb = 0;
3037 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
3038 double div = t_n_diff_hvec(j, j);
3039 t_data(i) += t_dof(i) * div;
3040 if constexpr (CoordSys == CYLINDRICAL) {
3041 t_data(i) += t_base(0) * (t_dof(i) / t_coords(0));
3042 }
3043 ++t_n_diff_hvec;
3044 ++t_dof;
3045 ++t_base;
3046 }
3047 for (; bb < nb_base_functions; ++bb) {
3048 ++t_base;
3049 ++t_n_diff_hvec;
3050 }
3051 ++t_data;
3052 ++t_coords;
3053 }
3054
3055 if (dataVec.use_count()) {
3056 data.getFieldData().swap(dotVector);
3057 }
3058 }
3060 }
3061
3062private:
3063 boost::shared_ptr<MatrixDouble> dataPtr;
3066 const int zeroSide;
3067
3068 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
3069};
3070
3071/**
3072 * @brief Calculate divergence of tonsorial field using vectorial base
3073 * \ingroup mofem_forces_and_sources_user_data_operators
3074 *
3075 * \warning This operator is not tested
3076 *
3077 * @tparam Tensor_Dim0 rank of the field
3078 * @tparam Tensor_Dim1 dimension of space
3079 */
3080template <int Tensor_Dim0, int Tensor_Dim1,
3081 CoordinateTypes CoordSys = CARTESIAN>
3084
3086
3088 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3089 SmartPetscObj<Vec> data_vec, EntityType broken_type,
3090 boost::shared_ptr<Range> broken_range_ptr = nullptr,
3091 boost::shared_ptr<double> scale_ptr = nullptr,
3092 const EntityType zero_type = MBEDGE)
3094 dataPtr(data_ptr), dataVec(data_vec), brokenType(broken_type),
3095 brokenRangePtr(broken_range_ptr), scalePtr(scale_ptr),
3096 zeroType(zero_type) {
3097 if (!dataPtr)
3098 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
3099 }
3100
3101 MoFEMErrorCode doWork(int side, EntityType type,
3104 const size_t nb_integration_points = getGaussPts().size2();
3105 if (type == zeroType && side == 0) {
3106 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
3107 dataPtr->clear();
3108 }
3109
3110 const size_t nb_dofs = data.getFieldData().size();
3111 if (nb_dofs) {
3112
3113 if (dataVec.use_count()) {
3114 dotVector.resize(nb_dofs, false);
3115 const double *array;
3116 CHKERR VecGetArrayRead(dataVec, &array);
3117 const auto &local_indices = data.getLocalIndices();
3118 for (int i = 0; i != local_indices.size(); ++i)
3119 if (local_indices[i] != -1)
3120 dotVector[i] = array[local_indices[i]];
3121 else
3122 dotVector[i] = 0;
3123 CHKERR VecRestoreArrayRead(dataVec, &array);
3124 data.getFieldData().swap(dotVector);
3125 }
3126
3127 /**
3128 * @brief Get side face dofs
3129 *
3130 * Find which base functions on borken space have adjacent given entity
3131 * type and are in the range ptr if given.
3132 *
3133 */
3134 auto get_get_side_face_dofs = [&]() {
3135 auto fe_type = OP::getFEType();
3136
3137 BaseFunction::DofsSideMap &side_dof_map =
3138 data.getFieldEntities()[0]->getDofSideMap().at(fe_type);
3139 std::vector<int> side_face_dofs;
3140 side_face_dofs.reserve(data.getIndices().size() / Tensor_Dim0);
3141
3142 for (
3143
3144 auto it = side_dof_map.get<1>().begin();
3145 it != side_dof_map.get<1>().end(); ++it
3146
3147 ) {
3148 if ((Tensor_Dim0 * it->dof) >= data.getIndices().size()) {
3149 break;
3150 }
3151 if (it->type == brokenType) {
3152 if (brokenRangePtr) {
3153 auto ent = OP::getSideEntity(it->side, brokenType);
3154 if (brokenRangePtr->find(ent) != brokenRangePtr->end()) {
3155 side_face_dofs.push_back(it->dof);
3156 }
3157 } else {
3158 side_face_dofs.push_back(it->dof);
3159 }
3160 }
3161 }
3162
3163 return side_face_dofs;
3164 };
3165
3166 auto side_face_dofs = get_get_side_face_dofs();
3167
3168 FTensor::Index<'i', Tensor_Dim0> i;
3169 FTensor::Index<'j', Tensor_Dim1> j;
3171 auto t_coords = getFTensor1CoordsAtGaussPts();
3172 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3173 for (auto b : side_face_dofs) {
3175 data.getFieldData().data() + b * Tensor_Dim0);
3176 auto t_base = data.getFTensor1N<3>(gg, b);
3177 auto t_diff_base = data.getFTensor2DiffN<3, Tensor_Dim1>(gg, b);
3178 double div = t_diff_base(j, j);
3179 t_data(i) += t_dof(i) * div;
3180 if constexpr (CoordSys == CYLINDRICAL) {
3181 t_data(i) += t_base(0) * (t_dof(i) / t_coords(0));
3182 }
3183 }
3184 ++t_data;
3185 ++t_coords;
3186 }
3187 }
3188
3189 if (dataVec.use_count()) {
3190 data.getFieldData().swap(dotVector);
3191 }
3192
3194 }
3195
3196private:
3197 boost::shared_ptr<MatrixDouble> dataPtr;
3199 EntityType brokenType;
3200 boost::shared_ptr<Range> brokenRangePtr;
3201 boost::shared_ptr<double> scalePtr;
3204};
3205
3206/**
3207 * @brief Calculate trace of vector (Hdiv/Hcurl) space
3208 *
3209 * @tparam Tensor_Dim
3210 * @tparam OpBase
3211 */
3212template <int Tensor_Dim, typename OpBase>
3214
3216 boost::shared_ptr<MatrixDouble> data_ptr,
3217 boost::shared_ptr<double> scale_ptr,
3218 const EntityType zero_type = MBEDGE,
3219 const int zero_side = 0)
3220 : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
3221 scalePtr(scale_ptr), zeroType(zero_type), zeroSide(zero_side) {
3222 if (!dataPtr)
3223 THROW_MESSAGE("Pointer is not set");
3224 }
3225
3227 boost::shared_ptr<MatrixDouble> data_ptr,
3228 const EntityType zero_type = MBEDGE,
3229 const int zero_side = 0)
3230 : OpCalculateHVecTensorTrace(field_name, data_ptr, nullptr, zero_type,
3231 zero_side) {}
3232
3233 MoFEMErrorCode doWork(int side, EntityType type,
3236 const size_t nb_integration_points = OpBase::getGaussPts().size2();
3237 if (type == zeroType && side == 0) {
3238 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
3239 dataPtr->clear();
3240 }
3241 const size_t nb_dofs = data.getFieldData().size();
3242 if (nb_dofs) {
3243 double scale_val = (scalePtr) ? *scalePtr : 1.0;
3244 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
3245 const size_t nb_base_functions = data.getN().size2() / 3;
3246 auto t_base = data.getFTensor1N<3>();
3248 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3249 FTensor::Tensor1<double, Tensor_Dim> t_normalized_normal;
3250 t_normalized_normal(j) = t_normal(j);
3251 t_normalized_normal.normalize();
3252 auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
3253 size_t bb = 0;
3254 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
3255 t_data(i) +=
3256 (scale_val * t_dof(i)) * (t_base(j) * t_normalized_normal(j));
3257 ++t_base;
3258 ++t_dof;
3259 }
3260 for (; bb < nb_base_functions; ++bb) {
3261 ++t_base;
3262 }
3263 ++t_data;
3264 ++t_normal;
3265 }
3266 }
3268 }
3269
3270private:
3271 boost::shared_ptr<MatrixDouble> dataPtr;
3272 boost::shared_ptr<double> scalePtr;
3274 const int zeroSide;
3275 FTensor::Index<'i', Tensor_Dim> i;
3276 FTensor::Index<'j', Tensor_Dim> j;
3277};
3278
3279/**@}*/
3280
3281/** \name Other operators */
3282
3283/**@{*/
3284
3285/**@}*/
3286
3287/** \name Operators for faces */
3288
3289/**@{*/
3290
3291/** \brief Transform local reference derivatives of shape functions to global
3292derivatives
3293
3294\ingroup mofem_forces_and_sources_tri_element
3295
3296*/
3297template <int DIM, int DERIVATIVE = 1> struct OpSetInvJacSpaceForFaceImpl;
3298
3301
3303 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3304
3305protected:
3306
3307 /**
3308 * @brief Apply transformation to the input matrix
3309 *
3310 * @tparam D1 dimension of the directive of base functions in input
3311 * @tparam D2 dimension of the directive of base functions in output
3312 * @tparam J1 nb of rows in jacobian (= dimension of space)
3313 * @tparam J2 nb of columns in jacobian (= dimension of reference element)
3314 * @param diff_n
3315 * @return MoFEMErrorCode
3316 */
3317 template <int D1, int D2, int J1, int J2>
3320
3321 static_assert(D2 == J2, "Dimension of jacobian and dimension of <out> "
3322 "directive does not match");
3323
3324 size_t nb_functions = diff_n.size2() / D1;
3325 if (nb_functions) {
3326 size_t nb_gauss_pts = diff_n.size1();
3327
3328 #ifndef NDEBUG
3329 if (nb_gauss_pts != getGaussPts().size2())
3330 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3331 "Wrong number of Gauss Pts");
3332 if (diff_n.size2() % D1)
3333 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3334 "Number of directives of base functions and D1 dimension does "
3335 "not match");
3336 #endif
3337
3338 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions, false);
3339
3340 FTensor::Index<'i', D2> i;
3341 FTensor::Index<'k', D1> k;
3342 auto t_diff_n = getFTensor1FromPtr<D2>(&*diffNinvJac.data().begin());
3343 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
3344 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*invJacPtr);
3345 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
3346 for (size_t dd = 0; dd != nb_functions; ++dd) {
3347 t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
3348 ++t_diff_n;
3349 ++t_diff_n_ref;
3350 }
3351 }
3352
3353 diff_n.swap(diffNinvJac);
3354 }
3356 }
3357
3358 boost::shared_ptr<MatrixDouble> invJacPtr;
3360};
3361
3362template <>
3365
3367
3368 MoFEMErrorCode doWork(int side, EntityType type,
3370};
3371
3372template <>
3374 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
3375
3376 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
3377
3378 MoFEMErrorCode doWork(int side, EntityType type,
3380};
3381
3382template <>
3384 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
3385
3386 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
3387
3388 MoFEMErrorCode doWork(int side, EntityType type,
3390};
3391
3392template <int DERIVARIVE = 1>
3394 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
3395 OpSetInvJacH1ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3396 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(H1, inv_jac_ptr) {}
3397};
3398
3399template <int DERIVARIVE = 1>
3401 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
3402 OpSetInvJacL2ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3403 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(L2, inv_jac_ptr) {}
3404};
3405
3406template <int DERIVARIVE = 1>
3408 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
3410 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3411 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(H1, inv_jac_ptr) {}
3412};
3413
3414template <int DERIVARIVE = 1>
3416 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
3418 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3419 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(L2, inv_jac_ptr) {}
3420};
3421
3422/**
3423 * \brief Transform local reference derivatives of shape function to
3424 global derivatives for face
3425
3426 * \ingroup mofem_forces_and_sources_tri_element
3427 */
3428template <int DIM> struct OpSetInvJacHcurlFaceImpl;
3429
3430template <>
3433
3434 OpSetInvJacHcurlFaceImpl(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3436 invJacPtr(inv_jac_ptr) {}
3437
3438 MoFEMErrorCode doWork(int side, EntityType type,
3440
3441protected:
3442 boost::shared_ptr<MatrixDouble> invJacPtr;
3444};
3445
3446template <>
3448 using OpSetInvJacHcurlFaceImpl<2>::OpSetInvJacHcurlFaceImpl;
3449 MoFEMErrorCode doWork(int side, EntityType type,
3451};
3452
3455
3456/**
3457 * @brief Make Hdiv space from Hcurl space in 2d
3458 * @ingroup mofem_forces_and_sources_tri_element
3459 */
3469
3470/** \brief Transform Hcurl base fluxes from reference element to physical
3471 * triangle
3472 */
3474
3475/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3476 *
3477 * Covariant Piola transformation
3478 \f[
3479 \psi_i|_t = J^{-1}_{ij}\hat{\psi}_j\\
3480 \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3481 = J^{-1}_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3482 \f]
3483
3484 */
3485template <>
3488
3490 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3491
3492 MoFEMErrorCode doWork(int side, EntityType type,
3494
3495private:
3496 boost::shared_ptr<MatrixDouble> invJacPtr;
3497
3500};
3501
3504
3505/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3506 *
3507 * \note Hdiv space is generated by Hcurl space in 2d.
3508 *
3509 * Contravariant Piola transformation
3510 * \f[
3511 * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
3512 * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3513 * =
3514 * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3515 * \f]
3516 *
3517 * \ingroup mofem_forces_and_sources
3518 *
3519 */
3521
3522template <>
3525
3527 boost::shared_ptr<MatrixDouble> jac_ptr)
3529 jacPtr(jac_ptr) {}
3530
3531 MoFEMErrorCode doWork(int side, EntityType type,
3533
3534protected:
3535 boost::shared_ptr<MatrixDouble> jacPtr;
3538};
3539
3540template <>
3544 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3545
3546 MoFEMErrorCode doWork(int side, EntityType type,
3548};
3549
3554
3555/**@}*/
3556
3557/** \name Operators for edges */
3558
3559/**@{*/
3560
3570
3571/**
3572 * @deprecated Name is deprecated and this is added for back compatibility
3573 */
3576
3577/**@}*/
3578
3579/** \name Operator for fat prisms */
3580
3581/**@{*/
3582
3583/**
3584 * @brief Operator for fat prism element updating integration weights in the
3585 * volume.
3586 *
3587 * Jacobian on the distorted element is nonconstant. This operator updates
3588 * integration weight on prism to take into account nonconstat jacobian.
3589 *
3590 * \f[
3591 * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
3592 * \pmb\xi} \right\| \right)
3593 * \f]
3594 * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
3595 * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
3596 * element coordinate.
3597 *
3598 */
3608
3609/** \brief Calculate inverse of jacobian for face element
3610
3611 It is assumed that face element is XY plane. Applied
3612 only for 2d problems.
3613
3614 FIXME Generalize function for arbitrary face orientation in 3d space
3615 FIXME Calculate to Jacobins for two faces
3616
3617 \ingroup mofem_forces_and_sources_prism_element
3618
3619*/
3622
3623 OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3625 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3626
3630 MoFEMErrorCode doWork(int side, EntityType type,
3632
3633private:
3634 const boost::shared_ptr<MatrixDouble> invJacPtr;
3636};
3637
3638/** \brief Transform local reference derivatives of shape functions to global
3639derivatives
3640
3641FIXME Generalize to curved shapes
3642FIXME Generalize to case that top and bottom face has different shape
3643
3644\ingroup mofem_forces_and_sources_prism_element
3645
3646*/
3649
3650 OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3652 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3653
3657
3658 MoFEMErrorCode doWork(int side, EntityType type,
3660
3661private:
3662 const boost::shared_ptr<MatrixDouble> invJacPtr;
3665};
3666
3667// Flat prism
3668
3669/** \brief Calculate inverse of jacobian for face element
3670
3671 It is assumed that face element is XY plane. Applied
3672 only for 2d problems.
3673
3674 FIXME Generalize function for arbitrary face orientation in 3d space
3675 FIXME Calculate to Jacobins for two faces
3676
3677 \ingroup mofem_forces_and_sources_prism_element
3678
3679*/
3692
3693/** \brief Transform local reference derivatives of shape functions to global
3694derivatives
3695
3696FIXME Generalize to curved shapes
3697FIXME Generalize to case that top and bottom face has different shape
3698
3699\ingroup mofem_forces_and_sources_prism_element
3700
3701*/
3716
3717/**@}*/
3718
3719/** \name Operation on matrices at integration points */
3720
3721/**@{*/
3722
3723template <int DIM>
3725
3726 OpInvertMatrix(boost::shared_ptr<MatrixDouble> in_ptr,
3727 boost::shared_ptr<VectorDouble> det_ptr,
3728 boost::shared_ptr<MatrixDouble> out_ptr)
3730 outPtr(out_ptr), detPtr(det_ptr) {}
3731
3732 MoFEMErrorCode doWork(int side, EntityType type,
3734
3735private:
3736 boost::shared_ptr<MatrixDouble> inPtr;
3737 boost::shared_ptr<MatrixDouble> outPtr;
3738 boost::shared_ptr<VectorDouble> detPtr;
3739};
3740
3741template <int DIM>
3745
3746 if (!inPtr)
3747 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3748 "Pointer for inPtr matrix not allocated");
3749 if (!detPtr)
3750 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3751 "Pointer for detPtr matrix not allocated");
3752
3753 const auto nb_rows = inPtr->size1();
3754 const auto nb_integration_pts = inPtr->size2();
3755
3756 // Calculate determinant
3757 {
3758 detPtr->resize(nb_integration_pts, false);
3759 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3760 auto t_det = getFTensor0FromVec(*detPtr);
3761 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3762 t_det = determinantTensor(t_in);
3763 ++t_in;
3764 ++t_det;
3765 }
3766 }
3767
3768 // Invert jacobian
3769 if (outPtr) {
3770 outPtr->resize(nb_rows, nb_integration_pts, false);
3771 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3772 auto t_out = getFTensor2FromMat<DIM, DIM>(*outPtr);
3773 auto t_det = getFTensor0FromVec(*detPtr);
3774 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3775 CHKERR invertTensor(t_in, t_det, t_out);
3776 ++t_in;
3777 ++t_out;
3778 ++t_det;
3779 }
3780 }
3781
3783}
3784
3785/**@}*/
3786
3787/** \brief Calculates the trace of an input matrix
3788
3789\ingroup mofem_forces_and_sources
3790
3791*/
3792
3793template <int DIM>
3795
3796 OpCalculateTraceFromMat(boost::shared_ptr<MatrixDouble> in_ptr,
3797 boost::shared_ptr<VectorDouble> out_ptr)
3799 outPtr(out_ptr) {}
3800
3801 MoFEMErrorCode doWork(int side, EntityType type,
3803
3804private:
3806 boost::shared_ptr<MatrixDouble> inPtr;
3807 boost::shared_ptr<VectorDouble> outPtr;
3808};
3809
3810template <int DIM>
3815
3816 if (!inPtr)
3817 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3818 "Pointer for inPtr matrix not allocated");
3819
3820 const auto nb_integration_pts = inPtr->size2();
3821 // Invert jacobian
3822 if (outPtr) {
3823 outPtr->resize(nb_integration_pts, false);
3824 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3825 auto t_out = getFTensor0FromVec(*outPtr);
3826
3827 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3828 t_out = t_in(i, i);
3829 ++t_in;
3830 ++t_out;
3831 }
3832 }
3833
3835}
3836
3837/**@}*/
3838
3839/** \brief Calculates the trace of an input matrix
3840
3841\ingroup mofem_forces_and_sources
3842
3843*/
3844
3845template <int DIM>
3848
3849 OpCalculateTraceFromSymmMat(boost::shared_ptr<MatrixDouble> in_ptr,
3850 boost::shared_ptr<VectorDouble> out_ptr)
3852 outPtr(out_ptr) {}
3853
3854 MoFEMErrorCode doWork(int side, EntityType type,
3856
3857private:
3859 boost::shared_ptr<MatrixDouble> inPtr;
3860 boost::shared_ptr<VectorDouble> outPtr;
3861};
3862
3863template <int DIM>
3868
3869 if (!inPtr)
3870 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3871 "Pointer for inPtr matrix not allocated");
3872
3873 const auto nb_integration_pts = inPtr->size2();
3874 // Invert jacobian
3875 if (outPtr) {
3876 outPtr->resize(nb_integration_pts, false);
3877 auto t_in = getFTensor2SymmetricFromMat<DIM>(*inPtr);
3878 auto t_out = getFTensor0FromVec(*outPtr);
3879
3880 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3881 t_out = t_in(i, i);
3882 ++t_in;
3883 ++t_out;
3884 }
3885 }
3886
3888}
3889
3890} // namespace MoFEM
3891
3892#endif // __USER_DATA_OPERATORS_HPP__
3893
3894/**
3895 * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
3896 *
3897 * \brief Classes and functions used to evaluate fields at integration pts,
3898 *jacobians, etc..
3899 *
3900 * \ingroup mofem_forces_and_sources
3901 **/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
Tensor1< T, Tensor_Dim > normalize()
static const char *const CoordinateTypesNames[]
Coordinate system names.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
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 ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
CoordinateTypes
Coodinate system.
@ CYLINDRICAL
@ POLAR
@ CARTESIAN
@ SPHERICAL
#define DEPRECATED
Definition definitions.h:17
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
constexpr int BASE_DIM
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VecAllocator< double > DoubleAllocator
Definition Types.hpp:62
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static FTensor::Tensor3< FTensor::PackPtr< T *, 1 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 3 (non symmetries) form data matrix.
static FTensor::Dg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim2 > getFTensor3DgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 3 on the first two indices from form data matrix.
auto getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
constexpr IntegrationType I
double scale
Definition plastic.cpp:123
constexpr auto field_name
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
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
const VectorFieldEntities & getFieldEntities() const
get field entities
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 field data coefficients.
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 field data coefficients.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from field data coefficients.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
EntityType getFEType() const
Get dimension of finite element.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
EntityHandle getSideEntity(const int side_number, const EntityType type)
Get the side entity.
@ 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 divergence of tonsorial field using vectorial base.
OpCalculateBrokenHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, EntityType broken_type, boost::shared_ptr< Range > broken_range_ptr=nullptr, boost::shared_ptr< double > scale_ptr=nullptr, const EntityType zero_type=MBEDGE)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get tensor field for H-div approximation.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateBrokenHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, EntityType broken_type, boost::shared_ptr< Range > broken_range_ptr=nullptr, boost::shared_ptr< double > scale_ptr=nullptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate values of vector field at integration points.
Calculate field values (template specialization) for tensor field rank 1, i.e. vector field.
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.
VectorDouble dotVector
Keeps temporary values of time derivatives.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, SmartPetscObj< Vec > data_vec=SmartPetscObj< Vec >(), const EntityType zero_type=MBEDGE, const int zero_side=0)
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
VectorDouble dotVector
Keeps temporary values of time derivatives.
OpCalculateHVecTensorDivergence(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)
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.
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, SmartPetscObj< Vec > data_vec=SmartPetscObj< Vec >(), const EntityType zero_type=MBEDGE, const int zero_side=0)
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)
VectorDouble dotVector
Keeps temporary values of time derivatives.
boost::shared_ptr< double > scalePtr
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
boost::shared_ptr< double > scalePtr
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
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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.
OpCalculateHVecVectorField_General(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)
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.
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.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate gradient values of scalar field at integration points
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector)
Get field gradients at integration pts for scalar field 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
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)
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 value at integration points for scalar field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar field at integration points
Get time direvarive values at integration pts for tensor field 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 temporary values of time derivatives.
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
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 values at integration pts for tensor field rank 2, i.e. matrix field.
OpCalculateTensor2SymmetricFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
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
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > outPtr
boost::shared_ptr< MatrixDouble > inPtr
OpCalculateTraceFromSymmMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Get field gradients time derivative at integration pts for scalar field 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 temporary values of time derivatives.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
Evaluate field gradient values for vector field, i.e. gradient is tensor rank 2.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
Approximate field values for given petsc vector.
OpCalculateVectorFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX, bool throw_error=true)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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
Get values at integration pts for tensor field 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.
boost::shared_ptr< MatrixDouble > outMat
DEPRECATED OpScaleMatrix(const std::string field_name, const double scale, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< double > scalePtr
OpScaleMatrix(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
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)
Apply transformation to the input matrix.
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
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > dMat
OpTensorTimesSymmetricTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
boost::shared_ptr< MatrixDouble > inMat
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
static constexpr Switches CtxSetX_TT
std::bitset< 8 > Switches
static constexpr Switches CtxSetX_T
intrusive_ptr for managing petsc objects