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