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 
41  const std::string field_name,
42  boost::shared_ptr<ublas::vector<T, A>> data_ptr,
43  const EntityType zero_type = MBVERTEX)
45  dataPtr(data_ptr), zeroType(zero_type) {
46  if (!dataPtr)
47  THROW_MESSAGE("Pointer is not set");
48  }
49 
50  /**
51  * \brief calculate values of scalar field at integration points
52  * @param side side entity number
53  * @param type side entity type
54  * @param data entity data
55  * @return error code
56  */
57  MoFEMErrorCode doWork(int side, EntityType type,
59 
60 protected:
61  boost::shared_ptr<ublas::vector<T, A>> dataPtr;
62  const EntityHandle zeroType;
63 };
64 
65 /**
66  * \brief Specialization of member function
67  *
68  */
69 template <class T, class A>
71  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
73  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented for T = %s",
74  typeid(T).name() // boost::core::demangle(typeid(T).name()).c_str()
75  );
77 }
78 
79 /**
80  * \brief Get value at integration points for scalar field
81  *
82  * \ingroup mofem_forces_and_sources_user_data_operators
83  */
85  : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
86 
89 
90  /**
91  * \brief calculate values of scalar field at integration points
92  * @param side side entity number
93  * @param type side entity type
94  * @param data entity data
95  * @return error code
96  */
97  MoFEMErrorCode doWork(int side, EntityType type,
100  VectorDouble &vec = *dataPtr;
101  const size_t nb_gauss_pts = getGaussPts().size2();
102  if (type == zeroType) {
103  vec.resize(nb_gauss_pts, false);
104  vec.clear();
105  }
106  const size_t nb_dofs = data.getFieldData().size();
107  if (nb_dofs) {
108  const size_t nb_base_functions = data.getN().size2();
109  auto base_function = data.getFTensor0N();
110  auto values_at_gauss_pts = getFTensor0FromVec(vec);
111  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
112  auto field_data = data.getFTensor0FieldData();
113  size_t 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  OpCalculateScalarFieldValuesDot(const std::string field_name,
138  boost::shared_ptr<VectorDouble> data_ptr,
139  const EntityType zero_at_type = MBVERTEX)
142  dataPtr(data_ptr), zeroAtType(zero_at_type) {
143  if (!dataPtr)
144  THROW_MESSAGE("Pointer is not set");
145  }
146 
147  MoFEMErrorCode doWork(int side, EntityType type,
150 
151  const size_t nb_gauss_pts = getGaussPts().size2();
152  VectorDouble &vec = *dataPtr;
153  if (type == zeroAtType) {
154  vec.resize(nb_gauss_pts, false);
155  vec.clear();
156  }
157 
158  auto &local_indices = data.getLocalIndices();
159  const size_t nb_dofs = local_indices.size();
160  if (nb_dofs) {
161 
162  std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
163  const double *array;
164  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
165  for (size_t i = 0; i != local_indices.size(); ++i)
166  if (local_indices[i] != -1)
167  dot_dofs_vector[i] = array[local_indices[i]];
168  else
169  dot_dofs_vector[i] = 0;
170  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
171 
172  const size_t nb_base_functions = data.getN().size2();
173  auto base_function = data.getFTensor0N();
174  auto values_at_gauss_pts = getFTensor0FromVec(vec);
175 
176  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
177  size_t bb = 0;
178  for (; bb != nb_dofs; ++bb) {
179  values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
180  ++base_function;
181  }
182  // Number of dofs can be smaller than number of Tensor_Dim x base
183  // functions
184  for (; bb != nb_base_functions; ++bb)
185  ++base_function;
186  ++values_at_gauss_pts;
187  }
188  }
190  }
191 
192 private:
193  boost::shared_ptr<VectorDouble> dataPtr;
194  const EntityHandle zeroAtType;
195 };
196 
197 /**
198  * \depreacted Name inconstent with other operators
199  *
200  */
202 
203 /**@}*/
204 
205 /** \name Vector field values at integration points */
206 
207 /**@{*/
208 
209 /** \brief Calculate field values for tenor field rank 1, i.e. vector field
210  *
211  */
212 template <int Tensor_Dim, class T, class L, class A>
215 
217  const std::string field_name,
218  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
219  const EntityType zero_type = MBVERTEX)
222  dataPtr(data_ptr), zeroType(zero_type) {
223  if (!dataPtr)
224  THROW_MESSAGE("Pointer is not set");
225  }
226 
227  /**
228  * \brief calculate values of vector field at integration points
229  * @param side side entity number
230  * @param type side entity type
231  * @param data entity data
232  * @return error code
233  */
234  MoFEMErrorCode doWork(int side, EntityType type,
236 
237 protected:
238  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
239  const EntityHandle zeroType;
240 };
241 
242 template <int Tensor_Dim, class T, class L, class A>
245  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
247  SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
248  "Not implemented for T = %s and dim = %d",
249  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
250  Tensor_Dim);
252 }
253 
254 /** \brief Calculate field values (template specialization) for tensor field
255  * rank 1, i.e. vector field
256  *
257  */
258 template <int Tensor_Dim>
259 struct OpCalculateVectorFieldValues_General<Tensor_Dim, double,
260  ublas::row_major, DoubleAllocator>
262 
263  OpCalculateVectorFieldValues_General(const std::string field_name,
264  boost::shared_ptr<MatrixDouble> data_ptr,
265  const EntityType zero_type = MBVERTEX)
268  dataPtr(data_ptr), zeroType(zero_type) {
269  if (!dataPtr)
270  THROW_MESSAGE("Pointer is not set");
271  }
272 
273  MoFEMErrorCode doWork(int side, EntityType type,
275 
276 protected:
277  boost::shared_ptr<MatrixDouble> dataPtr;
278  const EntityHandle zeroType;
279 };
280 
281 /**
282  * \brief Member function specialization calculating values for tenor field rank
283  *
284  */
285 template <int Tensor_Dim>
287  Tensor_Dim, double, ublas::row_major,
288  DoubleAllocator>::doWork(int side, EntityType type,
291 
292  const size_t nb_gauss_pts = getGaussPts().size2();
293  auto &mat = *dataPtr;
294  if (type == zeroType) {
295  mat.resize(Tensor_Dim, nb_gauss_pts, false);
296  mat.clear();
297  }
298 
299  const size_t nb_dofs = data.getFieldData().size();
300  if (nb_dofs) {
301 
302  if (nb_gauss_pts) {
303  const size_t nb_base_functions = data.getN().size2();
304  auto base_function = data.getFTensor0N();
305  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
307  const size_t size = nb_dofs / Tensor_Dim;
308  if (nb_dofs % Tensor_Dim) {
309  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
310  "Data inconsistency");
311  }
312  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
313  auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
314  size_t bb = 0;
315  for (; bb != size; ++bb) {
316  values_at_gauss_pts(I) += field_data(I) * base_function;
317  ++field_data;
318  ++base_function;
319  }
320  // Number of dofs can be smaller than number of Tensor_Dim x base
321  // functions
322  for (; bb != nb_base_functions; ++bb)
323  ++base_function;
324  ++values_at_gauss_pts;
325  }
326  }
327  }
329 }
330 
331 /** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
332  * field
333  *
334  * \ingroup mofem_forces_and_sources_user_data_operators
335  */
336 template <int Tensor_Dim>
339  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
340 
342  Tensor_Dim, double, ublas::row_major,
344 };
345 
346 /** \brief Approximate field valuse for given petsc vector
347  *
348  * \note Look at PetscData to see what vectors could be extarcted with that user
349  * data opetaor.
350  *
351  * \ingroup mofem_forces_and_sources_user_data_operators
352  */
353 template <int Tensor_Dim, PetscData::DataContext CTX>
356 
358  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
359  const EntityType zero_at_type = MBVERTEX)
362  dataPtr(data_ptr), zeroAtType(zero_at_type) {
363  if (!dataPtr)
364  THROW_MESSAGE("Pointer is not set");
365  }
366 
367  MoFEMErrorCode doWork(int side, EntityType type,
370  auto &local_indices = data.getLocalIndices();
371  const int nb_dofs = local_indices.size();
372  if (!nb_dofs && type == zeroAtType) {
373  dataPtr->resize(Tensor_Dim, 0, false);
375  }
376  if (!nb_dofs)
378 
379  const double *array;
380 
381  auto get_array = [&](const auto ctx, auto vec) {
383  if ((getFEMethod()->data_ctx & ctx).none())
384  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
385  CHKERR VecGetArrayRead(vec, &array);
387  };
388 
389  auto restore_array = [&](auto vec) {
390  return VecRestoreArrayRead(vec, &array);
391  };
392 
393  switch (CTX) {
395  CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
396  break;
398  CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
399  break;
401  CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
402  break;
403  default:
404  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
405  "That case is not implemented");
406  }
407 
408  dotVector.resize(local_indices.size());
409  for (int i = 0; i != local_indices.size(); ++i)
410  if (local_indices[i] != -1)
411  dotVector[i] = array[local_indices[i]];
412  else
413  dotVector[i] = 0;
414 
415  switch (CTX) {
417  CHKERR restore_array(getFEMethod()->ts_u);
418  break;
420  CHKERR restore_array(getFEMethod()->ts_u_t);
421  break;
423  CHKERR restore_array(getFEMethod()->ts_u_tt);
424  break;
425  default:
426  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
427  "That case is not implemented");
428  }
429 
430  const size_t nb_gauss_pts = data.getN().size1();
431  const size_t nb_base_functions = data.getN().size2();
432  MatrixDouble &mat = *dataPtr;
433  if (type == zeroAtType) {
434  mat.resize(Tensor_Dim, nb_gauss_pts, false);
435  mat.clear();
436  }
437  auto base_function = data.getFTensor0N();
438  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
440  const size_t size = nb_dofs / Tensor_Dim;
441  if (nb_dofs % Tensor_Dim) {
442  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
443  }
444  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
445  auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(dotVector);
446  size_t bb = 0;
447  for (; bb != size; ++bb) {
448  values_at_gauss_pts(I) += field_data(I) * base_function;
449  ++field_data;
450  ++base_function;
451  }
452  // Number of dofs can be smaller than number of Tensor_Dim x base
453  // functions
454  for (; bb != nb_base_functions; ++bb)
455  ++base_function;
456  ++values_at_gauss_pts;
457  }
459  }
460 
461 protected:
462  boost::shared_ptr<MatrixDouble> dataPtr;
463  const EntityHandle zeroAtType;
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>
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 
501  const std::string field_name,
502  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
503  const EntityType zero_type = MBVERTEX)
506  dataPtr(data_ptr), zeroType(zero_type) {
507  if (!dataPtr)
508  THROW_MESSAGE("Pointer is not set");
509  }
510 
511  MoFEMErrorCode doWork(int side, EntityType type,
513 
514 protected:
515  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
516  const EntityHandle zeroType;
517 };
518 
519 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
522  doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
524  SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
525  "Not implemented for T = %s, dim0 = %d and dim1 = %d",
526  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
527  Tensor_Dim0, Tensor_Dim1);
529 }
530 
531 template <int Tensor_Dim0, int Tensor_Dim1>
532 struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
533  ublas::row_major, DoubleAllocator>
535 
537  const std::string field_name,
538  boost::shared_ptr<
539  ublas::matrix<double, ublas::row_major, DoubleAllocator>> &data_ptr,
540  const EntityType zero_type = MBVERTEX)
543  dataPtr(data_ptr), zeroType(zero_type) {
544  if (!dataPtr)
545  THROW_MESSAGE("Pointer is not set");
546  }
547 
548  MoFEMErrorCode doWork(int side, EntityType type,
550 
551 protected:
552  boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
554  const EntityHandle zeroType;
555 };
556 
557 template <int Tensor_Dim0, int Tensor_Dim1>
559  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
560  DoubleAllocator>::doWork(int side, EntityType type,
563 
564  MatrixDouble &mat = *dataPtr;
565  const size_t nb_gauss_pts = data.getN().size1();
566  if (type == zeroType) {
567  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
568  mat.clear();
569  }
570 
571  const size_t nb_dofs = data.getFieldData().size();
572  if (nb_dofs) {
573  const size_t nb_base_functions = data.getN().size2();
574  auto base_function = data.getFTensor0N();
575  auto values_at_gauss_pts =
576  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
579  const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
580  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
581  auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
582  size_t bb = 0;
583  for (; bb != size; ++bb) {
584  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
585  ++field_data;
586  ++base_function;
587  }
588  for (; bb != nb_base_functions; ++bb)
589  ++base_function;
590  ++values_at_gauss_pts;
591  }
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  OpCalculateTensor2FieldValuesDot(const std::string field_name,
625  boost::shared_ptr<MatrixDouble> data_ptr,
626  const EntityType zero_at_type = MBVERTEX)
629  dataPtr(data_ptr), zeroAtType(zero_at_type) {
630  if (!dataPtr)
631  THROW_MESSAGE("Pointer is not set");
632  }
633 
634  MoFEMErrorCode doWork(int side, EntityType type,
637 
638  const size_t nb_gauss_pts = getGaussPts().size2();
639  MatrixDouble &mat = *dataPtr;
640  if (type == zeroAtType) {
641  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
642  mat.clear();
643  }
644  const auto &local_indices = data.getLocalIndices();
645  const size_t nb_dofs = local_indices.size();
646  if (nb_dofs) {
647  dotVector.resize(nb_dofs, false);
648  const double *array;
649  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
650  for (size_t i = 0; i != local_indices.size(); ++i)
651  if (local_indices[i] != -1)
652  dotVector[i] = array[local_indices[i]];
653  else
654  dotVector[i] = 0;
655  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
656 
657  const size_t nb_base_functions = data.getN().size2();
658 
659  auto base_function = data.getFTensor0N();
660  auto values_at_gauss_pts =
661  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
664  const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
665  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
666  auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
667  size_t bb = 0;
668  for (; bb != size; ++bb) {
669  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
670  ++field_data;
671  ++base_function;
672  }
673  for (; bb != nb_base_functions; ++bb)
674  ++base_function;
675  ++values_at_gauss_pts;
676  }
677  }
679  }
680 
681 protected:
682  boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
683  EntityType zeroAtType; ///< Zero values at Gauss point at this type
684  VectorDouble dotVector; ///< Keeps temoorary values of time directives
685 
686  template <int Dim0, int Dim1> auto getFTensorDotData() {
687  static_assert(Dim0 || !Dim0 || Dim1 || !Dim1, "not implemented");
688  }
689 };
690 
691 template <>
692 template <>
695  &dotVector[0], &dotVector[1], &dotVector[2],
696 
697  &dotVector[3], &dotVector[4], &dotVector[5],
698 
699  &dotVector[6], &dotVector[7], &dotVector[8]);
700 }
701 
702 /**
703  * @brief Calculate symmetric tensor field values at integration pts.
704  *
705  * @tparam Tensor_Dim
706 
707  * \ingroup mofem_forces_and_sources_user_data_operators
708  */
709 template <int Tensor_Dim>
712 
714  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
715  const EntityType zero_type = MBEDGE, const int zero_side = 0)
718  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
719  if (!dataPtr)
720  THROW_MESSAGE("Pointer is not set");
721  }
722 
723  MoFEMErrorCode doWork(int side, EntityType type,
726  MatrixDouble &mat = *dataPtr;
727  const int nb_gauss_pts = getGaussPts().size2();
728  if (type == this->zeroType && side == zeroSide) {
729  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
730  mat.clear();
731  }
732  const int nb_dofs = data.getFieldData().size();
733  if (!nb_dofs)
735 
736  const int nb_base_functions = data.getN().size2();
737  auto base_function = data.getFTensor0N();
738  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
741  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
742  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
743  auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
744  int bb = 0;
745  for (; bb != size; ++bb) {
746  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
747  ++field_data;
748  ++base_function;
749  }
750  for (; bb != nb_base_functions; ++bb)
751  ++base_function;
752  ++values_at_gauss_pts;
753  }
754 
756  }
757 
758 protected:
759  boost::shared_ptr<MatrixDouble> dataPtr;
760  const EntityHandle zeroType;
761  const int zeroSide;
762 };
763 
764 /**
765  * @brief Calculate symmetric tensor field rates ant integratio pts.
766  *
767  * @tparam Tensor_Dim
768  *
769  * \ingroup mofem_forces_and_sources_user_data_operators
770  */
771 template <int Tensor_Dim>
774 
776  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
777  const EntityType zero_type = MBEDGE, const int zero_side = 0)
780  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
781  if (!dataPtr)
782  THROW_MESSAGE("Pointer is not set");
783  }
784 
785  MoFEMErrorCode doWork(int side, EntityType type,
788  const int nb_gauss_pts = getGaussPts().size2();
789  MatrixDouble &mat = *dataPtr;
790  if (type == zeroType && side == zeroSide) {
791  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
792  mat.clear();
793  }
794  auto &local_indices = data.getLocalIndices();
795  const int nb_dofs = local_indices.size();
796  if (!nb_dofs)
798 
799  dotVector.resize(nb_dofs, false);
800  const double *array;
801  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
802  for (int i = 0; i != local_indices.size(); ++i)
803  if (local_indices[i] != -1)
804  dotVector[i] = array[local_indices[i]];
805  else
806  dotVector[i] = 0;
807  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
808 
809  const int nb_base_functions = data.getN().size2();
810 
811  auto base_function = data.getFTensor0N();
812  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
815  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
816  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
817  auto field_data = getFTensorDotData<Tensor_Dim>();
818  int bb = 0;
819  for (; bb != size; ++bb) {
820  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
821  ++field_data;
822  ++base_function;
823  }
824  for (; bb != nb_base_functions; ++bb)
825  ++base_function;
826  ++values_at_gauss_pts;
827  }
828 
830  }
831 
832 protected:
833  boost::shared_ptr<MatrixDouble> dataPtr;
834  const EntityHandle zeroType;
835  const int zeroSide;
837 
838  template <int Dim> inline auto getFTensorDotData() {
839  static_assert(Dim || !Dim, "not implemented");
840  }
841 };
842 
843 template <>
844 template <>
845 inline auto
848  &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
849  &dotVector[5]);
850 }
851 
852 template <>
853 template <>
854 inline auto
857  &dotVector[0], &dotVector[1], &dotVector[2]);
858 }
859 
860 /**@}*/
861 
862 /** \name Gradients of scalar fields at integration points */
863 
864 /**@{*/
865 
866 /**
867  * \brief Evaluate field gradient values for scalar field, i.e. gradient is
868  * tensor rank 1 (vector)
869  *
870  */
871 template <int Tensor_Dim, class T, class L, class A>
873  : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
874 
876  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
877  const EntityType zero_type = MBVERTEX)
878  : OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A>(
879  field_name, data_ptr, zero_type) {}
880 };
881 
882 /** \brief Evaluate field gradient values for scalar field, i.e. gradient is
883  * tensor rank 1 (vector), specialization
884  *
885  */
886 template <int Tensor_Dim>
887 struct OpCalculateScalarFieldGradient_General<Tensor_Dim, double,
888  ublas::row_major, DoubleAllocator>
890  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
891 
893  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
894  const EntityType zero_type = MBVERTEX)
895  : OpCalculateVectorFieldValues_General<Tensor_Dim, double,
896  ublas::row_major, DoubleAllocator>(
897  field_name, data_ptr, zero_type) {
898  if (!this->dataPtr)
899  THROW_MESSAGE("Data pointer not allocated");
900  }
901 
902  /**
903  * \brief calculate gradient values of scalar field at integration points
904  * @param side side entity number
905  * @param type side entity type
906  * @param data entity data
907  * @return error code
908  */
909  MoFEMErrorCode doWork(int side, EntityType type,
911 };
912 
913 /**
914  * \brief Member function specialization calculating scalar field gradients for
915  * tenor field rank 1
916  *
917  */
918 template <int Tensor_Dim>
919 MoFEMErrorCode OpCalculateScalarFieldGradient_General<
920  Tensor_Dim, double, ublas::row_major,
921  DoubleAllocator>::doWork(int side, EntityType type,
924 
925  const size_t nb_gauss_pts = this->getGaussPts().size2();
926  auto &mat = *this->dataPtr;
927  if (type == this->zeroType) {
928  mat.resize(Tensor_Dim, nb_gauss_pts, false);
929  mat.clear();
930  }
931 
932  const int nb_dofs = data.getFieldData().size();
933  if (nb_dofs) {
934 
935  if (nb_gauss_pts) {
936  const int nb_base_functions = data.getN().size2();
937  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
938  auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
939 
941  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
942  auto field_data = data.getFTensor0FieldData();
943  int bb = 0;
944  for (; bb != nb_dofs; ++bb) {
945  gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
946  ++field_data;
947  ++diff_base_function;
948  }
949  // Number of dofs can be smaller than number of base functions
950  for (; bb < nb_base_functions; ++bb)
951  ++diff_base_function;
952  ++gradients_at_gauss_pts;
953  }
954  }
955  }
956 
958 }
959 
960 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
961  * vector field
962  *
963  * \ingroup mofem_forces_and_sources_user_data_operators
964  */
965 template <int Tensor_Dim>
968  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
970  Tensor_Dim, double, ublas::row_major,
972 };
973 
974 /**}*/
975 
976 /** \name Gradients of tensor fields at integration points */
977 
978 /**@{*/
979 
980 /**
981  * \brief Evaluate field gradient values for vector field, i.e. gradient is
982  * tensor rank 2
983  *
984  */
985 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
987  : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
988  L, A> {
989 
991  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
992  const EntityType zero_type = MBVERTEX)
993  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
994  A>(field_name, data_ptr,
995  zero_type) {}
996 };
997 
998 template <int Tensor_Dim0, int Tensor_Dim1>
999 struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1000  ublas::row_major, DoubleAllocator>
1002  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1003 
1005  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1006  const EntityType zero_type = MBVERTEX)
1007  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1008  ublas::row_major,
1009  DoubleAllocator>(
1010  field_name, data_ptr, zero_type) {}
1011 
1012  /**
1013  * \brief calculate values of vector field at integration points
1014  * @param side side entity number
1015  * @param type side entity type
1016  * @param data entity data
1017  * @return error code
1018  */
1019  MoFEMErrorCode doWork(int side, EntityType type,
1021 };
1022 
1023 /**
1024  * \brief Member function specialization calculating vector field gradients for
1025  * tenor field rank 2
1026  *
1027  */
1028 template <int Tensor_Dim0, int Tensor_Dim1>
1029 MoFEMErrorCode OpCalculateVectorFieldGradient_General<
1030  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1031  DoubleAllocator>::doWork(int side, EntityType type,
1034  if (!this->dataPtr)
1035  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1036  "Data pointer not allocated");
1037 
1038  const size_t nb_gauss_pts = this->getGaussPts().size2();
1039  auto &mat = *this->dataPtr;
1040  if (type == this->zeroType) {
1041  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1042  mat.clear();
1043  }
1044 
1045  if (nb_gauss_pts) {
1046  const size_t nb_dofs = data.getFieldData().size();
1047 
1048  if (nb_dofs) {
1049 
1050  const int nb_base_functions = data.getN().size2();
1051  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1052  auto gradients_at_gauss_pts =
1053  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1056  int size = nb_dofs / Tensor_Dim0;
1057  if (nb_dofs % Tensor_Dim0) {
1058  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1059  "Data inconsistency");
1060  }
1061  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1062  auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1063  int bb = 0;
1064  for (; bb < size; ++bb) {
1065  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1066  ++field_data;
1067  ++diff_base_function;
1068  }
1069  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1070  // functions
1071  for (; bb != nb_base_functions; ++bb)
1072  ++diff_base_function;
1073  ++gradients_at_gauss_pts;
1074  }
1075  }
1076  }
1078 }
1079 
1080 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1081  * vector field
1082  *
1083  * \ingroup mofem_forces_and_sources_user_data_operators
1084  */
1085 template <int Tensor_Dim0, int Tensor_Dim1>
1088  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1089 
1090  OpCalculateVectorFieldGradient(const std::string field_name,
1091  boost::shared_ptr<MatrixDouble> data_ptr,
1092  const EntityType zero_type = MBVERTEX)
1093  : OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1094  ublas::row_major,
1095  DoubleAllocator>(
1096  field_name, data_ptr, zero_type) {}
1097 };
1098 
1099 /** \brief Get field gradients time derivative at integration pts for scalar
1100  * filed rank 0, i.e. vector field
1101  *
1102  * \ingroup mofem_forces_and_sources_user_data_operators
1103  */
1104 template <int Tensor_Dim0, int Tensor_Dim1>
1107 
1108  OpCalculateVectorFieldGradientDot(const std::string field_name,
1109  boost::shared_ptr<MatrixDouble> data_ptr,
1110  const EntityType zero_at_type = MBVERTEX)
1113  dataPtr(data_ptr), zeroAtType(zero_at_type) {
1114  if (!dataPtr)
1115  THROW_MESSAGE("Pointer is not set");
1116  }
1117 
1118  MoFEMErrorCode doWork(int side, EntityType type,
1121 
1122  const auto &local_indices = data.getLocalIndices();
1123  const int nb_dofs = local_indices.size();
1124  if (!nb_dofs && type == zeroAtType) {
1125  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
1127  }
1128  if (!nb_dofs)
1130 
1131  dotVector.resize(nb_dofs, false);
1132  const double *array;
1133  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1134  for (int i = 0; i != local_indices.size(); ++i)
1135  if (local_indices[i] != -1)
1136  dotVector[i] = array[local_indices[i]];
1137  else
1138  dotVector[i] = 0;
1139  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1140 
1141  const int nb_gauss_pts = data.getN().size1();
1142  const int nb_base_functions = data.getN().size2();
1143  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat = *dataPtr;
1144  if (type == zeroAtType) {
1145  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1146  mat.clear();
1147  }
1148  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1149  auto gradients_at_gauss_pts =
1150  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1153  int size = nb_dofs / Tensor_Dim0;
1154  if (nb_dofs % Tensor_Dim0) {
1155  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1156  }
1157  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1158  auto field_data = getFTensorDotData<Tensor_Dim0>();
1159  int bb = 0;
1160  for (; bb < size; ++bb) {
1161  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1162  ++field_data;
1163  ++diff_base_function;
1164  }
1165  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1166  // functions
1167  for (; bb != nb_base_functions; ++bb)
1168  ++diff_base_function;
1169  ++gradients_at_gauss_pts;
1170  }
1172  }
1173 
1174 private:
1175  boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1176  EntityType zeroAtType; ///< Zero values at Gauss point at this type
1177  VectorDouble dotVector; ///< Keeps temoorary values of time directives
1178 
1179  template <int Dim> inline auto getFTensorDotData() {
1180  static_assert(Dim || !Dim, "not implemented");
1181  }
1182 };
1183 
1184 template <>
1185 template <>
1188  &dotVector[0], &dotVector[1], &dotVector[2]);
1189 }
1190 
1191 template <>
1192 template <>
1194  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(&dotVector[0],
1195  &dotVector[1]);
1196 }
1197 
1198 /**@}*/
1199 
1200 /** \name Transform tensors and vectors */
1201 
1202 /**@{*/
1203 
1204 /**
1205  * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
1206  * \f$
1207  *
1208  * @tparam DIM
1209  *
1210  * \ingroup mofem_forces_and_sources_user_data_operators
1211  */
1212 template <int DIM_01, int DIM_23, int S = 0>
1215 
1218 
1219  OpTensorTimesSymmetricTensor(const std::string field_name,
1220  boost::shared_ptr<MatrixDouble> in_mat,
1221  boost::shared_ptr<MatrixDouble> out_mat,
1222  boost::shared_ptr<MatrixDouble> d_mat)
1223  : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
1224  // Only is run for vertices
1225  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1226  if (!inMat)
1227  THROW_MESSAGE("Pointer for in mat is null");
1228  if (!outMat)
1229  THROW_MESSAGE("Pointer for out mat is null");
1230  if (!dMat)
1231  THROW_MESSAGE("Pointer for tensor mat is null");
1232  }
1233 
1234  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1236  const size_t nb_gauss_pts = getGaussPts().size2();
1237  auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
1238  auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
1239  outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
1240  auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
1241  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1242  t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
1243  ++t_in;
1244  ++t_out;
1245  }
1247  }
1248 
1249 private:
1254 
1255  boost::shared_ptr<MatrixDouble> inMat;
1256  boost::shared_ptr<MatrixDouble> outMat;
1257  boost::shared_ptr<MatrixDouble> dMat;
1258 };
1259 
1260 template <int DIM>
1262 
1265 
1266  OpSymmetrizeTensor(const std::string field_name,
1267  boost::shared_ptr<MatrixDouble> in_mat,
1268  boost::shared_ptr<MatrixDouble> out_mat)
1269  : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
1270  // Only is run for vertices
1271  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1272  if (!inMat)
1273  THROW_MESSAGE("Pointer not set for in matrix");
1274  if (!outMat)
1275  THROW_MESSAGE("Pointer not set for in matrix");
1276  }
1277 
1278  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1280  const size_t nb_gauss_pts = getGaussPts().size2();
1281  auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
1282  outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
1283  auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
1284  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1285  t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
1286  ++t_in;
1287  ++t_out;
1288  }
1290  }
1291 
1292 private:
1295  boost::shared_ptr<MatrixDouble> inMat;
1296  boost::shared_ptr<MatrixDouble> outMat;
1297 };
1298 
1300 
1303 
1304  OpScaleMatrix(const std::string field_name, const double scale,
1305  boost::shared_ptr<MatrixDouble> in_mat,
1306  boost::shared_ptr<MatrixDouble> out_mat)
1307  : UserOp(field_name, OPROW), scale(scale), inMat(in_mat),
1308  outMat(out_mat) {
1309  // Only is run for vertices
1310  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1311  if (!inMat)
1312  THROW_MESSAGE("Pointer not set for in matrix");
1313  if (!outMat)
1314  THROW_MESSAGE("Pointer not set for in matrix");
1315  }
1316 
1317  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1319  outMat->resize(inMat->size1(), inMat->size2(), false);
1320  noalias(*outMat) = scale * (*inMat);
1322  }
1323 
1324 private:
1325  const double scale;
1326  boost::shared_ptr<MatrixDouble> inMat;
1327  boost::shared_ptr<MatrixDouble> outMat;
1328 };
1329 
1330 /**@}*/
1331 
1332 /** \name H-div/H-curls (Vectorial bases) values at integration points */
1333 
1334 /**@{*/
1335 
1336 /** \brief Get vector field for H-div approximation
1337  * \ingroup mofem_forces_and_sources_user_data_operators
1338  */
1339 template <int Tensor_Dim0, class T, class L, class A>
1342 
1344  const std::string field_name,
1345  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
1346  const EntityType zero_type = MBEDGE, const int zero_side = 0)
1349  dataPtr(data_ptr), zeroType(zero_type), zeroSide(0) {
1350  if (!dataPtr)
1351  THROW_MESSAGE("Pointer is not set");
1352  }
1353 
1354  /**
1355  * \brief calculate values of vector field at integration points
1356  * @param side side entity number
1357  * @param type side entity type
1358  * @param data entity data
1359  * @return error code
1360  */
1361  MoFEMErrorCode doWork(int side, EntityType type,
1363 
1364 private:
1365  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
1366  const EntityHandle zeroType;
1367  const int zeroSide;
1368 };
1369 
1370 template <int Tensor_Dim, class T, class L, class A>
1372  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1374  SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1375  "Not implemented for T = %s and dim = %d",
1376  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
1377  Tensor_Dim);
1379 }
1380 
1381 /** \brief Get vector field for H-div approximation
1382  * \ingroup mofem_forces_and_sources_user_data_operators
1383  */
1384 template <int Tensor_Dim>
1385 struct OpCalculateHdivVectorField_General<Tensor_Dim, double, ublas::row_major,
1388 
1389  OpCalculateHdivVectorField_General(const std::string field_name,
1390  boost::shared_ptr<MatrixDouble> data_ptr,
1391  const EntityType zero_type = MBEDGE,
1392  const int zero_side = 0)
1395  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1396  if (!dataPtr)
1397  THROW_MESSAGE("Pointer is not set");
1398  }
1399 
1400  /**
1401  * \brief Calculate values of vector field at integration points
1402  * @param side side entity number
1403  * @param type side entity type
1404  * @param data entity data
1405  * @return error code
1406  */
1407  MoFEMErrorCode doWork(int side, EntityType type,
1409 
1410 private:
1411  boost::shared_ptr<MatrixDouble> dataPtr;
1412  const EntityHandle zeroType;
1413  const int zeroSide;
1414 };
1415 
1416 template <int Tensor_Dim>
1418  Tensor_Dim, double, ublas::row_major,
1419  DoubleAllocator>::doWork(int side, EntityType type,
1422  const size_t nb_integration_points = this->getGaussPts().size2();
1423  if (type == zeroType && side == zeroSide) {
1424  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1425  dataPtr->clear();
1426  }
1427  const size_t nb_dofs = data.getFieldData().size();
1428  if (!nb_dofs)
1430  const size_t nb_base_functions = data.getN().size2() / Tensor_Dim;
1432  auto t_n_hdiv = data.getFTensor1N<Tensor_Dim>();
1433  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1434  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
1435  auto t_dof = data.getFTensor0FieldData();
1436  int bb = 0;
1437  for (; bb != nb_dofs; ++bb) {
1438  t_data(i) += t_n_hdiv(i) * t_dof;
1439  ++t_n_hdiv;
1440  ++t_dof;
1441  }
1442  for (; bb != nb_base_functions; ++bb)
1443  ++t_n_hdiv;
1444  ++t_data;
1445  }
1447 }
1448 
1449 /** \brief Get vector field for H-div approximation
1450  *
1451  * \ingroup mofem_forces_and_sources_user_data_operators
1452  */
1453 template <int Tensor_Dim>
1456  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1458  Tensor_Dim, double, ublas::row_major,
1460 };
1461 
1462 /** \brief Get vector field for H-div approximation
1463  * \ingroup mofem_forces_and_sources_user_data_operators
1464  */
1465 template <int Tensor_Dim>
1468 
1469  OpCalculateHdivVectorFieldDot(const std::string field_name,
1470  boost::shared_ptr<MatrixDouble> data_ptr,
1471  const EntityType zero_type = MBEDGE,
1472  const int zero_side = 0)
1475  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1476  if (!dataPtr)
1477  THROW_MESSAGE("Pointer is not set");
1478  }
1479 
1480  /**
1481  * \brief Calculate values of vector field at integration points
1482  * @param side side entity number
1483  * @param type side entity type
1484  * @param data entity data
1485  * @return error code
1486  */
1487  MoFEMErrorCode doWork(int side, EntityType type,
1489 
1490 private:
1491  boost::shared_ptr<MatrixDouble> dataPtr;
1492  const EntityHandle zeroType;
1493  const int zeroSide;
1494 };
1495 
1496 template <int Tensor_Dim>
1498  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1500  const size_t nb_integration_points = this->getGaussPts().size2();
1501  if (type == zeroType && side == zeroSide) {
1502  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1503  dataPtr->clear();
1504  }
1505 
1506  auto &local_indices = data.getIndices();
1507  const size_t nb_dofs = local_indices.size();
1508  if (nb_dofs) {
1509 
1510  std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
1511  const double *array;
1512  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1513  for (size_t i = 0; i != nb_dofs; ++i)
1514  if (local_indices[i] != -1)
1515  dot_dofs_vector[i] = array[local_indices[i]];
1516  else
1517  dot_dofs_vector[i] = 0;
1518  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1519 
1520  const size_t nb_base_functions = data.getN().size2() / 3;
1522  auto t_n_hdiv = data.getFTensor1N<Tensor_Dim>();
1523  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1524  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
1525  int bb = 0;
1526  for (; bb != nb_dofs; ++bb) {
1527  t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
1528  ++t_n_hdiv;
1529  }
1530  for (; bb != nb_base_functions; ++bb)
1531  ++t_n_hdiv;
1532  ++t_data;
1533  }
1534  }
1536 }
1537 
1538 /**
1539  * @brief Calculate divergence of vector field
1540  * @ingroup mofem_forces_and_sources_user_data_operators
1541  *
1542  * @tparam BASE_DIM
1543  * @tparam SPACE_DIM
1544  */
1545 template <int BASE_DIM, int SPACE_DIM>
1548 
1549  OpCalculateHdivVectorDivergence(const std::string field_name,
1550  boost::shared_ptr<VectorDouble> data_ptr,
1551  const EntityType zero_type = MBEDGE,
1552  const int zero_side = 0)
1555  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1556  if (!dataPtr)
1557  THROW_MESSAGE("Pointer is not set");
1558  }
1559 
1560  MoFEMErrorCode doWork(int side, EntityType type,
1563  const size_t nb_integration_points = getGaussPts().size2();
1564  if (type == zeroType && side == zeroSide) {
1565  dataPtr->resize(nb_integration_points, false);
1566  dataPtr->clear();
1567  }
1568  const size_t nb_dofs = data.getFieldData().size();
1569  if (!nb_dofs)
1571  const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
1574  auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
1575  auto t_data = getFTensor0FromVec(*dataPtr);
1576  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
1577  auto t_dof = data.getFTensor0FieldData();
1578  int bb = 0;
1579  for (; bb != nb_dofs; ++bb) {
1580  t_data += t_dof * t_n_diff_hdiv(j, j);
1581  ++t_n_diff_hdiv;
1582  ++t_dof;
1583  }
1584  for (; bb != nb_base_functions; ++bb)
1585  ++t_n_diff_hdiv;
1586  ++t_data;
1587  }
1589  }
1590 
1591 private:
1592  boost::shared_ptr<VectorDouble> dataPtr;
1593  const EntityHandle zeroType;
1594  const int zeroSide;
1595 };
1596 
1597 /**
1598  * @brief Calculate gradient of vector field
1599  * @ingroup mofem_forces_and_sources_user_data_operators
1600  *
1601  * @tparam BASE_DIM
1602  * @tparam SPACE_DIM
1603  */
1604 template <int BASE_DIM, int SPACE_DIM>
1607 
1608  OpCalculateHVecVectorGradient(const std::string field_name,
1609  boost::shared_ptr<MatrixDouble> data_ptr,
1610  const EntityType zero_type = MBEDGE,
1611  const int zero_side = 0)
1614  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1615  if (!dataPtr)
1616  THROW_MESSAGE("Pointer is not set");
1617  }
1618 
1619  MoFEMErrorCode doWork(int side, EntityType type,
1622  const size_t nb_integration_points = getGaussPts().size2();
1623  if (type == zeroType && side == zeroSide) {
1624  dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
1625  dataPtr->clear();
1626  }
1627  const size_t nb_dofs = data.getFieldData().size();
1628  if (!nb_dofs)
1630  const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
1633  auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
1634  auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*dataPtr);
1635  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
1636  auto t_dof = data.getFTensor0FieldData();
1637  int bb = 0;
1638  for (; bb != nb_dofs; ++bb) {
1639  t_data(i, j) += t_dof * t_base_diff(i, j);
1640  ++t_base_diff;
1641  ++t_dof;
1642  }
1643  for (; bb != nb_base_functions; ++bb)
1644  ++t_base_diff;
1645  ++t_data;
1646  }
1648  }
1649 
1650 private:
1651  boost::shared_ptr<MatrixDouble> dataPtr;
1652  const EntityHandle zeroType;
1653  const int zeroSide;
1654 };
1655 
1656 /**
1657  * @brief Calculate divergence of vector field dot
1658  * @ingroup mofem_forces_and_sources_user_data_operators
1659  *
1660  * @tparam Tensor_Dim dimension of space
1661  */
1662 template <int Tensor_Dim1, int Tensor_Dim2>
1665 
1666  OpCalculateHdivVectorDivergenceDot(const std::string field_name,
1667  boost::shared_ptr<VectorDouble> data_ptr,
1668  const EntityType zero_type = MBEDGE,
1669  const int zero_side = 0)
1672  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1673  if (!dataPtr)
1674  THROW_MESSAGE("Pointer is not set");
1675  }
1676 
1677  MoFEMErrorCode doWork(int side, EntityType type,
1680  const size_t nb_integration_points = getGaussPts().size2();
1681  if (type == zeroType && side == zeroSide) {
1682  dataPtr->resize(nb_integration_points, false);
1683  dataPtr->clear();
1684  }
1685 
1686  const auto &local_indices = data.getLocalIndices();
1687  const int nb_dofs = local_indices.size();
1688  if (nb_dofs) {
1689 
1690  std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
1691  const double *array;
1692  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1693  for (size_t i = 0; i != local_indices.size(); ++i)
1694  if (local_indices[i] != -1)
1695  dot_dofs_vector[i] = array[local_indices[i]];
1696  else
1697  dot_dofs_vector[i] = 0;
1698  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1699 
1700  const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
1702  auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
1703  auto t_data = getFTensor0FromVec(*dataPtr);
1704  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
1705  int bb = 0;
1706  for (; bb != nb_dofs; ++bb) {
1707  double div = 0;
1708  for (auto ii = 0; ii != Tensor_Dim2; ++ii)
1709  div += t_n_diff_hdiv(ii, ii);
1710  t_data += dot_dofs_vector[bb] * div;
1711  ++t_n_diff_hdiv;
1712  }
1713  for (; bb != nb_base_functions; ++bb)
1714  ++t_n_diff_hdiv;
1715  ++t_data;
1716  }
1717  }
1719  }
1720 
1721 private:
1722  boost::shared_ptr<VectorDouble> dataPtr;
1723  const EntityHandle zeroType;
1724  const int zeroSide;
1725 };
1726 
1727 /**
1728  * @brief Calculate curl of vector field
1729  * @ingroup mofem_forces_and_sources_user_data_operators
1730  *
1731  * @tparam Tensor_Dim dimension of space
1732  */
1733 template <int Tensor_Dim>
1736 
1737  OpCalculateHcurlVectorCurl(const std::string field_name,
1738  boost::shared_ptr<MatrixDouble> data_ptr,
1739  const EntityType zero_type = MBEDGE,
1740  const int zero_side = 0)
1743  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1744  if (!dataPtr)
1745  THROW_MESSAGE("Pointer is not set");
1746  }
1747 
1748  MoFEMErrorCode doWork(int side, EntityType type,
1751  const size_t nb_integration_points = getGaussPts().size2();
1752  if (type == zeroType && side == zeroSide) {
1753  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1754  dataPtr->clear();
1755  }
1756  const size_t nb_dofs = data.getFieldData().size();
1757  if (!nb_dofs)
1759 
1760  MatrixDouble curl_mat;
1761 
1762  const size_t nb_base_functions = data.getN().size2() / Tensor_Dim;
1766  auto t_n_diff_hcurl = data.getFTensor2DiffN<Tensor_Dim, Tensor_Dim>();
1767  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1768  for (int gg = 0; gg != nb_integration_points; ++gg) {
1769 
1770  auto t_dof = data.getFTensor0FieldData();
1771  int bb = 0;
1772  for (; bb != nb_dofs; ++bb) {
1773 
1774  t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
1775  ++t_n_diff_hcurl;
1776  ++t_dof;
1777  }
1778 
1779  for (; bb != nb_base_functions; ++bb)
1780  ++t_n_diff_hcurl;
1781  ++t_data;
1782  }
1783 
1785  }
1786 
1787 private:
1788  boost::shared_ptr<MatrixDouble> dataPtr;
1789  const EntityHandle zeroType;
1790  const int zeroSide;
1791 };
1792 
1793 /**
1794  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1795  * \ingroup mofem_forces_and_sources_user_data_operators
1796  *
1797  * @tparam Tensor_Dim0 rank of the filed
1798  * @tparam Tensor_Dim1 dimension of space
1799  */
1800 template <int Tensor_Dim0, int Tensor_Dim1>
1803 
1804  OpCalculateHVecTensorField(const std::string field_name,
1805  boost::shared_ptr<MatrixDouble> data_ptr,
1806  const EntityType zero_type = MBEDGE,
1807  const int zero_side = 0)
1810  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1811  if (!dataPtr)
1812  THROW_MESSAGE("Pointer is not set");
1813  }
1814 
1815  MoFEMErrorCode doWork(int side, EntityType type,
1818  const size_t nb_integration_points = getGaussPts().size2();
1819  if (type == zeroType && side == zeroSide) {
1820  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1821  dataPtr->clear();
1822  }
1823  const size_t nb_dofs = data.getFieldData().size();
1824  if (nb_dofs) {
1825  const size_t nb_base_functions = data.getN().size2() / 3;
1828  auto t_n_hvec = data.getFTensor1N<3>();
1829  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1830  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
1831  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1832  size_t bb = 0;
1833  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1834  t_data(i, j) += t_dof(i) * t_n_hvec(j);
1835  ++t_n_hvec;
1836  ++t_dof;
1837  }
1838  for (; bb != nb_base_functions; ++bb)
1839  ++t_n_hvec;
1840  ++t_data;
1841  }
1842  }
1844  }
1845 
1846 private:
1847  boost::shared_ptr<MatrixDouble> dataPtr;
1848  const EntityHandle zeroType;
1849  const int zeroSide;
1850 };
1851 
1852 /**
1853  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1854  * \ingroup mofem_forces_and_sources_user_data_operators
1855  *
1856  * @tparam Tensor_Dim0 rank of the filed
1857  * @tparam Tensor_Dim1 dimension of space
1858  */
1859 template <int Tensor_Dim0, int Tensor_Dim1>
1862 
1863  OpCalculateHTensorTensorField(const std::string field_name,
1864  boost::shared_ptr<MatrixDouble> data_ptr,
1865  const EntityType zero_type = MBEDGE,
1866  const int zero_side = 0)
1869  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1870  if (!dataPtr)
1871  THROW_MESSAGE("Pointer is not set");
1872  }
1873 
1874  MoFEMErrorCode doWork(int side, EntityType type,
1877  const size_t nb_integration_points = getGaussPts().size2();
1878  if (type == zeroType && side == zeroSide) {
1879  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1880  dataPtr->clear();
1881  }
1882  const size_t nb_dofs = data.getFieldData().size();
1883  if (!nb_dofs)
1885  const size_t nb_base_functions =
1886  data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
1889  auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
1890  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1891  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
1892  auto t_dof = data.getFTensor0FieldData();
1893  size_t bb = 0;
1894  for (; bb != nb_dofs; ++bb) {
1895  t_data(i, j) += t_dof * t_n_hten(i, j);
1896  ++t_n_hten;
1897  ++t_dof;
1898  }
1899  for (; bb != nb_base_functions; ++bb)
1900  ++t_n_hten;
1901  ++t_data;
1902  }
1904  }
1905 
1906 private:
1907  boost::shared_ptr<MatrixDouble> dataPtr;
1908  const EntityHandle zeroType;
1909  const int zeroSide;
1910 };
1911 
1912 /**
1913  * @brief Calculate divergence of tonsorial field using vectorial base
1914  * \ingroup mofem_forces_and_sources_user_data_operators
1915  *
1916  * @tparam Tensor_Dim0 rank of the field
1917  * @tparam Tensor_Dim1 dimension of space
1918  */
1919 template <int Tensor_Dim0, int Tensor_Dim1>
1922 
1923  OpCalculateHVecTensorDivergence(const std::string field_name,
1924  boost::shared_ptr<MatrixDouble> data_ptr,
1925  const EntityType zero_type = MBEDGE,
1926  const int zero_side = 0)
1929  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1930  if (!dataPtr)
1931  THROW_MESSAGE("Pointer is not set");
1932  }
1933 
1934  MoFEMErrorCode doWork(int side, EntityType type,
1937  const size_t nb_integration_points = getGaussPts().size2();
1938  if (type == zeroType && side == 0) {
1939  dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
1940  dataPtr->clear();
1941  }
1942  const size_t nb_dofs = data.getFieldData().size();
1943  if (nb_dofs) {
1944  const size_t nb_base_functions = data.getN().size2() / 3;
1947  auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
1948  auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
1949  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
1950  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1951  size_t bb = 0;
1952  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1953  double div = t_n_diff_hvec(j, j);
1954  t_data(i) += t_dof(i) * div;
1955  ++t_n_diff_hvec;
1956  ++t_dof;
1957  }
1958  for (; bb < nb_base_functions; ++bb)
1959  ++t_n_diff_hvec;
1960  ++t_data;
1961  }
1962  }
1964  }
1965 
1966 private:
1967  boost::shared_ptr<MatrixDouble> dataPtr;
1968  const EntityHandle zeroType;
1969  const int zeroSide;
1970 };
1971 
1972 /**
1973  * @brief Calculate trace of vector (Hdiv/Hcurl) space
1974  *
1975  * @tparam Tensor_Dim
1976  * @tparam OpBase
1977  */
1978 template <int Tensor_Dim, typename OpBase>
1980 
1981  OpCalculateHVecTensorTrace(const std::string field_name,
1982  boost::shared_ptr<MatrixDouble> data_ptr,
1983  const EntityType zero_type = MBEDGE,
1984  const int zero_side = 0)
1985  : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
1986  zeroType(zero_type), zeroSide(zero_side) {
1987  if (!dataPtr)
1988  THROW_MESSAGE("Pointer is not set");
1989  }
1990 
1991  MoFEMErrorCode doWork(int side, EntityType type,
1994  const size_t nb_integration_points = OpBase::getGaussPts().size2();
1995  if (type == zeroType && side == 0) {
1996  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1997  dataPtr->clear();
1998  }
1999  const size_t nb_dofs = data.getFieldData().size();
2000  if (nb_dofs) {
2001  auto t_normal = OpBase::getFTensor1Normal();
2002  t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
2003  const size_t nb_base_functions = data.getN().size2() / 3;
2006  auto t_base = data.getFTensor1N<3>();
2007  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
2008  for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2009  auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
2010  size_t bb = 0;
2011  for (; bb != nb_dofs / Tensor_Dim; ++bb) {
2012  t_data(i) += t_dof(i) * (t_base(j) * t_normal(j));
2013  ++t_base;
2014  ++t_dof;
2015  }
2016  for (; bb < nb_base_functions; ++bb)
2017  ++t_base;
2018  ++t_data;
2019  }
2020  }
2022  }
2023 
2024 private:
2025  boost::shared_ptr<MatrixDouble> dataPtr;
2026  const EntityHandle zeroType;
2027  const int zeroSide;
2030 };
2031 
2032 /**@}*/
2033 
2034 /**@}*/
2035 
2036 inline auto getFaceJac(MatrixDouble &jac, const FTensor::Number<2> &) {
2038  &jac(0, 0), &jac(1, 0), &jac(2, 0), &jac(3, 0)};
2039 }
2040 
2041 inline auto getFaceJac(MatrixDouble &jac, const FTensor::Number<3> &) {
2043  &jac(0, 0), &jac(1, 0), &jac(4, 0), &jac(2, 0), &jac(3, 0),
2044  &jac(5, 0), &jac(6, 0), &jac(7, 0), &jac(8, 0)};
2045 }
2046 
2047 /** \name Operators for faces */
2048 
2049 /**@{*/
2050 
2051 /** \brief Calculate jacobian for face element
2052 
2053  It is assumed that face element is XY plane. Applied
2054  only for 2d problems.
2055 
2056  \todo Generalize function for arbitrary face orientation in 3d space
2057 
2058  \ingroup mofem_forces_and_sources_tri_element
2059 
2060 */
2061 template <int DIM> struct OpCalculateJacForFaceImpl;
2062 
2063 template <>
2066 
2069  jac(jac) {}
2070 
2071  MoFEMErrorCode doWork(int side, EntityType type,
2073 
2074 protected:
2076 };
2077 
2078 template <> struct OpCalculateJacForFaceImpl<3> : public OpCalculateJacForFaceImpl<2> {
2079 
2081 
2082  MoFEMErrorCode doWork(int side, EntityType type,
2084 };
2085 
2087 
2089 
2090 /** \brief Calculate inverse of jacobian for face element
2091 
2092  It is assumed that face element is XY plane. Applied
2093  only for 2d problems.
2094 
2095  \todo Generalize function for arbitrary face orientation in 3d space
2096 
2097  \ingroup mofem_forces_and_sources_tri_element
2098 
2099 */
2100 template <int DIM> struct OpCalculateInvJacForFaceImpl;
2101 
2102 template <>
2105 
2108  invJac(inv_jac) {}
2109 
2110  MoFEMErrorCode doWork(int side, EntityType type,
2112 
2113 protected:
2115 
2116 };
2117 
2118 template <>
2120  : public OpCalculateInvJacForFaceImpl<2>{
2121 
2123 
2124  MoFEMErrorCode doWork(int side, EntityType type,
2126 };
2127 
2129 
2132 
2133 /** \brief Transform local reference derivatives of shape functions to global
2134 derivatives
2135 
2136 \ingroup mofem_forces_and_sources_tri_element
2137 
2138 */
2139 template <int DIM> struct OpSetInvJacSpaceForFaceImpl;
2140 
2141 template <>
2144 
2147  invJac(inv_jac) {}
2148 
2149  MoFEMErrorCode doWork(int side, EntityType type,
2151 
2152 protected:
2155 };
2156 
2157 template <>
2159  : public OpSetInvJacSpaceForFaceImpl<2> {
2160 
2162 
2163  MoFEMErrorCode doWork(int side, EntityType type,
2165 
2166 };
2167 
2170  : OpSetInvJacSpaceForFaceImpl(inv_jac, H1) {}
2171 };
2172 
2175  : OpSetInvJacSpaceForFaceImpl(inv_jac, L2) {}
2176 };
2177 
2179  : public OpSetInvJacSpaceForFaceImpl<3> {
2181  : OpSetInvJacSpaceForFaceImpl(inv_jac, H1) {}
2182 };
2183 
2185  : public OpSetInvJacSpaceForFaceImpl<3> {
2187  : OpSetInvJacSpaceForFaceImpl(inv_jac, L2) {}
2188 };
2189 
2190 /**
2191  * \brief brief Transform local reference derivatives of shape function to
2192  global derivatives for face
2193 
2194  * \ingroup mofem_forces_and_sources_tri_element
2195  */
2196 template<int DIM>
2198 
2199 template<>
2202 
2205  invJac(inv_jac) {}
2206 
2207  MoFEMErrorCode doWork(int side, EntityType type,
2209 
2210 protected:
2213 };
2214 
2215 template<>
2218  MoFEMErrorCode doWork(int side, EntityType type,
2220 };
2221 
2224 
2225 /**
2226  * @brief Make Hdiv space from Hcurl space in 2d
2227  * @ingroup mofem_forces_and_sources_tri_element
2228  */
2231 
2234 
2235  MoFEMErrorCode doWork(int side, EntityType type,
2237 };
2238 
2243  MoFEMErrorCode doWork(int side, EntityType type,
2245 };
2246 
2247 /** \brief Apply contravariant (Piola) transfer to Hdiv space on face
2248  *
2249  * \note Hdiv space is generated by Hcurl space in 2d.
2250  *
2251  * Contravariant Piola transformation
2252  * \f[
2253  * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
2254  * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
2255  * =
2256  * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
2257  * \f]
2258  *
2259  * \ingroup mofem_forces_and_sources
2260  *
2261  */
2262 template<int DIM>
2264 
2265 template <>
2268 
2271  }
2272 
2273  MoFEMErrorCode doWork(int side, EntityType type,
2275 
2276 protected:
2280 };
2281 
2282 template <>
2287 
2288  MoFEMErrorCode doWork(int side, EntityType type,
2290 };
2291 
2296 
2297 /**@*/
2298 
2299 /** \name Operators for edges */
2300 
2301 /**@{*/
2302 
2305 
2308 
2309  MoFEMErrorCode doWork(int side, EntityType type,
2311 };
2312 
2313 /**
2314  * @deprecated Name is deprecated and this is added for back compatibility
2315  */
2318 
2319 /**@}*/
2320 
2321 /** \name Operator for fat prisms */
2322 
2323 /**@{*/
2324 
2325 /**
2326  * @brief Operator for fat prism element updating integration weights in the
2327  * volume.
2328  *
2329  * Jacobian on the distorted element is nonconstant. This operator updates
2330  * integration weight on prism to take into account nonconstat jacobian.
2331  *
2332  * \f[
2333  * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
2334  * \pmb\xi} \right\| \right)
2335  * \f]
2336  * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
2337  * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
2338  * element coordinate.
2339  *
2340  */
2343 
2346 
2347  MoFEMErrorCode doWork(int side, EntityType type,
2349 };
2350 
2351 /** \brief Calculate inverse of jacobian for face element
2352 
2353  It is assumed that face element is XY plane. Applied
2354  only for 2d problems.
2355 
2356  FIXME Generalize function for arbitrary face orientation in 3d space
2357  FIXME Calculate to Jacobins for two faces
2358 
2359  \ingroup mofem_forces_and_sources_prism_element
2360 
2361 */
2364 
2365  OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2367  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
2368 
2371  invJac(inv_jac) {}
2372  MoFEMErrorCode doWork(int side, EntityType type,
2374 
2375 private:
2376  const boost::shared_ptr<MatrixDouble> invJacPtr;
2378 };
2379 
2380 /** \brief Transform local reference derivatives of shape functions to global
2381 derivatives
2382 
2383 FIXME Generalize to curved shapes
2384 FIXME Generalize to case that top and bottom face has different shape
2385 
2386 \ingroup mofem_forces_and_sources_prism_element
2387 
2388 */
2391 
2392  OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
2394  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
2395 
2398  invJac(inv_jac) {}
2399 
2400  MoFEMErrorCode doWork(int side, EntityType type,
2402 
2403 private:
2404  const boost::shared_ptr<MatrixDouble> invJacPtr;
2407 };
2408 
2409 // Flat prism
2410 
2411 /** \brief Calculate inverse of jacobian for face element
2412 
2413  It is assumed that face element is XY plane. Applied
2414  only for 2d problems.
2415 
2416  FIXME Generalize function for arbitrary face orientation in 3d space
2417  FIXME Calculate to Jacobins for two faces
2418 
2419  \ingroup mofem_forces_and_sources_prism_element
2420 
2421 */
2424 
2427  invJacF3(inv_jac_f3) {}
2428  MoFEMErrorCode doWork(int side, EntityType type,
2430 
2431 private:
2433 };
2434 
2435 /** \brief Transform local reference derivatives of shape functions to global
2436 derivatives
2437 
2438 FIXME Generalize to curved shapes
2439 FIXME Generalize to case that top and bottom face has different shape
2440 
2441 \ingroup mofem_forces_and_sources_prism_element
2442 
2443 */
2446 
2449  invJacF3(inv_jac_f3) {}
2450 
2451  MoFEMErrorCode doWork(int side, EntityType type,
2453 
2454 private:
2457 };
2458 
2459 /**@}*/
2460 
2461 } // namespace MoFEM
2462 
2463 #endif // __USER_DATA_OPERATORS_HPP__
2464 
2465 /**
2466  * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
2467  *
2468  * \brief Classes and functions used to evaluate fields at integration pts,
2469  *jacobians, etc..
2470  *
2471  * \ingroup mofem_forces_and_sources
2472  **/
static Index< 'L', 3 > L
static Index< 'J', 3 > J
static Index< 'I', 3 > I
MatrixDouble invJac
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
FieldSpace
approximation spaces
Definition: definitions.h:174
@ L2
field with C-1 continuity
Definition: definitions.h:180
@ H1
continuous field
Definition: definitions.h:177
@ NOSPACE
Definition: definitions.h:175
@ HCURL
field with continuous tangents
Definition: definitions.h:178
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
const double T
@ row_major
Definition: Layout.hpp:13
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
std::vector< double, std::allocator< double > > DoubleAllocator
Definition: Types.hpp:70
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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
auto getFaceJac(MatrixDouble &jac, const FTensor::Number< 2 > &)
constexpr int SPACE_DIM
Definition: plastic.cpp:50
DataForcesAndSourcesCore::EntData EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from filed data coeffinects.
const VectorDouble & getFieldData() const
get dofs values
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(Tensor_Dim *(Tensor_Dim+1))/2 >, Tensor_Dim > getFTensor2SymmetricFieldData()
Return symmetric FTensor rank 2, i.e. matrix from filed data coeffinects.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getLocalIndices() const
get local indices of dofs on entity
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data operator to do calculations at integration points.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
structure to get information form mofem into DataForcesAndSourcesCore
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.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorDivergence(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.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate trace of vector (Hdiv/Hcurl) space.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Index< 'j', Tensor_Dim > j
FTensor::Index< 'i', Tensor_Dim > i
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecVectorGradient(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.
Calculate curl of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHcurlVectorCurl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate divergence of vector field dot.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergenceDot(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate divergence of vector field.
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
OpCalculateHdivVectorDivergence(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
OpCalculateHdivVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate values of vector field at integration points
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
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)
Get vector field for H-div approximation.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Calculate values of vector field at integration points.
OpCalculateHdivVectorFieldDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< MatrixDouble > dataPtr
Get vector field for H-div approximation.
Calculate inverse of jacobian for face element.
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFatPrism(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateInvJacForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
const boost::shared_ptr< MatrixDouble > invJacPtr
Calculate inverse of jacobian for face element.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateInvJacForFlatPrism(MatrixDouble &inv_jac_f3)
Calculate jacobian for face element.
OpCalculateScalarFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector)
OpCalculateScalarFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Scalar field values at integration points.
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
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A >> data_ptr, const EntityType zero_type=MBVERTEX)
Get rate of scalar field at integration points.
OpCalculateScalarFieldValuesDot(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
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
Get value at integration points for scalar field.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate values of scalar field at integration points
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator >> &data_ptr, const EntityType zero_type=MBVERTEX)
Calculate field values for tenor field rank 2.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A >> data_ptr, const EntityType zero_type=MBVERTEX)
Get time direvarive 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)
EntityType zeroAtType
Zero values at Gauss point at this type.
VectorDouble dotVector
Keeps temoorary values of time directives.
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
OpCalculateTensor2FieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate symmetric tensor field rates ant integratio pts.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateTensor2SymmetricFieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate symmetric tensor field values at integration pts.
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< MatrixDouble > dataPtr
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)
Evaluate field gradient values for vector field, i.e. gradient is tensor rank 2.
OpCalculateVectorFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateVectorFieldGradientDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
EntityType zeroAtType
Zero values at Gauss point at this type.
VectorDouble dotVector
Keeps temoorary values of time directives.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
OpCalculateVectorFieldGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate field values for tenor field rank 1, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate values of vector field at integration points
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A >> data_ptr, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Approximate field valuse for given petsc vector.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateVectorFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Make Hdiv space from Hcurl space in 2d.
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.
Operator for fat prism element updating integration weights in the volume.
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.
OpScaleMatrix(const std::string field_name, const double scale, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< MatrixDouble > outMat
boost::shared_ptr< MatrixDouble > inMat
Apply contravariant (Piola) transfer to Hdiv space on face.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1ForFace(MatrixDouble &inv_jac)
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const boost::shared_ptr< MatrixDouble > invJacPtr
OpSetInvJacH1ForFatPrism(MatrixDouble &inv_jac)
OpSetInvJacH1ForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
OpSetInvJacH1ForFlatPrism(MatrixDouble &inv_jac_f3)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacHcurlFaceImpl(MatrixDouble &inv_jac)
brief Transform local reference derivatives of shape function to global derivatives for face
OpSetInvJacL2ForFace(MatrixDouble &inv_jac)
OpSetInvJacSpaceForFaceImpl(MatrixDouble &inv_jac, FieldSpace space)
Transform local reference derivatives of shape functions to global derivatives.
FTensor::Index< 'i', DIM > i
boost::shared_ptr< MatrixDouble > outMat
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSymmetrizeTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
FTensor::Index< 'j', DIM > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'i', DIM_01 > i
OpTensorTimesSymmetricTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'k', DIM_23 > k
boost::shared_ptr< MatrixDouble > dMat
FTensor::Index< 'j', DIM_01 > j
boost::shared_ptr< MatrixDouble > inMat
FTensor::Index< 'l', DIM_23 > l
boost::shared_ptr< MatrixDouble > outMat
static constexpr Switches CtxSetX
Definition: LoopMethods.hpp:65
static constexpr Switches CtxSetX_TT
Definition: LoopMethods.hpp:67
static constexpr Switches CtxSetX_T
Definition: LoopMethods.hpp:66