v0.10.0
UserDataOperators.hpp
Go to the documentation of this file.
1 /** \file UserDataOperators.hpp
2  * \brief User data Operators
3 
4 */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 #ifndef __USER_DATA_OPERATORS_HPP__
21 #define __USER_DATA_OPERATORS_HPP__
22 
23 namespace MoFEM {
24 
25 /** \name Get values at Gauss pts */
26 
27 /**@{*/
28 
29 /** \name Scalar values */
30 
31 /**@{*/
32 
33 /** \brief Scalar field values at integration points
34  *
35  */
36 template <class T, class A>
39 
40  boost::shared_ptr<ublas::vector<T, A>> dataPtr;
42 
44  const std::string field_name,
45  boost::shared_ptr<ublas::vector<T, A>> data_ptr,
46  const EntityType zero_type = MBVERTEX)
48  dataPtr(data_ptr), zeroType(zero_type) {
49  if (!dataPtr)
50  THROW_MESSAGE("Pointer is not set");
51  }
52 
53  /**
54  * \brief calculate values of scalar field at integration points
55  * @param side side entity number
56  * @param type side entity type
57  * @param data entity data
58  * @return error code
59  */
60  MoFEMErrorCode doWork(int side, EntityType type,
62 };
63 
64 /**
65  * \brief Specialization of member function
66  *
67  */
68 template <class T, class A>
70  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
72  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented for T = %s",
73  typeid(T).name() // boost::core::demangle(typeid(T).name()).c_str()
74  );
76 }
77 
78 /**
79  * \brief Get value at integration points for scalar field
80  *
81  * \ingroup mofem_forces_and_sources_user_data_operators
82  */
84  : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
85 
88 
89  /**
90  * \brief calculate values of scalar field at integration points
91  * @param side side entity number
92  * @param type side entity type
93  * @param data entity data
94  * @return error code
95  */
96  MoFEMErrorCode doWork(int side, EntityType type,
99  VectorDouble &vec = *dataPtr;
100  const int nb_gauss_pts = getGaussPts().size2();
101  if (type == zeroType) {
102  vec.resize(nb_gauss_pts, false);
103  vec.clear();
104  }
105  const int nb_dofs = data.getFieldData().size();
106  if (nb_dofs) {
107  const int nb_gauss_pts = getGaussPts().size2();
108  const int nb_base_functions = data.getN().size2();
109  auto base_function = data.getFTensor0N();
110  auto values_at_gauss_pts = getFTensor0FromVec(vec);
111  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
112  auto field_data = data.getFTensor0FieldData();
113  int bb = 0;
114  for (; bb != nb_dofs; ++bb) {
115  values_at_gauss_pts += field_data * base_function;
116  ++field_data;
117  ++base_function;
118  }
119  // It is possible to have more base functions than dofs
120  for (; bb != nb_base_functions; ++bb)
121  ++base_function;
122  ++values_at_gauss_pts;
123  }
124  }
126  }
127 };
128 
129 /**
130  * @brief Get rate of scalar field at integration points
131  *
132  * \ingroup mofem_forces_and_sources_user_data_operators
133  */
136 
137  boost::shared_ptr<VectorDouble> dataPtr;
140 
141  OpCalculateScalarFieldValuesDot(const std::string field_name,
142  boost::shared_ptr<VectorDouble> data_ptr,
143  const EntityType zero_at_type = MBVERTEX)
146  dataPtr(data_ptr), zeroAtType(zero_at_type) {
147  if (!dataPtr)
148  THROW_MESSAGE("Pointer is not set");
149  }
150 
151  MoFEMErrorCode doWork(int side, EntityType type,
154 
155  const int nb_gauss_pts = getGaussPts().size2();
156  VectorDouble &vec = *dataPtr;
157  if (type == zeroAtType) {
158  vec.resize(nb_gauss_pts, false);
159  vec.clear();
160  }
161 
162  auto &local_indices = data.getLocalIndices();
163  const int nb_dofs = local_indices.size();
164  if (nb_dofs) {
165 
166  dotVector.resize(nb_dofs, false);
167  const double *array;
168  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
169  for (int i = 0; i != local_indices.size(); ++i)
170  if (local_indices[i] != -1)
171  dotVector[i] = array[local_indices[i]];
172  else
173  dotVector[i] = 0;
174  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
175 
176  const int nb_base_functions = data.getN().size2();
177  auto base_function = data.getFTensor0N();
178  auto values_at_gauss_pts = getFTensor0FromVec(vec);
179 
180  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
181  auto field_data = getFTensor0FromVec(dotVector);
182  int bb = 0;
183  for (; bb != nb_dofs; ++bb) {
184  values_at_gauss_pts += field_data * base_function;
185  ++field_data;
186  ++base_function;
187  }
188  // Number of dofs can be smaller than number of Tensor_Dim x base
189  // functions
190  for (; bb != nb_base_functions; ++bb)
191  ++base_function;
192  ++values_at_gauss_pts;
193  }
194  }
196  }
197 };
198 
199 /**
200  * \depreacted Name inconstent with other operators
201  *
202  */
204 
205 /**@}*/
206 
207 /** \name Vector field values at integration points */
208 
209 /**@{*/
210 
211 /** \brief Calculate field values for tenor field rank 1, i.e. vector field
212  *
213  */
214 template <int Tensor_Dim, class T, class L, class A>
217 
218  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
220 
222  const std::string field_name,
223  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
224  const EntityType zero_type = MBVERTEX)
227  dataPtr(data_ptr), zeroType(zero_type) {
228  if (!dataPtr)
229  THROW_MESSAGE("Pointer is not set");
230  }
231 
232  /**
233  * \brief calculate values of vector field at integration points
234  * @param side side entity number
235  * @param type side entity type
236  * @param data entity data
237  * @return error code
238  */
239  MoFEMErrorCode doWork(int side, EntityType type,
241 };
242 
243 template <int Tensor_Dim, class T, class L, class A>
246  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
248  SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
249  "Not implemented for T = %s and dim = %d",
250  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
251  Tensor_Dim);
253 }
254 
255 /** \brief Calculate field values (template specialization) for tensor field
256  * rank 1, i.e. vector field
257  *
258  */
259 template <int Tensor_Dim>
261  ublas::row_major, DoubleAllocator>
263 
264  boost::shared_ptr<MatrixDouble> dataPtr;
266 
267  OpCalculateVectorFieldValues_General(const std::string field_name,
268  boost::shared_ptr<MatrixDouble> data_ptr,
269  const EntityType zero_type = MBVERTEX)
272  dataPtr(data_ptr), zeroType(zero_type) {
273  if (!dataPtr)
274  THROW_MESSAGE("Pointer is not set");
275  }
276 
277  MoFEMErrorCode doWork(int side, EntityType type,
279 };
280 
281 /**
282  * \brief Member function specialization calculating values for tenor field rank
283  *
284  */
285 template <int Tensor_Dim>
286 MoFEMErrorCode OpCalculateVectorFieldValues_General<
287  Tensor_Dim, double, ublas::row_major,
288  DoubleAllocator>::doWork(int side, EntityType type,
291 
292  const int nb_dofs = data.getFieldData().size();
293 
294  if (nb_dofs) {
295  const int nb_gauss_pts = getGaussPts().size2();
296  auto &mat = *dataPtr;
297  if (type == zeroType) {
298  mat.resize(Tensor_Dim, nb_gauss_pts, false);
299  mat.clear();
300  }
301  if (nb_gauss_pts) {
302  const int nb_base_functions = data.getN().size2();
303  auto base_function = data.getFTensor0N();
304  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
306  const int size = nb_dofs / Tensor_Dim;
307  if (nb_dofs % Tensor_Dim) {
308  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
309  "Data inconsistency");
310  }
311  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
312  auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
313  int bb = 0;
314  for (; bb != size; ++bb) {
315  values_at_gauss_pts(I) += field_data(I) * base_function;
316  ++field_data;
317  ++base_function;
318  }
319  // Number of dofs can be smaller than number of Tensor_Dim x base
320  // functions
321  for (; bb != nb_base_functions; ++bb)
322  ++base_function;
323  ++values_at_gauss_pts;
324  }
325  }
326  } else if (type == this->zeroType) {
327  dataPtr->resize(Tensor_Dim, 0, false);
328  }
330 }
331 
332 /** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
333  * field
334  *
335  * \ingroup mofem_forces_and_sources_user_data_operators
336  */
337 template <int Tensor_Dim>
340  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
341 
343  Tensor_Dim, double, ublas::row_major,
345 };
346 
347 /** \brief Approximate field valuse for given petsc vector
348  *
349  * \note Look at PetscData to see what vectors could be extarcted with that user
350  * data opetaor.
351  *
352  * \ingroup mofem_forces_and_sources_user_data_operators
353  */
354 template <int Tensor_Dim, PetscData::DataContext CTX>
357 
358  boost::shared_ptr<MatrixDouble> dataPtr;
361 
363  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
364  const EntityType zero_at_type = MBVERTEX)
367  dataPtr(data_ptr), zeroAtType(zero_at_type) {
368  if (!dataPtr)
369  THROW_MESSAGE("Pointer is not set");
370  }
371 
372  MoFEMErrorCode doWork(int side, EntityType type,
375  auto &local_indices = data.getLocalIndices();
376  const int nb_dofs = local_indices.size();
377  if (!nb_dofs && type == zeroAtType) {
378  dataPtr->resize(Tensor_Dim, 0, false);
380  }
381  if (!nb_dofs)
383 
384  const double *array;
385 
386  auto get_array = [&](const auto ctx, auto vec) {
388  if ((getFEMethod()->data_ctx & ctx).none())
389  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
390  CHKERR VecGetArrayRead(vec, &array);
392  };
393 
394  auto restore_array = [&](auto vec) {
395  return VecRestoreArrayRead(vec, &array);
396  };
397 
398  switch (CTX) {
400  CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
401  break;
403  CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
404  break;
406  CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
407  break;
408  default:
409  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
410  "That case is not implemented");
411  }
412 
413  dotVector.resize(local_indices.size());
414  for (int i = 0; i != local_indices.size(); ++i)
415  if (local_indices[i] != -1)
416  dotVector[i] = array[local_indices[i]];
417  else
418  dotVector[i] = 0;
419 
420  switch (CTX) {
422  CHKERR restore_array(getFEMethod()->ts_u);
423  break;
425  CHKERR restore_array(getFEMethod()->ts_u_t);
426  break;
428  CHKERR restore_array(getFEMethod()->ts_u_tt);
429  break;
430  default:
431  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
432  "That case is not implemented");
433  }
434 
435  const int nb_gauss_pts = data.getN().size1();
436  const int nb_base_functions = data.getN().size2();
437  MatrixDouble &mat = *dataPtr;
438  if (type == zeroAtType) {
439  mat.resize(Tensor_Dim, nb_gauss_pts, false);
440  mat.clear();
441  }
442  auto base_function = data.getFTensor0N();
443  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
445  const int size = nb_dofs / Tensor_Dim;
446  if (nb_dofs % Tensor_Dim) {
447  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
448  }
449  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
450  auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(dotVector);
451  int bb = 0;
452  for (; bb != size; ++bb) {
453  values_at_gauss_pts(I) += field_data(I) * base_function;
454  ++field_data;
455  ++base_function;
456  }
457  // Number of dofs can be smaller than number of Tensor_Dim x base
458  // functions
459  for (; bb != nb_base_functions; ++bb)
460  ++base_function;
461  ++values_at_gauss_pts;
462  }
464  }
465 };
466 
467 /** \brief Get time direvatives of values at integration pts for tensor filed
468  * rank 1, i.e. vector field
469  *
470  * \ingroup mofem_forces_and_sources_user_data_operators
471  */
472 template <int Tensor_Dim>
474  OpCalculateVectorFieldValuesFromPetscVecImpl<Tensor_Dim,
476 
477 /** \brief Get second time direvatives of values at integration pts for tensor
478  * filed rank 1, i.e. vector field
479  *
480  * \ingroup mofem_forces_and_sources_user_data_operators
481  */
482 template <int Tensor_Dim>
486 
487 /**@}*/
488 
489 /** \name Tensor field values at integration points */
490 
491 /**@{*/
492 
493 /** \brief Calculate field values for tenor field rank 2.
494  *
495  */
496 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
499 
500  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
502 
504  const std::string field_name,
505  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
506  const EntityType zero_type = MBVERTEX)
509  dataPtr(data_ptr), zeroType(zero_type) {
510  if (!dataPtr)
511  THROW_MESSAGE("Pointer is not set");
512  }
513 
514  MoFEMErrorCode doWork(int side, EntityType type,
516 };
517 
518 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
521  doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
523  SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
524  "Not implemented for T = %s, dim0 = %d and dim1 = %d",
525  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
526  Tensor_Dim0, Tensor_Dim1);
528 }
529 
530 template <int Tensor_Dim0, int Tensor_Dim1>
531 struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
532  ublas::row_major, DoubleAllocator>
534 
535  boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
538 
540  const std::string field_name,
541  boost::shared_ptr<
542  ublas::matrix<double, ublas::row_major, DoubleAllocator>> &data_ptr,
543  const EntityType zero_type = MBVERTEX)
546  dataPtr(data_ptr), zeroType(zero_type) {
547  if (!dataPtr)
548  THROW_MESSAGE("Pointer is not set");
549  }
550 
551  MoFEMErrorCode doWork(int side, EntityType type,
553 };
554 
555 template <int Tensor_Dim0, int Tensor_Dim1>
556 MoFEMErrorCode OpCalculateTensor2FieldValues_General<
557  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
558  DoubleAllocator>::doWork(int side, EntityType type,
561  const int nb_dofs = data.getFieldData().size();
562  const int nb_gauss_pts = data.getN().size1();
563  if (!nb_dofs && type == this->zeroType) {
564  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
566  }
567  if (!nb_dofs) {
569  }
570  const int nb_base_functions = data.getN().size2();
571  MatrixDouble &mat = *dataPtr;
572  if (type == zeroType) {
573  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
574  mat.clear();
575  }
576  auto base_function = data.getFTensor0N();
577  auto values_at_gauss_pts = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
580  const int size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
581  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
582  auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
583  int bb = 0;
584  for (; bb != size; ++bb) {
585  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
586  ++field_data;
587  ++base_function;
588  }
589  for (; bb != nb_base_functions; ++bb)
590  ++base_function;
591  ++values_at_gauss_pts;
592  }
594 }
595 
596 /** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
597  * field
598  *
599  * \ingroup mofem_forces_and_sources_user_data_operators
600  */
601 template <int Tensor_Dim0, int Tensor_Dim1>
604  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
605 
606  OpCalculateTensor2FieldValues(const std::string field_name,
607  boost::shared_ptr<MatrixDouble> data_ptr,
608  const EntityType zero_type = MBVERTEX)
609  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
610  ublas::row_major,
612  field_name, data_ptr, zero_type) {}
613 };
614 
615 /** \brief Get time direvarive values at integration pts for tensor filed rank
616  * 2, i.e. matrix field
617  *
618  * \ingroup mofem_forces_and_sources_user_data_operators
619  */
620 template <int Tensor_Dim0, int Tensor_Dim1>
623 
624  boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
625  EntityType zeroAtType; ///< Zero values at Gauss point at this type
626  VectorDouble dotVector; ///< Keeps temoorary values of time directives
627 
628  template <int Dim0, int Dim1> auto getFTensorDotData() {
629  static_assert(Dim0 || !Dim0 || Dim1 || !Dim1, "not implemented");
630  }
631 
632  OpCalculateTensor2FieldValuesDot(const std::string field_name,
633  boost::shared_ptr<MatrixDouble> data_ptr,
634  const EntityType zero_at_type = MBVERTEX)
637  dataPtr(data_ptr), zeroAtType(zero_at_type) {
638  if (!dataPtr)
639  THROW_MESSAGE("Pointer is not set");
640  }
641 
642  MoFEMErrorCode doWork(int side, EntityType type,
645  const auto &local_indices = data.getLocalIndices();
646  const int nb_dofs = local_indices.size();
647  const int nb_gauss_pts = data.getN().size1();
648  if (!nb_dofs && type == zeroAtType) {
649  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
651  }
652  if (!nb_dofs) {
654  }
655 
656  dotVector.resize(nb_dofs, false);
657  const double *array;
658  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
659  for (int i = 0; i != local_indices.size(); ++i)
660  if (local_indices[i] != -1)
661  dotVector[i] = array[local_indices[i]];
662  else
663  dotVector[i] = 0;
664  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
665 
666  const int nb_base_functions = data.getN().size2();
667  MatrixDouble &mat = *dataPtr;
668  if (type == zeroAtType) {
669  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
670  mat.clear();
671  }
672  auto base_function = data.getFTensor0N();
673  auto values_at_gauss_pts =
674  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
677  const int size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
678  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
679  auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
680  int bb = 0;
681  for (; bb != size; ++bb) {
682  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
683  ++field_data;
684  ++base_function;
685  }
686  for (; bb != nb_base_functions; ++bb)
687  ++base_function;
688  ++values_at_gauss_pts;
689  }
691  }
692 };
693 
694 template <>
695 template <>
698  &dotVector[0], &dotVector[1], &dotVector[2],
699 
700  &dotVector[3], &dotVector[4], &dotVector[5],
701 
702  &dotVector[6], &dotVector[7], &dotVector[8]);
703 }
704 
705 /**
706  * @brief Calculate symmetric tensor field values at integration pts.
707  *
708  * @tparam Tensor_Dim
709 
710  * \ingroup mofem_forces_and_sources_user_data_operators
711  */
712 template <int Tensor_Dim>
715 
716  boost::shared_ptr<MatrixDouble> dataPtr;
718  const int zeroSide;
719 
721  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
722  const EntityType zero_type = MBEDGE, const int zero_side = 0)
725  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
726  if (!dataPtr)
727  THROW_MESSAGE("Pointer is not set");
728  }
729 
730  MoFEMErrorCode doWork(int side, EntityType type,
733  MatrixDouble &mat = *dataPtr;
734  const int nb_gauss_pts = getGaussPts().size2();
735  if (type == this->zeroType && side == zeroSide) {
736  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
737  mat.clear();
738  }
739  const int nb_dofs = data.getFieldData().size();
740  if (!nb_dofs)
742 
743  const int nb_base_functions = data.getN().size2();
744  auto base_function = data.getFTensor0N();
745  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
748  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
749  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
750  auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
751  int bb = 0;
752  for (; bb != size; ++bb) {
753  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
754  ++field_data;
755  ++base_function;
756  }
757  for (; bb != nb_base_functions; ++bb)
758  ++base_function;
759  ++values_at_gauss_pts;
760  }
761 
763  }
764 };
765 
766 /**
767  * @brief Calculate symmetric tensor field rates ant integratio pts.
768  *
769  * @tparam Tensor_Dim
770  *
771  * \ingroup mofem_forces_and_sources_user_data_operators
772  */
773 template <int Tensor_Dim>
776 
777  boost::shared_ptr<MatrixDouble> dataPtr;
779  const int zeroSide;
781 
782  template <int Dim> inline auto getFTensorDotData() {
783  static_assert(Dim || !Dim, "not implemented");
784  }
785 
787  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
788  const EntityType zero_type = MBEDGE, const int zero_side = 0)
791  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
792  if (!dataPtr)
793  THROW_MESSAGE("Pointer is not set");
794  }
795 
796  MoFEMErrorCode doWork(int side, EntityType type,
799  const int nb_gauss_pts = getGaussPts().size2();
800  MatrixDouble &mat = *dataPtr;
801  if (type == zeroType && side == zeroSide) {
802  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
803  mat.clear();
804  }
805  auto &local_indices = data.getLocalIndices();
806  const int nb_dofs = local_indices.size();
807  if (!nb_dofs)
809 
810  dotVector.resize(nb_dofs, false);
811  const double *array;
812  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
813  for (int i = 0; i != local_indices.size(); ++i)
814  if (local_indices[i] != -1)
815  dotVector[i] = array[local_indices[i]];
816  else
817  dotVector[i] = 0;
818  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
819 
820  const int nb_base_functions = data.getN().size2();
821 
822  auto base_function = data.getFTensor0N();
823  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
826  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
827  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
828  auto field_data = getFTensorDotData<Tensor_Dim>();
829  int bb = 0;
830  for (; bb != size; ++bb) {
831  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
832  ++field_data;
833  ++base_function;
834  }
835  for (; bb != nb_base_functions; ++bb)
836  ++base_function;
837  ++values_at_gauss_pts;
838  }
839 
841  }
842 };
843 
844 template <>
845 template <>
846 inline auto
849  &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
850  &dotVector[5]);
851 }
852 
853 template <>
854 template <>
855 inline auto
858  &dotVector[0], &dotVector[1], &dotVector[2]);
859 }
860 
861 /**@}*/
862 
863 /** \name Gradients of scalar fields at integration points */
864 
865 /**@{*/
866 
867 /**
868  * \brief Evaluate field gradient values for scalar field, i.e. gradient is
869  * tensor rank 1 (vector)
870  *
871  */
872 template <int Tensor_Dim, class T, class L, class A>
874  : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
875 
877  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
878  const EntityType zero_type = MBVERTEX)
879  : OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A>(
880  field_name, data_ptr, zero_type) {}
881 };
882 
883 /** \brief Evaluate field gradient values for scalar field, i.e. gradient is
884  * tensor rank 1 (vector), specialization
885  *
886  */
887 template <int Tensor_Dim>
889  ublas::row_major, DoubleAllocator>
891  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
892 
894  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
895  const EntityType zero_type = MBVERTEX)
897  ublas::row_major, DoubleAllocator>(
898  field_name, data_ptr, zero_type) {}
899 
900  /**
901  * \brief calculate gradient values of scalar field at integration points
902  * @param side side entity number
903  * @param type side entity type
904  * @param data entity data
905  * @return error code
906  */
907  MoFEMErrorCode doWork(int side, EntityType type,
909 };
910 
911 /**
912  * \brief Member function specialization calculating scalar field gradients for
913  * tenor field rank 1
914  *
915  */
916 template <int Tensor_Dim>
917 MoFEMErrorCode OpCalculateScalarFieldGradient_General<
918  Tensor_Dim, double, ublas::row_major,
919  DoubleAllocator>::doWork(int side, EntityType type,
922  if (!this->dataPtr)
923  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
924  "Data pointer not allocated");
925 
926  const int nb_dofs = data.getFieldData().size();
927  if (nb_dofs) {
928  const int nb_gauss_pts = this->getGaussPts().size2();
929  auto &mat = *this->dataPtr;
930  if (type == this->zeroType) {
931  mat.resize(Tensor_Dim, nb_gauss_pts, false);
932  mat.clear();
933  }
934  if (nb_gauss_pts) {
935  const int nb_base_functions = data.getN().size2();
936  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
937  auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
938 
940  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
941  auto field_data = data.getFTensor0FieldData();
942  int bb = 0;
943  for (; bb != nb_dofs; ++bb) {
944  gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
945  ++field_data;
946  ++diff_base_function;
947  }
948  // Number of dofs can be smaller than number of base functions
949  for (; bb < nb_base_functions; ++bb)
950  ++diff_base_function;
951  ++gradients_at_gauss_pts;
952  }
953  }
954 
955  } else if (type == this->zeroType) {
956  this->dataPtr->resize(Tensor_Dim, 0, false);
957  }
958 
960 }
961 
962 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
963  * vector field
964  *
965  * \ingroup mofem_forces_and_sources_user_data_operators
966  */
967 template <int Tensor_Dim>
970  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
972  Tensor_Dim, double, ublas::row_major,
974 };
975 
976 /**}*/
977 
978 /** \name Gradients of tensor fields at integration points */
979 
980 /**@{*/
981 
982 /**
983  * \brief Evaluate field gradient values for vector field, i.e. gradient is
984  * tensor rank 2
985  *
986  */
987 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
989  : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
990  L, A> {
991 
993  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
994  const EntityType zero_type = MBVERTEX)
995  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
996  A>(field_name, data_ptr,
997  zero_type) {}
998 };
999 
1000 template <int Tensor_Dim0, int Tensor_Dim1>
1001 struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1002  ublas::row_major, DoubleAllocator>
1004  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1005 
1007  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1008  const EntityType zero_type = MBVERTEX)
1009  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1010  ublas::row_major,
1011  DoubleAllocator>(
1012  field_name, data_ptr, zero_type) {}
1013 
1014  /**
1015  * \brief calculate values of vector field at integration points
1016  * @param side side entity number
1017  * @param type side entity type
1018  * @param data entity data
1019  * @return error code
1020  */
1021  MoFEMErrorCode doWork(int side, EntityType type,
1023 };
1024 
1025 /**
1026  * \brief Member function specialization calculating vector field gradients for
1027  * tenor field rank 2
1028  *
1029  */
1030 template <int Tensor_Dim0, int Tensor_Dim1>
1031 MoFEMErrorCode OpCalculateVectorFieldGradient_General<
1032  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1033  DoubleAllocator>::doWork(int side, EntityType type,
1036  if (!this->dataPtr)
1037  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1038  "Data pointer not allocated");
1039 
1040  const int nb_dofs = data.getFieldData().size();
1041 
1042  if (nb_dofs) {
1043  const int nb_gauss_pts = this->getGaussPts().size2();
1044  auto &mat = *this->dataPtr;
1045  if (type == this->zeroType) {
1046  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1047  mat.clear();
1048  }
1049 
1050  if (nb_gauss_pts) {
1051  const int nb_base_functions = data.getN().size2();
1052  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1053  auto gradients_at_gauss_pts =
1054  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1057  int size = nb_dofs / Tensor_Dim0;
1058  if (nb_dofs % Tensor_Dim0) {
1059  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1060  "Data inconsistency");
1061  }
1062  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1063  auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1064  int bb = 0;
1065  for (; bb < size; ++bb) {
1066  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1067  ++field_data;
1068  ++diff_base_function;
1069  }
1070  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1071  // functions
1072  for (; bb != nb_base_functions; ++bb)
1073  ++diff_base_function;
1074  ++gradients_at_gauss_pts;
1075  }
1076  }
1077 
1078  } else if (type == this->zeroType) {
1079  this->dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
1080  }
1082 }
1083 
1084 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1085  * vector field
1086  *
1087  * \ingroup mofem_forces_and_sources_user_data_operators
1088  */
1089 template <int Tensor_Dim0, int Tensor_Dim1>
1092  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1093 
1094  OpCalculateVectorFieldGradient(const std::string field_name,
1095  boost::shared_ptr<MatrixDouble> data_ptr,
1096  const EntityType zero_type = MBVERTEX)
1097  : OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1098  ublas::row_major,
1099  DoubleAllocator>(
1100  field_name, data_ptr, zero_type) {}
1101 };
1102 
1103 /** \brief Get field gradients time derivative at integration pts for scalar
1104  * filed rank 0, i.e. vector field
1105  *
1106  * \ingroup mofem_forces_and_sources_user_data_operators
1107  */
1108 template <int Tensor_Dim0, int Tensor_Dim1>
1111 
1112  boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1113  EntityType zeroAtType; ///< Zero values at Gauss point at this type
1114  VectorDouble dotVector; ///< Keeps temoorary values of time directives
1115 
1116  template <int Dim> inline auto getFTensorDotData() {
1117  static_assert(Dim || !Dim, "not implemented");
1118  }
1119 
1120  OpCalculateVectorFieldGradientDot(const std::string field_name,
1121  boost::shared_ptr<MatrixDouble> data_ptr,
1122  const EntityType zero_at_type = MBVERTEX)
1125  dataPtr(data_ptr), zeroAtType(zero_at_type) {
1126  if (!dataPtr)
1127  THROW_MESSAGE("Pointer is not set");
1128  }
1129 
1130  MoFEMErrorCode doWork(int side, EntityType type,
1133 
1134  const auto &local_indices = data.getLocalIndices();
1135  const int nb_dofs = local_indices.size();
1136  if (!nb_dofs && type == zeroAtType) {
1137  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
1139  }
1140  if (!nb_dofs)
1142 
1143  dotVector.resize(nb_dofs, false);
1144  const double *array;
1145  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1146  for (int i = 0; i != local_indices.size(); ++i)
1147  if (local_indices[i] != -1)
1148  dotVector[i] = array[local_indices[i]];
1149  else
1150  dotVector[i] = 0;
1151  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1152 
1153  const int nb_gauss_pts = data.getN().size1();
1154  const int nb_base_functions = data.getN().size2();
1155  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat = *dataPtr;
1156  if (type == zeroAtType) {
1157  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1158  mat.clear();
1159  }
1160  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1161  auto gradients_at_gauss_pts =
1162  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1165  int size = nb_dofs / Tensor_Dim0;
1166  if (nb_dofs % Tensor_Dim0) {
1167  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1168  }
1169  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1170  auto field_data = getFTensorDotData<Tensor_Dim0>();
1171  int bb = 0;
1172  for (; bb < size; ++bb) {
1173  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1174  ++field_data;
1175  ++diff_base_function;
1176  }
1177  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1178  // functions
1179  for (; bb != nb_base_functions; ++bb)
1180  ++diff_base_function;
1181  ++gradients_at_gauss_pts;
1182  }
1184  }
1185 };
1186 
1187 template <>
1188 template <>
1191  &dotVector[0], &dotVector[1], &dotVector[2]);
1192 }
1193 
1194 template <>
1195 template <>
1197  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(&dotVector[0],
1198  &dotVector[1]);
1199 }
1200 
1201 /**@}*/
1202 
1203 /** \name Transform tensors and vectors */
1204 
1205 /**@{*/
1206 
1207 /**
1208  * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1209  * \f$
1210  *
1211  * @tparam DIM
1212  *
1213  * \ingroup mofem_forces_and_sources_user_data_operators
1214  */
1215 template <int DIM_01, int DIM_23, int S = 0>
1218 
1221 
1222  OpTensorTimesSymmetricTensor(const std::string field_name,
1223  boost::shared_ptr<MatrixDouble> in_mat,
1224  boost::shared_ptr<MatrixDouble> out_mat,
1225  boost::shared_ptr<MatrixDouble> d_mat)
1226  : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1227  // Only is run for vertices
1228  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1229  if (!inMat)
1230  THROW_MESSAGE("Pointer for in mat is null");
1231  if (!outMat)
1232  THROW_MESSAGE("Pointer for out mat is null");
1233  if (!dMat)
1234  THROW_MESSAGE("Pointer for tensor mat is null");
1235  }
1236 
1237  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1239  const size_t nb_gauss_pts = getGaussPts().size2();
1240  auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
1241  auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
1242  outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1243  auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
1244  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1245  t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1246  ++t_in;
1247  ++t_out;
1248  }
1250  }
1251 
1252 private:
1257 
1258  boost::shared_ptr<MatrixDouble> inMat;
1259  boost::shared_ptr<MatrixDouble> outMat;
1260  boost::shared_ptr<MatrixDouble> dMat;
1261 };
1262 
1263 template <int DIM>
1265 
1268 
1269  OpSymmetrizeTensor(const std::string field_name,
1270  boost::shared_ptr<MatrixDouble> in_mat,
1271  boost::shared_ptr<MatrixDouble> out_mat)
1272  : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1273  // Only is run for vertices
1274  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1275  if (!inMat)
1276  THROW_MESSAGE("Pointer not set for in matrix");
1277  if (!outMat)
1278  THROW_MESSAGE("Pointer not set for in matrix");
1279  }
1280 
1281  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1283  const size_t nb_gauss_pts = getGaussPts().size2();
1284  auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
1285  outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
1286  auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
1287  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1288  t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
1289  ++t_in;
1290  ++t_out;
1291  }
1293  }
1294 
1295 private:
1298  boost::shared_ptr<MatrixDouble> inMat;
1299  boost::shared_ptr<MatrixDouble> outMat;
1300 };
1301 
1303 
1306 
1307  OpScaleMatrix(const std::string field_name, const double scale,
1308  boost::shared_ptr<MatrixDouble> in_mat,
1309  boost::shared_ptr<MatrixDouble> out_mat)
1310  : UserOp(field_name, OPROW), scale(scale), inMat(in_mat),
1311  outMat(out_mat) {
1312  // Only is run for vertices
1313  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1314  if (!inMat)
1315  THROW_MESSAGE("Pointer not set for in matrix");
1316  if (!outMat)
1317  THROW_MESSAGE("Pointer not set for in matrix");
1318  }
1319 
1320  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1322  outMat->resize(inMat->size1(), inMat->size2(), false);
1323  noalias(*outMat) = scale * (*inMat);
1325  }
1326 
1327 private:
1328  const double scale;
1329  boost::shared_ptr<MatrixDouble> inMat;
1330  boost::shared_ptr<MatrixDouble> outMat;
1331 };
1332 
1333 /**@}*/
1334 
1335 /** \name H-div/H-curls (Vectorial bases) values at integration points */
1336 
1337 /**@{*/
1338 
1339 /** \brief Get vector field for H-div approximation
1340  * \ingroup mofem_forces_and_sources_user_data_operators
1341  */
1342 template <int Tensor_Dim0, class T, class L, class A>
1345 
1346  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
1348  const int zeroSide;
1349 
1351  const std::string field_name,
1352  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
1353  const EntityType zero_type = MBEDGE, const int zero_side = 0)
1356  dataPtr(data_ptr), zeroType(zero_type), zeroSide(0) {
1357  if (!dataPtr)
1358  THROW_MESSAGE("Pointer is not set");
1359  }
1360 
1361  /**
1362  * \brief calculate values of vector field at integration points
1363  * @param side side entity number
1364  * @param type side entity type
1365  * @param data entity data
1366  * @return error code
1367  */
1368  MoFEMErrorCode doWork(int side, EntityType type,
1370 };
1371 
1372 template <int Tensor_Dim, class T, class L, class A>
1374  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1376  SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1377  "Not implemented for T = %s and dim = %d",
1378  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
1379  Tensor_Dim);
1381 }
1382 
1383 /** \brief Get vector field for H-div approximation
1384  * \ingroup mofem_forces_and_sources_user_data_operators
1385  */
1386 template <int Tensor_Dim>
1390 
1391  boost::shared_ptr<MatrixDouble> dataPtr;
1393  const int zeroSide;
1394 
1395  OpCalculateHdivVectorField_General(const std::string field_name,
1396  boost::shared_ptr<MatrixDouble> data_ptr,
1397  const EntityType zero_type = MBEDGE,
1398  const int zero_side = 0)
1401  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1402  if (!dataPtr)
1403  THROW_MESSAGE("Pointer is not set");
1404  }
1405 
1406  /**
1407  * \brief Calculate values of vector field at integration points
1408  * @param side side entity number
1409  * @param type side entity type
1410  * @param data entity data
1411  * @return error code
1412  */
1413  MoFEMErrorCode doWork(int side, EntityType type,
1415 };
1416 
1417 template <int Tensor_Dim>
1418 MoFEMErrorCode OpCalculateHdivVectorField_General<
1419  Tensor_Dim, double, ublas::row_major,
1420  DoubleAllocator>::doWork(int side, EntityType type,
1423  const int nb_integration_points = this->getGaussPts().size2();
1424  if (type == zeroType && side == zeroSide) {
1425  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1426  dataPtr->clear();
1427  }
1428  const int nb_dofs = data.getFieldData().size();
1429  if (!nb_dofs)
1431  const int nb_base_functions = data.getN().size2() / 3;
1433  auto t_n_hdiv = data.getFTensor1N<Tensor_Dim>();
1434  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1435  for (int gg = 0; gg != nb_integration_points; ++gg) {
1436  auto t_dof = data.getFTensor0FieldData();
1437  int bb = 0;
1438  for (; bb != nb_dofs; ++bb) {
1439  t_data(i) += t_n_hdiv(i) * t_dof;
1440  ++t_n_hdiv;
1441  ++t_dof;
1442  }
1443  for (; bb != nb_base_functions; ++bb)
1444  ++t_n_hdiv;
1445  ++t_data;
1446  }
1448 }
1449 
1450 /** \brief Get vector field for H-div approximation
1451  *
1452  * \ingroup mofem_forces_and_sources_user_data_operators
1453  */
1454 template <int Tensor_Dim>
1457  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1458 
1459  OpCalculateHdivVectorField(const std::string field_name,
1460  boost::shared_ptr<MatrixDouble> data_ptr,
1461  const EntityType zero_type = MBEDGE,
1462  const int zero_side = 0)
1463  : OpCalculateHdivVectorField_General<Tensor_Dim, double, ublas::row_major,
1464  DoubleAllocator>(
1465  field_name, data_ptr, zero_type, zero_side) {}
1466 };
1467 
1468 /**
1469  * @brief Calculate divergence of vector field
1470  * @ingroup mofem_forces_and_sources_user_data_operators
1471  *
1472  * @tparam Tensor_Dim dimension of space
1473  */
1474 template <int Tensor_Dim1, int Tensor_Dim2>
1477 
1478  boost::shared_ptr<VectorDouble> dataPtr;
1480  const int zeroSide;
1481 
1482  OpCalculateHdivVectorDivergence(const std::string field_name,
1483  boost::shared_ptr<VectorDouble> data_ptr,
1484  const EntityType zero_type = MBEDGE,
1485  const int zero_side = 0)
1488  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1489  if (!dataPtr)
1490  THROW_MESSAGE("Pointer is not set");
1491  }
1492 
1493  MoFEMErrorCode doWork(int side, EntityType type,
1496  const int nb_integration_points = getGaussPts().size2();
1497  if (type == zeroType && side == zeroSide) {
1498  dataPtr->resize(nb_integration_points, false);
1499  dataPtr->clear();
1500  }
1501  const int nb_dofs = data.getFieldData().size();
1502  if (!nb_dofs)
1504  const int nb_base_functions = data.getN().size2() / Tensor_Dim1;
1506  auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
1507  auto t_data = getFTensor0FromVec(*dataPtr);
1508  for (int gg = 0; gg != nb_integration_points; ++gg) {
1509  auto t_dof = data.getFTensor0FieldData();
1510  int bb = 0;
1511  for (; bb != nb_dofs; ++bb) {
1512  double div = 0;
1513  for (auto ii = 0; ii != Tensor_Dim2; ++ii)
1514  div += t_n_diff_hdiv(ii, ii);
1515  t_data += t_dof * div;
1516  ++t_n_diff_hdiv;
1517  ++t_dof;
1518  }
1519  for (; bb != nb_base_functions; ++bb)
1520  ++t_n_diff_hdiv;
1521  ++t_data;
1522  }
1524  }
1525 };
1526 
1527 /**
1528  * @brief Calculate curl of vector field
1529  * @ingroup mofem_forces_and_sources_user_data_operators
1530  *
1531  * @tparam Tensor_Dim dimension of space
1532  */
1533 template <int Tensor_Dim>
1536 
1537  boost::shared_ptr<MatrixDouble> dataPtr;
1539  const int zeroSide;
1540 
1541  OpCalculateHcurlVectorCurl(const std::string field_name,
1542  boost::shared_ptr<MatrixDouble> data_ptr,
1543  const EntityType zero_type = MBEDGE,
1544  const int zero_side = 0)
1547  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1548  if (!dataPtr)
1549  THROW_MESSAGE("Pointer is not set");
1550  }
1551 
1552  MoFEMErrorCode doWork(int side, EntityType type,
1555  const int nb_integration_points = getGaussPts().size2();
1556  if (type == zeroType && side == zeroSide) {
1557  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1558  dataPtr->clear();
1559  }
1560  const int nb_dofs = data.getFieldData().size();
1561  if (!nb_dofs)
1563 
1564  MatrixDouble curl_mat;
1565 
1566  const int nb_base_functions = data.getN().size2() / Tensor_Dim;
1570  auto t_n_diff_hcurl = data.getFTensor2DiffN<Tensor_Dim, Tensor_Dim>();
1571  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1572  for (int gg = 0; gg != nb_integration_points; ++gg) {
1573 
1574  auto t_dof = data.getFTensor0FieldData();
1575  int bb = 0;
1576  for (; bb != nb_dofs; ++bb) {
1577 
1578  t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
1579  ++t_n_diff_hcurl;
1580  ++t_dof;
1581  }
1582 
1583  for (; bb != nb_base_functions; ++bb)
1584  ++t_n_diff_hcurl;
1585  ++t_data;
1586  }
1587 
1589  }
1590 };
1591 
1592 /**
1593  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1594  * \ingroup mofem_forces_and_sources_user_data_operators
1595  *
1596  * @tparam Tensor_Dim0 rank of the filed
1597  * @tparam Tensor_Dim1 dimension of space
1598  */
1599 template <int Tensor_Dim0, int Tensor_Dim1>
1602 
1603  boost::shared_ptr<MatrixDouble> dataPtr;
1605  const int zeroSide;
1606 
1607  OpCalculateHVecTensorField(const std::string field_name,
1608  boost::shared_ptr<MatrixDouble> data_ptr,
1609  const EntityType zero_type = MBEDGE,
1610  const int zero_side = 0)
1613  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1614  if (!dataPtr)
1615  THROW_MESSAGE("Pointer is not set");
1616  }
1617 
1618  MoFEMErrorCode doWork(int side, EntityType type,
1621  const int nb_integration_points = getGaussPts().size2();
1622  if (type == zeroType && side == zeroSide) {
1623  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1624  dataPtr->clear();
1625  }
1626  const int nb_dofs = data.getFieldData().size();
1627  if (nb_dofs) {
1628  const int nb_base_functions = data.getN().size2() / 3;
1631  auto t_n_hvec = data.getFTensor1N<3>();
1632  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1633  for (int gg = 0; gg != nb_integration_points; ++gg) {
1634  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1635  int bb = 0;
1636  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1637  t_data(i, j) += t_dof(i) * t_n_hvec(j);
1638  ++t_n_hvec;
1639  ++t_dof;
1640  }
1641  for (; bb != nb_base_functions; ++bb)
1642  ++t_n_hvec;
1643  ++t_data;
1644  }
1645  }
1647  }
1648 };
1649 
1650 /**
1651  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1652  * \ingroup mofem_forces_and_sources_user_data_operators
1653  *
1654  * @tparam Tensor_Dim0 rank of the filed
1655  * @tparam Tensor_Dim1 dimension of space
1656  */
1657 template <int Tensor_Dim0, int Tensor_Dim1>
1660 
1661  boost::shared_ptr<MatrixDouble> dataPtr;
1663  const int zeroSide;
1664 
1665  OpCalculateHTensorTensorField(const std::string field_name,
1666  boost::shared_ptr<MatrixDouble> data_ptr,
1667  const EntityType zero_type = MBEDGE,
1668  const int zero_side = 0)
1671  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1672  if (!dataPtr)
1673  THROW_MESSAGE("Pointer is not set");
1674  }
1675 
1676  MoFEMErrorCode doWork(int side, EntityType type,
1679  const int nb_integration_points = getGaussPts().size2();
1680  if (type == zeroType && side == zeroSide) {
1681  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1682  dataPtr->clear();
1683  }
1684  const int nb_dofs = data.getFieldData().size();
1685  if (!nb_dofs)
1687  const int nb_base_functions =
1688  data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
1691  auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
1692  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1693  for (int gg = 0; gg != nb_integration_points; ++gg) {
1694  auto t_dof = data.getFTensor0FieldData();
1695  int bb = 0;
1696  for (; bb != nb_dofs; ++bb) {
1697  t_data(i, j) += t_dof * t_n_hten(i, j);
1698  ++t_n_hten;
1699  ++t_dof;
1700  }
1701  for (; bb != nb_base_functions; ++bb)
1702  ++t_n_hten;
1703  ++t_data;
1704  }
1706  }
1707 };
1708 
1709 /**
1710  * @brief Calculate divergence of tonsorial field using vectorial base
1711  * \ingroup mofem_forces_and_sources_user_data_operators
1712  *
1713  * @tparam Tensor_Dim0 rank of the field
1714  * @tparam Tensor_Dim1 dimension of space
1715  */
1716 template <int Tensor_Dim0, int Tensor_Dim1>
1719 
1720  boost::shared_ptr<MatrixDouble> dataPtr;
1722  const int zeroSide;
1723 
1724  OpCalculateHVecTensorDivergence(const std::string field_name,
1725  boost::shared_ptr<MatrixDouble> data_ptr,
1726  const EntityType zero_type = MBEDGE,
1727  const int zero_side = 0)
1730  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1731  if (!dataPtr)
1732  THROW_MESSAGE("Pointer is not set");
1733  }
1734 
1735  MoFEMErrorCode doWork(int side, EntityType type,
1738  const int nb_integration_points = getGaussPts().size2();
1739  if (type == zeroType && side == 0) {
1740  dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
1741  dataPtr->clear();
1742  }
1743  const int nb_dofs = data.getFieldData().size();
1744  if (nb_dofs) {
1745  const int nb_base_functions = data.getN().size2() / 3;
1748  auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
1749  auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
1750  for (int gg = 0; gg != nb_integration_points; ++gg) {
1751  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1752  int bb = 0;
1753  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1754  double div = t_n_diff_hvec(j, j);
1755  t_data(i) += t_dof(i) * div;
1756  ++t_n_diff_hvec;
1757  ++t_dof;
1758  }
1759  for (; bb < nb_base_functions; ++bb)
1760  ++t_n_diff_hvec;
1761  ++t_data;
1762  }
1763  }
1765  }
1766 };
1767 
1768 /**@}*/
1769 
1770 /**@}*/
1771 
1772 /** \name Operators for faces */
1773 
1774 /**@{*/
1775 
1776 /** \brief Calculate jacobian for face element
1777 
1778  It is assumed that face element is XY plane. Applied
1779  only for 2d problems.
1780 
1781  \todo Generalize function for arbitrary face orientation in 3d space
1782 
1783  \ingroup mofem_forces_and_sources_tri_element
1784 
1785 */
1789 
1792 
1793  MoFEMErrorCode doWork(int side, EntityType type,
1795 };
1796 
1797 /** \brief Calculate inverse of jacobian for face element
1798 
1799  It is assumed that face element is XY plane. Applied
1800  only for 2d problems.
1801 
1802  \todo Generalize function for arbitrary face orientation in 3d space
1803 
1804  \ingroup mofem_forces_and_sources_tri_element
1805 
1806 */
1810 
1813  invJac(inv_jac) {}
1814 
1815  MoFEMErrorCode doWork(int side, EntityType type,
1817 };
1818 
1819 /** \brief Transform local reference derivatives of shape functions to global
1820 derivatives
1821 
1822 \ingroup mofem_forces_and_sources_tri_element
1823 
1824 */
1827 
1830  invJac(inv_jac) {}
1831 
1832  MoFEMErrorCode doWork(int side, EntityType type,
1834 
1835 private:
1838 };
1839 
1840 /**
1841  * \brief brief Transform local reference derivatives of shape function to
1842  global derivatives for face
1843 
1844  * \ingroup mofem_forces_and_sources_tri_element
1845  */
1848 
1851  invJac(inv_jac) {}
1852 
1853  MoFEMErrorCode doWork(int side, EntityType type,
1855 
1856 private:
1859 };
1860 
1861 /**
1862  * @brief Make Hdiv space from Hcurl space in 2d
1863  * @ingroup mofem_forces_and_sources_tri_element
1864  */
1867 
1870 
1871  MoFEMErrorCode doWork(int side, EntityType type,
1873 };
1874 
1875 /** \brief Apply contravariant (Piola) transfer to Hdiv space on face
1876  *
1877  * \note Hdiv space is generated by Hcurl space in 2d.
1878  *
1879  * Contravariant Piola transformation
1880  * \f[
1881  * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
1882  * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
1883  * =
1884  * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
1885  * \f]
1886  *
1887  * \ingroup mofem_forces_and_sources
1888  *
1889  */
1892 
1895  }
1896 
1899 
1900  MoFEMErrorCode doWork(int side, EntityType type,
1902 
1903 private:
1905 };
1906 
1907 /**@*/
1908 
1909 /** \name Operators for edges */
1910 
1911 /**@{*/
1912 
1915 
1918 
1919  MoFEMErrorCode doWork(int side, EntityType type,
1921 };
1922 
1923 /**@}*/
1924 
1925 /** \name Operator for fat prisms */
1926 
1927 /**@{*/
1928 
1929 /**
1930  * @brief Operator for fat prism element updating integration weights in the
1931  * volume.
1932  *
1933  * Jacobian on the distorted element is nonconstant. This operator updates
1934  * integration weight on prism to take into account nonconstat jacobian.
1935  *
1936  * \f[
1937  * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
1938  * \pmb\xi} \right\| \right)
1939  * \f]
1940  * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
1941  * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
1942  * element coordinate.
1943  *
1944  */
1947 
1950 
1951  MoFEMErrorCode doWork(int side, EntityType type,
1953 };
1954 
1955 /** \brief Calculate inverse of jacobian for face element
1956 
1957  It is assumed that face element is XY plane. Applied
1958  only for 2d problems.
1959 
1960  FIXME Generalize function for arbitrary face orientation in 3d space
1961  FIXME Calculate to Jacobins for two faces
1962 
1963  \ingroup mofem_forces_and_sources_prism_element
1964 
1965 */
1968 
1969  OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
1971  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
1972 
1975  invJac(inv_jac) {}
1976  MoFEMErrorCode doWork(int side, EntityType type,
1978 
1979 private:
1980  const boost::shared_ptr<MatrixDouble> invJacPtr;
1982 };
1983 
1984 /** \brief Transform local reference derivatives of shape functions to global
1985 derivatives
1986 
1987 FIXME Generalize to curved shapes
1988 FIXME Generalize to case that top and bottom face has different shape
1989 
1990 \ingroup mofem_forces_and_sources_prism_element
1991 
1992 */
1995 
1996  OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
1998  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
1999 
2002  invJac(inv_jac) {}
2003 
2004  MoFEMErrorCode doWork(int side, EntityType type,
2006 
2007 private:
2008  const boost::shared_ptr<MatrixDouble> invJacPtr;
2011 };
2012 
2013 // Flat prism
2014 
2015 /** \brief Calculate inverse of jacobian for face element
2016 
2017  It is assumed that face element is XY plane. Applied
2018  only for 2d problems.
2019 
2020  FIXME Generalize function for arbitrary face orientation in 3d space
2021  FIXME Calculate to Jacobins for two faces
2022 
2023  \ingroup mofem_forces_and_sources_prism_element
2024 
2025 */
2028 
2032  invJacF3(inv_jac_f3) {}
2033  MoFEMErrorCode doWork(int side, EntityType type,
2035 };
2036 
2037 /** \brief Transform local reference derivatives of shape functions to global
2038 derivatives
2039 
2040 FIXME Generalize to curved shapes
2041 FIXME Generalize to case that top and bottom face has different shape
2042 
2043 \ingroup mofem_forces_and_sources_prism_element
2044 
2045 */
2051  invJacF3(inv_jac_f3) {}
2052 
2054  MoFEMErrorCode doWork(int side, EntityType type,
2056 };
2057 
2058 /**@}*/
2059 
2060 } // namespace MoFEM
2061 
2062 #endif // __USER_DATA_OPERATORS_HPP__
2063 
2064 /**
2065  * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
2066  *
2067  * \brief Classes and functions used to evaluate fields at integration pts,
2068  *jacobians, etc..
2069  *
2070  * \ingroup mofem_forces_and_sources
2071  **/
const double T
Calculate divergence of vector field.
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
FTensor::Index< 'i', DIM > i
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A >> data_ptr, const EntityType zero_type=MBVERTEX)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::shared_ptr< MatrixDouble > dataPtr
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:143
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static constexpr Switches CtxSetX_T
Definition: LoopMethods.hpp:66
OpCalculateTensor2SymmetricFieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(Tensor_Dim *(Tensor_Dim+1))/2 >, Tensor_Dim > getFTensor2SymmetricFieldData()
Return symmetric FTensor rank 2, i.e. matrix from filed data coeffinects.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpTensorTimesSymmetricTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
OpSymmetrizeTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
OpCalculateScalarFieldValuesDot(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
OpCalculateHcurlVectorCurl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateVectorFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate field values for tenor field rank 1, i.e. vector field.
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
Operator for fat prism element updating integration weights in the volume.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
const boost::shared_ptr< MatrixDouble > invJacPtr
OpCalculateInvJacForFlatPrism(MatrixDouble &inv_jac_f3)
OpCalculateHdivVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A >> data_ptr, const EntityType zero_type=MBVERTEX)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
EntityType zeroAtType
Zero values at Gauss point at this type.
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:71
Approximate field valuse for given petsc vector.
OpCalculateScalarFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1ForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
FTensor::Index< 'l', DIM_23 > l
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
boost::shared_ptr< MatrixDouble > dMat
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:70
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator >> &data_ptr, const EntityType zero_type=MBVERTEX)
Get values at integration pts for tensor filed rank 1, i.e. vector field.
boost::shared_ptr< MatrixDouble > inMat
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Evaluate field gradient values for vector field, i.e. gradient is tensor rank 2.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate values of vector field at integration points
FTensor::Index< 'i', DIM_01 > i
EntityType zeroAtType
Zero values at Gauss point at this type.
boost::shared_ptr< MatrixDouble > dataPtr
Calculate curl of vector field.
OpCalculateHdivVectorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate values of vector field at integration points
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
boost::shared_ptr< MatrixDouble > outMat
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
OpSetInvJacHcurlFace(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateVectorFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
FTensor::Index< 'j', DIM > j
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
Get vector field for H-div approximation.
OpSetInvJacH1ForFatPrism(MatrixDouble &inv_jac)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
FTensor::Index< 'j', DIM_01 > j
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateJacForFace(MatrixDouble &jac)
Calculate inverse of jacobian for face element.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
static Index< 'i', 3 > i
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
OpCalculateInvJacForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
OpCalculateHdivVectorDivergence(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
DataForcesAndSourcesCore::EntData EntData
boost::shared_ptr< MatrixDouble > outMat
boost::shared_ptr< MatrixDouble > dataPtr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
boost::shared_ptr< VectorDouble > dataPtr
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from filed data coeffinects.
boost::shared_ptr< MatrixDouble > inMat
VectorDouble dotVector
Keeps temoorary values of time directives.
static Index< 'j', 3 > j
Calculate inverse of jacobian for face element.
Calculate inverse of jacobian for face element.
OpSetInvJacH1ForFlatPrism(MatrixDouble &inv_jac_f3)
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FlatPrism finite elementUser is implementing own operator at Gauss points level, by own object derive...
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
field with continuous tangents
Definition: definitions.h:178
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A >> data_ptr, const EntityType zero_type=MBVERTEX)
static constexpr Switches CtxSetX_TT
Definition: LoopMethods.hpp:67
Get vector field for H-div approximation.
Calculate divergence of tonsorial field using vectorial base.
Calculate field values for tenor field rank 2.
OpCalculateVectorFieldValuesFromPetscVecImpl< Tensor_Dim, PetscData::CTX_SET_X_T > OpCalculateVectorFieldValuesDot
Get time direvatives of values at integration pts for tensor filed rank 1, i.e. vector field.
Transform local reference derivatives of shape functions to global derivatives.
OpCalculateTensor2FieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Transform local reference derivatives of shape functions to global derivatives.
Transform local reference derivatives of shape functions to global derivatives.
OpScaleMatrix(const std::string field_name, const double scale, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Index< 'L', 3 > L
FTensor::Index< 'k', DIM_23 > k
OpCalculateVectorFieldGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateInvJacForFatPrism(MatrixDouble &inv_jac)
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
OpCalculateVectorFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< MatrixDouble > inMat
OpSetInvJacH1ForFace(MatrixDouble &inv_jac)
brief Transform local reference derivatives of shape function to global derivatives for face
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateInvJacForFace(MatrixDouble &inv_jac)
static Index< 'k', 3 > k
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Scalar field values at integration points.
OpCalculateHdivVectorField_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A >> data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Data on single entity (This is passed as argument to DataOperator::doWork)
Apply contravariant (Piola) transfer to Hdiv space on face.
structure to get information form mofem into DataForcesAndSourcesCore
Get rate of scalar field at integration points.
FatPrism finite elementUser is implementing own operator at Gauss points level, by own object derived...
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > outMat
Calculate jacobian for face element.
continuous field
Definition: definitions.h:177
boost::shared_ptr< MatrixDouble > dataPtr
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
Calculate symmetric tensor field rates ant integratio pts.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Get time direvarive values at integration pts for tensor filed rank 2, i.e. matrix field.
const VectorDouble & getFieldData() const
get dofs values
OpCalculateScalarFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
std::vector< double, std::allocator< double > > DoubleAllocator
Definition: Types.hpp:70
Calculate symmetric tensor field values at integration pts.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate values of scalar field at integration points
boost::shared_ptr< ublas::vector< T, A > > dataPtr
OpCalculateVectorFieldGradientDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
Get value at integration points for scalar field.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static constexpr Switches CtxSetX
Definition: LoopMethods.hpp:65
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
OpCalculateTensor2FieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
VectorDouble dotVector
Keeps temoorary values of time directives.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const boost::shared_ptr< MatrixDouble > invJacPtr
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate values of scalar field at integration points
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Make Hdiv space from Hcurl space in 2d.
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.