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  "Nb. of DOFs is inconsistent with Tensor_Dim");
418  }
419  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
420  auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
421 
422  #ifndef NDEBUG
423  if (field_data.l2() != field_data.l2()) {
424  MOFEM_LOG("SELF", Sev::error) << "field data: " << field_data;
425  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
426  "Wrong number in coefficients");
427  }
428  #endif
429 
430  size_t bb = 0;
431  for (; bb != size; ++bb) {
432 
433  #ifndef SINGULARITY
434  #ifndef NDEBUG
435  if (base_function != base_function) {
436  MOFEM_LOG("SELF", Sev::error) << "base finction: " << base_function;
437  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
438  "Wrong number number in base functions");
439  }
440  #endif
441  #endif
442 
443  values_at_gauss_pts(I) += field_data(I) * base_function;
444  ++field_data;
445  ++base_function;
446  }
447  // Number of dofs can be smaller than number of Tensor_Dim x base
448  // functions
449  for (; bb < nb_base_functions; ++bb)
450  ++base_function;
451  ++values_at_gauss_pts;
452  }
453  }
454 
455  if (dataVec.use_count()) {
456  data.getFieldData().swap(dotVector);
457  }
458  }
460 }
461 
462 /** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
463  * field
464  *
465  * \ingroup mofem_forces_and_sources_user_data_operators
466  */
467 template <int Tensor_Dim>
470  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
471 
473  Tensor_Dim, double, ublas::row_major,
475 };
476 
477 /**@}*/
478 
479 /** \name Vector field values at integration points */
480 
481 /**@{*/
482 
483 /** \brief Calculate field values (template specialization) for tensor field
484  * rank 1, i.e. vector field
485  *
486  */
487 template <int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
490 
492  const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
493  const EntityType zero_type = MBVERTEX)
496  dataPtr(data_ptr), zeroType(zero_type) {
497  if (!dataPtr)
498  THROW_MESSAGE("Pointer is not set");
499  }
500 
501  MoFEMErrorCode doWork(int side, EntityType type,
504 
505  // When we move to C++17 add if constexpr()
506  if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
507  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
508  "%s coordiante not implemented",
509  CoordinateTypesNames[COORDINATE_SYSTEM]);
510 
511  const size_t nb_gauss_pts = getGaussPts().size2();
512  auto &vec = *dataPtr;
513  if (type == zeroType) {
514  vec.resize(nb_gauss_pts, false);
515  vec.clear();
516  }
517 
518  const size_t nb_dofs = data.getFieldData().size();
519  if (nb_dofs) {
520 
521  if (nb_gauss_pts) {
522  const size_t nb_base_functions = data.getN().size2();
523  auto values_at_gauss_pts = getFTensor0FromVec(vec);
525  const size_t size = nb_dofs / Tensor_Dim;
526  #ifndef NDEBUG
527  if (nb_dofs % Tensor_Dim) {
528  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
529  "Number of dofs should multiple of dimensions");
530  }
531  #endif
532 
533  // When we move to C++17 add if constexpr()
534  if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
535  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
536  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
537  auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
538  size_t bb = 0;
539  for (; bb != size; ++bb) {
540  values_at_gauss_pts += field_data(I) * diff_base_function(I);
541  ++field_data;
542  ++diff_base_function;
543  }
544  // Number of dofs can be smaller than number of Tensor_Dim x base
545  // functions
546  for (; bb < nb_base_functions; ++bb)
547  ++diff_base_function;
548  ++values_at_gauss_pts;
549  }
550  }
551 
552  // When we move to C++17 add if constexpr()
553  if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
554  auto t_coords = getFTensor1CoordsAtGaussPts();
555  auto values_at_gauss_pts = getFTensor0FromVec(vec);
556  auto base_function = data.getFTensor0N();
557  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
558  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
559  auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
560  size_t bb = 0;
561  for (; bb != size; ++bb) {
562  values_at_gauss_pts += field_data(I) * diff_base_function(I);
563  values_at_gauss_pts +=
564  base_function * (field_data(0) / t_coords(0));
565  ++field_data;
566  ++base_function;
567  ++diff_base_function;
568  }
569  // Number of dofs can be smaller than number of Tensor_Dim x base
570  // functions
571  for (; bb < nb_base_functions; ++bb) {
572  ++base_function;
573  ++diff_base_function;
574  }
575  ++values_at_gauss_pts;
576  ++t_coords;
577  }
578  }
579  }
580  }
582  }
583 
584 protected:
585  boost::shared_ptr<VectorDouble> dataPtr;
587 };
588 
589 /** \brief Approximate field values for given petsc vector
590  *
591  * \note Look at PetscData to see what vectors could be extracted with that user
592  * data operator.
593  *
594  * \ingroup mofem_forces_and_sources_user_data_operators
595  */
596 template <int Tensor_Dim, PetscData::DataContext CTX>
599 
601  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
602  const EntityType zero_at_type = MBVERTEX, bool throw_error = true)
605  dataPtr(data_ptr), zeroAtType(zero_at_type), throwError(throw_error) {
606  if (!dataPtr)
607  THROW_MESSAGE("Pointer is not set");
608  }
609 
610  MoFEMErrorCode doWork(int side, EntityType type,
613 
614  auto &local_indices = data.getLocalIndices();
615  const size_t nb_dofs = local_indices.size();
616  const size_t nb_gauss_pts = getGaussPts().size2();
617 
618  MatrixDouble &mat = *dataPtr;
619  if (type == zeroAtType) {
620  mat.resize(Tensor_Dim, nb_gauss_pts, false);
621  mat.clear();
622  }
623  if (!nb_dofs)
625 
626  if (!throwError) {
627  if ((getFEMethod()->data_ctx & PetscData::Switches(CTX)).none()) {
629  }
630  }
631 
632  const double *array;
633 
634  auto get_array = [&](const auto ctx, auto vec) {
636  #ifndef NDEBUG
637  if ((getFEMethod()->data_ctx & ctx).none()) {
638  MOFEM_LOG_CHANNEL("SELF");
639  MOFEM_LOG("SELF", Sev::error)
640  << "In this case filed degrees of freedom are read from vector. "
641  "That usually happens when time solver is used, and access to "
642  "first or second rates is needed. You probably not set ts_u, "
643  "ts_u_t, or ts_u_tt and associated data structure, i.e. "
644  "data_ctx to CTX_SET_X, CTX_SET_DX, CTX_SET_X_T, or "
645  "CTX_SET_X_TT respectively";
646  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
647  }
648  #endif
649  CHKERR VecGetArrayRead(vec, &array);
651  };
652 
653  auto restore_array = [&](auto vec) {
654  return VecRestoreArrayRead(vec, &array);
655  };
656 
657  switch (CTX) {
659  CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
660  break;
662  CHKERR get_array(PetscData::CtxSetX, getFEMethod()->dx);
663  break;
665  CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
666  break;
668  CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
669  break;
670  default:
671  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
672  "That case is not implemented");
673  }
674 
675  dotVector.resize(local_indices.size());
676  for (int i = 0; i != local_indices.size(); ++i)
677  if (local_indices[i] != -1)
678  dotVector[i] = array[local_indices[i]];
679  else
680  dotVector[i] = 0;
681 
682  switch (CTX) {
684  CHKERR restore_array(getFEMethod()->ts_u);
685  break;
687  CHKERR restore_array(getFEMethod()->dx);
688  break;
690  CHKERR restore_array(getFEMethod()->ts_u_t);
691  break;
693  CHKERR restore_array(getFEMethod()->ts_u_tt);
694  break;
695  default:
696  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
697  "That case is not implemented");
698  }
699 
700  const size_t nb_base_functions = data.getN().size2();
701  auto base_function = data.getFTensor0N();
702  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
703 
705  const size_t size = nb_dofs / Tensor_Dim;
706  if (nb_dofs % Tensor_Dim) {
707  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
708  }
709  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
710  auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(dotVector);
711  size_t bb = 0;
712  for (; bb != size; ++bb) {
713  values_at_gauss_pts(I) += field_data(I) * base_function;
714  ++field_data;
715  ++base_function;
716  }
717  // Number of dofs can be smaller than number of Tensor_Dim x base
718  // functions
719  for (; bb < nb_base_functions; ++bb)
720  ++base_function;
721  ++values_at_gauss_pts;
722  }
724  }
725 
726 protected:
727  boost::shared_ptr<MatrixDouble> dataPtr;
731 };
732 
733 /** \brief Get rate of values at integration pts for tensor filed
734  * rank 1, i.e. vector field
735  *
736  * \ingroup mofem_forces_and_sources_user_data_operators
737  */
738 template <int Tensor_Dim>
742 
743 /** \brief Get second rate of values at integration pts for tensor
744  * filed rank 1, i.e. vector field
745  *
746  * \ingroup mofem_forces_and_sources_user_data_operators
747  */
748 template <int Tensor_Dim>
752 
753 /** \brief Get second time second update vector at integration pts for tensor
754  * filed rank 1, i.e. vector field
755  *
756  * \ingroup mofem_forces_and_sources_user_data_operators
757  */
758 template <int Tensor_Dim>
762 
763 /**@}*/
764 
765 /** \name Tensor field values at integration points */
766 
767 /**@{*/
768 
769 /** \brief Calculate field values for tenor field rank 2.
770  *
771  */
772 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
775 
777  const std::string field_name,
778  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
779  const EntityType zero_type = MBVERTEX)
782  dataPtr(data_ptr), zeroType(zero_type) {
783  if (!dataPtr)
784  THROW_MESSAGE("Pointer is not set");
785  }
786 
787  MoFEMErrorCode doWork(int side, EntityType type,
789 
790 protected:
791  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
793 };
794 
795 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
798  doWork(int side, EntityType type, EntitiesFieldData::EntData &data) {
800  SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
801  "Not implemented for T = %s, dim0 = %d and dim1 = %d",
802  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
803  Tensor_Dim0, Tensor_Dim1);
805 }
806 
807 template <int Tensor_Dim0, int Tensor_Dim1>
808 struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
809  ublas::row_major, DoubleAllocator>
811 
813  const std::string field_name,
814  boost::shared_ptr<
815  ublas::matrix<double, ublas::row_major, DoubleAllocator>>
816  data_ptr,
817  const EntityType zero_type = MBVERTEX)
820  dataPtr(data_ptr), zeroType(zero_type) {
821  if (!dataPtr)
822  THROW_MESSAGE("Pointer is not set");
823  }
824 
826  const std::string field_name,
827  boost::shared_ptr<
828  ublas::matrix<double, ublas::row_major, DoubleAllocator>>
829  data_ptr,
830  SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
833  dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
834  if (!dataPtr)
835  THROW_MESSAGE("Pointer is not set");
836  }
837 
838  MoFEMErrorCode doWork(int side, EntityType type,
840 
841 protected:
842  boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
847 };
848 
849 template <int Tensor_Dim0, int Tensor_Dim1>
851  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
852  DoubleAllocator>::doWork(int side, EntityType type,
855 
856  MatrixDouble &mat = *dataPtr;
857  const size_t nb_gauss_pts = data.getN().size1();
858  if (type == zeroType) {
859  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
860  mat.clear();
861  }
862 
863  const size_t nb_dofs = data.getFieldData().size();
864 
865  if (dataVec.use_count()) {
866  dotVector.resize(nb_dofs, false);
867  const double *array;
868  CHKERR VecGetArrayRead(dataVec, &array);
869  const auto &local_indices = data.getLocalIndices();
870  for (int i = 0; i != local_indices.size(); ++i)
871  if (local_indices[i] != -1)
872  dotVector[i] = array[local_indices[i]];
873  else
874  dotVector[i] = 0;
875  CHKERR VecRestoreArrayRead(dataVec, &array);
876  data.getFieldData().swap(dotVector);
877  }
878 
879  if (nb_dofs) {
880  const size_t nb_base_functions = data.getN().size2();
881  auto base_function = data.getFTensor0N();
882  auto values_at_gauss_pts =
883  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
886  const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
887  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
888  auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
889  size_t bb = 0;
890  for (; bb != size; ++bb) {
891  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
892  ++field_data;
893  ++base_function;
894  }
895  for (; bb < nb_base_functions; ++bb)
896  ++base_function;
897  ++values_at_gauss_pts;
898  }
899 
900  if (dataVec.use_count()) {
901  data.getFieldData().swap(dotVector);
902  }
903  }
905 }
906 
907 /** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
908  * field
909  *
910  * \ingroup mofem_forces_and_sources_user_data_operators
911  */
912 template <int Tensor_Dim0, int Tensor_Dim1>
915  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
916 
918  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
920 };
921 
922 /** \brief Get time direvarive values at integration pts for tensor filed rank
923  * 2, i.e. matrix field
924  *
925  * \ingroup mofem_forces_and_sources_user_data_operators
926  */
927 template <int Tensor_Dim0, int Tensor_Dim1>
930 
932  boost::shared_ptr<MatrixDouble> data_ptr,
933  const EntityType zero_at_type = MBVERTEX)
936  dataPtr(data_ptr), zeroAtType(zero_at_type) {
937  if (!dataPtr)
938  THROW_MESSAGE("Pointer is not set");
939  }
940 
941  MoFEMErrorCode doWork(int side, EntityType type,
944 
945  const size_t nb_gauss_pts = getGaussPts().size2();
946  MatrixDouble &mat = *dataPtr;
947  if (type == zeroAtType) {
948  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
949  mat.clear();
950  }
951  const auto &local_indices = data.getLocalIndices();
952  const size_t nb_dofs = local_indices.size();
953  if (nb_dofs) {
954  dotVector.resize(nb_dofs, false);
955  const double *array;
956  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
957  for (size_t i = 0; i != local_indices.size(); ++i)
958  if (local_indices[i] != -1)
959  dotVector[i] = array[local_indices[i]];
960  else
961  dotVector[i] = 0;
962  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
963 
964  const size_t nb_base_functions = data.getN().size2();
965 
966  auto base_function = data.getFTensor0N();
967  auto values_at_gauss_pts =
968  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
971  const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
972  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
973  auto field_data =
974  getFTensor2FromPtr<Tensor_Dim0, Tensor_Dim1>(&*dotVector.begin());
975  size_t bb = 0;
976  for (; bb != size; ++bb) {
977  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
978  ++field_data;
979  ++base_function;
980  }
981  for (; bb < nb_base_functions; ++bb)
982  ++base_function;
983  ++values_at_gauss_pts;
984  }
985  }
987  }
988 
989 protected:
990  boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
991  EntityType zeroAtType; ///< Zero values at Gauss point at this type
992  VectorDouble dotVector; ///< Keeps temporary values of time derivatives
993 };
994 
995 /**
996  * @brief Calculate symmetric tensor field values at integration pts.
997  *
998  * @tparam Tensor_Dim
999 
1000  * \ingroup mofem_forces_and_sources_user_data_operators
1001  */
1002 template <int Tensor_Dim>
1005 
1007  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1008  const EntityType zero_type = MBEDGE, const int zero_side = 0)
1011  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1012  if (!dataPtr)
1013  THROW_MESSAGE("Pointer is not set");
1014  }
1015 
1017  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1018  SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBEDGE,
1019  const int zero_side = 0)
1022  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side),
1023  dataVec(data_vec) {
1024  if (!dataPtr)
1025  THROW_MESSAGE("Pointer is not set");
1026  }
1027 
1028  MoFEMErrorCode doWork(int side, EntityType type,
1031  MatrixDouble &mat = *dataPtr;
1032  const int nb_gauss_pts = getGaussPts().size2();
1033  if (type == this->zeroType && side == zeroSide) {
1034  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
1035  mat.clear();
1036  }
1037  const int nb_dofs = data.getFieldData().size();
1038  if (!nb_dofs)
1040 
1041  if (dataVec.use_count()) {
1042  dotVector.resize(nb_dofs, false);
1043  const double *array;
1044  CHKERR VecGetArrayRead(dataVec, &array);
1045  const auto &local_indices = data.getLocalIndices();
1046  for (int i = 0; i != local_indices.size(); ++i)
1047  if (local_indices[i] != -1)
1048  dotVector[i] = array[local_indices[i]];
1049  else
1050  dotVector[i] = 0;
1051  CHKERR VecRestoreArrayRead(dataVec, &array);
1052  data.getFieldData().swap(dotVector);
1053  }
1054 
1055  const int nb_base_functions = data.getN().size2();
1056  auto base_function = data.getFTensor0N();
1057  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1060  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1061  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1062  auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
1063  int bb = 0;
1064  for (; bb != size; ++bb) {
1065  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1066  ++field_data;
1067  ++base_function;
1068  }
1069  for (; bb < nb_base_functions; ++bb)
1070  ++base_function;
1071  ++values_at_gauss_pts;
1072  }
1073 
1074  if (dataVec.use_count()) {
1075  data.getFieldData().swap(dotVector);
1076  }
1077 
1079  }
1080 
1081 protected:
1082  boost::shared_ptr<MatrixDouble> dataPtr;
1084  const int zeroSide;
1087 };
1088 
1089 /**
1090  * @brief Calculate symmetric tensor field rates ant integratio pts.
1091  *
1092  * @tparam Tensor_Dim
1093  *
1094  * \ingroup mofem_forces_and_sources_user_data_operators
1095  */
1096 template <int Tensor_Dim>
1099 
1101  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1102  const EntityType zero_type = MBEDGE, const int zero_side = 0)
1105  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1106  if (!dataPtr)
1107  THROW_MESSAGE("Pointer is not set");
1108  }
1109 
1110  MoFEMErrorCode doWork(int side, EntityType type,
1113  const int nb_gauss_pts = getGaussPts().size2();
1114  MatrixDouble &mat = *dataPtr;
1115  constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1116  if (type == zeroType && side == zeroSide) {
1117  mat.resize(symm_size, nb_gauss_pts, false);
1118  mat.clear();
1119  }
1120  auto &local_indices = data.getLocalIndices();
1121  const int nb_dofs = local_indices.size();
1122  if (!nb_dofs)
1124 
1125  #ifndef NDEBUG
1126  if ((getFEMethod()->data_ctx & PetscData::CtxSetX_T).none()) {
1127  MOFEM_LOG_CHANNEL("SELF");
1128  MOFEM_LOG("SELF", Sev::error)
1129  << "In this case filed degrees of freedom are read from vector. "
1130  "That usually happens when time solver is used, and acces to "
1131  "first rates is needed. You probably not set "
1132  "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1133  "respectively";
1134  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
1135  }
1136  #endif
1137 
1138  dotVector.resize(nb_dofs, false);
1139  const double *array;
1140  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1141  for (int i = 0; i != local_indices.size(); ++i)
1142  if (local_indices[i] != -1)
1143  dotVector[i] = array[local_indices[i]];
1144  else
1145  dotVector[i] = 0;
1146  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1147 
1148  const int nb_base_functions = data.getN().size2();
1149 
1150  auto base_function = data.getFTensor0N();
1151  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1154  const int size = nb_dofs / symm_size;
1155  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1156  auto field_data = getFTensorDotData<Tensor_Dim>();
1157  int bb = 0;
1158  for (; bb != size; ++bb) {
1159  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1160  ++field_data;
1161  ++base_function;
1162  }
1163  for (; bb < nb_base_functions; ++bb)
1164  ++base_function;
1165  ++values_at_gauss_pts;
1166  }
1167 
1169  }
1170 
1171 protected:
1172  boost::shared_ptr<MatrixDouble> dataPtr;
1174  const int zeroSide;
1176 
1177  template <int Dim> inline auto getFTensorDotData() {
1178  static_assert(Dim || !Dim, "not implemented");
1179  }
1180 };
1181 
1182 template <>
1183 template <>
1184 inline auto
1187  &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1188  &dotVector[5]);
1189 }
1190 
1191 template <>
1192 template <>
1193 inline auto
1196  &dotVector[0], &dotVector[1], &dotVector[2]);
1197 }
1198 
1199 /**@}*/
1200 
1201 /** \name Gradients and Hessian of scalar fields at integration points */
1202 
1203 /**@{*/
1204 
1205 /**
1206  * \brief Evaluate field gradient values for scalar field, i.e. gradient is
1207  * tensor rank 1 (vector)
1208  *
1209  */
1210 template <int Tensor_Dim, class T, class L, class A>
1212  : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
1213 
1215  Tensor_Dim, T, L, A>::OpCalculateVectorFieldValues_General;
1216 };
1217 
1218 /** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1219  * tensor rank 1 (vector), specialization
1220  *
1221  */
1222 template <int Tensor_Dim>
1224  ublas::row_major, DoubleAllocator>
1226  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1227 
1229  Tensor_Dim, double, ublas::row_major,
1231 
1232  /**
1233  * \brief calculate gradient values of scalar field at integration points
1234  * @param side side entity number
1235  * @param type side entity type
1236  * @param data entity data
1237  * @return error code
1238  */
1239  MoFEMErrorCode doWork(int side, EntityType type,
1241 };
1242 
1243 /**
1244  * \brief Member function specialization calculating scalar field gradients for
1245  * tenor field rank 1
1246  *
1247  */
1248 template <int Tensor_Dim>
1250  Tensor_Dim, double, ublas::row_major,
1251  DoubleAllocator>::doWork(int side, EntityType type,
1254 
1255  const size_t nb_gauss_pts = this->getGaussPts().size2();
1256  auto &mat = *this->dataPtr;
1257  if (type == this->zeroType) {
1258  mat.resize(Tensor_Dim, nb_gauss_pts, false);
1259  mat.clear();
1260  }
1261 
1262  const int nb_dofs = data.getFieldData().size();
1263  if (nb_dofs) {
1264 
1265  if (nb_gauss_pts) {
1266  const int nb_base_functions = data.getN().size2();
1267  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
1268  auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1269 
1270  #ifndef NDEBUG
1271  if (nb_dofs > nb_base_functions)
1272  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1273  "Number of base functions inconsistent with number of DOFs "
1274  "(%d > %d)",
1275  nb_dofs, nb_base_functions);
1276 
1277  if (data.getDiffN().size2() != nb_base_functions * Tensor_Dim)
1278  SETERRQ2(
1279  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1280  "Number of base functions inconsistent with number of derivatives "
1281  "(%d != %d)",
1282  data.getDiffN().size2(), nb_base_functions);
1283 
1284  if (data.getDiffN().size1() != nb_gauss_pts)
1285  SETERRQ2(
1286  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1287  "Number of base functions inconsistent with number of integration "
1288  "pts (%d != %d)",
1289  data.getDiffN().size2(), nb_gauss_pts);
1290 
1291  #endif
1292 
1294  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1295  auto field_data = data.getFTensor0FieldData();
1296  int bb = 0;
1297  for (; bb != nb_dofs; ++bb) {
1298  gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
1299  ++field_data;
1300  ++diff_base_function;
1301  }
1302  // Number of dofs can be smaller than number of base functions
1303  for (; bb < nb_base_functions; ++bb)
1304  ++diff_base_function;
1305  ++gradients_at_gauss_pts;
1306  }
1307  }
1308  }
1309 
1311 }
1312 
1313 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1314  * vector field
1315  *
1316  * \ingroup mofem_forces_and_sources_user_data_operators
1317  */
1318 template <int Tensor_Dim>
1321  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1323  Tensor_Dim, double, ublas::row_major,
1325 };
1326 
1327 /** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1328  * tensor rank 1 (vector), specialization
1329  *
1330  */
1331 template <int Tensor_Dim>
1334  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1335 
1337  Tensor_Dim, double, ublas::row_major,
1339 
1340  /**
1341  * \brief calculate gradient values of scalar field at integration points
1342  * @param side side entity number
1343  * @param type side entity type
1344  * @param data entity data
1345  * @return error code
1346  */
1347  MoFEMErrorCode doWork(int side, EntityType type,
1349 };
1350 
1351 template <int Tensor_Dim>
1353  int side, EntityType type, EntitiesFieldData::EntData &data) {
1355 
1356  const size_t nb_gauss_pts = this->getGaussPts().size2();
1357  auto &mat = *this->dataPtr;
1358  if (type == this->zeroType) {
1359  mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts, false);
1360  mat.clear();
1361  }
1362 
1363  const int nb_dofs = data.getFieldData().size();
1364  if (nb_dofs) {
1365 
1366  if (nb_gauss_pts) {
1367  const int nb_base_functions = data.getN().size2();
1368 
1369  auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1370  #ifndef NDEBUG
1371  if (hessian_base.size1() != nb_gauss_pts) {
1372  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1373  "Wrong number of integration pts (%d != %d)",
1374  hessian_base.size1(), nb_gauss_pts);
1375  }
1376  if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1377  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1378  "Wrong number of base functions (%d != %d)",
1379  hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1380  nb_base_functions);
1381  }
1382  if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1383  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1384  "Wrong number of base functions (%d < %d)",
1385  hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1386  }
1387  #endif
1388 
1389  auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1390  &*hessian_base.data().begin());
1391 
1392  auto hessian_at_gauss_pts =
1393  getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1394 
1397  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1398  auto field_data = data.getFTensor0FieldData();
1399  int bb = 0;
1400  for (; bb != nb_dofs; ++bb) {
1401  hessian_at_gauss_pts(I, J) +=
1402  field_data * t_diff2_base_function(I, J);
1403  ++field_data;
1404  ++t_diff2_base_function;
1405  }
1406  // Number of dofs can be smaller than number of base functions
1407  for (; bb < nb_base_functions; ++bb) {
1408  ++t_diff2_base_function;
1409  }
1410 
1411  ++hessian_at_gauss_pts;
1412  }
1413  }
1414  }
1415 
1417 }
1418 
1419 /**}*/
1420 
1421 /** \name Gradients and hessian of tensor fields at integration points */
1422 
1423 /**@{*/
1424 
1425 /**
1426  * \brief Evaluate field gradient values for vector field, i.e. gradient is
1427  * tensor rank 2
1428  *
1429  */
1430 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1432  : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1433  L, A> {
1434 
1436  Tensor_Dim0, Tensor_Dim1, T, L, A>::OpCalculateTensor2FieldValues_General;
1437 };
1438 
1439 template <int Tensor_Dim0, int Tensor_Dim1>
1440 struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1441  ublas::row_major, DoubleAllocator>
1443  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1444 
1446  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1448 
1449  /**
1450  * \brief calculate values of vector field at integration points
1451  * @param side side entity number
1452  * @param type side entity type
1453  * @param data entity data
1454  * @return error code
1455  */
1456  MoFEMErrorCode doWork(int side, EntityType type,
1458 };
1459 
1460 /**
1461  * \brief Member function specialization calculating vector field gradients for
1462  * tenor field rank 2
1463  *
1464  */
1465 template <int Tensor_Dim0, int Tensor_Dim1>
1467  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1468  DoubleAllocator>::doWork(int side, EntityType type,
1471  if (!this->dataPtr)
1472  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1473  "Data pointer not allocated");
1474 
1475  const size_t nb_gauss_pts = this->getGaussPts().size2();
1476  auto &mat = *this->dataPtr;
1477  if (type == this->zeroType) {
1478  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1479  mat.clear();
1480  }
1481 
1482  if (nb_gauss_pts) {
1483  const size_t nb_dofs = data.getFieldData().size();
1484 
1485  if (nb_dofs) {
1486 
1487  if (this->dataVec.use_count()) {
1488  this->dotVector.resize(nb_dofs, false);
1489  const double *array;
1490  CHKERR VecGetArrayRead(this->dataVec, &array);
1491  const auto &local_indices = data.getLocalIndices();
1492  for (int i = 0; i != local_indices.size(); ++i)
1493  if (local_indices[i] != -1)
1494  this->dotVector[i] = array[local_indices[i]];
1495  else
1496  this->dotVector[i] = 0;
1497  CHKERR VecRestoreArrayRead(this->dataVec, &array);
1498  data.getFieldData().swap(this->dotVector);
1499  }
1500 
1501  const int nb_base_functions = data.getN().size2();
1502  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1503  auto gradients_at_gauss_pts =
1504  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1507  int size = nb_dofs / Tensor_Dim0;
1508  if (nb_dofs % Tensor_Dim0) {
1509  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1510  "Data inconsistency");
1511  }
1512  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1513  auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1514 
1515  #ifndef NDEBUG
1516  if (field_data.l2() != field_data.l2()) {
1517  MOFEM_LOG("SELF", Sev::error) << "field data " << field_data;
1518  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1519  "Wrong number in coefficients");
1520  }
1521  #endif
1522 
1523  int bb = 0;
1524  for (; bb < size; ++bb) {
1525  #ifndef SINGULARITY
1526  #ifndef NDEBUG
1527  if (diff_base_function.l2() != diff_base_function.l2()) {
1528  MOFEM_LOG("SELF", Sev::error)
1529  << "diff_base_function: " << diff_base_function;
1530  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1531  "Wrong number number in base functions");
1532  }
1533  #endif
1534  #endif
1535 
1536  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1537  ++field_data;
1538  ++diff_base_function;
1539  }
1540  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1541  // functions
1542  for (; bb != nb_base_functions; ++bb)
1543  ++diff_base_function;
1544  ++gradients_at_gauss_pts;
1545  }
1546 
1547  if (this->dataVec.use_count()) {
1548  data.getFieldData().swap(this->dotVector);
1549  }
1550  }
1551  }
1553 }
1554 
1555 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1556  * vector field
1557  *
1558  * \ingroup mofem_forces_and_sources_user_data_operators
1559  */
1560 template <int Tensor_Dim0, int Tensor_Dim1>
1563  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1564 
1566  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1568 };
1569 
1570 /** \brief Get field gradients time derivative at integration pts for scalar
1571  * filed rank 0, i.e. vector field
1572  *
1573  * \ingroup mofem_forces_and_sources_user_data_operators
1574  */
1575 template <int Tensor_Dim0, int Tensor_Dim1>
1578 
1580  boost::shared_ptr<MatrixDouble> data_ptr,
1581  const EntityType zero_at_type = MBVERTEX)
1584  dataPtr(data_ptr), zeroAtType(zero_at_type) {
1585  if (!dataPtr)
1586  THROW_MESSAGE("Pointer is not set");
1587  }
1588 
1589  MoFEMErrorCode doWork(int side, EntityType type,
1592 
1593  const auto &local_indices = data.getLocalIndices();
1594  const int nb_dofs = local_indices.size();
1595  const int nb_gauss_pts = this->getGaussPts().size2();
1596 
1597  auto &mat = *this->dataPtr;
1598  if (type == this->zeroAtType) {
1599  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1600  mat.clear();
1601  }
1602  if (!nb_dofs)
1604 
1605  dotVector.resize(nb_dofs, false);
1606  const double *array;
1607  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1608  for (int i = 0; i != local_indices.size(); ++i)
1609  if (local_indices[i] != -1)
1610  dotVector[i] = array[local_indices[i]];
1611  else
1612  dotVector[i] = 0;
1613  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1614 
1615  const int nb_base_functions = data.getN().size2();
1616  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1617  auto gradients_at_gauss_pts =
1618  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1621  int size = nb_dofs / Tensor_Dim0;
1622  if (nb_dofs % Tensor_Dim0) {
1623  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1624  }
1625 
1626  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1627  auto field_data = getFTensor1FromPtr<Tensor_Dim0>(&*dotVector.begin());
1628  int bb = 0;
1629  for (; bb < size; ++bb) {
1630  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1631  ++field_data;
1632  ++diff_base_function;
1633  }
1634  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1635  // functions
1636  for (; bb != nb_base_functions; ++bb)
1637  ++diff_base_function;
1638  ++gradients_at_gauss_pts;
1639  }
1641  }
1642 
1643 private:
1644  boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1645  EntityType zeroAtType; ///< Zero values at Gauss point at this type
1646  VectorDouble dotVector; ///< Keeps temporary values of time derivatives
1647 };
1648 
1649 /**
1650  * \brief Evaluate field gradient values for symmetric 2nd order tensor field,
1651  * i.e. gradient is tensor rank 3
1652  *
1653  */
1654 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1656  : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1657  L, A> {
1658 
1660  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1661  const EntityType zero_type = MBVERTEX)
1662  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
1663  A>(field_name, data_ptr,
1664  zero_type) {}
1665 };
1666 
1667 template <int Tensor_Dim0, int Tensor_Dim1>
1669  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator>
1671  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1672 
1674  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1675  const EntityType zero_type = MBVERTEX)
1676  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1677  ublas::row_major,
1678  DoubleAllocator>(
1679  field_name, data_ptr, zero_type) {}
1680 
1681  /**
1682  * \brief calculate values of vector field at integration points
1683  * @param side side entity number
1684  * @param type side entity type
1685  * @param data entity data
1686  * @return error code
1687  */
1688  MoFEMErrorCode doWork(int side, EntityType type,
1690 };
1691 
1692 /**
1693  * \brief Member function specialization calculating tensor field gradients for
1694  * symmetric tensor field rank 2
1695  *
1696  */
1697 template <int Tensor_Dim0, int Tensor_Dim1>
1698 MoFEMErrorCode OpCalculateTensor2SymmetricFieldGradient_General<
1699  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1700  DoubleAllocator>::doWork(int side, EntityType type,
1703  if (!this->dataPtr)
1704  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1705  "Data pointer not allocated");
1706 
1707  const size_t nb_gauss_pts = this->getGaussPts().size2();
1708  constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1709  auto &mat = *this->dataPtr;
1710  if (type == this->zeroType) {
1711  mat.resize(msize * Tensor_Dim1, nb_gauss_pts, false);
1712  mat.clear();
1713  }
1714 
1715  if (nb_gauss_pts) {
1716  const size_t nb_dofs = data.getFieldData().size();
1717 
1718  if (nb_dofs) {
1719 
1720  const int nb_base_functions = data.getN().size2();
1721  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1722  auto gradients_at_gauss_pts =
1723  getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1727  int size = nb_dofs / msize;
1728  if (nb_dofs % msize) {
1729  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1730  "Data inconsistency");
1731  }
1732  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1733  auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim0>();
1734  int bb = 0;
1735  for (; bb < size; ++bb) {
1736  gradients_at_gauss_pts(I, J, K) +=
1737  field_data(I, J) * diff_base_function(K);
1738  ++field_data;
1739  ++diff_base_function;
1740  }
1741  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1742  // functions
1743  for (; bb != nb_base_functions; ++bb)
1744  ++diff_base_function;
1745  ++gradients_at_gauss_pts;
1746  }
1747  }
1748  }
1750 }
1751 
1752 /** \brief Get field gradients at integration pts for symmetric tensorial field
1753  * rank 2
1754  *
1755  * \ingroup mofem_forces_and_sources_user_data_operators
1756  */
1757 template <int Tensor_Dim0, int Tensor_Dim1>
1760  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1761 
1763  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1764  const EntityType zero_type = MBVERTEX)
1766  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1767  DoubleAllocator>(field_name, data_ptr, zero_type) {}
1768 };
1769 
1770 template <int Tensor_Dim0, int Tensor_Dim1>
1773  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1774 
1776  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1778 
1779  /**
1780  * \brief calculate values of vector field at integration points
1781  * @param side side entity number
1782  * @param type side entity type
1783  * @param data entity data
1784  * @return error code
1785  */
1786  MoFEMErrorCode doWork(int side, EntityType type,
1788 };
1789 
1790 /**
1791  * \brief Member function specialization calculating vector field gradients for
1792  * tenor field rank 2
1793  *
1794  */
1795 template <int Tensor_Dim0, int Tensor_Dim1>
1797  int side, EntityType type, EntitiesFieldData::EntData &data) {
1799  if (!this->dataPtr)
1800  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1801  "Data pointer not allocated");
1802 
1803  const size_t nb_gauss_pts = this->getGaussPts().size2();
1804  auto &mat = *this->dataPtr;
1805  if (type == this->zeroType) {
1806  mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts, false);
1807  mat.clear();
1808  }
1809 
1810  if (nb_gauss_pts) {
1811  const size_t nb_dofs = data.getFieldData().size();
1812 
1813  if (nb_dofs) {
1814 
1815  if (this->dataVec.use_count()) {
1816  this->dotVector.resize(nb_dofs, false);
1817  const double *array;
1818  CHKERR VecGetArrayRead(this->dataVec, &array);
1819  const auto &local_indices = data.getLocalIndices();
1820  for (int i = 0; i != local_indices.size(); ++i)
1821  if (local_indices[i] != -1)
1822  this->dotVector[i] = array[local_indices[i]];
1823  else
1824  this->dotVector[i] = 0;
1825  CHKERR VecRestoreArrayRead(this->dataVec, &array);
1826  data.getFieldData().swap(this->dotVector);
1827  }
1828 
1829  const int nb_base_functions = data.getN().size2();
1830 
1831  auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1832  #ifndef NDEBUG
1833  if (hessian_base.size1() != nb_gauss_pts) {
1834  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1835  "Wrong number of integration pts (%d != %d)",
1836  hessian_base.size1(), nb_gauss_pts);
1837  }
1838  if (hessian_base.size2() !=
1839  nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
1840  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1841  "Wrong number of base functions (%d != %d)",
1842  hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
1843  nb_base_functions);
1844  }
1845  if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
1846  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1847  "Wrong number of base functions (%d < %d)",
1848  hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
1849  }
1850  #endif
1851 
1852  auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
1853  &*hessian_base.data().begin());
1854 
1855  auto t_hessian_at_gauss_pts =
1856  getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
1857 
1861 
1862  int size = nb_dofs / Tensor_Dim0;
1863  #ifndef NDEBUG
1864  if (nb_dofs % Tensor_Dim0) {
1865  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1866  "Data inconsistency");
1867  }
1868  #endif
1869 
1870  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1871  auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1872  int bb = 0;
1873  for (; bb < size; ++bb) {
1874  t_hessian_at_gauss_pts(I, J, K) +=
1875  field_data(I) * t_diff2_base_function(J, K);
1876  ++field_data;
1877  ++t_diff2_base_function;
1878  }
1879  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1880  // functions
1881  for (; bb != nb_base_functions; ++bb)
1882  ++t_diff2_base_function;
1883  ++t_hessian_at_gauss_pts;
1884  }
1885 
1886  if (this->dataVec.use_count()) {
1887  data.getFieldData().swap(this->dotVector);
1888  }
1889  }
1890  }
1892 }
1893 
1894 /**@}*/
1895 
1896 /** \name Transform tensors and vectors */
1897 
1898 /**@{*/
1899 
1900 /**
1901  * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1902  * \f$
1903  *
1904  * @tparam DIM
1905  *
1906  * \ingroup mofem_forces_and_sources_user_data_operators
1907  */
1908 template <int DIM_01, int DIM_23, int S = 0>
1911 
1914 
1915  /**
1916  * @deprecated Do not use this constructor
1917  */
1918  DEPRECATED
1920  boost::shared_ptr<MatrixDouble> in_mat,
1921  boost::shared_ptr<MatrixDouble> out_mat,
1922  boost::shared_ptr<MatrixDouble> d_mat)
1923  : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1924  // Only is run for vertices
1925  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1926  if (!inMat)
1927  THROW_MESSAGE("Pointer for in mat is null");
1928  if (!outMat)
1929  THROW_MESSAGE("Pointer for out mat is null");
1930  if (!dMat)
1931  THROW_MESSAGE("Pointer for tensor mat is null");
1932  }
1933 
1934  OpTensorTimesSymmetricTensor(boost::shared_ptr<MatrixDouble> in_mat,
1935  boost::shared_ptr<MatrixDouble> out_mat,
1936  boost::shared_ptr<MatrixDouble> d_mat)
1937  : UserOp(NOSPACE, OPSPACE), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1938  // Only is run for vertices
1939  if (!inMat)
1940  THROW_MESSAGE("Pointer for in mat is null");
1941  if (!outMat)
1942  THROW_MESSAGE("Pointer for out mat is null");
1943  if (!dMat)
1944  THROW_MESSAGE("Pointer for tensor mat is null");
1945  }
1946 
1947  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1949  const size_t nb_gauss_pts = getGaussPts().size2();
1950  auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
1951  auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
1952  outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1953  auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
1954  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1955  t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1956  ++t_in;
1957  ++t_out;
1958  }
1960  }
1961 
1962 private:
1967 
1968  boost::shared_ptr<MatrixDouble> inMat;
1969  boost::shared_ptr<MatrixDouble> outMat;
1970  boost::shared_ptr<MatrixDouble> dMat;
1971 };
1972 
1973 template <int DIM>
1975 
1978 
1979  /**
1980  * @deprecated Do not use this constructor
1981  */
1983  boost::shared_ptr<MatrixDouble> in_mat,
1984  boost::shared_ptr<MatrixDouble> out_mat)
1985  : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1986  // Only is run for vertices
1987  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1988  if (!inMat)
1989  THROW_MESSAGE("Pointer not set for in matrix");
1990  if (!outMat)
1991  THROW_MESSAGE("Pointer not set for in matrix");
1992  }
1993 
1994  OpSymmetrizeTensor(boost::shared_ptr<MatrixDouble> in_mat,
1995  boost::shared_ptr<MatrixDouble> out_mat)
1996  : UserOp(NOSPACE, OPSPACE), inMat(in_mat), outMat(out_mat) {
1997  // Only is run for vertices
1998  if (!inMat)
1999  THROW_MESSAGE("Pointer not set for in matrix");
2000  if (!outMat)
2001  THROW_MESSAGE("Pointer not set for in matrix");
2002  }
2003 
2004  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
2006  const size_t nb_gauss_pts = getGaussPts().size2();
2007  auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
2008  outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
2009  auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
2010  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
2011  t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
2012  ++t_in;
2013  ++t_out;
2014  }
2016  }
2017 
2018 private:
2021  boost::shared_ptr<MatrixDouble> inMat;
2022  boost::shared_ptr<MatrixDouble> outMat;
2023 };
2024 
2026 
2029 
2030  OpScaleMatrix(const std::string field_name, const double scale,
2031  boost::shared_ptr<MatrixDouble> in_mat,
2032  boost::shared_ptr<MatrixDouble> out_mat)
2033  : UserOp(field_name, OPROW), scale(scale), inMat(in_mat),
2034  outMat(out_mat) {
2035  // Only is run for vertices
2036  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
2037  if (!inMat)
2038  THROW_MESSAGE("Pointer not set for in matrix");
2039  if (!outMat)
2040  THROW_MESSAGE("Pointer not set for in matrix");
2041  }
2042 
2043  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
2045  outMat->resize(inMat->size1(), inMat->size2(), false);
2046  noalias(*outMat) = scale * (*inMat);
2048  }
2049 
2050 private:
2051  const double scale;
2052  boost::shared_ptr<MatrixDouble> inMat;
2053  boost::shared_ptr<MatrixDouble> outMat;
2054 };
2055 
2056 /**@}*/
2057 
2058 /** \name H-div/H-curls (Vectorial bases) values at integration points */
2059 
2060 /**@{*/
2061 
2062 /** \brief Get vector field for H-div approximation
2063  * \ingroup mofem_forces_and_sources_user_data_operators
2064  */
2065 template <int Base_Dim, int Field_Dim, class T, class L, class A>
2067 
2068 /** \brief Get vector field for H-div approximation
2069  * \ingroup mofem_forces_and_sources_user_data_operators
2070  */
2071 template <int Field_Dim>
2073  ublas::row_major, DoubleAllocator>
2075 
2077  boost::shared_ptr<MatrixDouble> data_ptr,
2078  SmartPetscObj<Vec> data_vec,
2079  const EntityType zero_type = MBEDGE,
2080  const int zero_side = 0)
2083  dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type),
2084  zeroSide(zero_side) {
2085  if (!dataPtr)
2086  THROW_MESSAGE("Pointer is not set");
2087  }
2088 
2090  boost::shared_ptr<MatrixDouble> data_ptr,
2091  const EntityType zero_type = MBEDGE,
2092  const int zero_side = 0)
2094  field_name, data_ptr, SmartPetscObj<Vec>(), zero_type, zero_side) {}
2095 
2096  /**
2097  * \brief Calculate values of vector field at integration points
2098  * @param side side entity number
2099  * @param type side entity type
2100  * @param data entity data
2101  * @return error code
2102  */
2103  MoFEMErrorCode doWork(int side, EntityType type,
2105 
2106 private:
2107  boost::shared_ptr<MatrixDouble> dataPtr;
2110  const int zeroSide;
2112 };
2113 
2114 template <int Field_Dim>
2116  3, Field_Dim, double, ublas::row_major,
2117  DoubleAllocator>::doWork(int side, EntityType type,
2120  const size_t nb_integration_points = this->getGaussPts().size2();
2121  if (type == zeroType && side == zeroSide) {
2122  dataPtr->resize(Field_Dim, nb_integration_points, false);
2123  dataPtr->clear();
2124  }
2125  const size_t nb_dofs = data.getFieldData().size();
2126  if (!nb_dofs)
2128 
2129  if (dataVec.use_count()) {
2130  dotVector.resize(nb_dofs, false);
2131  const double *array;
2132  CHKERR VecGetArrayRead(dataVec, &array);
2133  const auto &local_indices = data.getLocalIndices();
2134  for (int i = 0; i != local_indices.size(); ++i)
2135  if (local_indices[i] != -1)
2136  dotVector[i] = array[local_indices[i]];
2137  else
2138  dotVector[i] = 0;
2139  CHKERR VecRestoreArrayRead(dataVec, &array);
2140  data.getFieldData().swap(dotVector);
2141  }
2142 
2143  const size_t nb_base_functions = data.getN().size2() / 3;
2145  auto t_n_hdiv = data.getFTensor1N<3>();
2146  auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2147  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2148  auto t_dof = data.getFTensor0FieldData();
2149  int bb = 0;
2150  for (; bb != nb_dofs; ++bb) {
2151  t_data(i) += t_n_hdiv(i) * t_dof;
2152  ++t_n_hdiv;
2153  ++t_dof;
2154  }
2155  for (; bb != nb_base_functions; ++bb)
2156  ++t_n_hdiv;
2157  ++t_data;
2158  }
2159 
2160  if (dataVec.use_count()) {
2161  data.getFieldData().swap(dotVector);
2162  }
2164 }
2165 
2166 /** \brief Get vector field for H-div approximation
2167  *
2168  * \ingroup mofem_forces_and_sources_user_data_operators
2169  */
2170 template <int Base_Dim, int Field_Dim = Base_Dim>
2173  Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2175  Base_Dim, Field_Dim, double, ublas::row_major,
2177 };
2178 
2179 /** \brief Get vector field for H-div approximation
2180  * \ingroup mofem_forces_and_sources_user_data_operators
2181  */
2182 template <int Base_Dim, int Field_Dim = Base_Dim>
2184 
2185 template <int Field_Dim>
2186 struct OpCalculateHVecVectorFieldDot<3, Field_Dim>
2188 
2190  boost::shared_ptr<MatrixDouble> data_ptr,
2191  const EntityType zero_type = MBEDGE,
2192  const int zero_side = 0)
2195  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2196  if (!dataPtr)
2197  THROW_MESSAGE("Pointer is not set");
2198  }
2199 
2200  /**
2201  * \brief Calculate values of vector field at integration points
2202  * @param side side entity number
2203  * @param type side entity type
2204  * @param data entity data
2205  * @return error code
2206  */
2207  MoFEMErrorCode doWork(int side, EntityType type,
2209 
2210 private:
2211  boost::shared_ptr<MatrixDouble> dataPtr;
2213  const int zeroSide;
2214 };
2215 
2216 template <int Field_Dim>
2218  int side, EntityType type, EntitiesFieldData::EntData &data) {
2220 
2221  const size_t nb_integration_points = this->getGaussPts().size2();
2222  if (type == zeroType && side == zeroSide) {
2223  dataPtr->resize(Field_Dim, nb_integration_points, false);
2224  dataPtr->clear();
2225  }
2226 
2227  auto &local_indices = data.getIndices();
2228  const size_t nb_dofs = local_indices.size();
2229  if (nb_dofs) {
2230 
2231  std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2232  const double *array;
2233  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2234  for (size_t i = 0; i != nb_dofs; ++i)
2235  if (local_indices[i] != -1)
2236  dot_dofs_vector[i] = array[local_indices[i]];
2237  else
2238  dot_dofs_vector[i] = 0;
2239  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2240 
2241  const size_t nb_base_functions = data.getN().size2() / 3;
2243  auto t_n_hdiv = data.getFTensor1N<3>();
2244  auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2245  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2246  int bb = 0;
2247  for (; bb != nb_dofs; ++bb) {
2248  t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
2249  ++t_n_hdiv;
2250  }
2251  for (; bb != nb_base_functions; ++bb)
2252  ++t_n_hdiv;
2253  ++t_data;
2254  }
2255  }
2256 
2258 }
2259 
2260 /**
2261  * @brief Calculate divergence of vector field
2262  * @ingroup mofem_forces_and_sources_user_data_operators
2263  *
2264  * @tparam BASE_DIM
2265  * @tparam SPACE_DIM
2266  */
2267 template <int BASE_DIM, int SPACE_DIM>
2270 
2272  boost::shared_ptr<VectorDouble> data_ptr,
2273  const EntityType zero_type = MBEDGE,
2274  const int zero_side = 0)
2277  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2278  if (!dataPtr)
2279  THROW_MESSAGE("Pointer is not set");
2280  }
2281 
2282  MoFEMErrorCode doWork(int side, EntityType type,
2285  const size_t nb_integration_points = getGaussPts().size2();
2286  if (type == zeroType && side == zeroSide) {
2287  dataPtr->resize(nb_integration_points, false);
2288  dataPtr->clear();
2289  }
2290  const size_t nb_dofs = data.getFieldData().size();
2291  if (!nb_dofs)
2293  const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2296  auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2297  auto t_data = getFTensor0FromVec(*dataPtr);
2298  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2299  auto t_dof = data.getFTensor0FieldData();
2300  int bb = 0;
2301  for (; bb != nb_dofs; ++bb) {
2302  t_data += t_dof * t_n_diff_hdiv(j, j);
2303  ++t_n_diff_hdiv;
2304  ++t_dof;
2305  }
2306  for (; bb != nb_base_functions; ++bb)
2307  ++t_n_diff_hdiv;
2308  ++t_data;
2309  }
2311  }
2312 
2313 private:
2314  boost::shared_ptr<VectorDouble> dataPtr;
2316  const int zeroSide;
2317 };
2318 
2319 /**
2320  * @brief Calculate gradient of vector field
2321  * @ingroup mofem_forces_and_sources_user_data_operators
2322  *
2323  * @tparam BASE_DIM
2324  * @tparam SPACE_DIM
2325  */
2326 template <int BASE_DIM, int SPACE_DIM>
2329 
2331  boost::shared_ptr<MatrixDouble> data_ptr,
2332  const EntityType zero_type = MBEDGE,
2333  const int zero_side = 0)
2336  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2337  if (!dataPtr)
2338  THROW_MESSAGE("Pointer is not set");
2339  }
2340 
2341  MoFEMErrorCode doWork(int side, EntityType type,
2344  const size_t nb_integration_points = getGaussPts().size2();
2345  if (type == zeroType && side == zeroSide) {
2346  dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
2347  dataPtr->clear();
2348  }
2349  const size_t nb_dofs = data.getFieldData().size();
2350  if (!nb_dofs)
2352  const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2355  auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2356  auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*dataPtr);
2357  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2358  auto t_dof = data.getFTensor0FieldData();
2359  int bb = 0;
2360  for (; bb != nb_dofs; ++bb) {
2361  t_data(i, j) += t_dof * t_base_diff(i, j);
2362  ++t_base_diff;
2363  ++t_dof;
2364  }
2365  for (; bb != nb_base_functions; ++bb)
2366  ++t_base_diff;
2367  ++t_data;
2368  }
2370  }
2371 
2372 private:
2373  boost::shared_ptr<MatrixDouble> dataPtr;
2375  const int zeroSide;
2376 };
2377 
2378 /**
2379  * @brief Calculate gradient of vector field
2380  * @ingroup mofem_forces_and_sources_user_data_operators
2381  *
2382  * @tparam BASE_DIM
2383  * @tparam SPACE_DIM
2384  */
2385 template <int BASE_DIM, int SPACE_DIM>
2388 
2390  boost::shared_ptr<MatrixDouble> data_ptr,
2391  const EntityType zero_type = MBEDGE,
2392  const int zero_side = 0)
2395  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2396  if (!dataPtr)
2397  THROW_MESSAGE("Pointer is not set");
2398  }
2399 
2400  MoFEMErrorCode doWork(int side, EntityType type,
2403  const size_t nb_integration_points = getGaussPts().size2();
2404  if (type == zeroType && side == zeroSide) {
2405  dataPtr->resize(BASE_DIM * SPACE_DIM * SPACE_DIM, nb_integration_points,
2406  false);
2407  dataPtr->clear();
2408  }
2409  const size_t nb_dofs = data.getFieldData().size();
2410  if (!nb_dofs)
2412 
2413  const int nb_base_functions = data.getN().size2() / BASE_DIM;
2414 
2415  #ifndef NDEBUG
2416  auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2417  if (hessian_base.size1() != nb_integration_points) {
2418  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2419  "Wrong number of integration pts (%d != %d)",
2420  hessian_base.size1(), nb_integration_points);
2421  }
2422  if (hessian_base.size2() !=
2423  BASE_DIM * nb_base_functions * SPACE_DIM * SPACE_DIM) {
2424  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2425  "Wrong number of base functions (%d != %d)",
2426  hessian_base.size2() / (BASE_DIM * SPACE_DIM * SPACE_DIM),
2427  nb_base_functions);
2428  }
2429  if (hessian_base.size2() < BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM) {
2430  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2431  "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2432  BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM);
2433  }
2434  #endif
2435 
2439 
2440  auto t_base_diff2 =
2442  auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*dataPtr);
2443  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2444  auto t_dof = data.getFTensor0FieldData();
2445  int bb = 0;
2446  for (; bb != nb_dofs; ++bb) {
2447  t_data(i, j, k) += t_dof * t_base_diff2(i, j, k);
2448 
2449  ++t_base_diff2;
2450  ++t_dof;
2451  }
2452  for (; bb != nb_base_functions; ++bb)
2453  ++t_base_diff2;
2454  ++t_data;
2455  }
2457  }
2458 
2459 private:
2460  boost::shared_ptr<MatrixDouble> dataPtr;
2462  const int zeroSide;
2463 };
2464 
2465 /**
2466  * @brief Calculate divergence of vector field dot
2467  * @ingroup mofem_forces_and_sources_user_data_operators
2468  *
2469  * @tparam Tensor_Dim dimension of space
2470  */
2471 template <int Tensor_Dim1, int Tensor_Dim2>
2474 
2476  boost::shared_ptr<VectorDouble> data_ptr,
2477  const EntityType zero_type = MBEDGE,
2478  const int zero_side = 0)
2481  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2482  if (!dataPtr)
2483  THROW_MESSAGE("Pointer is not set");
2484  }
2485 
2486  MoFEMErrorCode doWork(int side, EntityType type,
2489  const size_t nb_integration_points = getGaussPts().size2();
2490  if (type == zeroType && side == zeroSide) {
2491  dataPtr->resize(nb_integration_points, false);
2492  dataPtr->clear();
2493  }
2494 
2495  const auto &local_indices = data.getLocalIndices();
2496  const int nb_dofs = local_indices.size();
2497  if (nb_dofs) {
2498 
2499  std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2500  const double *array;
2501  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2502  for (size_t i = 0; i != local_indices.size(); ++i)
2503  if (local_indices[i] != -1)
2504  dot_dofs_vector[i] = array[local_indices[i]];
2505  else
2506  dot_dofs_vector[i] = 0;
2507  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2508 
2509  const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
2511  auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
2512  auto t_data = getFTensor0FromVec(*dataPtr);
2513  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2514  int bb = 0;
2515  for (; bb != nb_dofs; ++bb) {
2516  double div = 0;
2517  for (auto ii = 0; ii != Tensor_Dim2; ++ii)
2518  div += t_n_diff_hdiv(ii, ii);
2519  t_data += dot_dofs_vector[bb] * div;
2520  ++t_n_diff_hdiv;
2521  }
2522  for (; bb != nb_base_functions; ++bb)
2523  ++t_n_diff_hdiv;
2524  ++t_data;
2525  }
2526  }
2528  }
2529 
2530 private:
2531  boost::shared_ptr<VectorDouble> dataPtr;
2533  const int zeroSide;
2534 };
2535 
2536 /**
2537  * @brief Calculate curl of vector field
2538  * @ingroup mofem_forces_and_sources_user_data_operators
2539  *
2540  * @tparam Base_Dim base function dimension
2541  * @tparam Space_Dim dimension of space
2542  * @tparam Hcurl field dimension
2543  */
2544 template <int Base_Dim, int Space_Dim> struct OpCalculateHcurlVectorCurl;
2545 
2546 /**
2547  * @brief Calculate curl of vector field
2548  * @ingroup mofem_forces_and_sources_user_data_operators
2549  *
2550  * @tparam Base_Dim base function dimension
2551  * @tparam Space_Dim dimension of space
2552  * @tparam Hcurl field dimension
2553  */
2554 template <>
2557  OpCalculateHcurlVectorCurl(const std::string field_name,
2558  boost::shared_ptr<MatrixDouble> data_ptr,
2559  const EntityType zero_type = MBEDGE,
2560  const int zero_side = 0);
2561  MoFEMErrorCode doWork(int side, EntityType type,
2563 
2564 private:
2565  boost::shared_ptr<MatrixDouble> dataPtr;
2567  const int zeroSide;
2568 };
2569 
2570 /**
2571  * @brief Calculate curl of vector field
2572  * @ingroup mofem_forces_and_sources_user_data_operators
2573  *
2574  * @tparam Field_Dim dimension of field
2575  * @tparam Space_Dim dimension of space
2576  */
2577 template <>
2580 
2581  OpCalculateHcurlVectorCurl(const std::string field_name,
2582  boost::shared_ptr<MatrixDouble> data_ptr,
2583  const EntityType zero_type = MBVERTEX,
2584  const int zero_side = 0);
2585 
2586  MoFEMErrorCode doWork(int side, EntityType type,
2588 
2589 private:
2590  boost::shared_ptr<MatrixDouble> dataPtr;
2592  const int zeroSide;
2593 };
2594 
2595 /**
2596  * @brief Calculate curl of vector field
2597  * @ingroup mofem_forces_and_sources_user_data_operators
2598  *
2599  * @tparam Field_Dim dimension of field
2600  * @tparam Space_Dim dimension of space
2601  */
2602 template <>
2605 
2606  OpCalculateHcurlVectorCurl(const std::string field_name,
2607  boost::shared_ptr<MatrixDouble> data_ptr,
2608  const EntityType zero_type = MBVERTEX,
2609  const int zero_side = 0);
2610 
2611  MoFEMErrorCode doWork(int side, EntityType type,
2613 
2614 private:
2615  boost::shared_ptr<MatrixDouble> dataPtr;
2617  const int zeroSide;
2618 };
2619 
2620 /**
2621  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
2622  * \ingroup mofem_forces_and_sources_user_data_operators
2623  *
2624  * @tparam Tensor_Dim0 rank of the filed
2625  * @tparam Tensor_Dim1 dimension of space
2626  */
2627 template <int Tensor_Dim0, int Tensor_Dim1>
2630 
2632  boost::shared_ptr<MatrixDouble> data_ptr,
2633  boost::shared_ptr<double> scale_ptr,
2635  const EntityType zero_type = MBEDGE,
2636  const int zero_side = 0)
2639  dataPtr(data_ptr), scalePtr(scale_ptr), dataVec(data_vec),
2640  zeroType(zero_type), zeroSide(zero_side) {
2641  if (!dataPtr)
2642  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
2643  }
2644 
2646  boost::shared_ptr<MatrixDouble> data_ptr,
2647  const EntityType zero_type = MBEDGE,
2648  const int zero_side = 0)
2649  : OpCalculateHVecTensorField(field_name, data_ptr, nullptr,
2650  SmartPetscObj<Vec>(), zero_type, zero_side) {
2651  }
2652 
2653  MoFEMErrorCode doWork(int side, EntityType type,
2656  const size_t nb_integration_points = getGaussPts().size2();
2657  if (type == zeroType && side == zeroSide) {
2658  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2659  dataPtr->clear();
2660  }
2661  const size_t nb_dofs = data.getFieldData().size();
2662  if (nb_dofs) {
2663 
2664  if (dataVec.use_count()) {
2665  dotVector.resize(nb_dofs, false);
2666  const double *array;
2667  CHKERR VecGetArrayRead(dataVec, &array);
2668  const auto &local_indices = data.getLocalIndices();
2669  for (int i = 0; i != local_indices.size(); ++i)
2670  if (local_indices[i] != -1)
2671  dotVector[i] = array[local_indices[i]];
2672  else
2673  dotVector[i] = 0;
2674  CHKERR VecRestoreArrayRead(dataVec, &array);
2675  data.getFieldData().swap(dotVector);
2676  }
2677 
2678  double scale = (scalePtr) ? *scalePtr : 1.0;
2679  const size_t nb_base_functions = data.getN().size2() / 3;
2682  auto t_n_hvec = data.getFTensor1N<3>();
2683  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2684  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2685  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
2686  size_t bb = 0;
2687  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
2688  t_data(i, j) += (scale * t_dof(i)) * t_n_hvec(j);
2689  ++t_n_hvec;
2690  ++t_dof;
2691  }
2692  for (; bb < nb_base_functions; ++bb)
2693  ++t_n_hvec;
2694  ++t_data;
2695  }
2696 
2697  if (dataVec.use_count()) {
2698  data.getFieldData().swap(dotVector);
2699  }
2700  }
2702  }
2703 
2704 private:
2705  boost::shared_ptr<MatrixDouble> dataPtr;
2706  boost::shared_ptr<double> scalePtr;
2709  const int zeroSide;
2710  VectorDouble dotVector; ///< Keeps temporary values of time derivatives
2711 };
2712 
2713 /** \brief Get tensor field for H-div approximation
2714  * \ingroup mofem_forces_and_sources_user_data_operators
2715  *
2716  * \warning This operator is not tested
2717  */
2718 template <int Tensor_Dim0, int Tensor_Dim1>
2721 
2723 
2725  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
2726  SmartPetscObj<Vec> data_vec, EntityType broken_type,
2727  boost::shared_ptr<Range> broken_range_ptr = nullptr,
2728  boost::shared_ptr<double> scale_ptr = nullptr,
2729  const EntityType zero_type = MBEDGE, const int zero_side = 0)
2731  dataPtr(data_ptr), dataVec(data_vec), brokenType(broken_type),
2732  brokenRangePtr(broken_range_ptr), zeroType(zero_type) {
2733  if (!dataPtr)
2734  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
2735  }
2736 
2737  /**
2738  * \brief Calculate values of vector field at integration points
2739  * @param side side entity number
2740  * @param type side entity type
2741  * @param data entity data
2742  * @return error code
2743  */
2744  MoFEMErrorCode doWork(int side, EntityType type,
2746 
2747 private:
2748  boost::shared_ptr<MatrixDouble> dataPtr;
2750  EntityType brokenType;
2751  boost::shared_ptr<Range> brokenRangePtr;
2752  boost::shared_ptr<double> scalePtr;
2755 };
2756 
2757 template <int Tensor_Dim0, int Tensor_Dim1>
2760  int side, EntityType type, EntitiesFieldData::EntData &data) {
2762  const size_t nb_integration_points = OP::getGaussPts().size2();
2763  if (type == zeroType) {
2764  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2765  dataPtr->clear();
2766  }
2767  const size_t nb_dofs = data.getFieldData().size();
2768  if (!nb_dofs)
2770 
2771  if (dataVec.use_count()) {
2772  dotVector.resize(nb_dofs, false);
2773  const double *array;
2774  CHKERR VecGetArrayRead(dataVec, &array);
2775  const auto &local_indices = data.getLocalIndices();
2776  for (int i = 0; i != local_indices.size(); ++i)
2777  if (local_indices[i] != -1)
2778  dotVector[i] = array[local_indices[i]];
2779  else
2780  dotVector[i] = 0;
2781  CHKERR VecRestoreArrayRead(dataVec, &array);
2782  data.getFieldData().swap(dotVector);
2783  }
2784 
2785  /**
2786  * @brief Get side face dofs
2787  *
2788  * Find which base functions on borken space have adjacent given entity type
2789  * and are in the range ptr if given.
2790  *
2791  */
2792  auto get_get_side_face_dofs = [&]() {
2793  auto fe_type = OP::getFEType();
2794 
2795  BaseFunction::DofsSideMap &side_dof_map =
2796  data.getFieldEntities()[0]->getDofSideMap().at(fe_type);
2797  std::vector<int> side_face_dofs;
2798  side_face_dofs.reserve(data.getIndices().size() / Tensor_Dim0);
2799 
2800  for (
2801 
2802  auto it = side_dof_map.get<1>().begin();
2803  it != side_dof_map.get<1>().end(); ++it
2804 
2805  ) {
2806  if ((Tensor_Dim0 * it->dof) >= data.getIndices().size()) {
2807  break;
2808  }
2809  if (it->type == brokenType) {
2810  if (brokenRangePtr) {
2811  auto ent = OP::getSideEntity(it->side, brokenType);
2812  if (brokenRangePtr->find(ent) != brokenRangePtr->end()) {
2813  side_face_dofs.push_back(it->dof);
2814  }
2815  } else {
2816  side_face_dofs.push_back(it->dof);
2817  }
2818  }
2819  }
2820 
2821  return side_face_dofs;
2822  };
2823 
2824  auto side_face_dofs = get_get_side_face_dofs();
2825 
2828  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2829  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2830  for (auto b : side_face_dofs) {
2831  auto t_row_base = data.getFTensor1N<3>(gg, b);
2832  auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(data.getFieldData().data() +
2833  b * Tensor_Dim0);
2834  t_data(i, j) += t_dof(i) * t_row_base(j);
2835  }
2836  ++t_data;
2837  }
2838  *dataPtr *= (scalePtr) ? *scalePtr : 1.0;
2839 
2840  if (dataVec.use_count()) {
2841  data.getFieldData().swap(dotVector);
2842  }
2843 
2845 }
2846 
2847 /**
2848  * @brief Calculate tenor field using tensor base, i.e. Hdiv/Hcurl
2849  * \ingroup mofem_forces_and_sources_user_data_operators
2850  *
2851  * @tparam Tensor_Dim0 rank of the filed
2852  * @tparam Tensor_Dim1 dimension of space
2853  */
2854 template <int Tensor_Dim0, int Tensor_Dim1>
2857 
2859  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
2860  boost::shared_ptr<double> scale_ptr,
2862  const EntityType zero_type = MBEDGE, const int zero_side = 0)
2865  dataPtr(data_ptr), scalePtr(scale_ptr), dataVec(data_vec),
2866  zeroType(zero_type), zeroSide(zero_side) {
2867  if (!dataPtr)
2868  THROW_MESSAGE("Pointer is not set");
2869  }
2870 
2872  boost::shared_ptr<MatrixDouble> data_ptr,
2873  const EntityType zero_type = MBEDGE,
2874  const int zero_side = 0)
2875  : OpCalculateHTensorTensorField(field_name, data_ptr, nullptr,
2876  SmartPetscObj<Vec>(), zero_type,
2877  zero_side) {}
2878 
2879  MoFEMErrorCode doWork(int side, EntityType type,
2882  const size_t nb_integration_points = getGaussPts().size2();
2883  if (type == zeroType && side == zeroSide) {
2884  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
2885  dataPtr->clear();
2886  }
2887  const size_t nb_dofs = data.getFieldData().size();
2888  if (!nb_dofs)
2890 
2891  if (dataVec.use_count()) {
2892  dotVector.resize(nb_dofs, false);
2893  const double *array;
2894  CHKERR VecGetArrayRead(dataVec, &array);
2895  const auto &local_indices = data.getLocalIndices();
2896  for (int i = 0; i != local_indices.size(); ++i)
2897  if (local_indices[i] != -1)
2898  dotVector[i] = array[local_indices[i]];
2899  else
2900  dotVector[i] = 0;
2901  CHKERR VecRestoreArrayRead(dataVec, &array);
2902  data.getFieldData().swap(dotVector);
2903  }
2904 
2905  double scale = (scalePtr) ? *scalePtr : 1.0;
2906  const size_t nb_base_functions =
2907  data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
2910  auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
2911  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
2912  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2913  auto t_dof = data.getFTensor0FieldData();
2914  size_t bb = 0;
2915  for (; bb != nb_dofs; ++bb) {
2916  t_data(i, j) += (scale * t_dof) * t_n_hten(i, j);
2917  ++t_n_hten;
2918  ++t_dof;
2919  }
2920  for (; bb < nb_base_functions; ++bb)
2921  ++t_n_hten;
2922  ++t_data;
2923  }
2924 
2925  if (dataVec.use_count()) {
2926  data.getFieldData().swap(dotVector);
2927  }
2928 
2930  }
2931 
2932 private:
2933  boost::shared_ptr<MatrixDouble> dataPtr;
2934  boost::shared_ptr<double> scalePtr;
2937  const int zeroSide;
2938  VectorDouble dotVector; ///< Keeps temporary values of time derivatives
2939 };
2940 
2941 /**
2942  * @brief Calculate divergence of tonsorial field using vectorial base
2943  * \ingroup mofem_forces_and_sources_user_data_operators
2944  *
2945  * @tparam Tensor_Dim0 rank of the field
2946  * @tparam Tensor_Dim1 dimension of space
2947  */
2948 template <int Tensor_Dim0, int Tensor_Dim1,
2949  CoordinateTypes CoordSys = CARTESIAN>
2952 
2954  boost::shared_ptr<MatrixDouble> data_ptr,
2955  SmartPetscObj<Vec> data_vec,
2956  const EntityType zero_type = MBEDGE,
2957  const int zero_side = 0)
2960  dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type),
2961  zeroSide(zero_side) {
2962  if (!dataPtr)
2963  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
2964  }
2965 
2967  boost::shared_ptr<MatrixDouble> data_ptr,
2968  const EntityType zero_type = MBEDGE,
2969  const int zero_side = 0)
2971  field_name, data_ptr, SmartPetscObj<Vec>(), zero_type, zero_side) {}
2972 
2973  MoFEMErrorCode doWork(int side, EntityType type,
2976  const size_t nb_integration_points = getGaussPts().size2();
2977  if (type == zeroType && side == zeroSide) {
2978  dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
2979  dataPtr->clear();
2980  }
2981  const size_t nb_dofs = data.getFieldData().size();
2982  if (nb_dofs) {
2983 
2984  if (dataVec.use_count()) {
2985  dotVector.resize(nb_dofs, false);
2986  const double *array;
2987  CHKERR VecGetArrayRead(dataVec, &array);
2988  const auto &local_indices = data.getLocalIndices();
2989  for (int i = 0; i != local_indices.size(); ++i)
2990  if (local_indices[i] != -1)
2991  dotVector[i] = array[local_indices[i]];
2992  else
2993  dotVector[i] = 0;
2994  CHKERR VecRestoreArrayRead(dataVec, &array);
2995  data.getFieldData().swap(dotVector);
2996  }
2997 
2998  const size_t nb_base_functions = data.getN().size2() / 3;
3001  auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
3002  auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
3003  auto t_base = data.getFTensor1N<3>();
3004  auto t_coords = getFTensor1CoordsAtGaussPts();
3005  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3006  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
3007  size_t bb = 0;
3008  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
3009  double div = t_n_diff_hvec(j, j);
3010  t_data(i) += t_dof(i) * div;
3011  if constexpr (CoordSys == CYLINDRICAL) {
3012  t_data(i) += t_base(0) * (t_dof(i) / t_coords(0));
3013  }
3014  ++t_n_diff_hvec;
3015  ++t_dof;
3016  ++t_base;
3017  }
3018  for (; bb < nb_base_functions; ++bb) {
3019  ++t_base;
3020  ++t_n_diff_hvec;
3021  }
3022  ++t_data;
3023  ++t_coords;
3024  }
3025 
3026  if (dataVec.use_count()) {
3027  data.getFieldData().swap(dotVector);
3028  }
3029  }
3031  }
3032 
3033 private:
3034  boost::shared_ptr<MatrixDouble> dataPtr;
3037  const int zeroSide;
3038 
3039  VectorDouble dotVector; ///< Keeps temporary values of time derivatives
3040 };
3041 
3042 /**
3043  * @brief Calculate divergence of tonsorial field using vectorial base
3044  * \ingroup mofem_forces_and_sources_user_data_operators
3045  *
3046  * \warning This operator is not tested
3047  *
3048  * @tparam Tensor_Dim0 rank of the field
3049  * @tparam Tensor_Dim1 dimension of space
3050  */
3051 template <int Tensor_Dim0, int Tensor_Dim1,
3052  CoordinateTypes CoordSys = CARTESIAN>
3055 
3057 
3059  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3060  SmartPetscObj<Vec> data_vec, EntityType broken_type,
3061  boost::shared_ptr<Range> broken_range_ptr = nullptr,
3062  boost::shared_ptr<double> scale_ptr = nullptr,
3063  const EntityType zero_type = MBEDGE)
3065  dataPtr(data_ptr), dataVec(data_vec), brokenType(broken_type),
3066  brokenRangePtr(broken_range_ptr), scalePtr(scale_ptr),
3067  zeroType(zero_type) {
3068  if (!dataPtr)
3069  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
3070  }
3071 
3072  MoFEMErrorCode doWork(int side, EntityType type,
3075  const size_t nb_integration_points = getGaussPts().size2();
3076  if (type == zeroType && side == 0) {
3077  dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
3078  dataPtr->clear();
3079  }
3080 
3081  const size_t nb_dofs = data.getFieldData().size();
3082  if (nb_dofs) {
3083 
3084  if (dataVec.use_count()) {
3085  dotVector.resize(nb_dofs, false);
3086  const double *array;
3087  CHKERR VecGetArrayRead(dataVec, &array);
3088  const auto &local_indices = data.getLocalIndices();
3089  for (int i = 0; i != local_indices.size(); ++i)
3090  if (local_indices[i] != -1)
3091  dotVector[i] = array[local_indices[i]];
3092  else
3093  dotVector[i] = 0;
3094  CHKERR VecRestoreArrayRead(dataVec, &array);
3095  data.getFieldData().swap(dotVector);
3096  }
3097 
3098  /**
3099  * @brief Get side face dofs
3100  *
3101  * Find which base functions on borken space have adjacent given entity
3102  * type and are in the range ptr if given.
3103  *
3104  */
3105  auto get_get_side_face_dofs = [&]() {
3106  auto fe_type = OP::getFEType();
3107 
3108  BaseFunction::DofsSideMap &side_dof_map =
3109  data.getFieldEntities()[0]->getDofSideMap().at(fe_type);
3110  std::vector<int> side_face_dofs;
3111  side_face_dofs.reserve(data.getIndices().size() / Tensor_Dim0);
3112 
3113  for (
3114 
3115  auto it = side_dof_map.get<1>().begin();
3116  it != side_dof_map.get<1>().end(); ++it
3117 
3118  ) {
3119  if ((Tensor_Dim0 * it->dof) >= data.getIndices().size()) {
3120  break;
3121  }
3122  if (it->type == brokenType) {
3123  if (brokenRangePtr) {
3124  auto ent = OP::getSideEntity(it->side, brokenType);
3125  if (brokenRangePtr->find(ent) != brokenRangePtr->end()) {
3126  side_face_dofs.push_back(it->dof);
3127  }
3128  } else {
3129  side_face_dofs.push_back(it->dof);
3130  }
3131  }
3132  }
3133 
3134  return side_face_dofs;
3135  };
3136 
3137  auto side_face_dofs = get_get_side_face_dofs();
3138 
3141  auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
3142  auto t_coords = getFTensor1CoordsAtGaussPts();
3143  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3144  for (auto b : side_face_dofs) {
3145  auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(
3146  data.getFieldData().data() + b * Tensor_Dim0);
3147  auto t_base = data.getFTensor1N<3>(gg, b);
3148  auto t_diff_base = data.getFTensor2DiffN<3, Tensor_Dim1>(gg, b);
3149  double div = t_diff_base(j, j);
3150  t_data(i) += t_dof(i) * div;
3151  if constexpr (CoordSys == CYLINDRICAL) {
3152  t_data(i) += t_base(0) * (t_dof(i) / t_coords(0));
3153  }
3154  }
3155  ++t_data;
3156  ++t_coords;
3157  }
3158  }
3159 
3160  if (dataVec.use_count()) {
3161  data.getFieldData().swap(dotVector);
3162  }
3163 
3165  }
3166 
3167 private:
3168  boost::shared_ptr<MatrixDouble> dataPtr;
3170  EntityType brokenType;
3171  boost::shared_ptr<Range> brokenRangePtr;
3172  boost::shared_ptr<double> scalePtr;
3175 };
3176 
3177 /**
3178  * @brief Calculate trace of vector (Hdiv/Hcurl) space
3179  *
3180  * @tparam Tensor_Dim
3181  * @tparam OpBase
3182  */
3183 template <int Tensor_Dim, typename OpBase>
3185 
3187  boost::shared_ptr<MatrixDouble> data_ptr,
3188  boost::shared_ptr<double> scale_ptr,
3189  const EntityType zero_type = MBEDGE,
3190  const int zero_side = 0)
3191  : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
3192  scalePtr(scale_ptr), zeroType(zero_type), zeroSide(zero_side) {
3193  if (!dataPtr)
3194  THROW_MESSAGE("Pointer is not set");
3195  }
3196 
3198  boost::shared_ptr<MatrixDouble> data_ptr,
3199  const EntityType zero_type = MBEDGE,
3200  const int zero_side = 0)
3201  : OpCalculateHVecTensorTrace(field_name, data_ptr, nullptr, zero_type,
3202  zero_side) {}
3203 
3204  MoFEMErrorCode doWork(int side, EntityType type,
3207  const size_t nb_integration_points = OpBase::getGaussPts().size2();
3208  if (type == zeroType && side == 0) {
3209  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
3210  dataPtr->clear();
3211  }
3212  const size_t nb_dofs = data.getFieldData().size();
3213  if (nb_dofs) {
3214  double scale_val = (scalePtr) ? *scalePtr : 1.0;
3215  auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
3216  const size_t nb_base_functions = data.getN().size2() / 3;
3217  auto t_base = data.getFTensor1N<3>();
3218  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
3219  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3220  FTensor::Tensor1<double, Tensor_Dim> t_normalized_normal;
3221  t_normalized_normal(j) = t_normal(j);
3222  t_normalized_normal.normalize();
3223  auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
3224  size_t bb = 0;
3225  for (; bb != nb_dofs / Tensor_Dim; ++bb) {
3226  t_data(i) +=
3227  (scale_val * t_dof(i)) * (t_base(j) * t_normalized_normal(j));
3228  ++t_base;
3229  ++t_dof;
3230  }
3231  for (; bb < nb_base_functions; ++bb) {
3232  ++t_base;
3233  }
3234  ++t_data;
3235  ++t_normal;
3236  }
3237  }
3239  }
3240 
3241 private:
3242  boost::shared_ptr<MatrixDouble> dataPtr;
3243  boost::shared_ptr<double> scalePtr;
3245  const int zeroSide;
3248 };
3249 
3250 /**@}*/
3251 
3252 /** \name Other operators */
3253 
3254 /**@{*/
3255 
3256 /**@}*/
3257 
3258 /** \name Operators for faces */
3259 
3260 /**@{*/
3261 
3262 /** \brief Transform local reference derivatives of shape functions to global
3263 derivatives
3264 
3265 \ingroup mofem_forces_and_sources_tri_element
3266 
3267 */
3268 template <int DIM, int DERIVATIVE = 1> struct OpSetInvJacSpaceForFaceImpl;
3269 
3272 
3274  boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3275 
3276 protected:
3277  template <int D1, int D2, int J1, int J2>
3280 
3281  static_assert(D2 == J2, "Dimension of jacobian and dimension of <out> "
3282  "directive does not match");
3283 
3284  size_t nb_functions = diff_n.size2() / D1;
3285  if (nb_functions) {
3286  size_t nb_gauss_pts = diff_n.size1();
3287 
3288  #ifndef NDEBUG
3289  if (nb_gauss_pts != getGaussPts().size2())
3290  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3291  "Wrong number of Gauss Pts");
3292  if (diff_n.size2() % D1)
3293  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3294  "Number of directives of base functions and D1 dimension does "
3295  "not match");
3296  #endif
3297 
3298  diffNinvJac.resize(diff_n.size1(), D2 * nb_functions, false);
3299 
3302  auto t_diff_n = getFTensor1FromPtr<D2>(&*diffNinvJac.data().begin());
3303  auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
3304  auto t_inv_jac = getFTensor2FromMat<J1, J2>(*invJacPtr);
3305  for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
3306  for (size_t dd = 0; dd != nb_functions; ++dd) {
3307  t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
3308  ++t_diff_n;
3309  ++t_diff_n_ref;
3310  }
3311  }
3312 
3313  diff_n.swap(diffNinvJac);
3314  }
3316  }
3317 
3318  boost::shared_ptr<MatrixDouble> invJacPtr;
3320 };
3321 
3322 template <>
3325 
3327 
3328  MoFEMErrorCode doWork(int side, EntityType type,
3330 };
3331 
3332 template <>
3334  : public OpSetInvJacSpaceForFaceImpl<2, 1> {
3335 
3337 
3338  MoFEMErrorCode doWork(int side, EntityType type,
3340 };
3341 
3342 template <>
3344  : public OpSetInvJacSpaceForFaceImpl<2, 1> {
3345 
3347 
3348  MoFEMErrorCode doWork(int side, EntityType type,
3350 };
3351 
3352 template <int DERIVARIVE = 1>
3354  : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
3355  OpSetInvJacH1ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3356  : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(H1, inv_jac_ptr) {}
3357 };
3358 
3359 template <int DERIVARIVE = 1>
3361  : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
3362  OpSetInvJacL2ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3363  : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(L2, inv_jac_ptr) {}
3364 };
3365 
3366 template <int DERIVARIVE = 1>
3368  : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
3370  boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3371  : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(H1, inv_jac_ptr) {}
3372 };
3373 
3374 template <int DERIVARIVE = 1>
3376  : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
3378  boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3379  : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(L2, inv_jac_ptr) {}
3380 };
3381 
3382 /**
3383  * \brief Transform local reference derivatives of shape function to
3384  global derivatives for face
3385 
3386  * \ingroup mofem_forces_and_sources_tri_element
3387  */
3388 template <int DIM> struct OpSetInvJacHcurlFaceImpl;
3389 
3390 template <>
3393 
3394  OpSetInvJacHcurlFaceImpl(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3396  invJacPtr(inv_jac_ptr) {}
3397 
3398  MoFEMErrorCode doWork(int side, EntityType type,
3400 
3401 protected:
3402  boost::shared_ptr<MatrixDouble> invJacPtr;
3404 };
3405 
3406 template <>
3409  MoFEMErrorCode doWork(int side, EntityType type,
3411 };
3412 
3415 
3416 /**
3417  * @brief Make Hdiv space from Hcurl space in 2d
3418  * @ingroup mofem_forces_and_sources_tri_element
3419  */
3422 
3425 
3426  MoFEMErrorCode doWork(int side, EntityType type,
3428 };
3429 
3430 /** \brief Transform Hcurl base fluxes from reference element to physical
3431  * triangle
3432  */
3434 
3435 /** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3436  *
3437  * Covariant Piola transformation
3438  \f[
3439  \psi_i|_t = J^{-1}_{ij}\hat{\psi}_j\\
3440  \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3441  = J^{-1}_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3442  \f]
3443 
3444  */
3445 template <>
3448 
3450  boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3451 
3452  MoFEMErrorCode doWork(int side, EntityType type,
3454 
3455 private:
3456  boost::shared_ptr<MatrixDouble> invJacPtr;
3457 
3460 };
3461 
3464 
3465 /** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3466  *
3467  * \note Hdiv space is generated by Hcurl space in 2d.
3468  *
3469  * Contravariant Piola transformation
3470  * \f[
3471  * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
3472  * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3473  * =
3474  * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3475  * \f]
3476  *
3477  * \ingroup mofem_forces_and_sources
3478  *
3479  */
3481 
3482 template <>
3485 
3487  boost::shared_ptr<MatrixDouble> jac_ptr)
3489  jacPtr(jac_ptr) {}
3490 
3491  MoFEMErrorCode doWork(int side, EntityType type,
3493 
3494 protected:
3495  boost::shared_ptr<MatrixDouble> jacPtr;
3498 };
3499 
3500 template <>
3505 
3506  MoFEMErrorCode doWork(int side, EntityType type,
3508 };
3509 
3514 
3515 /**@}*/
3516 
3517 /** \name Operators for edges */
3518 
3519 /**@{*/
3520 
3523 
3526 
3527  MoFEMErrorCode doWork(int side, EntityType type,
3529 };
3530 
3531 /**
3532  * @deprecated Name is deprecated and this is added for back compatibility
3533  */
3536 
3537 /**@}*/
3538 
3539 /** \name Operator for fat prisms */
3540 
3541 /**@{*/
3542 
3543 /**
3544  * @brief Operator for fat prism element updating integration weights in the
3545  * volume.
3546  *
3547  * Jacobian on the distorted element is nonconstant. This operator updates
3548  * integration weight on prism to take into account nonconstat jacobian.
3549  *
3550  * \f[
3551  * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
3552  * \pmb\xi} \right\| \right)
3553  * \f]
3554  * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
3555  * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
3556  * element coordinate.
3557  *
3558  */
3561 
3564 
3565  MoFEMErrorCode doWork(int side, EntityType type,
3567 };
3568 
3569 /** \brief Calculate inverse of jacobian for face element
3570 
3571  It is assumed that face element is XY plane. Applied
3572  only for 2d problems.
3573 
3574  FIXME Generalize function for arbitrary face orientation in 3d space
3575  FIXME Calculate to Jacobins for two faces
3576 
3577  \ingroup mofem_forces_and_sources_prism_element
3578 
3579 */
3582 
3583  OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3585  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3586 
3589  invJac(inv_jac) {}
3590  MoFEMErrorCode doWork(int side, EntityType type,
3592 
3593 private:
3594  const boost::shared_ptr<MatrixDouble> invJacPtr;
3596 };
3597 
3598 /** \brief Transform local reference derivatives of shape functions to global
3599 derivatives
3600 
3601 FIXME Generalize to curved shapes
3602 FIXME Generalize to case that top and bottom face has different shape
3603 
3604 \ingroup mofem_forces_and_sources_prism_element
3605 
3606 */
3609 
3610  OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3612  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3613 
3616  invJac(inv_jac) {}
3617 
3618  MoFEMErrorCode doWork(int side, EntityType type,
3620 
3621 private:
3622  const boost::shared_ptr<MatrixDouble> invJacPtr;
3625 };
3626 
3627 // Flat prism
3628 
3629 /** \brief Calculate inverse of jacobian for face element
3630 
3631  It is assumed that face element is XY plane. Applied
3632  only for 2d problems.
3633 
3634  FIXME Generalize function for arbitrary face orientation in 3d space
3635  FIXME Calculate to Jacobins for two faces
3636 
3637  \ingroup mofem_forces_and_sources_prism_element
3638 
3639 */
3642 
3645  invJacF3(inv_jac_f3) {}
3646  MoFEMErrorCode doWork(int side, EntityType type,
3648 
3649 private:
3651 };
3652 
3653 /** \brief Transform local reference derivatives of shape functions to global
3654 derivatives
3655 
3656 FIXME Generalize to curved shapes
3657 FIXME Generalize to case that top and bottom face has different shape
3658 
3659 \ingroup mofem_forces_and_sources_prism_element
3660 
3661 */
3664 
3667  invJacF3(inv_jac_f3) {}
3668 
3669  MoFEMErrorCode doWork(int side, EntityType type,
3671 
3672 private:
3675 };
3676 
3677 /**@}*/
3678 
3679 /** \name Operation on matrices at integration points */
3680 
3681 /**@{*/
3682 
3683 template <int DIM>
3685 
3686  OpInvertMatrix(boost::shared_ptr<MatrixDouble> in_ptr,
3687  boost::shared_ptr<VectorDouble> det_ptr,
3688  boost::shared_ptr<MatrixDouble> out_ptr)
3690  outPtr(out_ptr), detPtr(det_ptr) {}
3691 
3692  MoFEMErrorCode doWork(int side, EntityType type,
3694 
3695 private:
3696  boost::shared_ptr<MatrixDouble> inPtr;
3697  boost::shared_ptr<MatrixDouble> outPtr;
3698  boost::shared_ptr<VectorDouble> detPtr;
3699 };
3700 
3701 template <int DIM>
3705 
3706  if (!inPtr)
3707  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3708  "Pointer for inPtr matrix not allocated");
3709  if (!detPtr)
3710  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3711  "Pointer for detPtr matrix not allocated");
3712 
3713  const auto nb_rows = inPtr->size1();
3714  const auto nb_integration_pts = inPtr->size2();
3715 
3716  // Calculate determinant
3717  {
3718  detPtr->resize(nb_integration_pts, false);
3719  auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3720  auto t_det = getFTensor0FromVec(*detPtr);
3721  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3722  t_det = determinantTensor(t_in);
3723  ++t_in;
3724  ++t_det;
3725  }
3726  }
3727 
3728  // Invert jacobian
3729  if (outPtr) {
3730  outPtr->resize(nb_rows, nb_integration_pts, false);
3731  auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3732  auto t_out = getFTensor2FromMat<DIM, DIM>(*outPtr);
3733  auto t_det = getFTensor0FromVec(*detPtr);
3734  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3735  CHKERR invertTensor(t_in, t_det, t_out);
3736  ++t_in;
3737  ++t_out;
3738  ++t_det;
3739  }
3740  }
3741 
3743 }
3744 
3745 /**@}*/
3746 
3747 /** \brief Calculates the trace of an input matrix
3748 
3749 \ingroup mofem_forces_and_sources
3750 
3751 */
3752 
3753 template <int DIM>
3755 
3756  OpCalculateTraceFromMat(boost::shared_ptr<MatrixDouble> in_ptr,
3757  boost::shared_ptr<VectorDouble> out_ptr)
3759  outPtr(out_ptr) {}
3760 
3761  MoFEMErrorCode doWork(int side, EntityType type,
3763 
3764 private:
3766  boost::shared_ptr<MatrixDouble> inPtr;
3767  boost::shared_ptr<VectorDouble> outPtr;
3768 };
3769 
3770 template <int DIM>
3775 
3776  if (!inPtr)
3777  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3778  "Pointer for inPtr matrix not allocated");
3779 
3780  const auto nb_integration_pts = inPtr->size2();
3781  // Invert jacobian
3782  if (outPtr) {
3783  outPtr->resize(nb_integration_pts, false);
3784  auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
3785  auto t_out = getFTensor0FromVec(*outPtr);
3786 
3787  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3788  t_out = t_in(i, i);
3789  ++t_in;
3790  ++t_out;
3791  }
3792  }
3793 
3795 }
3796 
3797 /**@}*/
3798 
3799 /** \brief Calculates the trace of an input matrix
3800 
3801 \ingroup mofem_forces_and_sources
3802 
3803 */
3804 
3805 template <int DIM>
3808 
3809  OpCalculateTraceFromSymmMat(boost::shared_ptr<MatrixDouble> in_ptr,
3810  boost::shared_ptr<VectorDouble> out_ptr)
3812  outPtr(out_ptr) {}
3813 
3814  MoFEMErrorCode doWork(int side, EntityType type,
3816 
3817 private:
3819  boost::shared_ptr<MatrixDouble> inPtr;
3820  boost::shared_ptr<VectorDouble> outPtr;
3821 };
3822 
3823 template <int DIM>
3828 
3829  if (!inPtr)
3830  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3831  "Pointer for inPtr matrix not allocated");
3832 
3833  const auto nb_integration_pts = inPtr->size2();
3834  // Invert jacobian
3835  if (outPtr) {
3836  outPtr->resize(nb_integration_pts, false);
3837  auto t_in = getFTensor2SymmetricFromMat<DIM>(*inPtr);
3838  auto t_out = getFTensor0FromVec(*outPtr);
3839 
3840  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
3841  t_out = t_in(i, i);
3842  ++t_in;
3843  ++t_out;
3844  }
3845  }
3846 
3848 }
3849 
3850 } // namespace MoFEM
3851 
3852 #endif // __USER_DATA_OPERATORS_HPP__
3853 
3854 /**
3855  * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
3856  *
3857  * \brief Classes and functions used to evaluate fields at integration pts,
3858  *jacobians, etc..
3859  *
3860  * \ingroup mofem_forces_and_sources
3861  **/
UBlasMatrix< double >
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:460
MoFEM::OpCalculateVectorFieldGradientDot::dotVector
VectorDouble dotVector
Keeps temporary values of time derivatives.
Definition: UserDataOperators.hpp:1646
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:3686
MoFEM::PetscData::Switches
std::bitset< 8 > Switches
Definition: LoopMethods.hpp:34
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
MoFEM::ForcesAndSourcesCore::UserDataOperator::getSideEntity
EntityHandle getSideEntity(const int side_number, const EntityType type)
Get the side entity.
Definition: ForcesAndSourcesCore.hpp:1031
MoFEM::OpTensorTimesSymmetricTensor::i
FTensor::Index< 'i', DIM_01 > i
Definition: UserDataOperators.hpp:1963
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:1919
MoFEM::OpCalculateHVecTensorDivergence::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:3036
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::OpCalculateInvJacForFatPrism::OpCalculateInvJacForFatPrism
OpCalculateInvJacForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Definition: UserDataOperators.hpp:3583
MoFEM::OpSetContravariantPiolaTransformOnFace2DImpl< 2 >::piolaDiffN
MatrixDouble piolaDiffN
Definition: UserDataOperators.hpp:3497
MoFEM::EntitiesFieldData::EntData::getFieldEntities
const VectorFieldEntities & getFieldEntities() const
get field entities
Definition: EntitiesFieldData.hpp:1277
MoFEM::OpCalculateHcurlVectorCurl< 1, 3 >::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2617
MoFEM::OpSymmetrizeTensor::i
FTensor::Index< 'i', DIM > i
Definition: UserDataOperators.hpp:2019
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:2004
MoFEM::OpCalculateHcurlVectorCurl< 3, 3 >::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2566
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:610
MoFEM::OpSetCovariantPiolaTransformOnFace2DImpl
Transform Hcurl base fluxes from reference element to physical triangle.
Definition: UserDataOperators.hpp:3433
MoFEM::OpSetInvJacToScalarBasesBasic::diffNinvJac
MatrixDouble diffNinvJac
Definition: UserDataOperators.hpp:3319
MoFEM::OpTensorTimesSymmetricTensor::j
FTensor::Index< 'j', DIM_01 > j
Definition: UserDataOperators.hpp:1964
MoFEM::OpTensorTimesSymmetricTensor::inMat
boost::shared_ptr< MatrixDouble > inMat
Definition: UserDataOperators.hpp:1968
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:2189
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:1083
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
MoFEM::OpCalculateHTensorTensorField::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:2935
MoFEM::OpCalculateHTensorTensorField::dotVector
VectorDouble dotVector
Keeps temporary values of time derivatives.
Definition: UserDataOperators.hpp:2938
MoFEM::OpCalculateBrokenHVecTensorDivergence::brokenType
EntityType brokenType
Definition: UserDataOperators.hpp:3170
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MoFEM::OpCalculateScalarFieldValuesFromPetscVecImpl::zeroAtType
const EntityHandle zeroAtType
Definition: UserDataOperators.hpp:269
CARTESIAN
@ CARTESIAN
Definition: definitions.h:128
MoFEM::OpSetContravariantPiolaTransformOnFace2DImpl< 3 >
Definition: UserDataOperators.hpp:3501
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
MoFEM::OpMakeHdivFromHcurl::OpMakeHdivFromHcurl
OpMakeHdivFromHcurl()
Definition: UserDataOperators.hpp:3423
MoFEM::OpCalculateTensor2FieldValues
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:913
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:488
MoFEM::OpCalculateTraceFromMat::OpCalculateTraceFromMat
OpCalculateTraceFromMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Definition: UserDataOperators.hpp:3756
MoFEM::OpCalculateBrokenHVecTensorField::brokenType
EntityType brokenType
Definition: UserDataOperators.hpp:2750
MoFEM::OpCalculateTraceFromMat::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:3772
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:1006
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::OpCalculateHTensorTensorField::scalePtr
boost::shared_ptr< double > scalePtr
Definition: UserDataOperators.hpp:2934
MoFEM::BaseFunction::DofsSideMap
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > >>, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Definition: BaseFunction.hpp:73
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:2937
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:1097
MoFEM::OpSetContravariantPiolaTransformOnEdge2D::OpSetContravariantPiolaTransformOnEdge2D
OpSetContravariantPiolaTransformOnEdge2D(const FieldSpace space=HCURL)
Definition: UserDataOperators.hpp:3524
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:812
MoFEM::OpCalculateHVecTensorField::dotVector
VectorDouble dotVector
Keeps temporary values of time derivatives.
Definition: UserDataOperators.hpp:2710
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot::dotVector
VectorDouble dotVector
Definition: UserDataOperators.hpp:1175
MoFEM::OpSetInvJacH1ForFlatPrism::OpSetInvJacH1ForFlatPrism
OpSetInvJacH1ForFlatPrism(MatrixDouble &inv_jac_f3)
Definition: UserDataOperators.hpp:3665
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::ForcesAndSourcesCore::UserDataOperator::getFEType
EntityType getFEType() const
Get dimension of finite element.
Definition: ForcesAndSourcesCore.hpp:1013
MoFEM::EntitiesFieldData::EntData::getLocalIndices
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
Definition: EntitiesFieldData.hpp:1229
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::OpSymmetrizeTensor::inMat
boost::shared_ptr< MatrixDouble > inMat
Definition: UserDataOperators.hpp:2021
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::OpCalculateHVecTensorTrace::scalePtr
boost::shared_ptr< double > scalePtr
Definition: UserDataOperators.hpp:3243
MoFEM::OpCalculateHTensorTensorField
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2855
MoFEM::OpCalculateInvJacForFatPrism::invJacPtr
const boost::shared_ptr< MatrixDouble > invJacPtr
Definition: UserDataOperators.hpp:3594
MoFEM::OpSetInvJacHcurlFaceImpl< 2 >::OpSetInvJacHcurlFaceImpl
OpSetInvJacHcurlFaceImpl(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Definition: UserDataOperators.hpp:3394
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:1174
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:2879
MoFEM::OpSetInvJacH1ForFatPrism::OpSetInvJacH1ForFatPrism
OpSetInvJacH1ForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Definition: UserDataOperators.hpp:3610
MoFEM::OpCalculateHVecVectorGradient::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2374
MoFEM::OpCalculateVectorFieldValues_General::dataPtr
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Definition: UserDataOperators.hpp:318
MoFEM::OpCalculateBrokenHVecTensorField::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:2749
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateHVecTensorDivergence::OpCalculateHVecTensorDivergence
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBEDGE, const int zero_side=0)
Definition: UserDataOperators.hpp:2953
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:2475
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:468
MoFEM::OpCalculateBrokenHVecTensorDivergence::scalePtr
boost::shared_ptr< double > scalePtr
Definition: UserDataOperators.hpp:3172
MoFEM::OpCalculateHVecTensorTrace::OpCalculateHVecTensorTrace
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Definition: UserDataOperators.hpp:3186
MoFEM::OpSetContravariantPiolaTransformOnFace2DImpl< 2 >::piolaN
MatrixDouble piolaN
Definition: UserDataOperators.hpp:3496
MoFEM::OpSetCovariantPiolaTransformOnFace2DImpl< 2 >::piolaN
MatrixDouble piolaN
Definition: UserDataOperators.hpp:3458
MoFEM::OpCalculateHVecTensorField::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2709
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:1028
MoFEM::OpCalculateHVecVectorHessian::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2462
MoFEM::OpCalculateBrokenHVecTensorDivergence::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:3168
MoFEM::OpSetInvJacToScalarBasesBasic
Definition: UserDataOperators.hpp:3270
MoFEM::OpCalculateVectorFieldGradientDot::zeroAtType
EntityType zeroAtType
Zero values at Gauss point at this type.
Definition: UserDataOperators.hpp:1645
MoFEM::OpCalculateTensor2SymmetricFieldValues::dotVector
VectorDouble dotVector
Definition: UserDataOperators.hpp:1086
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:1254
MoFEM::OpCalculateScalarFieldHessian::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate gradient values of scalar field at integration points
Definition: UserDataOperators.hpp:1352
MoFEM::OpCalculateHVecTensorDivergence::dotVector
VectorDouble dotVector
Keeps temporary values of time derivatives.
Definition: UserDataOperators.hpp:3039
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:1762
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:2330
MoFEM::OpCalculateVectorFieldGradientDot::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
Definition: UserDataOperators.hpp:1644
MoFEM::OpCalculateTensor2FieldValues_General< Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator >::dotVector
VectorDouble dotVector
Definition: UserDataOperators.hpp:846
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::OpCalculateBrokenHVecTensorField::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2748
MoFEM::OpCalculateTensor2FieldValues_General
Calculate field values for tenor field rank 2.
Definition: UserDataOperators.hpp:773
MoFEM::OpCalculateScalarFieldValues_General::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:59
MoFEM::OpCalculateHcurlVectorCurl< 1, 3 >::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2615
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1329
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:2966
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl::throwError
bool throwError
Definition: UserDataOperators.hpp:730
MoFEM::OpCalculateVectorFieldHessian::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
Definition: UserDataOperators.hpp:1796
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:1970
MoFEM::OpSetInvJacH1ForFatPrism::diffNinvJac
MatrixDouble diffNinvJac
Definition: UserDataOperators.hpp:3624
MoFEM::OpCalculateHVecVectorHessian::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2460
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
scale
double scale
Definition: plastic.cpp:123
MoFEM::OpCalculateHdivVectorDivergence::dataPtr
boost::shared_ptr< VectorDouble > dataPtr
Definition: UserDataOperators.hpp:2314
MoFEM::OpSymmetrizeTensor::j
FTensor::Index< 'j', DIM > j
Definition: UserDataOperators.hpp:2020
MoFEM::OpCalculateTensor2FieldValues_General< Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator >::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:845
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:1947
CoordinateTypes
CoordinateTypes
Coodinate system.
Definition: definitions.h:127
MoFEM::OpSetContrariantPiolaTransformOnEdge
OpSetContravariantPiolaTransformOnEdge2D OpSetContrariantPiolaTransformOnEdge
Definition: UserDataOperators.hpp:3535
MoFEM::OpCalculateHVecTensorField::scalePtr
boost::shared_ptr< double > scalePtr
Definition: UserDataOperators.hpp:2706
MoFEM::OpCalculateHdivVectorDivergence::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2315
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
MoFEM::OpCalculateHVecTensorTrace::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:3245
MoFEM::OpCalculateHVecVectorField_General< 3, Field_Dim, double, ublas::row_major, DoubleAllocator >::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:2108
MoFEM::OpCalculateHVecVectorGradient
Calculate gradient of vector field.
Definition: UserDataOperators.hpp:2327
MoFEM::OpCalculateTensor2SymmetricFieldValues::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:1082
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3580
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:2645
MoFEM::OpCalculateInvJacForFlatPrism::OpCalculateInvJacForFlatPrism
OpCalculateInvJacForFlatPrism(MatrixDouble &inv_jac_f3)
Definition: UserDataOperators.hpp:3643
MoFEM::OpCalculateScalarFieldGradient_General
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector)
Definition: UserDataOperators.hpp:1211
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:1319
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
MoFEM::OpSetContravariantPiolaTransformOnFace2DImpl< 2 >
Definition: UserDataOperators.hpp:3483
MoFEM::OpCalculateHVecTensorTrace::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:3244
MoFEM::OpCalculateTensor2FieldValues_General< Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator >::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:844
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1237
MoFEM::OpCalculateTensor2FieldValuesDot::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
Definition: UserDataOperators.hpp:990
MoFEM::OpSetInvJacL2ForFace
Definition: UserDataOperators.hpp:3360
MoFEM::OpSetInvJacH1ForFatPrism::invJacPtr
const boost::shared_ptr< MatrixDouble > invJacPtr
Definition: UserDataOperators.hpp:3622
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:597
MoFEM::PetscData::CTX_SET_X_TT
@ CTX_SET_X_TT
Definition: LoopMethods.hpp:30
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
MoFEM::OpCalculateBrokenHVecTensorDivergence::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:3072
OpBase
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:1332
MoFEM::OpSetInvJacL2ForFaceEmbeddedIn3DSpace::OpSetInvJacL2ForFaceEmbeddedIn3DSpace
OpSetInvJacL2ForFaceEmbeddedIn3DSpace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Definition: UserDataOperators.hpp:3377
MoFEM::OpCalculateTraceFromSymmMat::OpCalculateTraceFromSymmMat
OpCalculateTraceFromSymmMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Definition: UserDataOperators.hpp:3809
MoFEM::OpCalculateBrokenHVecTensorField::dotVector
VectorDouble dotVector
Definition: UserDataOperators.hpp:2754
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:3362
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot::getFTensorDotData
auto getFTensorDotData()
Definition: UserDataOperators.hpp:1177
MoFEM::OpCalculateTraceFromSymmMat
Calculates the trace of an input matrix.
Definition: UserDataOperators.hpp:3806
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:727
MoFEM::OpCalculateBrokenHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:3053
CoordinateTypesNames
const static char *const CoordinateTypesNames[]
Coordinate system names.
Definition: definitions.h:121
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:2400
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::OpSetInvJacL2ForFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:3375
MoFEM::OpCalculateHVecTensorField::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2708
MoFEM::OpCalculateHcurlVectorCurl< 3, 3 >::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2567
MoFEM::OpInvertMatrix::inPtr
boost::shared_ptr< MatrixDouble > inPtr
Definition: UserDataOperators.hpp:3696
MoFEM::OpCalculateHVecVectorField_General< 3, Field_Dim, double, ublas::row_major, DoubleAllocator >::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2107
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:843
MoFEM::OpCalculateTensor2SymmetricFieldGradient
Get field gradients at integration pts for symmetric tensorial field rank 2.
Definition: UserDataOperators.hpp:1758
MoFEM::OpCalculateHVecTensorTrace::i
FTensor::Index< 'i', Tensor_Dim > i
Definition: UserDataOperators.hpp:3246
MoFEM::OpSetInvJacH1ForFaceEmbeddedIn3DSpace::OpSetInvJacH1ForFaceEmbeddedIn3DSpace
OpSetInvJacH1ForFaceEmbeddedIn3DSpace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Definition: UserDataOperators.hpp:3369
MoFEM::OpCalculateBrokenHVecTensorDivergence::OpCalculateBrokenHVecTensorDivergence
OpCalculateBrokenHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, EntityType broken_type, boost::shared_ptr< Range > broken_range_ptr=nullptr, boost::shared_ptr< double > scale_ptr=nullptr, const EntityType zero_type=MBEDGE)
Definition: UserDataOperators.hpp:3058
MoFEM::SmartPetscObj::use_count
int use_count() const
Definition: PetscSmartObj.hpp:109
MoFEM::OpCalculateHVecTensorField::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:2707
MoFEM::OpCalculateHdivVectorDivergenceDot::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2532
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl::dotVector
VectorDouble dotVector
Definition: UserDataOperators.hpp:729
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:2973
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
MoFEM::OpCalculateHVecVectorHessian
Calculate gradient of vector field.
Definition: UserDataOperators.hpp:2386
MoFEM::OpCalculateVectorFieldGradientDot::OpCalculateVectorFieldGradientDot
OpCalculateVectorFieldGradientDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
Definition: UserDataOperators.hpp:1579
MoFEM::OpCalculateBrokenHVecTensorDivergence::dotVector
VectorDouble dotVector
Definition: UserDataOperators.hpp:3174
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3607
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:2341
MoFEM::OpInvertMatrix::detPtr
boost::shared_ptr< VectorDouble > detPtr
Definition: UserDataOperators.hpp:3698
POLAR
@ POLAR
Definition: definitions.h:129
MoFEM::OpInvertMatrix::outPtr
boost::shared_ptr< MatrixDouble > outPtr
Definition: UserDataOperators.hpp:3697
MoFEM::OpCalculateHVecTensorTrace::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:3242
MoFEM::OpCalculateBrokenHVecTensorDivergence::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:3169
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:825
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:3767
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2171
MoFEM::OpCalculateTraceFromMat::i
FTensor::Index< 'i', DIM > i
Definition: UserDataOperators.hpp:3765
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, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBEDGE, const int zero_side=0)
Definition: UserDataOperators.hpp:2076
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:798
MoFEM::OpCalculateTensor2SymmetricFieldValues::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:1084
MoFEM::OpSetCovariantPiolaTransformOnFace2DImpl< 2 >
Apply contravariant (Piola) transfer to Hdiv space on face.
Definition: UserDataOperators.hpp:3446
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:2871
MoFEM::OpScaleMatrix::outMat
boost::shared_ptr< MatrixDouble > outMat
Definition: UserDataOperators.hpp:2053
convert.type
type
Definition: convert.py:64
MoFEM::OpCalculateBrokenHVecTensorField::scalePtr
boost::shared_ptr< double > scalePtr
Definition: UserDataOperators.hpp:2752
MoFEM::OpSetInvJacToScalarBasesBasic::applyTransform
MoFEMErrorCode applyTransform(MatrixDouble &diff_n)
Definition: UserDataOperators.hpp:3278
MoFEM::OpSetContravariantPiolaTransformOnFace2DImpl< 2 >::OpSetContravariantPiolaTransformOnFace2DImpl
OpSetContravariantPiolaTransformOnFace2DImpl(boost::shared_ptr< MatrixDouble > jac_ptr)
Definition: UserDataOperators.hpp:3486
MoFEM::OpCalculateTensor2FieldValuesDot::dotVector
VectorDouble dotVector
Keeps temporary values of time derivatives.
Definition: UserDataOperators.hpp:992
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:2933
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:1934
MoFEM::invertTensor
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
Definition: Templates.hpp:1771
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:2389
MoFEM::OpSetContravariantPiolaTransformOnEdge2D
Definition: UserDataOperators.hpp:3521
MoFEM::PetscData::CtxSetX_T
static constexpr Switches CtxSetX_T
Definition: LoopMethods.hpp:42
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::OpSetInvJacH1ForFlatPrism::invJacF3
MatrixDouble & invJacF3
Definition: UserDataOperators.hpp:3673
MoFEM::OpCalculateTraceFromSymmMat::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:3825
MoFEM::determinantTensor
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
Definition: Templates.hpp:1571
MoFEM::OpCalculateBrokenHVecTensorField::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate values of vector field at integration points.
Definition: UserDataOperators.hpp:2759
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
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:2271
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:501
MoFEM::OpCalculateHVecTensorField::OpCalculateHVecTensorField
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, SmartPetscObj< Vec > data_vec=SmartPetscObj< Vec >(), const EntityType zero_type=MBEDGE, const int zero_side=0)
Definition: UserDataOperators.hpp:2631
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2628
MoFEM::OpCalculateTraceFromSymmMat::inPtr
boost::shared_ptr< MatrixDouble > inPtr
Definition: UserDataOperators.hpp:3819
MoFEM::OpCalculateTensor2SymmetricFieldGradient_General
Evaluate field gradient values for symmetric 2nd order tensor field, i.e. gradient is tensor rank 3.
Definition: UserDataOperators.hpp:1655
MoFEM::OpCalculateHdivVectorDivergenceDot
Calculate divergence of vector field dot.
Definition: UserDataOperators.hpp:2472
SPHERICAL
@ SPHERICAL
Definition: definitions.h:131
MoFEM::OpCalculateHVecVectorHessian::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2461
MoFEM::OpCalculateInvJacForFatPrism::OpCalculateInvJacForFatPrism
OpCalculateInvJacForFatPrism(MatrixDouble &inv_jac)
Definition: UserDataOperators.hpp:3587
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:1016
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:3820
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl::zeroAtType
const EntityHandle zeroAtType
Definition: UserDataOperators.hpp:728
MoFEM::OpCalculateTensor2FieldValuesDot::zeroAtType
EntityType zeroAtType
Zero values at Gauss point at this type.
Definition: UserDataOperators.hpp:991
MoFEM::OpTensorTimesSymmetricTensor::l
FTensor::Index< 'l', DIM_23 > l
Definition: UserDataOperators.hpp:1966
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:1173
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::OpCalculateHVecTensorTrace::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: UserDataOperators.hpp:3204
MoFEM::OpSetCovariantPiolaTransformOnFace2DImpl< 2 >::diffPiolaN
MatrixDouble diffPiolaN
Definition: UserDataOperators.hpp:3459
MoFEM::OpSetInvJacHcurlFaceImpl
Transform local reference derivatives of shape function to global derivatives for face.
Definition: UserDataOperators.hpp:3388
MoFEM::OpSetInvJacHcurlFaceImpl< 2 >::diffHcurlInvJac
MatrixDouble diffHcurlInvJac
Definition: UserDataOperators.hpp:3403
MoFEM::OpScaleMatrix
Definition: UserDataOperators.hpp:2025
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:2486
MoFEM::OpCalculateHVecVectorGradient::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2375
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:3766
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:1909
MoFEM::OpCalculateScalarFieldValuesFromPetscVecImpl
Get rate of scalar field at integration points.
Definition: UserDataOperators.hpp:156
MoFEM::OpCalculateHVecVectorField_General< 3, Field_Dim, double, ublas::row_major, DoubleAllocator >::dotVector
VectorDouble dotVector
Definition: UserDataOperators.hpp:2111
MoFEM::OpSymmetrizeTensor::OpSymmetrizeTensor
OpSymmetrizeTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
Definition: UserDataOperators.hpp:1994
MoFEM::OpSetCovariantPiolaTransformOnFace2DImpl< 2 >::invJacPtr
boost::shared_ptr< MatrixDouble > invJacPtr
Definition: UserDataOperators.hpp:3456
MoFEM::OpCalculateHdivVectorDivergence::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2316
MoFEM::OpCalculateBrokenHVecTensorDivergence::brokenRangePtr
boost::shared_ptr< Range > brokenRangePtr
Definition: UserDataOperators.hpp:3171
MoFEM::OpCalculateTensor2FieldValuesDot
Get time direvarive values at integration pts for tensor filed rank 2, i.e. matrix field.
Definition: UserDataOperators.hpp:928
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:1965
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1043
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:1172
MoFEM::OpCalculateHVecTensorTrace
Calculate trace of vector (Hdiv/Hcurl) space.
Definition: UserDataOperators.hpp:3184
MoFEM::OpCalculateHVecVectorGradient::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2373
FTensor::Tensor1::normalize
Tensor1< T, Tensor_Dim > normalize()
Definition: Tensor1_value.hpp:77
MoFEM::OpCalculateTensor2SymmetricFieldValues::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:1085
MoFEM::OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms
Operator for fat prism element updating integration weights in the volume.
Definition: UserDataOperators.hpp:3559
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:1561
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:776
MoFEM::PetscData::CtxSetX_TT
static constexpr Switches CtxSetX_TT
Definition: LoopMethods.hpp:43
MoFEM::OpCalculateHVecTensorTrace::j
FTensor::Index< 'j', Tensor_Dim > j
Definition: UserDataOperators.hpp:3247
MoFEM::OpCalculateHTensorTensorField::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2936
MoFEM::OpCalculateHcurlVectorCurl< 3, 3 >::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2565
MoFEM::PetscData::CtxSetX
static constexpr Switches CtxSetX
Definition: LoopMethods.hpp:40
MoFEM::OpCalculateHcurlVectorCurl< 1, 2 >::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2591
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:941
MoFEM::OpSetInvJacH1ForFace
Definition: UserDataOperators.hpp:3353
MoFEM::OpCalculateHcurlVectorCurl< 1, 2 >::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2592
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:2950
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:791
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:3197
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:3402
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:130
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:1288
MoFEM::OpCalculateDivergenceVectorFieldValues::dataPtr
boost::shared_ptr< VectorDouble > dataPtr
Definition: UserDataOperators.hpp:585
MoFEM::OpSymmetrizeTensor
Definition: UserDataOperators.hpp:1974
MoFEM::PetscData::CTX_SET_DX
@ CTX_SET_DX
Definition: LoopMethods.hpp:28
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:2653
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:1318
MoFEM::OpCalculateHVecTensorDivergence::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:3034
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:2282
MoFEM::OpCalculateHcurlVectorCurl< 1, 2 >::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2590
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:3614
MoFEM::OpCalculateHVecTensorField::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2705
MoFEM::OpCalculateInvJacForFatPrism::invJac
MatrixDouble & invJac
Definition: UserDataOperators.hpp:3595
MoFEM::OpCalculateHVecVectorFieldDot< 3, Field_Dim >::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2213
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:3702
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:1982
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:1673
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:1110
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl::OpCalculateVectorFieldValuesFromPetscVecImpl
OpCalculateVectorFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX, bool throw_error=true)
Definition: UserDataOperators.hpp:600
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:3391
MoFEM::OpCalculateTensor2FieldValuesDot::OpCalculateTensor2FieldValuesDot
OpCalculateTensor2FieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
Definition: UserDataOperators.hpp:931
MoFEM::OpCalculateTraceFromMat
Calculates the trace of an input matrix.
Definition: UserDataOperators.hpp:3754
MoFEM::OpSetContravariantPiolaTransformOnFace2DImpl< 2 >::jacPtr
boost::shared_ptr< MatrixDouble > jacPtr
Definition: UserDataOperators.hpp:3495
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:2089
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:1003
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:64
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::OpSetInvJacH1ForFace::OpSetInvJacH1ForFace
OpSetInvJacH1ForFace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Definition: UserDataOperators.hpp:3355
MoFEM::OpCalculateHTensorTensorField::OpCalculateHTensorTensorField
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, SmartPetscObj< Vec > data_vec=SmartPetscObj< Vec >(), const EntityType zero_type=MBEDGE, const int zero_side=0)
Definition: UserDataOperators.hpp:2858
MoFEM::OpCalculateBrokenHVecTensorField::brokenRangePtr
boost::shared_ptr< Range > brokenRangePtr
Definition: UserDataOperators.hpp:2751
MoFEM::OpCalculateHdivVectorDivergence
Calculate divergence of vector field.
Definition: UserDataOperators.hpp:2268
MoFEM::OpCalculateBrokenHVecTensorDivergence::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:3173
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3684
MoFEM::OpCalculateInvJacForFlatPrism::invJacF3
MatrixDouble & invJacF3
Definition: UserDataOperators.hpp:3650
MoFEM::OpSetInvJacHcurlFaceImpl< 3 >
Definition: UserDataOperators.hpp:3407
MoFEM::PetscData::CTX_SET_X_T
@ CTX_SET_X_T
Definition: LoopMethods.hpp:29
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:2030
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:2066
MoFEM::OpCalculateHcurlVectorCurl< 1, 3 >::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2616
MoFEM::OpCalculateDivergenceVectorFieldValues::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:586
MoFEM::OpCalculateHVecVectorField_General< 3, Field_Dim, double, ublas::row_major, DoubleAllocator >::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2109
MoFEM::OpCalculateHdivVectorDivergenceDot::dataPtr
boost::shared_ptr< VectorDouble > dataPtr
Definition: UserDataOperators.hpp:2531
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:1589
MoFEM::OpSetInvJacToScalarBasesBasic::invJacPtr
boost::shared_ptr< MatrixDouble > invJacPtr
Definition: UserDataOperators.hpp:3318
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:1307
MoFEM::OpScaleMatrix::scale
const double scale
Definition: UserDataOperators.hpp:2051
MoFEM::OpCalculateDivergenceVectorFieldValues::OpCalculateDivergenceVectorFieldValues
OpCalculateDivergenceVectorFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Definition: UserDataOperators.hpp:491
MoFEM::OpSetInvJacH1ForFatPrism::invJac
MatrixDouble & invJac
Definition: UserDataOperators.hpp:3623
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::OpCalculateVectorFieldGradient_General
Evaluate field gradient values for vector field, i.e. gradient is tensor rank 2.
Definition: UserDataOperators.hpp:1431
MoFEM::OpSetInvJacH1ForFlatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3662
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:3420
MoFEM::OpCalculateTraceFromSymmMat::i
FTensor::Index< 'i', DIM > i
Definition: UserDataOperators.hpp:3818
MoFEM::OpCalculateHVecTensorDivergence::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:3037
MoFEM::OpCalculateTensor2FieldValues_General::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:792
MoFEM::OpCalculateHVecVectorFieldDot< 3, Field_Dim >::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2212
MoFEM::SmartPetscObj< Vec >
MoFEM::PetscData::CTX_SET_X
@ CTX_SET_X
Definition: LoopMethods.hpp:27
MoFEM::OpCalculateBrokenHVecTensorField::zeroType
const EntityHandle zeroType
Definition: UserDataOperators.hpp:2753
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:1659
MoFEM::OpTensorTimesSymmetricTensor::outMat
boost::shared_ptr< MatrixDouble > outMat
Definition: UserDataOperators.hpp:1969
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::OpSetInvJacH1ForFlatPrism::diffNinvJac
MatrixDouble diffNinvJac
Definition: UserDataOperators.hpp:3674
MoFEM::OpCalculateHVecTensorDivergence::dataVec
SmartPetscObj< Vec > dataVec
Definition: UserDataOperators.hpp:3035
MoFEM::OpCalculateHVecVectorFieldDot
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2183
MoFEM::OpSetInvJacH1ForFaceEmbeddedIn3DSpace
Definition: UserDataOperators.hpp:3367
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:1576
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
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:2043
MoFEM::OpSetContravariantPiolaTransformOnFace2DImpl
Apply contravariant (Piola) transfer to Hdiv space on face.
Definition: UserDataOperators.hpp:3480
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:1297
MoFEM::OpCalculateHcurlVectorCurl
Calculate curl of vector field.
Definition: UserDataOperators.hpp:2544
MoFEM::OpCalculateHVecVectorFieldDot< 3, Field_Dim >::dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Definition: UserDataOperators.hpp:2211
MoFEM::OpCalculateInvJacForFlatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3640
MoFEM::OpSetInvJacSpaceForFaceImpl
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3268
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:1100
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor1CoordsAtGaussPts
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
Definition: ForcesAndSourcesCore.hpp:1270
MoFEM::OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms::OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms
OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms()
Definition: UserDataOperators.hpp:3562
MoFEM::OpCalculateHdivVectorDivergenceDot::zeroSide
const int zeroSide
Definition: UserDataOperators.hpp:2533
MoFEM::OpCalculateBrokenHVecTensorField
Get tensor field for H-div approximation.
Definition: UserDataOperators.hpp:2719
MoFEM::OpScaleMatrix::inMat
boost::shared_ptr< MatrixDouble > inMat
Definition: UserDataOperators.hpp:2052
MoFEM::OpSymmetrizeTensor::outMat
boost::shared_ptr< MatrixDouble > outMat
Definition: UserDataOperators.hpp:2022
MoFEM::OpCalculateBrokenHVecTensorField::OpCalculateBrokenHVecTensorField
OpCalculateBrokenHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, EntityType broken_type, boost::shared_ptr< Range > broken_range_ptr=nullptr, boost::shared_ptr< double > scale_ptr=nullptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Definition: UserDataOperators.hpp:2724
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:2110
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:1771