v0.9.2
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 // GET VALUES AT GAUSS PTS
26 
27 // TENSOR0
28 
29 /** \brief Calculate field values for tenor field rank 0, i.e. scalar field
30  *
31  */
32 template <class T, class A>
35 
36  boost::shared_ptr<ublas::vector<T, A>> dataPtr;
38 
40  const std::string field_name,
41  boost::shared_ptr<ublas::vector<T, A>> data_ptr,
42  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 
61 /**
62  * \brief Specialization of member function
63  *
64  */
65 template <class T, class A>
67  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
69  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented for T = %s",
70  typeid(T).name() // boost::core::demangle(typeid(T).name()).c_str()
71  );
73 }
74 
75 /**
76  * \brief Get value at integration points for scalar field
77  *
78  * \ingroup mofem_forces_and_sources_user_data_operators
79  */
81  : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
82 
83  OpCalculateScalarFieldValues(const std::string field_name,
84  boost::shared_ptr<VectorDouble> data_ptr,
85  const EntityType zero_type = MBVERTEX)
87  field_name, data_ptr, zero_type) {
88  if (!dataPtr)
89  THROW_MESSAGE("Pointer is not set");
90  }
91 
92  /**
93  * \brief calculate values of scalar field at integration points
94  * @param side side entity number
95  * @param type side entity type
96  * @param data entity data
97  * @return error code
98  */
99  MoFEMErrorCode doWork(int side, EntityType type,
102  VectorDouble &vec = *dataPtr;
103  const int nb_gauss_pts = getGaussPts().size2();
104  if (type == zeroType) {
105  vec.resize(nb_gauss_pts, false);
106  vec.clear();
107  }
108  const int nb_dofs = data.getFieldData().size();
109  if (nb_dofs) {
110  const int nb_gauss_pts = getGaussPts().size2();
111  const int nb_base_functions = data.getN().size2();
112  auto base_function = data.getFTensor0N();
113  auto values_at_gauss_pts = getFTensor0FromVec(vec);
114  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
115  auto field_data = data.getFTensor0FieldData();
116  int bb = 0;
117  for (; bb != nb_dofs; ++bb) {
118  values_at_gauss_pts += field_data * base_function;
119  ++field_data;
120  ++base_function;
121  }
122  // It is possible to have more base functions than dofs
123  for (; bb != nb_base_functions; ++bb)
124  ++base_function;
125  ++values_at_gauss_pts;
126  }
127  }
129  }
130 };
131 
132 /**
133  * @brief Get rate of scalar field at integration points
134  *
135  * \ingroup mofem_forces_and_sources_user_data_operators
136  */
139 
140  boost::shared_ptr<VectorDouble> dataPtr;
143 
144  OpCalculateScalarFieldValuesDot(const std::string field_name,
145  boost::shared_ptr<VectorDouble> data_ptr,
146  const EntityType zero_at_type = MBVERTEX)
149  dataPtr(data_ptr), zeroAtType(zero_at_type) {
150  if (!dataPtr)
151  THROW_MESSAGE("Pointer is not set");
152  }
153 
154  MoFEMErrorCode doWork(int side, EntityType type,
157 
158  const int nb_gauss_pts = getGaussPts().size2();
159  VectorDouble &vec = *dataPtr;
160  if (type == zeroAtType) {
161  vec.resize(nb_gauss_pts, false);
162  vec.clear();
163  }
164 
165  auto &local_indices = data.getLocalIndices();
166  const int nb_dofs = local_indices.size();
167  if (nb_dofs) {
168 
169  dotVector.resize(nb_dofs, false);
170  const double *array;
171  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
172  for (int i = 0; i != local_indices.size(); ++i)
173  if (local_indices[i] != -1)
174  dotVector[i] = array[local_indices[i]];
175  else
176  dotVector[i] = 0;
177  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
178 
179  const int nb_base_functions = data.getN().size2();
180  auto base_function = data.getFTensor0N();
181  auto values_at_gauss_pts = getFTensor0FromVec(vec);
182 
183  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
184  auto field_data = getFTensor0FromVec(dotVector);
185  int bb = 0;
186  for (; bb != nb_dofs; ++bb) {
187  values_at_gauss_pts += field_data * base_function;
188  ++field_data;
189  ++base_function;
190  }
191  // Number of dofs can be smaller than number of Tensor_Dim x base
192  // functions
193  for (; bb != nb_base_functions; ++bb)
194  ++base_function;
195  ++values_at_gauss_pts;
196  }
197  }
199  }
200 };
201 
202 /**
203  * \depreacted Name inconstent with other operators
204  *
205  */
207 
208 // TENSOR1
209 
210 /** \brief Calculate field values for tenor field rank 1, i.e. vector field
211  *
212  */
213 template <int Tensor_Dim, class T, class L, class A>
216 
217  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
219 
221  const std::string field_name,
222  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
223  const EntityType zero_type = MBVERTEX)
226  dataPtr(data_ptr), zeroType(zero_type) {
227  if (!dataPtr)
228  THROW_MESSAGE("Pointer is not set");
229  }
230 
231  /**
232  * \brief calculate values of vector field at integration points
233  * @param side side entity number
234  * @param type side entity type
235  * @param data entity data
236  * @return error code
237  */
238  MoFEMErrorCode doWork(int side, EntityType type,
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>
260  ublas::row_major, DoubleAllocator>
262 
263  boost::shared_ptr<MatrixDouble> dataPtr;
265 
267  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
268  const EntityType zero_type = MBVERTEX)
271  dataPtr(data_ptr), zeroType(zero_type) {
272  if (!dataPtr)
273  THROW_MESSAGE("Pointer is not set");
274  }
275 
276  MoFEMErrorCode doWork(int side, EntityType type,
278 };
279 
280 /**
281  * \brief Member function specialization calculating values for tenor field rank
282  *
283  */
284 template <int Tensor_Dim>
285 MoFEMErrorCode OpCalculateVectorFieldValues_General<
286  Tensor_Dim, double, ublas::row_major,
287  DoubleAllocator>::doWork(int side, EntityType type,
290 
291  const int nb_dofs = data.getFieldData().size();
292 
293  if (nb_dofs) {
294  const int nb_gauss_pts = getGaussPts().size2();
295  auto &mat = *dataPtr;
296  if (type == zeroType) {
297  mat.resize(Tensor_Dim, nb_gauss_pts, false);
298  mat.clear();
299  }
300  if (nb_gauss_pts) {
301  const int nb_base_functions = data.getN().size2();
302  auto base_function = data.getFTensor0N();
303  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
305  const int size = nb_dofs / Tensor_Dim;
306  if (nb_dofs % Tensor_Dim) {
307  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
308  "Data inconsistency");
309  }
310  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
311  auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
312  int bb = 0;
313  for (; bb != size; ++bb) {
314  values_at_gauss_pts(I) += field_data(I) * base_function;
315  ++field_data;
316  ++base_function;
317  }
318  // Number of dofs can be smaller than number of Tensor_Dim x base
319  // functions
320  for (; bb != nb_base_functions; ++bb)
321  ++base_function;
322  ++values_at_gauss_pts;
323  }
324  }
325  } else if (type == this->zeroType) {
326  dataPtr->resize(Tensor_Dim, 0, false);
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 
341  OpCalculateVectorFieldValues(const std::string field_name,
342  boost::shared_ptr<MatrixDouble> data_ptr,
343  const EntityType zero_type = MBVERTEX)
345  ublas::row_major, DoubleAllocator>(
346  field_name, data_ptr, zero_type) {}
347 };
348 
349 /** \brief Get time direvatives of values at integration pts for tensor filed
350  * rank 1, i.e. vector field
351  *
352  * \ingroup mofem_forces_and_sources_user_data_operators
353  */
354 template <int Tensor_Dim>
357 
358  boost::shared_ptr<MatrixDouble> dataPtr;
361 
362  template <int Dim> inline auto getFTensorDotData() {
363  static_assert(Dim || !Dim, "not implemented");
364  }
365 
366  OpCalculateVectorFieldValuesDot(const std::string field_name,
367  boost::shared_ptr<MatrixDouble> data_ptr,
368  const EntityType zero_at_type = MBVERTEX)
371  dataPtr(data_ptr), zeroAtType(zero_at_type) {
372  if (!dataPtr)
373  THROW_MESSAGE("Pointer is not set");
374  }
375 
376  MoFEMErrorCode doWork(int side, EntityType type,
379  auto &local_indices = data.getLocalIndices();
380  const int nb_dofs = local_indices.size();
381  if (!nb_dofs && type == zeroAtType) {
382  dataPtr->resize(Tensor_Dim, 0, false);
384  }
385  if (!nb_dofs)
387 
388  dotVector.resize(nb_dofs, false);
389  const double *array;
390  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
391  for (int i = 0; i != local_indices.size(); ++i)
392  if (local_indices[i] != -1)
393  dotVector[i] = array[local_indices[i]];
394  else
395  dotVector[i] = 0;
396  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
397 
398  const int nb_gauss_pts = data.getN().size1();
399  const int nb_base_functions = data.getN().size2();
400  MatrixDouble &mat = *dataPtr;
401  if (type == zeroAtType) {
402  mat.resize(Tensor_Dim, nb_gauss_pts, false);
403  mat.clear();
404  }
405  auto base_function = data.getFTensor0N();
406  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
408  const int size = nb_dofs / Tensor_Dim;
409  if (nb_dofs % Tensor_Dim) {
410  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
411  }
412  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
413  auto field_data = getFTensorDotData<3>();
414  int bb = 0;
415  for (; bb != size; ++bb) {
416  values_at_gauss_pts(I) += field_data(I) * base_function;
417  ++field_data;
418  ++base_function;
419  }
420  // Number of dofs can be smaller than number of Tensor_Dim x base
421  // functions
422  for (; bb != nb_base_functions; ++bb)
423  ++base_function;
424  ++values_at_gauss_pts;
425  }
427  }
428 };
429 
430 template <>
431 template <>
434  &dotVector[0], &dotVector[1], &dotVector[2]);
435 }
436 
437 /** \brief Calculate field values for tenor field rank 2.
438  *
439  */
440 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
443 
444  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
446 
448  const std::string field_name,
449  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
450  const EntityType zero_type = MBVERTEX)
453  dataPtr(data_ptr), zeroType(zero_type) {
454  if (!dataPtr)
455  THROW_MESSAGE("Pointer is not set");
456  }
457 
458  MoFEMErrorCode doWork(int side, EntityType type,
460 };
461 
462 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
465  doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
467  SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
468  "Not implemented for T = %s, dim0 = %d and dim1 = %d",
469  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
470  Tensor_Dim0, Tensor_Dim1);
472 }
473 
474 template <int Tensor_Dim0, int Tensor_Dim1>
475 struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
476  ublas::row_major, DoubleAllocator>
478 
479  boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
482 
484  const std::string field_name,
485  boost::shared_ptr<
486  ublas::matrix<double, ublas::row_major, DoubleAllocator>> &data_ptr,
487  const EntityType zero_type = MBVERTEX)
490  dataPtr(data_ptr), zeroType(zero_type) {
491  if (!dataPtr)
492  THROW_MESSAGE("Pointer is not set");
493  }
494 
495  MoFEMErrorCode doWork(int side, EntityType type,
497 };
498 
499 template <int Tensor_Dim0, int Tensor_Dim1>
500 MoFEMErrorCode OpCalculateTensor2FieldValues_General<
501  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
502  DoubleAllocator>::doWork(int side, EntityType type,
505  const int nb_dofs = data.getFieldData().size();
506  const int nb_gauss_pts = data.getN().size1();
507  if (!nb_dofs && type == this->zeroType) {
508  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
510  }
511  if (!nb_dofs) {
513  }
514  const int nb_base_functions = data.getN().size2();
515  MatrixDouble &mat = *dataPtr;
516  if (type == zeroType) {
517  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
518  mat.clear();
519  }
520  auto base_function = data.getFTensor0N();
521  auto values_at_gauss_pts = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
524  const int size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
525  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
526  auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
527  int bb = 0;
528  for (; bb != size; ++bb) {
529  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
530  ++field_data;
531  ++base_function;
532  }
533  for (; bb != nb_base_functions; ++bb)
534  ++base_function;
535  ++values_at_gauss_pts;
536  }
538 }
539 
540 /** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
541  * field
542  *
543  * \ingroup mofem_forces_and_sources_user_data_operators
544  */
545 template <int Tensor_Dim0, int Tensor_Dim1>
548  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
549 
550  OpCalculateTensor2FieldValues(const std::string field_name,
551  boost::shared_ptr<MatrixDouble> data_ptr,
552  const EntityType zero_type = MBVERTEX)
553  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
554  ublas::row_major,
556  field_name, data_ptr, zero_type) {}
557 };
558 
559 /** \brief Get time direvarive values at integration pts for tensor filed rank
560  * 2, i.e. matrix field
561  *
562  * \ingroup mofem_forces_and_sources_user_data_operators
563  */
564 template <int Tensor_Dim0, int Tensor_Dim1>
567 
568  boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
569  EntityType zeroAtType; ///< Zero values at Gauss point at this type
570  VectorDouble dotVector; ///< Keeps temoorary values of time directives
571 
572  template <int Dim0, int Dim1> auto getFTensorDotData() {
573  static_assert(Dim0 || !Dim0 || Dim1 || !Dim1, "not implemented");
574  }
575 
576  OpCalculateTensor2FieldValuesDot(const std::string field_name,
577  boost::shared_ptr<MatrixDouble> data_ptr,
578  const EntityType zero_at_type = MBVERTEX)
581  dataPtr(data_ptr), zeroAtType(zero_at_type) {
582  if (!dataPtr)
583  THROW_MESSAGE("Pointer is not set");
584  }
585 
586  MoFEMErrorCode doWork(int side, EntityType type,
589  const auto &local_indices = data.getLocalIndices();
590  const int nb_dofs = local_indices.size();
591  const int nb_gauss_pts = data.getN().size1();
592  if (!nb_dofs && type == zeroAtType) {
593  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
595  }
596  if (!nb_dofs) {
598  }
599 
600  dotVector.resize(nb_dofs, false);
601  const double *array;
602  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
603  for (int i = 0; i != local_indices.size(); ++i)
604  if (local_indices[i] != -1)
605  dotVector[i] = array[local_indices[i]];
606  else
607  dotVector[i] = 0;
608  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
609 
610  const int nb_base_functions = data.getN().size2();
611  MatrixDouble &mat = *dataPtr;
612  if (type == zeroAtType) {
613  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
614  mat.clear();
615  }
616  auto base_function = data.getFTensor0N();
617  auto values_at_gauss_pts =
618  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
621  const int size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
622  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
623  auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
624  int bb = 0;
625  for (; bb != size; ++bb) {
626  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
627  ++field_data;
628  ++base_function;
629  }
630  for (; bb != nb_base_functions; ++bb)
631  ++base_function;
632  ++values_at_gauss_pts;
633  }
635  }
636 };
637 
638 template <>
639 template <>
642  &dotVector[0], &dotVector[1], &dotVector[2],
643 
644  &dotVector[3], &dotVector[4], &dotVector[5],
645 
646  &dotVector[6], &dotVector[7], &dotVector[8]);
647 }
648 
649 /**
650  * @brief Calculate symmetric tensor field values at integration pts.
651  *
652  * @tparam Tensor_Dim
653 
654  * \ingroup mofem_forces_and_sources_user_data_operators
655  */
656 template <int Tensor_Dim>
659 
660  boost::shared_ptr<MatrixDouble> dataPtr;
662  const int zeroSide;
663 
665  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
666  const EntityType zero_type = MBEDGE, const int zero_side = 0)
669  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
670  if (!dataPtr)
671  THROW_MESSAGE("Pointer is not set");
672  }
673 
674  MoFEMErrorCode doWork(int side, EntityType type,
677  MatrixDouble &mat = *dataPtr;
678  const int nb_gauss_pts = getGaussPts().size2();
679  if (type == this->zeroType && side == zeroSide) {
680  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
681  mat.clear();
682  }
683  const int nb_dofs = data.getFieldData().size();
684  if (!nb_dofs)
686 
687  const int nb_base_functions = data.getN().size2();
688  auto base_function = data.getFTensor0N();
689  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
692  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
693  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
694  auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
695  int bb = 0;
696  for (; bb != size; ++bb) {
697  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
698  ++field_data;
699  ++base_function;
700  }
701  for (; bb != nb_base_functions; ++bb)
702  ++base_function;
703  ++values_at_gauss_pts;
704  }
705 
707  }
708 };
709 
710 /**
711  * @brief Calculate symmetric tensor field rates ant integratio pts.
712  *
713  * @tparam Tensor_Dim
714  *
715  * \ingroup mofem_forces_and_sources_user_data_operators
716  */
717 template <int Tensor_Dim>
720 
721  boost::shared_ptr<MatrixDouble> dataPtr;
723  const int zeroSide;
725 
726  template <int Dim> inline auto getFTensorDotData() {
727  static_assert(Dim || !Dim, "not implemented");
728  }
729 
731  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
732  const EntityType zero_type = MBEDGE, const int zero_side = 0)
735  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
736  if (!dataPtr)
737  THROW_MESSAGE("Pointer is not set");
738  }
739 
740  MoFEMErrorCode doWork(int side, EntityType type,
743  const int nb_gauss_pts = getGaussPts().size2();
744  MatrixDouble &mat = *dataPtr;
745  if (type == zeroType && side == zeroSide) {
746  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
747  mat.clear();
748  }
749  auto &local_indices = data.getLocalIndices();
750  const int nb_dofs = local_indices.size();
751  if (!nb_dofs)
753 
754  dotVector.resize(nb_dofs, false);
755  const double *array;
756  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
757  for (int i = 0; i != local_indices.size(); ++i)
758  if (local_indices[i] != -1)
759  dotVector[i] = array[local_indices[i]];
760  else
761  dotVector[i] = 0;
762  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
763 
764  const int nb_base_functions = data.getN().size2();
765 
766  auto base_function = data.getFTensor0N();
767  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
770  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
771  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
772  auto field_data = getFTensorDotData<Tensor_Dim>();
773  int bb = 0;
774  for (; bb != size; ++bb) {
775  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
776  ++field_data;
777  ++base_function;
778  }
779  for (; bb != nb_base_functions; ++bb)
780  ++base_function;
781  ++values_at_gauss_pts;
782  }
783 
785  }
786 };
787 
788 template <>
789 template <>
790 inline auto
793  &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
794  &dotVector[5]);
795 }
796 
797 template <>
798 template <>
799 inline auto
802  &dotVector[0], &dotVector[1], &dotVector[2]);
803 }
804 
805 // GET GRADIENTS AT GAUSS POINTS
806 
807 /**
808  * \brief Evaluate field gradient values for scalar field, i.e. gradient is
809  * tensor rank 1 (vector)
810  *
811  */
812 template <int Tensor_Dim, class T, class L, class A>
814  : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
815 
817  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
818  const EntityType zero_type = MBVERTEX)
819  : OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A>(
820  field_name, data_ptr, zero_type) {}
821 };
822 
823 /** \brief Evaluate field gradient values for scalar field, i.e. gradient is
824  * tensor rank 1 (vector), specialization
825  *
826  */
827 template <int Tensor_Dim>
829  ublas::row_major, DoubleAllocator>
831  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
832 
834  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
835  const EntityType zero_type = MBVERTEX)
837  ublas::row_major, DoubleAllocator>(
838  field_name, data_ptr, zero_type) {}
839 
840  /**
841  * \brief calculate gradient values of scalar field at integration points
842  * @param side side entity number
843  * @param type side entity type
844  * @param data entity data
845  * @return error code
846  */
847  MoFEMErrorCode doWork(int side, EntityType type,
849 };
850 
851 /**
852  * \brief Member function specialization calculating scalar field gradients for
853  * tenor field rank 1
854  *
855  */
856 template <int Tensor_Dim>
857 MoFEMErrorCode OpCalculateScalarFieldGradient_General<
858  Tensor_Dim, double, ublas::row_major,
859  DoubleAllocator>::doWork(int side, EntityType type,
862  if (!this->dataPtr)
863  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
864  "Data pointer not allocated");
865 
866  const int nb_dofs = data.getFieldData().size();
867  if (nb_dofs) {
868  const int nb_gauss_pts = this->getGaussPts().size2();
869  auto &mat = *this->dataPtr;
870  if (type == this->zeroType) {
871  mat.resize(Tensor_Dim, nb_gauss_pts, false);
872  mat.clear();
873  }
874  if (nb_gauss_pts) {
875  const int nb_base_functions = data.getN().size2();
876  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
877  auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
878 
880  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
881  auto field_data = data.getFTensor0FieldData();
882  int bb = 0;
883  for (; bb != nb_dofs; ++bb) {
884  gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
885  ++field_data;
886  ++diff_base_function;
887  }
888  // Number of dofs can be smaller than number of base functions
889  for (; bb < nb_base_functions; ++bb)
890  ++diff_base_function;
891  ++gradients_at_gauss_pts;
892  }
893  }
894 
895  } else if(type == this->zeroType) {
896  this->dataPtr->resize(Tensor_Dim, 0, false);
897  }
898 
900 }
901 
902 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
903  * vector field
904  *
905  * \ingroup mofem_forces_and_sources_user_data_operators
906  */
907 template <int Tensor_Dim>
910  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
911 
912  OpCalculateScalarFieldGradient(const std::string field_name,
913  boost::shared_ptr<MatrixDouble> data_ptr,
914  const EntityType zero_type = MBVERTEX)
916  Tensor_Dim, double, ublas::row_major, DoubleAllocator>(
917  field_name, data_ptr, zero_type) {}
918 };
919 
920 /**
921  * \brief Evaluate field gradient values for vector field, i.e. gradient is
922  * tensor rank 2
923  *
924  */
925 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
927  : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
928  L, A> {
929 
931  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
932  const EntityType zero_type = MBVERTEX)
933  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
934  A>(field_name, data_ptr,
935  zero_type) {}
936 };
937 
938 template <int Tensor_Dim0, int Tensor_Dim1>
939 struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
940  ublas::row_major, DoubleAllocator>
942  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
943 
945  const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
946  const EntityType zero_type = MBVERTEX)
947  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
948  ublas::row_major,
950  field_name, data_ptr, zero_type) {}
951 
952  /**
953  * \brief calculate values of vector field at integration points
954  * @param side side entity number
955  * @param type side entity type
956  * @param data entity data
957  * @return error code
958  */
959  MoFEMErrorCode doWork(int side, EntityType type,
961 };
962 
963 /**
964  * \brief Member function specialization calculating vector field gradients for
965  * tenor field rank 2
966  *
967  */
968 template <int Tensor_Dim0, int Tensor_Dim1>
969 MoFEMErrorCode OpCalculateVectorFieldGradient_General<
970  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
971  DoubleAllocator>::doWork(int side, EntityType type,
974  if (!this->dataPtr)
975  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
976  "Data pointer not allocated");
977 
978  const int nb_dofs = data.getFieldData().size();
979 
980  if (nb_dofs) {
981  const int nb_gauss_pts = this->getGaussPts().size2();
982  auto &mat = *this->dataPtr;
983  if (type == this->zeroType) {
984  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
985  mat.clear();
986  }
987 
988  if (nb_gauss_pts) {
989  const int nb_base_functions = data.getN().size2();
990  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
991  auto gradients_at_gauss_pts =
992  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
995  int size = nb_dofs / Tensor_Dim0;
996  if (nb_dofs % Tensor_Dim0) {
997  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
998  "Data inconsistency");
999  }
1000  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1001  auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
1002  int bb = 0;
1003  for (; bb < size; ++bb) {
1004  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1005  ++field_data;
1006  ++diff_base_function;
1007  }
1008  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1009  // functions
1010  for (; bb != nb_base_functions; ++bb)
1011  ++diff_base_function;
1012  ++gradients_at_gauss_pts;
1013  }
1014  }
1015 
1016  } else if (type == this->zeroType) {
1017  this->dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
1018  }
1020 }
1021 
1022 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1023  * vector field
1024  *
1025  * \ingroup mofem_forces_and_sources_user_data_operators
1026  */
1027 template <int Tensor_Dim0, int Tensor_Dim1>
1030  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1031 
1032  OpCalculateVectorFieldGradient(const std::string field_name,
1033  boost::shared_ptr<MatrixDouble> data_ptr,
1034  const EntityType zero_type = MBVERTEX)
1035  : OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1036  ublas::row_major,
1037  DoubleAllocator>(
1038  field_name, data_ptr, zero_type) {}
1039 };
1040 
1041 /** \brief Get field gradients time derivative at integration pts for scalar
1042  * filed rank 0, i.e. vector field
1043  *
1044  * \ingroup mofem_forces_and_sources_user_data_operators
1045  */
1046 template <int Tensor_Dim0, int Tensor_Dim1>
1049 
1050  boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1051  EntityType zeroAtType; ///< Zero values at Gauss point at this type
1052  VectorDouble dotVector; ///< Keeps temoorary values of time directives
1053 
1054  template <int Dim> inline auto getFTensorDotData() {
1055  static_assert(Dim || !Dim, "not implemented");
1056  }
1057 
1058  OpCalculateVectorFieldGradientDot(const std::string field_name,
1059  boost::shared_ptr<MatrixDouble> data_ptr,
1060  const EntityType zero_at_type = MBVERTEX)
1063  dataPtr(data_ptr), zeroAtType(zero_at_type) {
1064  if (!dataPtr)
1065  THROW_MESSAGE("Pointer is not set");
1066  }
1067 
1068  MoFEMErrorCode doWork(int side, EntityType type,
1071 
1072  const auto &local_indices = data.getLocalIndices();
1073  const int nb_dofs = local_indices.size();
1074  if (!nb_dofs && type == zeroAtType) {
1075  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
1077  }
1078  if (!nb_dofs)
1080 
1081  dotVector.resize(nb_dofs, false);
1082  const double *array;
1083  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1084  for (int i = 0; i != local_indices.size(); ++i)
1085  if (local_indices[i] != -1)
1086  dotVector[i] = array[local_indices[i]];
1087  else
1088  dotVector[i] = 0;
1089  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1090 
1091  const int nb_gauss_pts = data.getN().size1();
1092  const int nb_base_functions = data.getN().size2();
1093  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat = *dataPtr;
1094  if (type == zeroAtType) {
1095  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1096  mat.clear();
1097  }
1098  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1099  auto gradients_at_gauss_pts =
1100  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1103  int size = nb_dofs / Tensor_Dim0;
1104  if (nb_dofs % Tensor_Dim0) {
1105  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1106  }
1107  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1108  auto field_data = getFTensorDotData<Tensor_Dim0>();
1109  int bb = 0;
1110  for (; bb < size; ++bb) {
1111  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1112  ++field_data;
1113  ++diff_base_function;
1114  }
1115  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1116  // functions
1117  for (; bb != nb_base_functions; ++bb)
1118  ++diff_base_function;
1119  ++gradients_at_gauss_pts;
1120  }
1122  }
1123 };
1124 
1125 template <>
1126 template <>
1129  &dotVector[0], &dotVector[1], &dotVector[2]);
1130 }
1131 
1132 template <>
1133 template <>
1135  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(&dotVector[0],
1136  &dotVector[1]);
1137 }
1138 
1139 /** \brief Get vector field for H-div approximation
1140  *
1141  */
1142 template <int Tensor_Dim0, class T, class L, class A>
1145 
1146  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
1148  const int zeroSide;
1149 
1151  const std::string field_name,
1152  boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
1153  const EntityType zero_type = MBEDGE, const int zero_side = 0)
1156  dataPtr(data_ptr), zeroType(zero_type), zeroSide(0) {
1157  if (!dataPtr)
1158  THROW_MESSAGE("Pointer is not set");
1159  }
1160 
1161  /**
1162  * \brief calculate values of vector field at integration points
1163  * @param side side entity number
1164  * @param type side entity type
1165  * @param data entity data
1166  * @return error code
1167  */
1168  MoFEMErrorCode doWork(int side, EntityType type,
1170 };
1171 
1172 template <int Tensor_Dim, class T, class L, class A>
1174  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1176  SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1177  "Not implemented for T = %s and dim = %d",
1178  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
1179  Tensor_Dim);
1181 }
1182 
1183 /** \brief Get vector field for H-div approximation
1184  *
1185  */
1186 template <int Tensor_Dim>
1190 
1191  boost::shared_ptr<MatrixDouble> dataPtr;
1193  const int zeroSide;
1194 
1195  OpCalculateHdivVectorField_General(const std::string field_name,
1196  boost::shared_ptr<MatrixDouble> data_ptr,
1197  const EntityType zero_type = MBEDGE,
1198  const int zero_side = 0)
1201  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1202  if (!dataPtr)
1203  THROW_MESSAGE("Pointer is not set");
1204  }
1205 
1206  /**
1207  * \brief Calculate values of vector field at integration points
1208  * @param side side entity number
1209  * @param type side entity type
1210  * @param data entity data
1211  * @return error code
1212  */
1213  MoFEMErrorCode doWork(int side, EntityType type,
1215 };
1216 
1217 template <int Tensor_Dim>
1218 MoFEMErrorCode OpCalculateHdivVectorField_General<
1219  Tensor_Dim, double, ublas::row_major,
1220  DoubleAllocator>::doWork(int side, EntityType type,
1223  const int nb_integration_points = this->getGaussPts().size2();
1224  if (type == zeroType && side == zeroSide) {
1225  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1226  dataPtr->clear();
1227  }
1228  const int nb_dofs = data.getFieldData().size();
1229  if (!nb_dofs)
1231  const int nb_base_functions = data.getN().size2() / 3;
1233  auto t_n_hdiv = data.getFTensor1N<Tensor_Dim>();
1234  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1235  for (int gg = 0; gg != nb_integration_points; ++gg) {
1236  auto t_dof = data.getFTensor0FieldData();
1237  int bb = 0;
1238  for (; bb != nb_dofs; ++bb) {
1239  t_data(i) += t_n_hdiv(i) * t_dof;
1240  ++t_n_hdiv;
1241  ++t_dof;
1242  }
1243  for (; bb != nb_base_functions; ++bb)
1244  ++t_n_hdiv;
1245  ++t_data;
1246  }
1248 }
1249 
1250 /** \brief Get vector field for H-div approximation
1251  *
1252  * \ingroup mofem_forces_and_sources_user_data_operators
1253  */
1254 template <int Tensor_Dim>
1257  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1258 
1259  OpCalculateHdivVectorField(const std::string field_name,
1260  boost::shared_ptr<MatrixDouble> data_ptr,
1261  const EntityType zero_type = MBEDGE,
1262  const int zero_side = 0)
1263  : OpCalculateHdivVectorField_General<Tensor_Dim, double, ublas::row_major,
1264  DoubleAllocator>(
1265  field_name, data_ptr, zero_type, zero_side) {}
1266 };
1267 
1268 /**
1269  * @brief Calculate divergence of vector field
1270  * @ingroup mofem_forces_and_sources_user_data_operators
1271  *
1272  * @tparam Tensor_Dim dimension of space
1273  */
1274 template <int Tensor_Dim1, int Tensor_Dim2>
1277 
1278  boost::shared_ptr<VectorDouble> dataPtr;
1280  const int zeroSide;
1281 
1282  OpCalculateHdivVectorDivergence(const std::string field_name,
1283  boost::shared_ptr<VectorDouble> data_ptr,
1284  const EntityType zero_type = MBEDGE,
1285  const int zero_side = 0)
1288  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1289  if (!dataPtr)
1290  THROW_MESSAGE("Pointer is not set");
1291  }
1292 
1293  MoFEMErrorCode doWork(int side, EntityType type,
1296  const int nb_integration_points = getGaussPts().size2();
1297  if (type == zeroType && side == zeroSide) {
1298  dataPtr->resize(nb_integration_points, false);
1299  dataPtr->clear();
1300  }
1301  const int nb_dofs = data.getFieldData().size();
1302  if (!nb_dofs)
1304  const int nb_base_functions = data.getN().size2() / Tensor_Dim1;
1306  auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
1307  auto t_data = getFTensor0FromVec(*dataPtr);
1308  for (int gg = 0; gg != nb_integration_points; ++gg) {
1309  auto t_dof = data.getFTensor0FieldData();
1310  int bb = 0;
1311  for (; bb != nb_dofs; ++bb) {
1312  double div = 0;
1313  for (auto ii = 0; ii != Tensor_Dim2; ++ii)
1314  div += t_n_diff_hdiv(ii, ii);
1315  t_data += t_dof * div;
1316  ++t_n_diff_hdiv;
1317  ++t_dof;
1318  }
1319  for (; bb != nb_base_functions; ++bb)
1320  ++t_n_diff_hdiv;
1321  ++t_data;
1322  }
1324  }
1325 };
1326 
1327 /**
1328  * @brief Calculate curl of vector field
1329  * @ingroup mofem_forces_and_sources_user_data_operators
1330  *
1331  * @tparam Tensor_Dim dimension of space
1332  */
1333 template <int Tensor_Dim>
1336 
1337  boost::shared_ptr<MatrixDouble> dataPtr;
1339  const int zeroSide;
1340 
1341  OpCalculateHcurlVectorCurl(const std::string field_name,
1342  boost::shared_ptr<MatrixDouble> data_ptr,
1343  const EntityType zero_type = MBEDGE,
1344  const int zero_side = 0)
1347  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1348  if (!dataPtr)
1349  THROW_MESSAGE("Pointer is not set");
1350  }
1351 
1352  MoFEMErrorCode doWork(int side, EntityType type,
1355  const int nb_integration_points = getGaussPts().size2();
1356  if (type == zeroType && side == zeroSide) {
1357  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1358  dataPtr->clear();
1359  }
1360  const int nb_dofs = data.getFieldData().size();
1361  if (!nb_dofs)
1363 
1364  MatrixDouble curl_mat;
1365 
1366  const int nb_base_functions = data.getN().size2() / Tensor_Dim;
1370  auto t_n_diff_hcurl = data.getFTensor2DiffN<Tensor_Dim, Tensor_Dim>();
1371  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1372  for (int gg = 0; gg != nb_integration_points; ++gg) {
1373 
1374  auto t_dof = data.getFTensor0FieldData();
1375  int bb = 0;
1376  for (; bb != nb_dofs; ++bb) {
1377 
1378  t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
1379  ++t_n_diff_hcurl;
1380  ++t_dof;
1381  }
1382 
1383  for (; bb != nb_base_functions; ++bb)
1384  ++t_n_diff_hcurl;
1385  ++t_data;
1386  }
1387 
1389  }
1390 };
1391 
1392 /**
1393  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1394  * \ingroup mofem_forces_and_sources_user_data_operators
1395  *
1396  * @tparam Tensor_Dim0 rank of the filed
1397  * @tparam Tensor_Dim1 dimension of space
1398  */
1399 template <int Tensor_Dim0, int Tensor_Dim1>
1402 
1403  boost::shared_ptr<MatrixDouble> dataPtr;
1405  const int zeroSide;
1406 
1407  OpCalculateHVecTensorField(const std::string field_name,
1408  boost::shared_ptr<MatrixDouble> data_ptr,
1409  const EntityType zero_type = MBEDGE,
1410  const int zero_side = 0)
1413  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1414  if (!dataPtr)
1415  THROW_MESSAGE("Pointer is not set");
1416  }
1417 
1418  MoFEMErrorCode doWork(int side, EntityType type,
1421  const int nb_integration_points = getGaussPts().size2();
1422  if (type == zeroType && side == zeroSide) {
1423  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1424  dataPtr->clear();
1425  }
1426  const int nb_dofs = data.getFieldData().size();
1427  if (nb_dofs) {
1428  const int nb_base_functions = data.getN().size2() / 3;
1431  auto t_n_hvec = data.getFTensor1N<3>();
1432  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1433  for (int gg = 0; gg != nb_integration_points; ++gg) {
1434  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1435  int bb = 0;
1436  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1437  t_data(i, j) += t_dof(i) * t_n_hvec(j);
1438  ++t_n_hvec;
1439  ++t_dof;
1440  }
1441  for (; bb != nb_base_functions; ++bb)
1442  ++t_n_hvec;
1443  ++t_data;
1444  }
1445  }
1447  }
1448 };
1449 
1450 /**
1451  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1452  * \ingroup mofem_forces_and_sources_user_data_operators
1453  *
1454  * @tparam Tensor_Dim0 rank of the filed
1455  * @tparam Tensor_Dim1 dimension of space
1456  */
1457 template <int Tensor_Dim0, int Tensor_Dim1>
1460 
1461  boost::shared_ptr<MatrixDouble> dataPtr;
1463  const int zeroSide;
1464 
1465  OpCalculateHTensorTensorField(const std::string field_name,
1466  boost::shared_ptr<MatrixDouble> data_ptr,
1467  const EntityType zero_type = MBEDGE,
1468  const int zero_side = 0)
1471  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1472  if (!dataPtr)
1473  THROW_MESSAGE("Pointer is not set");
1474  }
1475 
1476  MoFEMErrorCode doWork(int side, EntityType type,
1479  const int nb_integration_points = getGaussPts().size2();
1480  if (type == zeroType && side == zeroSide) {
1481  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1482  dataPtr->clear();
1483  }
1484  const int nb_dofs = data.getFieldData().size();
1485  if (!nb_dofs)
1487  const int nb_base_functions =
1488  data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
1491  auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
1492  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1493  for (int gg = 0; gg != nb_integration_points; ++gg) {
1494  auto t_dof = data.getFTensor0FieldData();
1495  int bb = 0;
1496  for (; bb != nb_dofs; ++bb) {
1497  t_data(i, j) += t_dof * t_n_hten(i, j);
1498  ++t_n_hten;
1499  ++t_dof;
1500  }
1501  for (; bb != nb_base_functions; ++bb)
1502  ++t_n_hten;
1503  ++t_data;
1504  }
1506  }
1507 };
1508 
1509 /**
1510  * @brief Calculate divergence of tonsorial field using vectorial base
1511  * \ingroup mofem_forces_and_sources_user_data_operators
1512  *
1513  * @tparam Tensor_Dim0 rank of the field
1514  * @tparam Tensor_Dim1 dimension of space
1515  */
1516 template <int Tensor_Dim0, int Tensor_Dim1>
1519 
1520  boost::shared_ptr<MatrixDouble> dataPtr;
1522  const int zeroSide;
1523 
1524  OpCalculateHVecTensorDivergence(const std::string field_name,
1525  boost::shared_ptr<MatrixDouble> data_ptr,
1526  const EntityType zero_type = MBEDGE,
1527  const int zero_side = 0)
1530  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1531  if (!dataPtr)
1532  THROW_MESSAGE("Pointer is not set");
1533  }
1534 
1535  MoFEMErrorCode doWork(int side, EntityType type,
1538  const int nb_integration_points = getGaussPts().size2();
1539  if (type == zeroType && side == 0) {
1540  dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
1541  dataPtr->clear();
1542  }
1543  const int nb_dofs = data.getFieldData().size();
1544  if (nb_dofs) {
1545  const int nb_base_functions = data.getN().size2() / 3;
1548  auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
1549  auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
1550  for (int gg = 0; gg != nb_integration_points; ++gg) {
1551  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1552  int bb = 0;
1553  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1554  double div = t_n_diff_hvec(j, j);
1555  t_data(i) += t_dof(i) * div;
1556  ++t_n_diff_hvec;
1557  ++t_dof;
1558  }
1559  for (; bb < nb_base_functions; ++bb)
1560  ++t_n_diff_hvec;
1561  ++t_data;
1562  }
1563  }
1565  }
1566 };
1567 
1568 // Face opeartors
1569 
1570 /** \brief Calculate jacobian for face element
1571 
1572  It is assumed that face element is XY plane. Applied
1573  only for 2d problems.
1574 
1575  \todo Generalize function for arbitrary face orientation in 3d space
1576 
1577  \ingroup mofem_forces_and_sources_tri_element
1578 
1579 */
1583 
1586 
1587  MoFEMErrorCode doWork(int side, EntityType type,
1589 };
1590 
1591 /** \brief Calculate inverse of jacobian for face element
1592 
1593  It is assumed that face element is XY plane. Applied
1594  only for 2d problems.
1595 
1596  \todo Generalize function for arbitrary face orientation in 3d space
1597 
1598  \ingroup mofem_forces_and_sources_tri_element
1599 
1600 */
1604 
1607  invJac(inv_jac) {}
1608 
1609  MoFEMErrorCode doWork(int side, EntityType type,
1611 };
1612 
1613 /** \brief Transform local reference derivatives of shape functions to global
1614 derivatives
1615 
1616 \ingroup mofem_forces_and_sources_tri_element
1617 
1618 */
1621 
1624  invJac(inv_jac) {}
1625 
1626  MoFEMErrorCode doWork(int side, EntityType type,
1628 
1629 private:
1632 };
1633 
1634 /**
1635  * \brief brief Transform local reference derivatives of shape function to
1636  global derivatives for face
1637 
1638  * \ingroup mofem_forces_and_sources_tri_element
1639  */
1642 
1645  invJac(inv_jac) {}
1646 
1647  MoFEMErrorCode doWork(int side, EntityType type,
1649 
1650 private:
1653 };
1654 
1655 /**
1656  * @brief Make Hdiv space from Hcurl space in 2d
1657  * @ingroup mofem_forces_and_sources_tri_element
1658  */
1661 
1664 
1665  MoFEMErrorCode doWork(int side, EntityType type,
1667 };
1668 
1669 /** \brief Apply contravariant (Piola) transfer to Hdiv space on face
1670  *
1671  * \note Hdiv space is generated by Hcurl space in 2d.
1672  *
1673  * Contravariant Piola transformation
1674  * \f[
1675  * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
1676  * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
1677  * =
1678  * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
1679  * \f]
1680  *
1681  * \ingroup mofem_forces_and_sources
1682  *
1683  */
1686 
1689  }
1690 
1693 
1694  MoFEMErrorCode doWork(int side, EntityType type,
1696 
1697  private:
1699 };
1700 
1701 // Edge
1702 
1705 
1708 
1709  MoFEMErrorCode doWork(int side, EntityType type,
1711 };
1712 
1713 // Fat prims
1714 
1715 /**
1716  * @brief Operator for fat prism element updating integration weights in the
1717  * volume.
1718  *
1719  * Jacobian on the distorted element is nonconstant. This operator updates
1720  * integration weight on prism to take into account nonconstat jacobian.
1721  *
1722  * \f[
1723  * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
1724  * \pmb\xi} \right\| \right)
1725  * \f]
1726  * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
1727  * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
1728  * element coordinate.
1729  *
1730  */
1733 
1736 
1737  MoFEMErrorCode doWork(int side, EntityType type,
1739 };
1740 
1741 /** \brief Calculate inverse of jacobian for face element
1742 
1743  It is assumed that face element is XY plane. Applied
1744  only for 2d problems.
1745 
1746  FIXME Generalize function for arbitrary face orientation in 3d space
1747  FIXME Calculate to Jacobins for two faces
1748 
1749  \ingroup mofem_forces_and_sources_prism_element
1750 
1751 */
1754 
1755  OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
1757  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
1758 
1761  invJac(inv_jac) {}
1762  MoFEMErrorCode doWork(int side, EntityType type,
1764 
1765  private:
1766  const boost::shared_ptr<MatrixDouble> invJacPtr;
1768 };
1769 
1770 /** \brief Transform local reference derivatives of shape functions to global
1771 derivatives
1772 
1773 FIXME Generalize to curved shapes
1774 FIXME Generalize to case that top and bottom face has different shape
1775 
1776 \ingroup mofem_forces_and_sources_prism_element
1777 
1778 */
1781 
1782  OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
1784  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
1785 
1788  invJac(inv_jac) {}
1789 
1790  MoFEMErrorCode doWork(int side, EntityType type,
1792 
1793  private:
1794  const boost::shared_ptr<MatrixDouble> invJacPtr;
1797 
1798 };
1799 
1800 // Flat prism
1801 
1802 /** \brief Calculate inverse of jacobian for face element
1803 
1804  It is assumed that face element is XY plane. Applied
1805  only for 2d problems.
1806 
1807  FIXME Generalize function for arbitrary face orientation in 3d space
1808  FIXME Calculate to Jacobins for two faces
1809 
1810  \ingroup mofem_forces_and_sources_prism_element
1811 
1812 */
1815 
1819  invJacF3(inv_jac_f3) {}
1820  MoFEMErrorCode doWork(int side, EntityType type,
1822 };
1823 
1824 /** \brief Transform local reference derivatives of shape functions to global
1825 derivatives
1826 
1827 FIXME Generalize to curved shapes
1828 FIXME Generalize to case that top and bottom face has different shape
1829 
1830 \ingroup mofem_forces_and_sources_prism_element
1831 
1832 */
1838  invJacF3(inv_jac_f3) {}
1839 
1841  MoFEMErrorCode doWork(int side, EntityType type,
1843 };
1844 
1845 } // namespace MoFEM
1846 
1847 #endif // __USER_DATA_OPERATORS_HPP__
1848 
1849 /**
1850  * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
1851  *
1852  * \brief Classes and functions used to evaluate fields at integration pts,
1853  *jacobians, etc..
1854  *
1855  * \ingroup mofem_forces_and_sources
1856  **/
const double T
Calculate divergence of vector field.
Get field gradients time derivative at integration pts for scalar filed rank 0, 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
Get time direvatives of values at integration pts for tensor filed rank 1, i.e. vector field.
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
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)
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.
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)
OpCalculateVectorFieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_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:507
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:75
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
EntityType zeroAtType
Zero values at Gauss point at this type.
OpCalculateScalarFieldGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
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.
OpCalculateScalarFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:70
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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:626
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
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.
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.
FTensor::Index< 'J', 3 > J
Definition: PlasticOps.hpp:71
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.
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
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...
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
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.
OpCalculateVectorFieldValues(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)
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)
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)
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.
VectorDouble dotVector
Keeps temoorary values of time directives.
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)
Get vector field for H-div approximation.
Calculate divergence of tonsorial field using vectorial base.
Calculate field values for tenor field rank 2.
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:602
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.
Transform local reference derivatives of shape functions to global derivatives.
Transform local reference derivatives of shape functions to global derivatives.
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
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)
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)
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
Calculate field values for tenor field rank 0, i.e. scalar field.
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.
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
structure to get information form mofem into DataForcesAndSourcesCore
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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:413
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.
Calculate jacobian for face element.
continuous field
Definition: definitions.h:177
boost::shared_ptr< MatrixDouble > dataPtr
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
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
OpCalculateScalarFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBVERTEX)
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:69
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.
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.
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
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.