v0.9.0
UserDataOperators.hpp
Go to the documentation of this file.
1 /** \file UserDataOperators.hpp
2  * \brief User data Operators
3 
4 */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 #ifndef __USER_DATA_OPERATORS_HPP__
21 #define __USER_DATA_OPERATORS_HPP__
22 
23 namespace MoFEM {
24 
25 // GET VALUES AT GAUSS PTS
26 
27 // TENSOR0
28 
29 /** \brief Calculate field values for tenor field rank 0, i.e. scalar field
30  * \ingroup mofem_forces_and_sources_user_data_operators
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  * \ingroup mofem_forces_and_sources_user_data_operators
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  * \ingroup mofem_forces_and_sources_user_data_operators
78  */
80  : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
81 
82  OpCalculateScalarFieldValues(const std::string &field_name,
83  boost::shared_ptr<VectorDouble> &data_ptr,
84  const EntityType zero_type = MBVERTEX)
86  field_name, data_ptr, zero_type) {
87  if (!dataPtr)
88  THROW_MESSAGE("Pointer is not set");
89  }
90 
91  /**
92  * \brief calculate values of scalar field at integration points
93  * @param side side entity number
94  * @param type side entity type
95  * @param data entity data
96  * @return error code
97  */
98  MoFEMErrorCode doWork(int side, EntityType type,
101  VectorDouble &vec = *dataPtr;
102  const int nb_gauss_pts = getGaussPts().size2();
103  if (type == zeroType) {
104  vec.resize(nb_gauss_pts, false);
105  vec.clear();
106  }
107  const int nb_dofs = data.getFieldData().size();
108  if (nb_dofs) {
109  const int nb_gauss_pts = getGaussPts().size2();
110  const int nb_base_functions = data.getN().size2();
111  auto base_function = data.getFTensor0N();
112  auto values_at_gauss_pts = getFTensor0FromVec(vec);
113  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
114  auto field_data = data.getFTensor0FieldData();
115  int bb = 0;
116  for (; bb != nb_dofs; ++bb) {
117  values_at_gauss_pts += field_data * base_function;
118  ++field_data;
119  ++base_function;
120  }
121  // It is possible to have more base functions than dofs
122  for (; bb != nb_base_functions; ++bb)
123  ++base_function;
124  ++values_at_gauss_pts;
125  }
126  }
128  }
129 };
130 
133 
134  boost::shared_ptr<VectorDouble> dataPtr;
137 
138  OpCalculateScalarValuesDot(const std::string field_name,
139  boost::shared_ptr<VectorDouble> &data_ptr,
140  const EntityType zero_at_type = MBVERTEX)
143  dataPtr(data_ptr), zeroAtType(zero_at_type) {
144  if (!dataPtr)
145  THROW_MESSAGE("Pointer is not set");
146  }
147 
148  MoFEMErrorCode doWork(int side, EntityType type,
151 
152  const int nb_gauss_pts = getGaussPts().size2();
153  VectorDouble &vec = *dataPtr;
154  if (type == zeroAtType) {
155  vec.resize(nb_gauss_pts, false);
156  vec.clear();
157  }
158 
159  auto &local_indices = data.getLocalIndices();
160  const int nb_dofs = local_indices.size();
161  if (nb_dofs) {
162 
163  dotVector.resize(nb_dofs, false);
164  const double *array;
165  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
166  for (int i = 0; i != local_indices.size(); ++i)
167  if (local_indices[i] != -1)
168  dotVector[i] = array[local_indices[i]];
169  else
170  dotVector[i] = 0;
171  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
172 
173  const int nb_base_functions = data.getN().size2();
174  auto base_function = data.getFTensor0N();
175  auto values_at_gauss_pts = getFTensor0FromVec(vec);
176 
177  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
178  auto field_data = getFTensor0FromVec(dotVector);
179  int bb = 0;
180  for (; bb != nb_dofs; ++bb) {
181  values_at_gauss_pts += field_data * base_function;
182  ++field_data;
183  ++base_function;
184  }
185  // Number of dofs can be smaller than number of Tensor_Dim x base
186  // functions
187  for (; bb != nb_base_functions; ++bb)
188  ++base_function;
189  ++values_at_gauss_pts;
190  }
191  }
193  }
194 };
195 
196 // TENSOR1
197 
198 /** \brief Calculate field values for tenor field rank 1, i.e. vector field
199  * \ingroup mofem_forces_and_sources_user_data_operators
200  */
201 template <int Tensor_Dim, class T, class L, class A>
204 
205  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
207 
209  const std::string &field_name,
210  boost::shared_ptr<ublas::matrix<T, L, A>> &data_ptr,
211  const EntityType zero_type = MBVERTEX)
214  dataPtr(data_ptr), zeroType(zero_type) {
215  if (!dataPtr)
216  THROW_MESSAGE("Pointer is not set");
217  }
218 
219  /**
220  * \brief calculate values of vector field at integration points
221  * @param side side entity number
222  * @param type side entity type
223  * @param data entity data
224  * @return error code
225  */
226  MoFEMErrorCode doWork(int side, EntityType type,
228 };
229 
230 template <int Tensor_Dim, class T, class L, class A>
233  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
235  SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
236  "Not implemented for T = %s and dim = %d",
237  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
238  Tensor_Dim);
240 }
241 
242 /** \brief Calculate field values (template specialization) for tensor field
243  * rank 1, i.e. vector field \ingroup
244  * mofem_forces_and_sources_user_data_operators
245  */
246 template <int Tensor_Dim>
248  ublas::row_major, DoubleAllocator>
250 
251  boost::shared_ptr<MatrixDouble> dataPtr;
253 
255  const std::string &field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
256  const EntityType zero_type = MBVERTEX)
259  dataPtr(data_ptr), zeroType(zero_type) {
260  if (!dataPtr)
261  THROW_MESSAGE("Pointer is not set");
262  }
263 
264  MoFEMErrorCode doWork(int side, EntityType type,
266 };
267 
268 /**
269  * \brief Member function specialization calculating values for tenor field rank
270  * 1 \ingroup mofem_forces_and_sources_user_data_operators
271  */
272 template <int Tensor_Dim>
273 MoFEMErrorCode OpCalculateVectorFieldValues_General<
274  Tensor_Dim, double, ublas::row_major,
275  DoubleAllocator>::doWork(int side, EntityType type,
278  const int nb_dofs = data.getFieldData().size();
279  if (!nb_dofs && type == this->zeroType) {
280  dataPtr->resize(Tensor_Dim, 0, false);
282  }
283  if (!nb_dofs) {
285  }
286  const int nb_gauss_pts = data.getN().size1();
287  const int nb_base_functions = data.getN().size2();
288  MatrixDouble &mat = *dataPtr;
289  if (type == zeroType) {
290  mat.resize(Tensor_Dim, nb_gauss_pts, false);
291  mat.clear();
292  }
293  auto base_function = data.getFTensor0N();
294  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
296  const int size = nb_dofs / Tensor_Dim;
297  if (nb_dofs % Tensor_Dim) {
298  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
299  }
300  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
301  auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
302  int bb = 0;
303  for (; bb != size; ++bb) {
304  values_at_gauss_pts(I) += field_data(I) * base_function;
305  ++field_data;
306  ++base_function;
307  }
308  // Number of dofs can be smaller than number of Tensor_Dim x base functions
309  for (; bb != nb_base_functions; ++bb)
310  ++base_function;
311  ++values_at_gauss_pts;
312  }
314 }
315 
316 /** \brief Get values at integration pts for tensor filed rank 1, i.e. vector
317  * field \ingroup mofem_forces_and_sources_user_data_operators
318  */
319 template <int Tensor_Dim>
322  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
323 
324  OpCalculateVectorFieldValues(const std::string &field_name,
325  boost::shared_ptr<MatrixDouble> &data_ptr,
326  const EntityType zero_type = MBVERTEX)
328  ublas::row_major, DoubleAllocator>(
329  field_name, data_ptr, zero_type) {}
330 };
331 
332 /** \brief Get time direvatives of values at integration pts for tensor filed
333  * rank 1, i.e. vector field \ingroup
334  * mofem_forces_and_sources_user_data_operators
335  */
336 template <int Tensor_Dim>
339 
340  boost::shared_ptr<MatrixDouble> dataPtr;
343 
344  template <int Dim> inline auto getFTensorDotData() {
345  static_assert(Dim || !Dim, "not implemented");
346  }
347 
348  OpCalculateVectorFieldValuesDot(const std::string field_name,
349  boost::shared_ptr<MatrixDouble> &data_ptr,
350  const EntityType zero_at_type = MBVERTEX)
353  dataPtr(data_ptr), zeroAtType(zero_at_type) {
354  if (!dataPtr)
355  THROW_MESSAGE("Pointer is not set");
356  }
357 
358  MoFEMErrorCode doWork(int side, EntityType type,
361  auto &local_indices = data.getLocalIndices();
362  const int nb_dofs = local_indices.size();
363  if (!nb_dofs && type == zeroAtType) {
364  dataPtr->resize(Tensor_Dim, 0, false);
366  }
367  if (!nb_dofs)
369 
370  dotVector.resize(nb_dofs, false);
371  const double *array;
372  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
373  for (int i = 0; i != local_indices.size(); ++i)
374  if (local_indices[i] != -1)
375  dotVector[i] = array[local_indices[i]];
376  else
377  dotVector[i] = 0;
378  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
379 
380  const int nb_gauss_pts = data.getN().size1();
381  const int nb_base_functions = data.getN().size2();
382  MatrixDouble &mat = *dataPtr;
383  if (type == zeroAtType) {
384  mat.resize(Tensor_Dim, nb_gauss_pts, false);
385  mat.clear();
386  }
387  auto base_function = data.getFTensor0N();
388  auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
390  const int size = nb_dofs / Tensor_Dim;
391  if (nb_dofs % Tensor_Dim) {
392  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
393  }
394  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
395  auto field_data = getFTensorDotData<3>();
396  int bb = 0;
397  for (; bb != size; ++bb) {
398  values_at_gauss_pts(I) += field_data(I) * base_function;
399  ++field_data;
400  ++base_function;
401  }
402  // Number of dofs can be smaller than number of Tensor_Dim x base
403  // functions
404  for (; bb != nb_base_functions; ++bb)
405  ++base_function;
406  ++values_at_gauss_pts;
407  }
409  }
410 };
411 
412 template <>
413 template <>
416  &dotVector[0], &dotVector[1], &dotVector[2]);
417 }
418 
419 /** \brief Calculate field values for tenor field rank 2.
420  * \ingroup mofem_forces_and_sources_user_data_operators
421  */
422 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
425 
426  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
428 
430  const std::string &field_name,
431  boost::shared_ptr<ublas::matrix<T, L, A>> &data_ptr,
432  const EntityType zero_type = MBVERTEX)
435  dataPtr(data_ptr), zeroType(zero_type) {
436  if (!dataPtr)
437  THROW_MESSAGE("Pointer is not set");
438  }
439 
440  MoFEMErrorCode doWork(int side, EntityType type,
442 };
443 
444 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
447  doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
449  SETERRQ3(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
450  "Not implemented for T = %s, dim0 = %d and dim1 = %d",
451  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
452  Tensor_Dim0, Tensor_Dim1);
454 }
455 
456 template <int Tensor_Dim0, int Tensor_Dim1>
457 struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
458  ublas::row_major, DoubleAllocator>
460 
461  boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
464 
466  const std::string &field_name,
467  boost::shared_ptr<
468  ublas::matrix<double, ublas::row_major, DoubleAllocator>> &data_ptr,
469  const EntityType zero_type = MBVERTEX)
472  dataPtr(data_ptr), zeroType(zero_type) {
473  if (!dataPtr)
474  THROW_MESSAGE("Pointer is not set");
475  }
476 
477  MoFEMErrorCode doWork(int side, EntityType type,
479 };
480 
481 template <int Tensor_Dim0, int Tensor_Dim1>
482 MoFEMErrorCode OpCalculateTensor2FieldValues_General<
483  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
484  DoubleAllocator>::doWork(int side, EntityType type,
487  const int nb_dofs = data.getFieldData().size();
488  const int nb_gauss_pts = data.getN().size1();
489  if (!nb_dofs && type == this->zeroType) {
490  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
492  }
493  if (!nb_dofs) {
495  }
496  const int nb_base_functions = data.getN().size2();
497  MatrixDouble &mat = *dataPtr;
498  if (type == zeroType) {
499  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
500  mat.clear();
501  }
502  auto base_function = data.getFTensor0N();
503  auto values_at_gauss_pts = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
506  const int size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
507  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
508  auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
509  int bb = 0;
510  for (; bb != size; ++bb) {
511  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
512  ++field_data;
513  ++base_function;
514  }
515  for (; bb != nb_base_functions; ++bb)
516  ++base_function;
517  ++values_at_gauss_pts;
518  }
520 }
521 
522 /** \brief Get values at integration pts for tensor filed rank 2, i.e. matrix
523  * field \ingroup mofem_forces_and_sources_user_data_operators
524  */
525 template <int Tensor_Dim0, int Tensor_Dim1>
528  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
529 
530  OpCalculateTensor2FieldValues(const std::string &field_name,
531  boost::shared_ptr<MatrixDouble> &data_ptr,
532  const EntityType zero_type = MBVERTEX)
533  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
534  ublas::row_major,
536  field_name, data_ptr, zero_type) {}
537 };
538 
539 /** \brief Get time direvarive values at integration pts for tensor filed rank
540  * 2, i.e. matrix field
541  * \ingroup mofem_forces_and_sources_user_data_operators
542  */
543 template <int Tensor_Dim0, int Tensor_Dim1>
546 
547  boost::shared_ptr<MatrixDouble> &dataPtr; ///< Data computed into this matrix
548  EntityType zeroAtType; ///< Zero values at Gauss point at this type
549  VectorDouble dotVector; ///< Keeps temoorary values of time directives
550 
551  template <int Dim0, int Dim1> auto getFTensorDotData() {
552  static_assert(Dim0 || !Dim0 || Dim1 || !Dim1, "not implemented");
553  }
554 
555  OpCalculateTensor2FieldValuesDot(const std::string field_name,
556  boost::shared_ptr<MatrixDouble> &data_ptr,
557  const EntityType zero_at_type = MBVERTEX)
560  dataPtr(data_ptr), zeroAtType(zero_at_type) {
561  if (!dataPtr)
562  THROW_MESSAGE("Pointer is not set");
563  }
564 
565  MoFEMErrorCode doWork(int side, EntityType type,
568  const auto &local_indices = data.getLocalIndices();
569  const int nb_dofs = local_indices.size();
570  const int nb_gauss_pts = data.getN().size1();
571  if (!nb_dofs && type == zeroAtType) {
572  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
574  }
575  if (!nb_dofs) {
577  }
578 
579  dotVector.resize(nb_dofs, false);
580  const double *array;
581  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
582  for (int i = 0; i != local_indices.size(); ++i)
583  if (local_indices[i] != -1)
584  dotVector[i] = array[local_indices[i]];
585  else
586  dotVector[i] = 0;
587  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
588 
589  const int nb_base_functions = data.getN().size2();
590  MatrixDouble &mat = *dataPtr;
591  if (type == zeroAtType) {
592  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
593  mat.clear();
594  }
595  auto base_function = data.getFTensor0N();
596  auto values_at_gauss_pts =
597  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
600  const int size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
601  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
602  auto field_data = getFTensorDotData<Tensor_Dim0, Tensor_Dim1>();
603  int bb = 0;
604  for (; bb != size; ++bb) {
605  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
606  ++field_data;
607  ++base_function;
608  }
609  for (; bb != nb_base_functions; ++bb)
610  ++base_function;
611  ++values_at_gauss_pts;
612  }
614  }
615 };
616 
617 template <>
618 template <>
621  &dotVector[0], &dotVector[1], &dotVector[2],
622 
623  &dotVector[3], &dotVector[4], &dotVector[5],
624 
625  &dotVector[6], &dotVector[7], &dotVector[8]);
626 }
627 
628 template <int Tensor_Dim>
631 
632  boost::shared_ptr<MatrixDouble> dataPtr;
634  const int zeroSide;
635 
637  const std::string &field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
638  const EntityType zero_type = MBEDGE, const int zero_side = 0)
641  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
642  if (!dataPtr)
643  THROW_MESSAGE("Pointer is not set");
644  }
645 
646  MoFEMErrorCode doWork(int side, EntityType type,
649  MatrixDouble &mat = *dataPtr;
650  const int nb_gauss_pts = getGaussPts().size2();
651  if (type == this->zeroType && side == zeroSide) {
652  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
653  mat.clear();
654  }
655  const int nb_dofs = data.getFieldData().size();
656  if (!nb_dofs)
658 
659  const int nb_base_functions = data.getN().size2();
660  auto base_function = data.getFTensor0N();
661  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
664  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
665  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
666  auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
667  int bb = 0;
668  for (; bb != size; ++bb) {
669  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
670  ++field_data;
671  ++base_function;
672  }
673  for (; bb != nb_base_functions; ++bb)
674  ++base_function;
675  ++values_at_gauss_pts;
676  }
677 
679  }
680 };
681 
682 template <int Tensor_Dim>
685 
686  boost::shared_ptr<MatrixDouble> dataPtr;
688  const int zeroSide;
690 
691  template <int Dim> inline auto getFTensorDotData() {
692  static_assert(Dim || !Dim, "not implemented");
693  }
694 
696  const std::string &field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
697  const EntityType zero_type = MBEDGE, const int zero_side = 0)
700  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
701  if (!dataPtr)
702  THROW_MESSAGE("Pointer is not set");
703  }
704 
705  MoFEMErrorCode doWork(int side, EntityType type,
708  const int nb_gauss_pts = getGaussPts().size2();
709  MatrixDouble &mat = *dataPtr;
710  if (type == zeroType && side == zeroSide) {
711  mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
712  mat.clear();
713  }
714  auto &local_indices = data.getLocalIndices();
715  const int nb_dofs = local_indices.size();
716  if (!nb_dofs)
718 
719  dotVector.resize(nb_dofs, false);
720  const double *array;
721  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
722  for (int i = 0; i != local_indices.size(); ++i)
723  if (local_indices[i] != -1)
724  dotVector[i] = array[local_indices[i]];
725  else
726  dotVector[i] = 0;
727  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
728 
729  const int nb_base_functions = data.getN().size2();
730 
731  auto base_function = data.getFTensor0N();
732  auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
735  const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
736  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
737  auto field_data = getFTensorDotData<Tensor_Dim>();
738  int bb = 0;
739  for (; bb != size; ++bb) {
740  values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
741  ++field_data;
742  ++base_function;
743  }
744  for (; bb != nb_base_functions; ++bb)
745  ++base_function;
746  ++values_at_gauss_pts;
747  }
748 
750  }
751 };
752 
753 template <>
754 template <>
755 inline auto
758  &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
759  &dotVector[5]);
760 }
761 
762 // GET GRADIENTS AT GAUSS POINTS
763 
764 /**
765  * \brief Evaluate field gradient values for scalar field, i.e. gradient is
766  * tensor rank 1 (vector) \ingroup mofem_forces_and_sources_user_data_operators
767  */
768 template <int Tensor_Dim, class T, class L, class A>
770  : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
771 
773  const std::string &field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
774  const EntityType zero_type = MBVERTEX)
775  : OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A>(
776  field_name, data_ptr, zero_type) {}
777 };
778 
779 /** \brief Evaluate field gradient values for scalar field, i.e. gradient is
780  * tensor rank 1 (vector), specialization \ingroup
781  * mofem_forces_and_sources_user_data_operators
782  */
783 template <int Tensor_Dim>
785  ublas::row_major, DoubleAllocator>
787  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
788 
790  const std::string &field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
791  const EntityType zero_type = MBVERTEX)
793  ublas::row_major, DoubleAllocator>(
794  field_name, data_ptr, zero_type) {}
795 
796  /**
797  * \brief calculate gradient values of scalar field at integration points
798  * @param side side entity number
799  * @param type side entity type
800  * @param data entity data
801  * @return error code
802  */
803  MoFEMErrorCode doWork(int side, EntityType type,
805 };
806 
807 /**
808  * \brief Member function specialization calculating scalar field gradients for
809  * tenor field rank 1 \ingroup mofem_forces_and_sources_user_data_operators
810  */
811 template <int Tensor_Dim>
812 MoFEMErrorCode OpCalculateScalarFieldGradient_General<
813  Tensor_Dim, double, ublas::row_major,
814  DoubleAllocator>::doWork(int side, EntityType type,
817  const int nb_dofs = data.getFieldData().size();
818  if (!this->dataPtr) {
819  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
820  "Data pointer not allocated");
821  }
822  if (!nb_dofs && type == this->zeroType) {
823  this->dataPtr->resize(Tensor_Dim, 0, false);
825  }
826  if (!nb_dofs) {
828  }
829  const int nb_gauss_pts = data.getN().size1();
830  const int nb_base_functions = data.getN().size2();
831  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat =
832  *this->dataPtr;
833  if (type == this->zeroType) {
834  mat.resize(Tensor_Dim, nb_gauss_pts, false);
835  mat.clear();
836  }
837  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
838  auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
840  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
841  auto field_data = data.getFTensor0FieldData();
842  int bb = 0;
843  for (; bb < nb_dofs; ++bb) {
844  gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
845  ++field_data;
846  ++diff_base_function;
847  }
848  // Number of dofs can be smaller than number of base functions
849  for (; bb != nb_base_functions; ++bb)
850  ++diff_base_function;
851  ++gradients_at_gauss_pts;
852  }
854 }
855 
856 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
857  * vector field \ingroup mofem_forces_and_sources_user_data_operators
858  */
859 template <int Tensor_Dim>
862  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
863 
864  OpCalculateScalarFieldGradient(const std::string &field_name,
865  boost::shared_ptr<MatrixDouble> &data_ptr,
866  const EntityType zero_type = MBVERTEX)
868  Tensor_Dim, double, ublas::row_major, DoubleAllocator>(
869  field_name, data_ptr, zero_type) {}
870 };
871 
872 /**
873  * \brief Evaluate field gradient values for vector field, i.e. gradient is
874  * tensor rank 2 \ingroup mofem_forces_and_sources_user_data_operators
875  */
876 template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
878  : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
879  L, A> {
880 
882  const std::string &field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
883  const EntityType zero_type = MBVERTEX)
884  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
885  A>(field_name, data_ptr,
886  zero_type) {}
887 };
888 
889 template <int Tensor_Dim0, int Tensor_Dim1>
890 struct OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
891  ublas::row_major, DoubleAllocator>
893  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
894 
896  const std::string &field_name, boost::shared_ptr<MatrixDouble> &data_ptr,
897  const EntityType zero_type = MBVERTEX)
898  : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
899  ublas::row_major,
901  field_name, data_ptr, zero_type) {}
902 
903  /**
904  * \brief calculate values of vector field at integration points
905  * @param side side entity number
906  * @param type side entity type
907  * @param data entity data
908  * @return error code
909  */
910  MoFEMErrorCode doWork(int side, EntityType type,
912 };
913 
914 /**
915  * \brief Member function specialization calculating vector field gradients for
916  * tenor field rank 2 \ingroup mofem_forces_and_sources_user_data_operators
917  */
918 template <int Tensor_Dim0, int Tensor_Dim1>
919 MoFEMErrorCode OpCalculateVectorFieldGradient_General<
920  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
921  DoubleAllocator>::doWork(int side, EntityType type,
924  const int nb_dofs = data.getFieldData().size();
925  if (!nb_dofs && type == this->zeroType) {
926  this->dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
928  }
929  if (!nb_dofs) {
931  }
932  const int nb_gauss_pts = data.getN().size1();
933  const int nb_base_functions = data.getN().size2();
934  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat =
935  *this->dataPtr;
936  if (type == this->zeroType) {
937  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
938  mat.clear();
939  }
940  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
941  auto gradients_at_gauss_pts =
942  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
945  int size = nb_dofs / Tensor_Dim0;
946  if (nb_dofs % Tensor_Dim0) {
947  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
948  }
949  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
950  auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
951  int bb = 0;
952  for (; bb < size; ++bb) {
953  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
954  ++field_data;
955  ++diff_base_function;
956  }
957  // Number of dofs can be smaller than number of Tensor_Dim0 x base functions
958  for (; bb != nb_base_functions; ++bb)
959  ++diff_base_function;
960  ++gradients_at_gauss_pts;
961  }
963 }
964 
965 /** \brief Get field gradients at integration pts for scalar filed rank 0, i.e.
966  * vector field \ingroup mofem_forces_and_sources_user_data_operators
967  */
968 template <int Tensor_Dim0, int Tensor_Dim1>
971  Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
972 
973  OpCalculateVectorFieldGradient(const std::string &field_name,
974  boost::shared_ptr<MatrixDouble> &data_ptr,
975  const EntityType zero_type = MBVERTEX)
976  : OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, double,
977  ublas::row_major,
979  field_name, data_ptr, zero_type) {}
980 };
981 
982 /** \brief Get field gradients time derivative at integration pts for scalar
983  * filed rank 0, i.e. vector field \ingroup
984  * mofem_forces_and_sources_user_data_operators
985  */
986 template <int Tensor_Dim0, int Tensor_Dim1>
989 
990  boost::shared_ptr<MatrixDouble> &dataPtr; ///< Data computed into this matrix
991  EntityType zeroAtType; ///< Zero values at Gauss point at this type
992  VectorDouble dotVector; ///< Keeps temoorary values of time directives
993 
994  template <int Dim> inline auto getFTensorDotData() {
995  static_assert(Dim || !Dim, "not implemented");
996  }
997 
998  OpCalculateVectorFieldGradientDot(const std::string field_name,
999  boost::shared_ptr<MatrixDouble> &data_ptr,
1000  const EntityType zero_at_type = MBVERTEX)
1003  dataPtr(data_ptr), zeroAtType(zero_at_type) {
1004  if (!dataPtr)
1005  THROW_MESSAGE("Pointer is not set");
1006  }
1007 
1008  MoFEMErrorCode doWork(int side, EntityType type,
1011 
1012  const auto &local_indices = data.getLocalIndices();
1013  const int nb_dofs = local_indices.size();
1014  if (!nb_dofs && type == zeroAtType) {
1015  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, 0, false);
1017  }
1018  if (!nb_dofs)
1020 
1021  dotVector.resize(nb_dofs, false);
1022  const double *array;
1023  CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1024  for (int i = 0; i != local_indices.size(); ++i)
1025  if (local_indices[i] != -1)
1026  dotVector[i] = array[local_indices[i]];
1027  else
1028  dotVector[i] = 0;
1029  CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1030 
1031  const int nb_gauss_pts = data.getN().size1();
1032  const int nb_base_functions = data.getN().size2();
1033  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat = *dataPtr;
1034  if (type == zeroAtType) {
1035  mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1036  mat.clear();
1037  }
1038  auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1039  auto gradients_at_gauss_pts =
1040  getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1043  int size = nb_dofs / Tensor_Dim0;
1044  if (nb_dofs % Tensor_Dim0) {
1045  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1046  }
1047  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1048  auto field_data = getFTensorDotData<3>();
1049  int bb = 0;
1050  for (; bb < size; ++bb) {
1051  gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1052  ++field_data;
1053  ++diff_base_function;
1054  }
1055  // Number of dofs can be smaller than number of Tensor_Dim0 x base
1056  // functions
1057  for (; bb != nb_base_functions; ++bb)
1058  ++diff_base_function;
1059  ++gradients_at_gauss_pts;
1060  }
1062  }
1063 };
1064 
1065 template <>
1066 template <>
1069  &dotVector[0], &dotVector[1], &dotVector[2]);
1070 }
1071 
1072 /** \brief Get vector field for H-div approximation
1073  * \ingroup mofem_forces_and_sources_user_data_operators
1074  */
1075 template <int Tensor_Dim0, class T, class L, class A>
1078 
1079  boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
1081  const int zeroSide;
1082 
1084  const std::string &field_name,
1085  boost::shared_ptr<ublas::matrix<T, L, A>> &data_ptr,
1086  const EntityType zero_type = MBEDGE, const int zero_side = 0)
1089  dataPtr(data_ptr), zeroType(zero_type), zeroSide(0) {
1090  if (!dataPtr)
1091  THROW_MESSAGE("Pointer is not set");
1092  }
1093 
1094  /**
1095  * \brief calculate values of vector field at integration points
1096  * @param side side entity number
1097  * @param type side entity type
1098  * @param data entity data
1099  * @return error code
1100  */
1101  MoFEMErrorCode doWork(int side, EntityType type,
1103 };
1104 
1105 template <int Tensor_Dim, class T, class L, class A>
1107  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1109  SETERRQ2(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1110  "Not implemented for T = %s and dim = %d",
1111  typeid(T).name(), // boost::core::demangle(typeid(T).name()),
1112  Tensor_Dim);
1114 }
1115 
1116 /** \brief Get vector field for H-div approximation
1117  * \ingroup mofem_forces_and_sources_user_data_operators
1118  */
1119 template <int Tensor_Dim>
1123 
1124  boost::shared_ptr<MatrixDouble> dataPtr;
1126  const int zeroSide;
1127 
1128  OpCalculateHdivVectorField_General(const std::string &field_name,
1129  boost::shared_ptr<MatrixDouble> &data_ptr,
1130  const EntityType zero_type = MBEDGE,
1131  const int zero_side = 0)
1134  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1135  if (!dataPtr)
1136  THROW_MESSAGE("Pointer is not set");
1137  }
1138 
1139  /**
1140  * \brief Calculate values of vector field at integration points
1141  * @param side side entity number
1142  * @param type side entity type
1143  * @param data entity data
1144  * @return error code
1145  */
1146  MoFEMErrorCode doWork(int side, EntityType type,
1148 };
1149 
1150 template <int Tensor_Dim>
1151 MoFEMErrorCode OpCalculateHdivVectorField_General<
1152  Tensor_Dim, double, ublas::row_major,
1153  DoubleAllocator>::doWork(int side, EntityType type,
1156  const int nb_integration_points = getGaussPts().size2();
1157  if (type == zeroType && side == zeroSide) {
1158  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1159  dataPtr->clear();
1160  }
1161  const int nb_dofs = data.getFieldData().size();
1162  if (!nb_dofs)
1164  const int nb_base_functions = data.getN().size2() / Tensor_Dim;
1166  auto t_n_hdiv = data.getFTensor1N<Tensor_Dim>();
1167  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1168  for (int gg = 0; gg != nb_integration_points; ++gg) {
1169  auto t_dof = data.getFTensor0FieldData();
1170  int bb = 0;
1171  for (; bb != nb_dofs; ++bb) {
1172  t_data(i) += t_n_hdiv(i) * t_dof;
1173  ++t_n_hdiv;
1174  ++t_dof;
1175  }
1176  for (; bb != nb_base_functions; ++bb)
1177  ++t_n_hdiv;
1178  ++t_data;
1179  }
1181 }
1182 
1183 /** \brief Get vector field for H-div approximation
1184  * \ingroup mofem_forces_and_sources_user_data_operators
1185  */
1186 template <int Tensor_Dim>
1189  Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1190 
1191  OpCalculateHdivVectorField(const std::string &field_name,
1192  boost::shared_ptr<MatrixDouble> &data_ptr,
1193  const EntityType zero_type = MBEDGE,
1194  const int zero_side = 0)
1195  : OpCalculateHdivVectorField_General<Tensor_Dim, double, ublas::row_major,
1196  DoubleAllocator>(
1197  field_name, data_ptr, zero_type, zero_side) {}
1198 };
1199 
1200 /**
1201  * @brief Calculate divergence of vector field
1202  * @ingroup mofem_forces_and_sources_user_data_operators
1203  *
1204  * @tparam Tensor_Dim dimension of space
1205  */
1206 template <int Tensor_Dim1, int Tensor_Dim2>
1209 
1210  boost::shared_ptr<VectorDouble> dataPtr;
1212  const int zeroSide;
1213 
1214  OpCalculateHdivVectorDivergence(const std::string &field_name,
1215  boost::shared_ptr<VectorDouble> &data_ptr,
1216  const EntityType zero_type = MBEDGE,
1217  const int zero_side = 0)
1220  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1221  if (!dataPtr)
1222  THROW_MESSAGE("Pointer is not set");
1223  }
1224 
1225  MoFEMErrorCode doWork(int side, EntityType type,
1228  const int nb_integration_points = getGaussPts().size2();
1229  if (type == zeroType && side == zeroSide) {
1230  dataPtr->resize(nb_integration_points, false);
1231  dataPtr->clear();
1232  }
1233  const int nb_dofs = data.getFieldData().size();
1234  if (!nb_dofs)
1236  const int nb_base_functions = data.getN().size2() / Tensor_Dim1;
1238  auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
1239  auto t_data = getFTensor0FromVec(*dataPtr);
1240  for (int gg = 0; gg != nb_integration_points; ++gg) {
1241  auto t_dof = data.getFTensor0FieldData();
1242  int bb = 0;
1243  for (; bb != nb_dofs; ++bb) {
1244  double div = 0;
1245  for (auto ii = 0; ii != Tensor_Dim2; ++ii)
1246  div += t_n_diff_hdiv(ii, ii);
1247  t_data += t_dof * div;
1248  ++t_n_diff_hdiv;
1249  ++t_dof;
1250  }
1251  for (; bb != nb_base_functions; ++bb)
1252  ++t_n_diff_hdiv;
1253  ++t_data;
1254  }
1256  }
1257 };
1258 
1259 /**
1260  * @brief Calculate curl of vector field
1261  * @ingroup mofem_forces_and_sources_user_data_operators
1262  *
1263  * @tparam Tensor_Dim dimension of space
1264  */
1265 template <int Tensor_Dim>
1268 
1269  boost::shared_ptr<MatrixDouble> dataPtr;
1271  const int zeroSide;
1272 
1273  OpCalculateHcurlVectorCurl(const std::string &field_name,
1274  boost::shared_ptr<MatrixDouble> &data_ptr,
1275  const EntityType zero_type = MBEDGE,
1276  const int zero_side = 0)
1279  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1280  if (!dataPtr)
1281  THROW_MESSAGE("Pointer is not set");
1282  }
1283 
1284  MoFEMErrorCode doWork(int side, EntityType type,
1287  const int nb_integration_points = getGaussPts().size2();
1288  if (type == zeroType && side == zeroSide) {
1289  dataPtr->resize(Tensor_Dim, nb_integration_points, false);
1290  dataPtr->clear();
1291  }
1292  const int nb_dofs = data.getFieldData().size();
1293  if (!nb_dofs)
1295 
1296  MatrixDouble curl_mat;
1297 
1298  const int nb_base_functions = data.getN().size2() / Tensor_Dim;
1302  auto t_n_diff_hcurl = data.getFTensor2DiffN<Tensor_Dim, Tensor_Dim>();
1303  auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
1304  for (int gg = 0; gg != nb_integration_points; ++gg) {
1305 
1306  // // get curl of base functions
1307  // CHKERR getCurlOfHCurlBaseFunctions(side, type, data, gg, curl_mat);
1308  // FTensor::Tensor1<double *, 3> t_base_curl(
1309  // &curl_mat(0, HVEC0), &curl_mat(0, HVEC1), &curl_mat(0, HVEC2), 3);
1310 
1311  auto t_dof = data.getFTensor0FieldData();
1312  int bb = 0;
1313  for (; bb != nb_dofs; ++bb) {
1314 
1315  // FTensor::Tensor1<double, 3> t_test;
1316  // t_test(k) = levi_civita(j, i, k) * t_n_diff_hcurl(i, j);
1317  // t_base_curl(k) -= t_test(k);
1318  // std::cerr << "error: " << t_base_curl(0) << " " << t_base_curl(1) <<
1319  // " "
1320  // << t_base_curl(2) << std::endl;
1321  // ++t_base_curl;
1322 
1323  t_data(k) += t_dof * (levi_civita(j, i, k) * t_n_diff_hcurl(i, j));
1324  ++t_n_diff_hcurl;
1325  ++t_dof;
1326  }
1327 
1328  for (; bb != nb_base_functions; ++bb)
1329  ++t_n_diff_hcurl;
1330  ++t_data;
1331  }
1332 
1334  }
1335 };
1336 
1337 /**
1338  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1339  * \ingroup mofem_forces_and_sources_user_data_operators
1340  *
1341  * @tparam Tensor_Dim0 rank of the filed
1342  * @tparam Tensor_Dim1 dimension of space
1343  */
1344 template <int Tensor_Dim0, int Tensor_Dim1>
1347 
1348  boost::shared_ptr<MatrixDouble> dataPtr;
1350  const int zeroSide;
1351 
1352  OpCalculateHVecTensorField(const std::string &field_name,
1353  boost::shared_ptr<MatrixDouble> &data_ptr,
1354  const EntityType zero_type = MBEDGE,
1355  const int zero_side = 0)
1358  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1359  if (!dataPtr)
1360  THROW_MESSAGE("Pointer is not set");
1361  }
1362 
1363  MoFEMErrorCode doWork(int side, EntityType type,
1366  const int nb_integration_points = getGaussPts().size2();
1367  if (type == zeroType && side == zeroSide) {
1368  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1369  dataPtr->clear();
1370  }
1371  const int nb_dofs = data.getFieldData().size();
1372  if (!nb_dofs)
1374  const int nb_base_functions = data.getN().size2() / Tensor_Dim1;
1377  auto t_n_hvec = data.getFTensor1N<Tensor_Dim1>();
1378  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1379  for (int gg = 0; gg != nb_integration_points; ++gg) {
1380  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1381  int bb = 0;
1382  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1383  t_data(i, j) += t_dof(i) * t_n_hvec(j);
1384  ++t_n_hvec;
1385  ++t_dof;
1386  }
1387  for (; bb != nb_base_functions; ++bb)
1388  ++t_n_hvec;
1389  ++t_data;
1390  }
1392  }
1393 };
1394 
1395 /**
1396  * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
1397  * \ingroup mofem_forces_and_sources_user_data_operators
1398  *
1399  * @tparam Tensor_Dim0 rank of the filed
1400  * @tparam Tensor_Dim1 dimension of space
1401  */
1402 template <int Tensor_Dim0, int Tensor_Dim1>
1405 
1406  boost::shared_ptr<MatrixDouble> dataPtr;
1408  const int zeroSide;
1409 
1410  OpCalculateHTensorTensorField(const std::string &field_name,
1411  boost::shared_ptr<MatrixDouble> &data_ptr,
1412  const EntityType zero_type = MBEDGE,
1413  const int zero_side = 0)
1416  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1417  if (!dataPtr)
1418  THROW_MESSAGE("Pointer is not set");
1419  }
1420 
1421  MoFEMErrorCode doWork(int side, EntityType type,
1424  const int nb_integration_points = getGaussPts().size2();
1425  if (type == zeroType && side == zeroSide) {
1426  dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
1427  dataPtr->clear();
1428  }
1429  const int nb_dofs = data.getFieldData().size();
1430  if (!nb_dofs)
1432  const int nb_base_functions =
1433  data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
1436  auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
1437  auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
1438  for (int gg = 0; gg != nb_integration_points; ++gg) {
1439  auto t_dof = data.getFTensor0FieldData();
1440  int bb = 0;
1441  for (; bb != nb_dofs; ++bb) {
1442  t_data(i, j) += t_dof * t_n_hten(i, j);
1443  ++t_n_hten;
1444  ++t_dof;
1445  }
1446  for (; bb != nb_base_functions; ++bb)
1447  ++t_n_hten;
1448  ++t_data;
1449  }
1451  }
1452 };
1453 
1454 /**
1455  * @brief Calculate divergence of tonsorial field using vectorial base
1456  * \ingroup mofem_forces_and_sources_user_data_operators
1457  *
1458  * @tparam Tensor_Dim0 rank of the field
1459  * @tparam Tensor_Dim1 dimension of space
1460  */
1461 template <int Tensor_Dim0, int Tensor_Dim1>
1464 
1465  boost::shared_ptr<MatrixDouble> dataPtr;
1467  const int zeroSide;
1468 
1469  OpCalculateHVecTensorDivergence(const std::string &field_name,
1470  boost::shared_ptr<MatrixDouble> &data_ptr,
1471  const EntityType zero_type = MBEDGE,
1472  const int zero_side = 0)
1475  dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1476  if (!dataPtr)
1477  THROW_MESSAGE("Pointer is not set");
1478  }
1479 
1480  MoFEMErrorCode doWork(int side, EntityType type,
1483  const int nb_integration_points = getGaussPts().size2();
1484  if (type == zeroType && side == 0) {
1485  dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
1486  dataPtr->clear();
1487  }
1488  const int nb_dofs = data.getFieldData().size();
1489  if (!nb_dofs)
1491  const int nb_base_functions = data.getN().size2() / Tensor_Dim1;
1494  auto t_n_diff_hvec = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim1>();
1495  auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
1496  for (int gg = 0; gg != nb_integration_points; ++gg) {
1497  auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
1498  int bb = 0;
1499  for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
1500  double div = t_n_diff_hvec(j, j);
1501  t_data(i) += t_dof(i) * div;
1502  ++t_n_diff_hvec;
1503  ++t_dof;
1504  }
1505  for (; bb != nb_base_functions; ++bb)
1506  ++t_n_diff_hvec;
1507  ++t_data;
1508  }
1510  }
1511 };
1512 
1513 // Face opeartors
1514 
1515 /** \brief Calculate jacobian for face element
1516 
1517  It is assumed that face element is XY plane. Applied
1518  only for 2d problems.
1519 
1520  \todo Generalize function for arbitrary face orientation in 3d space
1521 
1522  \ingroup mofem_forces_and_sources_tri_element
1523 
1524 */
1528 
1531 
1532  MoFEMErrorCode doWork(int side, EntityType type,
1534 };
1535 
1536 /** \brief Calculate inverse of jacobian for face element
1537 
1538  It is assumed that face element is XY plane. Applied
1539  only for 2d problems.
1540 
1541  \todo Generalize function for arbitrary face orientation in 3d space
1542 
1543  \ingroup mofem_forces_and_sources_tri_element
1544 
1545 */
1549 
1552  invJac(inv_jac) {}
1553 
1554  MoFEMErrorCode doWork(int side, EntityType type,
1556 };
1557 
1558 /** \brief Transform local reference derivatives of shape functions to global
1559 derivatives
1560 
1561 \ingroup mofem_forces_and_sources_tri_element
1562 
1563 */
1566 
1569  invJac(inv_jac) {}
1570 
1571  MoFEMErrorCode doWork(int side, EntityType type,
1573 
1574 private:
1577 };
1578 
1579 /**
1580  * \brief brief Transform local reference derivatives of shape function to
1581  global derivatives for face
1582 
1583  * \ingroup mofem_forces_and_sources_tri_element
1584  */
1587 
1590  invJac(inv_jac) {}
1591 
1592  MoFEMErrorCode doWork(int side, EntityType type,
1594 
1595 private:
1598 };
1599 
1600 /**
1601  * @brief Make Hdiv space from Hcurl space in 2d
1602  * @ingroup mofem_forces_and_sources_tri_element
1603  */
1606 
1609 
1610  MoFEMErrorCode doWork(int side, EntityType type,
1612 };
1613 
1614 /** \brief Apply contravariant (Piola) transfer to Hdiv space on face
1615  *
1616  * \note Hdiv space is generated by Hcurl space in 2d.
1617  *
1618  * Contravariant Piola transformation
1619  * \f[
1620  * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
1621  * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
1622  * =
1623  * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
1624  * \f]
1625  *
1626  * \ingroup mofem_forces_and_sources
1627  *
1628  */
1631 
1634  }
1635 
1638 
1639  MoFEMErrorCode doWork(int side, EntityType type,
1641 
1642  private:
1644 };
1645 
1646 // Edge
1647 
1650 
1653 
1654  MoFEMErrorCode doWork(int side, EntityType type,
1656 };
1657 
1658 // Fat prims
1659 
1660 /** \brief Calculate inverse of jacobian for face element
1661 
1662  It is assumed that face element is XY plane. Applied
1663  only for 2d problems.
1664 
1665  FIXME Generalize function for arbitrary face orientation in 3d space
1666  FIXME Calculate to Jacobins for two faces
1667 
1668  \ingroup mofem_forces_and_sources_prism_element
1669 
1670 */
1673 
1674  OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
1676  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
1677 
1680  invJac(inv_jac) {}
1681  MoFEMErrorCode doWork(int side, EntityType type,
1683 
1684  private:
1685  const boost::shared_ptr<MatrixDouble> invJacPtr;
1687 };
1688 
1689 /** \brief Transform local reference derivatives of shape functions to global
1690 derivatives
1691 
1692 FIXME Generalize to curved shapes
1693 FIXME Generalize to case that top and bottom face has different shape
1694 
1695 \ingroup mofem_forces_and_sources_prism_element
1696 
1697 */
1700 
1701  OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
1703  invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
1704 
1707  invJac(inv_jac) {}
1708 
1709  MoFEMErrorCode doWork(int side, EntityType type,
1711 
1712  private:
1713  const boost::shared_ptr<MatrixDouble> invJacPtr;
1716 
1717 };
1718 
1719 // Flat prism
1720 
1721 /** \brief Calculate inverse of jacobian for face element
1722 
1723  It is assumed that face element is XY plane. Applied
1724  only for 2d problems.
1725 
1726  FIXME Generalize function for arbitrary face orientation in 3d space
1727  FIXME Calculate to Jacobins for two faces
1728 
1729  \ingroup mofem_forces_and_sources_prism_element
1730 
1731 */
1734 
1738  invJacF3(inv_jac_f3) {}
1739  MoFEMErrorCode doWork(int side, EntityType type,
1741 };
1742 
1743 /** \brief Transform local reference derivatives of shape functions to global
1744 derivatives
1745 
1746 FIXME Generalize to curved shapes
1747 FIXME Generalize to case that top and bottom face has different shape
1748 
1749 \ingroup mofem_forces_and_sources_prism_element
1750 
1751 */
1757  invJacF3(inv_jac_f3) {}
1758 
1760  MoFEMErrorCode doWork(int side, EntityType type,
1762 };
1763 
1764 } // namespace MoFEM
1765 
1766 #endif // __USER_DATA_OPERATORS_HPP__
1767 
1768 /**
1769  * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
1770  *
1771  * \brief Classes and functions used to evaluate fields at integration pts,
1772  *jacobians, etc..
1773  *
1774  * \ingroup mofem_forces_and_sources
1775  **/
Calculate divergence of vector field.
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
OpCalculateTensor2FieldValues_General(const std::string &field_name, boost::shared_ptr< ublas::matrix< T, L, A >> &data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateScalarFieldValues(const std::string &field_name, boost::shared_ptr< VectorDouble > &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 mofem_fo...
OpCalculateScalarFieldValues_General(const std::string &field_name, boost::shared_ptr< ublas::vector< T, A >> &data_ptr, const EntityType zero_type=MBVERTEX)
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:142
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateScalarFieldGradient(const std::string &field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
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.
OpCalculateHVecTensorField(const std::string &field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHdivVectorField_General(const std::string &field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
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:477
EntityType zeroAtType
Zero values at Gauss point at this type.
OpCalculateVectorFieldGradient_General(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.
OpCalculateScalarValuesDot(const std::string field_name, boost::shared_ptr< VectorDouble > &data_ptr, const EntityType zero_at_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::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
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
Data computed into this matrix.
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:620
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)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
OpCalculateScalarFieldGradient_General(const std::string &field_name, boost::shared_ptr< MatrixDouble > &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.
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...
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
OpSetInvJacHcurlFace(MatrixDouble &inv_jac)
OpCalculateTensor2SymmetricFieldValues(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.
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
OpCalculateHdivVectorDivergence(const std::string &field_name, boost::shared_ptr< VectorDouble > &data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
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)
OpCalculateScalarFieldGradient_General(const std::string &field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< MatrixDouble > dataPtr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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)
OpCalculateVectorFieldGradient_General(const std::string &field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateHVecTensorDivergence(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:172
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.
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:596
OpCalculateVectorFieldGradient(const std::string &field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateTensor2FieldValues_General(const std::string &field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator >> &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.
Transform local reference derivatives of shape functions to global derivatives.
OpCalculateHcurlVectorCurl(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.
OpCalculateVectorFieldValues_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.
OpCalculateTensor2SymmetricFieldValuesDot(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 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.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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.
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:407
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:171
boost::shared_ptr< MatrixDouble > dataPtr
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
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
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)
calculate values of scalar field at integration points
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.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorField(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.
OpCalculateTensor2FieldValues(const std::string &field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
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...
Make Hdiv space from Hcurl space in 2d.
OpCalculateVectorFieldValues(const std::string &field_name, boost::shared_ptr< MatrixDouble > &data_ptr, const EntityType zero_type=MBVERTEX)
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.