v0.9.1
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  const int nb_dofs = data.getFieldData().size();
291  if (!nb_dofs && type == this->zeroType) {
292  dataPtr->resize(Tensor_Dim, 0, false);
294  }
295  if (!nb_dofs) {
297  }
298  const int nb_gauss_pts = data.getN().size1();
299  const int nb_base_functions = data.getN().size2();
300  MatrixDouble &mat = *dataPtr;
301  if (type == zeroType) {
302  mat.resize(Tensor_Dim, nb_gauss_pts, false);
303  mat.clear();
304  }
305  auto base_function = data.getFTensor0N();
306  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
308  const int size = nb_dofs / Tensor_Dim;
309  if (nb_dofs % Tensor_Dim) {
310  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
311  }
312  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
313  auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
314  int 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 functions
321  for (; bb != nb_base_functions; ++bb)
322  ++base_function;
323  ++values_at_gauss_pts;
324  }
326 }
327 
328 /** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
329  * field
330  *
331  * \ingroup mofem_forces_and_sources_user_data_operators
332  */
333 template <int Tensor_Dim>
336  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
337 
338  OpCalculateVectorFieldValues(const std::string field_name,
339  boost::shared_ptr<MatrixDouble> &data_ptr,
340  const EntityType zero_type = MBVERTEX)
342  ublas::row_major, DoubleAllocator>(
343  field_name, data_ptr, zero_type) {}
344 };
345 
346 /** \brief Get time direvatives of values at integration pts for tensor filed
347  * rank 1, i.e. vector field
348  *
349  * \ingroup mofem_forces_and_sources_user_data_operators
350  */
351 template <int Tensor_Dim>
354 
355  boost::shared_ptr<MatrixDouble> dataPtr;
358 
359  template <int Dim> inline auto getFTensorDotData() {
360  static_assert(Dim || !Dim, "not implemented");
361  }
362 
363  OpCalculateVectorFieldValuesDot(const std::string field_name,
364  boost::shared_ptr<MatrixDouble> &data_ptr,
365  const EntityType zero_at_type = MBVERTEX)
368  dataPtr(data_ptr), zeroAtType(zero_at_type) {
369  if (!dataPtr)
370  THROW_MESSAGE("Pointer is not set");
371  }
372 
373  MoFEMErrorCode doWork(int side, EntityType type,
376  auto &local_indices = data.getLocalIndices();
377  const int nb_dofs = local_indices.size();
378  if (!nb_dofs && type == zeroAtType) {
379  dataPtr->resize(Tensor_Dim, 0, false);
381  }
382  if (!nb_dofs)
384 
385  dotVector.resize(nb_dofs, false);
386  const double *array;
387  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
388  for (int i = 0; i != local_indices.size(); ++i)
389  if (local_indices[i] != -1)
390  dotVector[i] = array[local_indices[i]];
391  else
392  dotVector[i] = 0;
393  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
394 
395  const int nb_gauss_pts = data.getN().size1();
396  const int nb_base_functions = data.getN().size2();
397  MatrixDouble &mat = *dataPtr;
398  if (type == zeroAtType) {
399  mat.resize(Tensor_Dim, nb_gauss_pts, false);
400  mat.clear();
401  }
402  auto base_function = data.getFTensor0N();
403  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
405  const int size = nb_dofs / Tensor_Dim;
406  if (nb_dofs % Tensor_Dim) {
407  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
408  }
409  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
410  auto field_data = getFTensorDotData<3>();
411  int bb = 0;
412  for (; bb != size; ++bb) {
413  values_at_gauss_pts(I) += field_data(I) * base_function;
414  ++field_data;
415  ++base_function;
416  }
417  // Number of dofs can be smaller than number of Tensor_Dim x base
418  // functions
419  for (; bb != nb_base_functions; ++bb)
420  ++base_function;
421  ++values_at_gauss_pts;
422  }
424  }
425 };
426 
427 template <>
428 template <>
431  &dotVector[0], &dotVector[1], &dotVector[2]);
432 }
433 
434 /** \brief Calculate field values for tenor field rank 2.
435  *
436  */
437 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
440 
441  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
443 
445  const std::string field_name,
446  boost::shared_ptr<ublas::matrix<T, L, A>> &data_ptr,
447  const EntityType zero_type = MBVERTEX)
450  dataPtr(data_ptr), zeroType(zero_type) {
451  if (!dataPtr)
452  THROW_MESSAGE("Pointer is not set");
453  }
454 
455  MoFEMErrorCode doWork(int side, EntityType type,
457 };
458 
459 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
462  doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
464  SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
465  "Not implemented for T = %s, dim0 = %d and dim1 = %d",
466  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
467  Tensor_Dim0, Tensor_Dim1);
469 }
470 
471 template <int Tensor_Dim0, int Tensor_Dim1>
472 struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
473  ublas::row_major, DoubleAllocator>
475 
476  boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
479 
481  const std::string field_name,
482  boost::shared_ptr<
483  ublas::matrix<double, ublas::row_major, DoubleAllocator>> &data_ptr,
484  const EntityType zero_type = MBVERTEX)
487  dataPtr(data_ptr), zeroType(zero_type) {
488  if (!dataPtr)
489  THROW_MESSAGE("Pointer is not set");
490  }
491 
492  MoFEMErrorCode doWork(int side, EntityType type,
494 };
495 
496 template <int Tensor_Dim0, int Tensor_Dim1>
497 MoFEMErrorCode OpCalculateTensor2FieldValues_General<
498  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
499  DoubleAllocator>::doWork(int side, EntityType type,
502  const int nb_dofs = data.getFieldData().size();
503  const int nb_gauss_pts = data.getN().size1();
504  if (!nb_dofs && type == this->zeroType) {
505  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
507  }
508  if (!nb_dofs) {
510  }
511  const int nb_base_functions = data.getN().size2();
512  MatrixDouble &mat = *dataPtr;
513  if (type == zeroType) {
514  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
515  mat.clear();
516  }
517  auto base_function = data.getFTensor0N();
518  auto values_at_gauss_pts = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
521  const int size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
522  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
523  auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
524  int bb = 0;
525  for (; bb != size; ++bb) {
526  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
527  ++field_data;
528  ++base_function;
529  }
530  for (; bb != nb_base_functions; ++bb)
531  ++base_function;
532  ++values_at_gauss_pts;
533  }
535 }
536 
537 /** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
538  * field
539  *
540  * \ingroup mofem_forces_and_sources_user_data_operators
541  */
542 template <int Tensor_Dim0, int Tensor_Dim1>
545  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
546 
547  OpCalculateTensor2FieldValues(const std::string field_name,
548  boost::shared_ptr<MatrixDouble> &data_ptr,
549  const EntityType zero_type = MBVERTEX)
550  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
551  ublas::row_major,
553  field_name, data_ptr, zero_type) {}
554 };
555 
556 /** \brief Get time direvarive values at integration pts for tensor filed rank
557  * 2, i.e. matrix field
558  *
559  * \ingroup mofem_forces_and_sources_user_data_operators
560  */
561 template <int Tensor_Dim0, int Tensor_Dim1>
564 
565  boost::shared_ptr<MatrixDouble> &dataPtr; ///< Data computed into this matrix
566  EntityType zeroAtType; ///< Zero values at Gauss point at this type
567  VectorDouble dotVector; ///< Keeps temoorary values of time directives
568 
569  template <int Dim0, int Dim1> auto getFTensorDotData() {
570  static_assert(Dim0 || !Dim0 || Dim1 || !Dim1, "not implemented");
571  }
572 
573  OpCalculateTensor2FieldValuesDot(const std::string field_name,
574  boost::shared_ptr<MatrixDouble> &data_ptr,
575  const EntityType zero_at_type = MBVERTEX)
578  dataPtr(data_ptr), zeroAtType(zero_at_type) {
579  if (!dataPtr)
580  THROW_MESSAGE("Pointer is not set");
581  }
582 
583  MoFEMErrorCode doWork(int side, EntityType type,
586  const auto &local_indices = data.getLocalIndices();
587  const int nb_dofs = local_indices.size();
588  const int nb_gauss_pts = data.getN().size1();
589  if (!nb_dofs && type == zeroAtType) {
590  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
592  }
593  if (!nb_dofs) {
595  }
596 
597  dotVector.resize(nb_dofs, false);
598  const double *array;
599  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
600  for (int i = 0; i != local_indices.size(); ++i)
601  if (local_indices[i] != -1)
602  dotVector[i] = array[local_indices[i]];
603  else
604  dotVector[i] = 0;
605  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
606 
607  const int nb_base_functions = data.getN().size2();
608  MatrixDouble &mat = *dataPtr;
609  if (type == zeroAtType) {
610  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
611  mat.clear();
612  }
613  auto base_function = data.getFTensor0N();
614  auto values_at_gauss_pts =
615  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
618  const int size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
619  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
620  auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
621  int bb = 0;
622  for (; bb != size; ++bb) {
623  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
624  ++field_data;
625  ++base_function;
626  }
627  for (; bb != nb_base_functions; ++bb)
628  ++base_function;
629  ++values_at_gauss_pts;
630  }
632  }
633 };
634 
635 template <>
636 template <>
639  &dotVector[0], &dotVector[1], &dotVector[2],
640 
641  &dotVector[3], &dotVector[4], &dotVector[5],
642 
643  &dotVector[6], &dotVector[7], &dotVector[8]);
644 }
645 
646 /**
647  * @brief Calculate symmetric tensor field values at integration pts.
648  *
649  * @tparam Tensor_Dim
650 
651  * \ingroup mofem_forces_and_sources_user_data_operators
652  */
653 template <int Tensor_Dim>
656 
657  boost::shared_ptr<MatrixDouble> dataPtr;
659  const int zeroSide;
660 
662  const std::string field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
663  const EntityType zero_type = MBEDGE, const int zero_side = 0)
666  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
667  if (!dataPtr)
668  THROW_MESSAGE("Pointer is not set");
669  }
670 
671  MoFEMErrorCode doWork(int side, EntityType type,
674  MatrixDouble &mat = *dataPtr;
675  const int nb_gauss_pts = getGaussPts().size2();
676  if (type == this->zeroType && side == zeroSide) {
677  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
678  mat.clear();
679  }
680  const int nb_dofs = data.getFieldData().size();
681  if (!nb_dofs)
683 
684  const int nb_base_functions = data.getN().size2();
685  auto base_function = data.getFTensor0N();
686  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
689  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
690  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
691  auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
692  int bb = 0;
693  for (; bb != size; ++bb) {
694  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
695  ++field_data;
696  ++base_function;
697  }
698  for (; bb != nb_base_functions; ++bb)
699  ++base_function;
700  ++values_at_gauss_pts;
701  }
702 
704  }
705 };
706 
707 /**
708  * @brief Calculate symmetric tensor field rates ant integratio pts.
709  *
710  * @tparam Tensor_Dim
711  *
712  * \ingroup mofem_forces_and_sources_user_data_operators
713  */
714 template <int Tensor_Dim>
717 
718  boost::shared_ptr<MatrixDouble> dataPtr;
720  const int zeroSide;
722 
723  template <int Dim> inline auto getFTensorDotData() {
724  static_assert(Dim || !Dim, "not implemented");
725  }
726 
728  const std::string field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
729  const EntityType zero_type = MBEDGE, const int zero_side = 0)
732  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
733  if (!dataPtr)
734  THROW_MESSAGE("Pointer is not set");
735  }
736 
737  MoFEMErrorCode doWork(int side, EntityType type,
740  const int nb_gauss_pts = getGaussPts().size2();
741  MatrixDouble &mat = *dataPtr;
742  if (type == zeroType && side == zeroSide) {
743  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
744  mat.clear();
745  }
746  auto &local_indices = data.getLocalIndices();
747  const int nb_dofs = local_indices.size();
748  if (!nb_dofs)
750 
751  dotVector.resize(nb_dofs, false);
752  const double *array;
753  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
754  for (int i = 0; i != local_indices.size(); ++i)
755  if (local_indices[i] != -1)
756  dotVector[i] = array[local_indices[i]];
757  else
758  dotVector[i] = 0;
759  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
760 
761  const int nb_base_functions = data.getN().size2();
762 
763  auto base_function = data.getFTensor0N();
764  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
767  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
768  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
769  auto field_data = getFTensorDotData<Tensor_Dim>();
770  int bb = 0;
771  for (; bb != size; ++bb) {
772  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
773  ++field_data;
774  ++base_function;
775  }
776  for (; bb != nb_base_functions; ++bb)
777  ++base_function;
778  ++values_at_gauss_pts;
779  }
780 
782  }
783 };
784 
785 template <>
786 template <>
787 inline auto
790  &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
791  &dotVector[5]);
792 }
793 
794 template <>
795 template <>
796 inline auto
799  &dotVector[0], &dotVector[1], &dotVector[2]);
800 }
801 
802 // GET GRADIENTS AT GAUSS POINTS
803 
804 /**
805  * \brief Evaluate field gradient values for scalar field, i.e. gradient is
806  * tensor rank 1 (vector)
807  *
808  */
809 template <int Tensor_Dim, class T, class L, class A>
811  : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
812 
814  const std::string field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
815  const EntityType zero_type = MBVERTEX)
816  : OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A>(
817  field_name, data_ptr, zero_type) {}
818 };
819 
820 /** \brief Evaluate field gradient values for scalar field, i.e. gradient is
821  * tensor rank 1 (vector), specialization
822  *
823  */
824 template <int Tensor_Dim>
826  ublas::row_major, DoubleAllocator>
828  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
829 
831  const std::string field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
832  const EntityType zero_type = MBVERTEX)
834  ublas::row_major, DoubleAllocator>(
835  field_name, data_ptr, zero_type) {}
836 
837  /**
838  * \brief calculate gradient values of scalar field at integration points
839  * @param side side entity number
840  * @param type side entity type
841  * @param data entity data
842  * @return error code
843  */
844  MoFEMErrorCode doWork(int side, EntityType type,
846 };
847 
848 /**
849  * \brief Member function specialization calculating scalar field gradients for
850  * tenor field rank 1
851  *
852  */
853 template <int Tensor_Dim>
854 MoFEMErrorCode OpCalculateScalarFieldGradient_General<
855  Tensor_Dim, double, ublas::row_major,
856  DoubleAllocator>::doWork(int side, EntityType type,
859  const int nb_dofs = data.getFieldData().size();
860  if (!this->dataPtr) {
861  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
862  "Data pointer not allocated");
863  }
864  if (!nb_dofs && type == this->zeroType) {
865  this->dataPtr->resize(Tensor_Dim, 0, false);
867  }
868  if (!nb_dofs) {
870  }
871  const int nb_gauss_pts = data.getN().size1();
872  const int nb_base_functions = data.getN().size2();
873  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat =
874  *this->dataPtr;
875  if (type == this->zeroType) {
876  mat.resize(Tensor_Dim, nb_gauss_pts, false);
877  mat.clear();
878  }
879  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
880  auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
882  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
883  auto field_data = data.getFTensor0FieldData();
884  int bb = 0;
885  for (; bb < nb_dofs; ++bb) {
886  gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
887  ++field_data;
888  ++diff_base_function;
889  }
890  // Number of dofs can be smaller than number of base functions
891  for (; bb != nb_base_functions; ++bb)
892  ++diff_base_function;
893  ++gradients_at_gauss_pts;
894  }
896 }
897 
898 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
899  * vector field
900  *
901  * \ingroup mofem_forces_and_sources_user_data_operators
902  */
903 template <int Tensor_Dim>
906  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
907 
908  OpCalculateScalarFieldGradient(const std::string field_name,
909  boost::shared_ptr<MatrixDouble> &data_ptr,
910  const EntityType zero_type = MBVERTEX)
912  Tensor_Dim, double, ublas::row_major, DoubleAllocator>(
913  field_name, data_ptr, zero_type) {}
914 };
915 
916 /**
917  * \brief Evaluate field gradient values for vector field, i.e. gradient is
918  * tensor rank 2
919  *
920  */
921 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
923  : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
924  L, A> {
925 
927  const std::string field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
928  const EntityType zero_type = MBVERTEX)
929  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
930  A>(field_name, data_ptr,
931  zero_type) {}
932 };
933 
934 template <int Tensor_Dim0, int Tensor_Dim1>
935 struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
936  ublas::row_major, DoubleAllocator>
938  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
939 
941  const std::string field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
942  const EntityType zero_type = MBVERTEX)
943  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
944  ublas::row_major,
946  field_name, data_ptr, zero_type) {}
947 
948  /**
949  * \brief calculate values of vector field at integration points
950  * @param side side entity number
951  * @param type side entity type
952  * @param data entity data
953  * @return error code
954  */
955  MoFEMErrorCode doWork(int side, EntityType type,
957 };
958 
959 /**
960  * \brief Member function specialization calculating vector field gradients for
961  * tenor field rank 2
962  *
963  */
964 template <int Tensor_Dim0, int Tensor_Dim1>
965 MoFEMErrorCode OpCalculateVectorFieldGradient_General<
966  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
967  DoubleAllocator>::doWork(int side, EntityType type,
970  const int nb_dofs = data.getFieldData().size();
971  if (!nb_dofs && type == this->zeroType) {
972  this->dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
974  }
975  if (!nb_dofs) {
977  }
978  const int nb_gauss_pts = data.getN().size1();
979  const int nb_base_functions = data.getN().size2();
980  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat =
981  *this->dataPtr;
982  if (type == this->zeroType) {
983  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
984  mat.clear();
985  }
986  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
987  auto gradients_at_gauss_pts =
988  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
991  int size = nb_dofs / Tensor_Dim0;
992  if (nb_dofs % Tensor_Dim0) {
993  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
994  }
995  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
996  auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
997  int bb = 0;
998  for (; bb < size; ++bb) {
999  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1000  ++field_data;
1001  ++diff_base_function;
1002  }
1003  // Number of dofs can be smaller than number of Tensor_Dim0 x base functions
1004  for (; bb != nb_base_functions; ++bb)
1005  ++diff_base_function;
1006  ++gradients_at_gauss_pts;
1007  }
1009 }
1010 
1011 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
1012  * vector field
1013  *
1014  * \ingroup mofem_forces_and_sources_user_data_operators
1015  */
1016 template <int Tensor_Dim0, int Tensor_Dim1>
1019  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1020 
1021  OpCalculateVectorFieldGradient(const std::string field_name,
1022  boost::shared_ptr<MatrixDouble> &data_ptr,
1023  const EntityType zero_type = MBVERTEX)
1024  : OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
1025  ublas::row_major,
1026  DoubleAllocator>(
1027  field_name, data_ptr, zero_type) {}
1028 };
1029 
1030 /** \brief Get field gradients time derivative at integration pts for scalar
1031  * filed rank 0, i.e. vector field
1032  *
1033  * \ingroup mofem_forces_and_sources_user_data_operators
1034  */
1035 template <int Tensor_Dim0, int Tensor_Dim1>
1038 
1039  boost::shared_ptr<MatrixDouble> &dataPtr; ///< Data computed into this matrix
1040  EntityType zeroAtType; ///< Zero values at Gauss point at this type
1041  VectorDouble dotVector; ///< Keeps temoorary values of time directives
1042 
1043  template <int Dim> inline auto getFTensorDotData() {
1044  static_assert(Dim || !Dim, "not implemented");
1045  }
1046 
1047  OpCalculateVectorFieldGradientDot(const std::string field_name,
1048  boost::shared_ptr<MatrixDouble> &data_ptr,
1049  const EntityType zero_at_type = MBVERTEX)
1052  dataPtr(data_ptr), zeroAtType(zero_at_type) {
1053  if (!dataPtr)
1054  THROW_MESSAGE("Pointer is not set");
1055  }
1056 
1057  MoFEMErrorCode doWork(int side, EntityType type,
1060 
1061  const auto &local_indices = data.getLocalIndices();
1062  const int nb_dofs = local_indices.size();
1063  if (!nb_dofs && type == zeroAtType) {
1064  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
1066  }
1067  if (!nb_dofs)
1069 
1070  dotVector.resize(nb_dofs, false);
1071  const double *array;
1072  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1073  for (int i = 0; i != local_indices.size(); ++i)
1074  if (local_indices[i] != -1)
1075  dotVector[i] = array[local_indices[i]];
1076  else
1077  dotVector[i] = 0;
1078  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1079 
1080  const int nb_gauss_pts = data.getN().size1();
1081  const int nb_base_functions = data.getN().size2();
1082  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat = *dataPtr;
1083  if (type == zeroAtType) {
1084  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1085  mat.clear();
1086  }
1087  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1088  auto gradients_at_gauss_pts =
1089  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1092  int size = nb_dofs / Tensor_Dim0;
1093  if (nb_dofs % Tensor_Dim0) {
1094  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1095  }
1096  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1097  auto field_data = getFTensorDotData<Tensor_Dim0>();
1098  int bb = 0;
1099  for (; bb < size; ++bb) {
1100  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1101  ++field_data;
1102  ++diff_base_function;
1103  }
1104  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1105  // functions
1106  for (; bb != nb_base_functions; ++bb)
1107  ++diff_base_function;
1108  ++gradients_at_gauss_pts;
1109  }
1111  }
1112 };
1113 
1114 template <>
1115 template <>
1118  &dotVector[0], &dotVector[1], &dotVector[2]);
1119 }
1120 
1121 template <>
1122 template <>
1124  return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>(&dotVector[0],
1125  &dotVector[1]);
1126 }
1127 
1128 /** \brief Get vector field for H-div approximation
1129  *
1130  */
1131 template <int Tensor_Dim0, class T, class L, class A>
1134 
1135  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
1137  const int zeroSide;
1138 
1140  const std::string field_name,
1141  boost::shared_ptr<ublas::matrix<T, L, A>> &data_ptr,
1142  const EntityType zero_type = MBEDGE, const int zero_side = 0)
1145  dataPtr(data_ptr), zeroType(zero_type), zeroSide(0) {
1146  if (!dataPtr)
1147  THROW_MESSAGE("Pointer is not set");
1148  }
1149 
1150  /**
1151  * \brief calculate values of vector field at integration points
1152  * @param side side entity number
1153  * @param type side entity type
1154  * @param data entity data
1155  * @return error code
1156  */
1157  MoFEMErrorCode doWork(int side, EntityType type,
1159 };
1160 
1161 template <int Tensor_Dim, class T, class L, class A>
1163  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1165  SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1166  "Not implemented for T = %s and dim = %d",
1167  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
1168  Tensor_Dim);
1170 }
1171 
1172 /** \brief Get vector field for H-div approximation
1173  *
1174  */
1175 template <int Tensor_Dim>
1179 
1180  boost::shared_ptr<MatrixDouble> dataPtr;
1182  const int zeroSide;
1183 
1184  OpCalculateHdivVectorField_General(const std::string field_name,
1185  boost::shared_ptr<MatrixDouble> &data_ptr,
1186  const EntityType zero_type = MBEDGE,
1187  const int zero_side = 0)
1190  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1191  if (!dataPtr)
1192  THROW_MESSAGE("Pointer is not set");
1193  }
1194 
1195  /**
1196  * \brief Calculate values of vector field at integration points
1197  * @param side side entity number
1198  * @param type side entity type
1199  * @param data entity data
1200  * @return error code
1201  */
1202  MoFEMErrorCode doWork(int side, EntityType type,
1204 };
1205 
1206 template <int Tensor_Dim>
1207 MoFEMErrorCode OpCalculateHdivVectorField_General<
1208  Tensor_Dim, double, ublas::row_major,
1209  DoubleAllocator>::doWork(int side, EntityType type,
1212  const int nb_integration_points = getGaussPts().size2();
1213  if (type == zeroType && side == zeroSide) {
1214  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1215  dataPtr->clear();
1216  }
1217  const int nb_dofs = data.getFieldData().size();
1218  if (!nb_dofs)
1220  const int nb_base_functions = data.getN().size2() / Tensor_Dim;
1222  auto t_n_hdiv = data.getFTensor1N<Tensor_Dim>();
1223  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1224  for (int gg = 0; gg != nb_integration_points; ++gg) {
1225  auto t_dof = data.getFTensor0FieldData();
1226  int bb = 0;
1227  for (; bb != nb_dofs; ++bb) {
1228  t_data(i) += t_n_hdiv(i) * t_dof;
1229  ++t_n_hdiv;
1230  ++t_dof;
1231  }
1232  for (; bb != nb_base_functions; ++bb)
1233  ++t_n_hdiv;
1234  ++t_data;
1235  }
1237 }
1238 
1239 /** \brief Get vector field for H-div approximation
1240  *
1241  * \ingroup mofem_forces_and_sources_user_data_operators
1242  */
1243 template <int Tensor_Dim>
1246  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1247 
1248  OpCalculateHdivVectorField(const std::string field_name,
1249  boost::shared_ptr<MatrixDouble> &data_ptr,
1250  const EntityType zero_type = MBEDGE,
1251  const int zero_side = 0)
1252  : OpCalculateHdivVectorField_General<Tensor_Dim, double, ublas::row_major,
1253  DoubleAllocator>(
1254  field_name, data_ptr, zero_type, zero_side) {}
1255 };
1256 
1257 /**
1258  * @brief Calculate divergence of vector field
1259  * @ingroup mofem_forces_and_sources_user_data_operators
1260  *
1261  * @tparam Tensor_Dim dimension of space
1262  */
1263 template <int Tensor_Dim1, int Tensor_Dim2>
1266 
1267  boost::shared_ptr<VectorDouble> dataPtr;
1269  const int zeroSide;
1270 
1271  OpCalculateHdivVectorDivergence(const std::string field_name,
1272  boost::shared_ptr<VectorDouble> &data_ptr,
1273  const EntityType zero_type = MBEDGE,
1274  const int zero_side = 0)
1277  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1278  if (!dataPtr)
1279  THROW_MESSAGE("Pointer is not set");
1280  }
1281 
1282  MoFEMErrorCode doWork(int side, EntityType type,
1285  const int nb_integration_points = getGaussPts().size2();
1286  if (type == zeroType && side == zeroSide) {
1287  dataPtr->resize(nb_integration_points, false);
1288  dataPtr->clear();
1289  }
1290  const int nb_dofs = data.getFieldData().size();
1291  if (!nb_dofs)
1293  const int nb_base_functions = data.getN().size2() / Tensor_Dim1;
1295  auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
1296  auto t_data = getFTensor0FromVec(*dataPtr);
1297  for (int gg = 0; gg != nb_integration_points; ++gg) {
1298  auto t_dof = data.getFTensor0FieldData();
1299  int bb = 0;
1300  for (; bb != nb_dofs; ++bb) {
1301  double div = 0;
1302  for (auto ii = 0; ii != Tensor_Dim2; ++ii)
1303  div += t_n_diff_hdiv(ii, ii);
1304  t_data += t_dof * div;
1305  ++t_n_diff_hdiv;
1306  ++t_dof;
1307  }
1308  for (; bb != nb_base_functions; ++bb)
1309  ++t_n_diff_hdiv;
1310  ++t_data;
1311  }
1313  }
1314 };
1315 
1316 /**
1317  * @brief Calculate curl of vector field
1318  * @ingroup mofem_forces_and_sources_user_data_operators
1319  *
1320  * @tparam Tensor_Dim dimension of space
1321  */
1322 template <int Tensor_Dim>
1325 
1326  boost::shared_ptr<MatrixDouble> dataPtr;
1328  const int zeroSide;
1329 
1330  OpCalculateHcurlVectorCurl(const std::string field_name,
1331  boost::shared_ptr<MatrixDouble> &data_ptr,
1332  const EntityType zero_type = MBEDGE,
1333  const int zero_side = 0)
1336  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1337  if (!dataPtr)
1338  THROW_MESSAGE("Pointer is not set");
1339  }
1340 
1341  MoFEMErrorCode doWork(int side, EntityType type,
1344  const int nb_integration_points = getGaussPts().size2();
1345  if (type == zeroType && side == zeroSide) {
1346  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1347  dataPtr->clear();
1348  }
1349  const int nb_dofs = data.getFieldData().size();
1350  if (!nb_dofs)
1352 
1353  MatrixDouble curl_mat;
1354 
1355  const int nb_base_functions = data.getN().size2() / Tensor_Dim;
1359  auto t_n_diff_hcurl = data.getFTensor2DiffN<Tensor_Dim, Tensor_Dim>();
1360  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1361  for (int gg = 0; gg != nb_integration_points; ++gg) {
1362 
1363  // // get curl of base functions
1364  // CHKERR getCurlOfHCurlBaseFunctions(side, type, data, gg, curl_mat);
1365  // FTensor::Tensor1<double *, 3> t_base_curl(
1366  // &curl_mat(0, HVEC0), &curl_mat(0, HVEC1), &curl_mat(0, HVEC2), 3);
1367 
1368  auto t_dof = data.getFTensor0FieldData();
1369  int bb = 0;
1370  for (; bb != nb_dofs; ++bb) {
1371 
1372  // FTensor::Tensor1<double, 3> t_test;
1373  // t_test(k) = levi_civita(j, i, k) * t_n_diff_hcurl(i, j);
1374  // t_base_curl(k) -= t_test(k);
1375  // std::cerr << "error: " << t_base_curl(0) << " " << t_base_curl(1) <<
1376  // " "
1377  // << t_base_curl(2) << std::endl;
1378  // ++t_base_curl;
1379 
1380  t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
1381  ++t_n_diff_hcurl;
1382  ++t_dof;
1383  }
1384 
1385  for (; bb != nb_base_functions; ++bb)
1386  ++t_n_diff_hcurl;
1387  ++t_data;
1388  }
1389 
1391  }
1392 };
1393 
1394 /**
1395  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1396  * \ingroup mofem_forces_and_sources_user_data_operators
1397  *
1398  * @tparam Tensor_Dim0 rank of the filed
1399  * @tparam Tensor_Dim1 dimension of space
1400  */
1401 template <int Tensor_Dim0, int Tensor_Dim1>
1404 
1405  boost::shared_ptr<MatrixDouble> dataPtr;
1407  const int zeroSide;
1408 
1409  OpCalculateHVecTensorField(const std::string field_name,
1410  boost::shared_ptr<MatrixDouble> &data_ptr,
1411  const EntityType zero_type = MBEDGE,
1412  const int zero_side = 0)
1415  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1416  if (!dataPtr)
1417  THROW_MESSAGE("Pointer is not set");
1418  }
1419 
1420  MoFEMErrorCode doWork(int side, EntityType type,
1423  const int nb_integration_points = getGaussPts().size2();
1424  if (type == zeroType && side == zeroSide) {
1425  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1426  dataPtr->clear();
1427  }
1428  const int nb_dofs = data.getFieldData().size();
1429  if (!nb_dofs)
1431  const int nb_base_functions = data.getN().size2() / Tensor_Dim1;
1434  auto t_n_hvec = data.getFTensor1N<Tensor_Dim1>();
1435  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1436  for (int gg = 0; gg != nb_integration_points; ++gg) {
1437  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1438  int bb = 0;
1439  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1440  t_data(i, j) += t_dof(i) * t_n_hvec(j);
1441  ++t_n_hvec;
1442  ++t_dof;
1443  }
1444  for (; bb != nb_base_functions; ++bb)
1445  ++t_n_hvec;
1446  ++t_data;
1447  }
1449  }
1450 };
1451 
1452 /**
1453  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1454  * \ingroup mofem_forces_and_sources_user_data_operators
1455  *
1456  * @tparam Tensor_Dim0 rank of the filed
1457  * @tparam Tensor_Dim1 dimension of space
1458  */
1459 template <int Tensor_Dim0, int Tensor_Dim1>
1462 
1463  boost::shared_ptr<MatrixDouble> dataPtr;
1465  const int zeroSide;
1466 
1467  OpCalculateHTensorTensorField(const std::string field_name,
1468  boost::shared_ptr<MatrixDouble> &data_ptr,
1469  const EntityType zero_type = MBEDGE,
1470  const int zero_side = 0)
1473  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1474  if (!dataPtr)
1475  THROW_MESSAGE("Pointer is not set");
1476  }
1477 
1478  MoFEMErrorCode doWork(int side, EntityType type,
1481  const int nb_integration_points = getGaussPts().size2();
1482  if (type == zeroType && side == zeroSide) {
1483  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1484  dataPtr->clear();
1485  }
1486  const int nb_dofs = data.getFieldData().size();
1487  if (!nb_dofs)
1489  const int nb_base_functions =
1490  data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
1493  auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
1494  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1495  for (int gg = 0; gg != nb_integration_points; ++gg) {
1496  auto t_dof = data.getFTensor0FieldData();
1497  int bb = 0;
1498  for (; bb != nb_dofs; ++bb) {
1499  t_data(i, j) += t_dof * t_n_hten(i, j);
1500  ++t_n_hten;
1501  ++t_dof;
1502  }
1503  for (; bb != nb_base_functions; ++bb)
1504  ++t_n_hten;
1505  ++t_data;
1506  }
1508  }
1509 };
1510 
1511 /**
1512  * @brief Calculate divergence of tonsorial field using vectorial base
1513  * \ingroup mofem_forces_and_sources_user_data_operators
1514  *
1515  * @tparam Tensor_Dim0 rank of the field
1516  * @tparam Tensor_Dim1 dimension of space
1517  */
1518 template <int Tensor_Dim0, int Tensor_Dim1>
1521 
1522  boost::shared_ptr<MatrixDouble> dataPtr;
1524  const int zeroSide;
1525 
1526  OpCalculateHVecTensorDivergence(const std::string field_name,
1527  boost::shared_ptr<MatrixDouble> &data_ptr,
1528  const EntityType zero_type = MBEDGE,
1529  const int zero_side = 0)
1532  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1533  if (!dataPtr)
1534  THROW_MESSAGE("Pointer is not set");
1535  }
1536 
1537  MoFEMErrorCode doWork(int side, EntityType type,
1540  const int nb_integration_points = getGaussPts().size2();
1541  if (type == zeroType && side == 0) {
1542  dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
1543  dataPtr->clear();
1544  }
1545  const int nb_dofs = data.getFieldData().size();
1546  if (!nb_dofs)
1548  const int nb_base_functions = data.getN().size2() / Tensor_Dim1;
1551  auto t_n_diff_hvec = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim1>();
1552  auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
1553  for (int gg = 0; gg != nb_integration_points; ++gg) {
1554  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1555  int bb = 0;
1556  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1557  double div = t_n_diff_hvec(j, j);
1558  t_data(i) += t_dof(i) * div;
1559  ++t_n_diff_hvec;
1560  ++t_dof;
1561  }
1562  for (; bb != nb_base_functions; ++bb)
1563  ++t_n_diff_hvec;
1564  ++t_data;
1565  }
1567  }
1568 };
1569 
1570 // Face opeartors
1571 
1572 /** \brief Calculate jacobian for face element
1573 
1574  It is assumed that face element is XY plane. Applied
1575  only for 2d problems.
1576 
1577  \todo Generalize function for arbitrary face orientation in 3d space
1578 
1579  \ingroup mofem_forces_and_sources_tri_element
1580 
1581 */
1585 
1588 
1589  MoFEMErrorCode doWork(int side, EntityType type,
1591 };
1592 
1593 /** \brief Calculate inverse of jacobian for face element
1594 
1595  It is assumed that face element is XY plane. Applied
1596  only for 2d problems.
1597 
1598  \todo Generalize function for arbitrary face orientation in 3d space
1599 
1600  \ingroup mofem_forces_and_sources_tri_element
1601 
1602 */
1606 
1609  invJac(inv_jac) {}
1610 
1611  MoFEMErrorCode doWork(int side, EntityType type,
1613 };
1614 
1615 /** \brief Transform local reference derivatives of shape functions to global
1616 derivatives
1617 
1618 \ingroup mofem_forces_and_sources_tri_element
1619 
1620 */
1623 
1626  invJac(inv_jac) {}
1627 
1628  MoFEMErrorCode doWork(int side, EntityType type,
1630 
1631 private:
1634 };
1635 
1636 /**
1637  * \brief brief Transform local reference derivatives of shape function to
1638  global derivatives for face
1639 
1640  * \ingroup mofem_forces_and_sources_tri_element
1641  */
1644 
1647  invJac(inv_jac) {}
1648 
1649  MoFEMErrorCode doWork(int side, EntityType type,
1651 
1652 private:
1655 };
1656 
1657 /**
1658  * @brief Make Hdiv space from Hcurl space in 2d
1659  * @ingroup mofem_forces_and_sources_tri_element
1660  */
1663 
1666 
1667  MoFEMErrorCode doWork(int side, EntityType type,
1669 };
1670 
1671 /** \brief Apply contravariant (Piola) transfer to Hdiv space on face
1672  *
1673  * \note Hdiv space is generated by Hcurl space in 2d.
1674  *
1675  * Contravariant Piola transformation
1676  * \f[
1677  * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
1678  * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
1679  * =
1680  * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
1681  * \f]
1682  *
1683  * \ingroup mofem_forces_and_sources
1684  *
1685  */
1688 
1691  }
1692 
1695 
1696  MoFEMErrorCode doWork(int side, EntityType type,
1698 
1699  private:
1701 };
1702 
1703 // Edge
1704 
1707 
1710 
1711  MoFEMErrorCode doWork(int side, EntityType type,
1713 };
1714 
1715 // Fat prims
1716 
1717 /** \brief Calculate inverse of jacobian for face element
1718 
1719  It is assumed that face element is XY plane. Applied
1720  only for 2d problems.
1721 
1722  FIXME Generalize function for arbitrary face orientation in 3d space
1723  FIXME Calculate to Jacobins for two faces
1724 
1725  \ingroup mofem_forces_and_sources_prism_element
1726 
1727 */
1730 
1731  OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
1733  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
1734 
1737  invJac(inv_jac) {}
1738  MoFEMErrorCode doWork(int side, EntityType type,
1740 
1741  private:
1742  const boost::shared_ptr<MatrixDouble> invJacPtr;
1744 };
1745 
1746 /** \brief Transform local reference derivatives of shape functions to global
1747 derivatives
1748 
1749 FIXME Generalize to curved shapes
1750 FIXME Generalize to case that top and bottom face has different shape
1751 
1752 \ingroup mofem_forces_and_sources_prism_element
1753 
1754 */
1757 
1758  OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
1760  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
1761 
1764  invJac(inv_jac) {}
1765 
1766  MoFEMErrorCode doWork(int side, EntityType type,
1768 
1769  private:
1770  const boost::shared_ptr<MatrixDouble> invJacPtr;
1773 
1774 };
1775 
1776 // Flat prism
1777 
1778 /** \brief Calculate inverse of jacobian for face element
1779 
1780  It is assumed that face element is XY plane. Applied
1781  only for 2d problems.
1782 
1783  FIXME Generalize function for arbitrary face orientation in 3d space
1784  FIXME Calculate to Jacobins for two faces
1785 
1786  \ingroup mofem_forces_and_sources_prism_element
1787 
1788 */
1791 
1795  invJacF3(inv_jac_f3) {}
1796  MoFEMErrorCode doWork(int side, EntityType type,
1798 };
1799 
1800 /** \brief Transform local reference derivatives of shape functions to global
1801 derivatives
1802 
1803 FIXME Generalize to curved shapes
1804 FIXME Generalize to case that top and bottom face has different shape
1805 
1806 \ingroup mofem_forces_and_sources_prism_element
1807 
1808 */
1814  invJacF3(inv_jac_f3) {}
1815 
1817  MoFEMErrorCode doWork(int side, EntityType type,
1819 };
1820 
1821 } // namespace MoFEM
1822 
1823 #endif // __USER_DATA_OPERATORS_HPP__
1824 
1825 /**
1826  * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
1827  *
1828  * \brief Classes and functions used to evaluate fields at integration pts,
1829  *jacobians, etc..
1830  *
1831  * \ingroup mofem_forces_and_sources
1832  **/
OpCalculateScalarFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateScalarFieldValuesDot(const std::string field_name, boost::shared_ptr< VectorDouble > &data_ptr, const EntityType zero_at_type=MBVERTEX)
OpCalculateScalarFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
Calculate divergence of vector field.
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
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.
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.
Calculate field values for tenor field rank 1, i.e. vector field.
OpCalculateTensor2FieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_at_type=MBVERTEX)
OpCalculateVectorFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateTensor2FieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:506
const boost::shared_ptr< MatrixDouble > invJacPtr
OpCalculateInvJacForFlatPrism(MatrixDouble &inv_jac_f3)
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
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:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
EntityType zeroAtType
Zero values at Gauss point at this type.
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'I', 3 > I
Definition: PlasticOps.hpp:35
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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
OpCalculateHVecTensorDivergence(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
Data computed into this matrix.
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:625
OpCalculateHcurlVectorCurl(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
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:36
OpCalculateVectorFieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_at_type=MBVERTEX)
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.
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...
OpCalculateHdivVectorField(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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
OpCalculateVectorFieldGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
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.
OpCalculateJacForFace(MatrixDouble &jac)
Calculate inverse of jacobian for face element.
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)
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateVectorFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Index< 'k', 2 > k
Definition: ElasticOps.hpp:28
boost::shared_ptr< MatrixDouble > dataPtr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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.
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A >> &data_ptr, const EntityType zero_type=MBVERTEX)
OpSetInvJacH1ForFlatPrism(MatrixDouble &inv_jac_f3)
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:177
Get vector field for H-div approximation.
Calculate divergence of tonsorial field using vectorial base.
OpCalculateScalarFieldGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > &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)
Calculate field values for tenor field rank 2.
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A >> &data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateHdivVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Transform local reference derivatives of shape functions to global derivatives.
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:601
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< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:24
OpCalculateVectorFieldGradient_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.
OpCalculateInvJacForFatPrism(MatrixDouble &inv_jac)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
Apply contravariant (Piola) transfer to Hdiv space on face.
structure to get information form mofem into DataForcesAndSourcesCore
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
Data computed into this matrix.
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:412
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:176
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< MatrixDouble > dataPtr
OpCalculateHdivVectorDivergence(const std::string field_name, boost::shared_ptr< VectorDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
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
OpCalculateVectorFieldGradientDot(const std::string field_name, boost::shared_ptr< MatrixDouble > &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.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
std::vector< double, std::allocator< double > > DoubleAllocator
Definition: Types.hpp:70
Calculate symmetric tensor field values at integration pts.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
calculate values of scalar field at integration points
boost::shared_ptr< ublas::vector< T, A > > dataPtr
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.
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:25
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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
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)
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...
Make Hdiv space from Hcurl space in 2d.
OpCalculateTensor2SymmetricFieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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.