v0.10.0
DataOperators.hpp
Go to the documentation of this file.
1 /** \file DataOperators.hpp
2 
3  Number of structures preforming operations on integration points.
4 
5  For example:
6  - calculate Jacobian
7  - calculate Piola-Transform on shape functions
8  - calculate normals and tangent vectors at integration points on triangle
9 
10  */
11 
12 /* This file is part of MoFEM.
13  * MoFEM is free software: you can redistribute it and/or modify it under
14  * the terms of the GNU Lesser General Public License as published by the
15  * Free Software Foundation, either version 3 of the License, or (at your
16  * option) any later version.
17  *
18  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
19  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21  * License for more details.
22  *
23  * You should have received a copy of the GNU Lesser General Public
24  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
25 
26 #ifndef __DATAOPERATORS_HPP
27 #define __DATAOPERATORS_HPP
28 
29 using namespace boost::numeric;
30 
31 namespace MoFEM {
32 
33 /** \brief base operator to do operations at Gauss Pt. level
34  * \ingroup mofem_forces_and_sources
35  */
36 struct DataOperator {
37 
38  DataOperator(const bool symm = true);
39 
40  virtual ~DataOperator() = default;
41 
42  /** \brief Operator for bi-linear form, usually to calculate values on left
43  * hand side
44  */
45  virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
46  EntityType col_type,
50  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
51  "doWork function is not implemented for this operator");
53  }
54 
55  virtual MoFEMErrorCode opLhs(DataForcesAndSourcesCore &row_data,
56  DataForcesAndSourcesCore &col_data);
57 
58  /** \brief Operator for linear form, usually to calculate values on right hand
59  * side
60  */
61  virtual MoFEMErrorCode doWork(int side, EntityType type,
64  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
65  "doWork function is not implemented for this operator");
67  }
68 
69  virtual MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data,
70  const bool error_if_no_base = false);
71 
72  bool sYmm; ///< If true assume that matrix is symmetric structure
73 
74  std::array<bool, MBMAXTYPE>
75  doEntities; ///< If true operator is executed for entity.
76 
77  // Deprecated variables. Use doEntities instead. I keep them for back
78  // compatibility with some older modules. It will be removed in some future.
79 
80  bool &doVertices; ///< \deprectaed If false skip vertices
81  bool &doEdges; ///< \deprectaed If false skip edges
82  bool &doQuads; ///< \deprectaed
83  bool &doTris; ///< \deprectaed
84  bool &doTets; ///< \deprectaed
85  bool &doPrisms; ///< \deprectaed
86 
87  /**
88  * \brief Get if operator uses symmetry of DOFs or not
89  *
90  * If symmetry is used, only not repeating combinations of entities are
91  * looped. For an example pair of (Vertex, Edge_0) and (Edge_0, Vertex) will
92  * calculate the same matrices only transposed. Implementing that this can be
93  * exploited by integrating only one pair.
94  *
95  * @return true if symmetry
96  */
97  inline bool getSymm() const { return sYmm; }
98 
99  /// set if operator is executed taking in account symmetry
100  inline void setSymm() { sYmm = true; }
101 
102  /// unset if operator is executed for non symmetric problem
103  inline void unSetSymm() { sYmm = false; }
104 
105 private:
106  template <bool Symm>
107  inline MoFEMErrorCode opLhs(DataForcesAndSourcesCore &row_data,
108  DataForcesAndSourcesCore &col_data);
109 
110  template <bool ErrorIfNoBase>
111  inline MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data,
112  const std::array<bool, MBMAXTYPE> &do_entities);
113 };
114 
115 /**
116  * \brief Transform local reference derivatives of shape function to global
117  derivatives
118 
119  * \ingroup mofem_forces_and_sources
120  */
121 struct OpSetInvJacH1 : public DataOperator {
122 
126 
128  : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
129  &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
130  &inv_jac(2, 2)) {}
131 
133  MoFEMErrorCode doWork(int side, EntityType type,
135 };
136 
137 /**
138  * \brief brief Transform local reference derivatives of shape function to
139  global derivatives
140 
141  * \ingroup mofem_forces_and_sources
142  */
144 
149 
151  : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
152  &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
153  &inv_jac(2, 2)) {}
154 
156  MoFEMErrorCode doWork(int side, EntityType type,
158 };
159 
160 /**
161  * \brief transform local reference derivatives of shape function to global
162  derivatives if higher order geometry is given
163 
164  * \ingroup mofem_forces_and_sources
165 */
166 struct OpSetHoInvJacH1 : public DataOperator {
167 
171  OpSetHoInvJacH1(MatrixDouble &inv_ho_jac) : invHoJac(inv_ho_jac) {}
172 
174  MoFEMErrorCode doWork(int side, EntityType type,
176 };
177 
178 /**
179  * \brief transform local reference derivatives of shape function to global
180  derivatives if higher order geometry is given
181  *
182 
183  * \ingroup mofem_forces_and_sources
184 */
186 
191 
192  OpSetHoInvJacHdivAndHcurl(MatrixDouble &inv_ho_jac) : invHoJac(inv_ho_jac) {}
193 
195  MoFEMErrorCode doWork(int side, EntityType type,
197 };
198 
199 /** \brief apply contravariant (Piola) transfer to Hdiv space
200 
201 Contravariant Piola transformation
202 \f[
203 \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
204 \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
205 =
206 \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
207 \f]
208 
209 * \ingroup mofem_forces_and_sources
210 
211 */
213 
214  double &vOlume;
215 
220 
222  : vOlume(volume),
223  // jAc(jac),
224  tJac(&jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
225  &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2)) {}
226 
229 
230  MoFEMErrorCode doWork(int side, EntityType type,
232 };
233 
234 /** \brief Apply contravariant (Piola) transfer to Hdiv space for HO geometr
235 
236 * \ingroup mofem_forces_and_sources
237 */
239 
245 
247  : detHoJac(det_jac), hoJac(jac) {}
248 
251  MoFEMErrorCode doWork(int side, EntityType type,
253 };
254 
255 /** \brief Apply covariant (Piola) transfer to Hcurl space for HO geometry
256  */
258 
263 
264  OpSetHoCovariantPiolaTransform(MatrixDouble &inv_jac) : hoInvJac(inv_jac) {}
265 
268  MoFEMErrorCode doWork(int side, EntityType type,
270 };
271 
272 /** \brief apply covariant transfer to Hcurl space
273 
274 Contravariant Piola transformation
275 \f[
276 \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
277 \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
278 =
279 \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
280 \f]
281 
282 
283 * \ingroup mofem_forces_and_sources
284 */
286 
291 
293  : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
294  &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
295  &inv_jac(2, 2)) {}
296 
299 
300  MoFEMErrorCode doWork(int side, EntityType type,
302 };
303 
304 /** \brief Get field values and gradients at Gauss points
305  * \ingroup mofem_forces_and_sources
306  */
307 template <int RANK, int DIM> struct OpGetDataAndGradient : public DataOperator {
308 
311 
313  MatrixDouble &data_grad_at_gauss_pt)
314  : dataAtGaussPts(data_at_gauss_pt),
315  dataGradAtGaussPts(data_grad_at_gauss_pt) {}
316 
317  /**
318  * Return tensor associated with matrix storing values
319  */
320  template <int R>
322  THROW_MESSAGE("Not implemented");
323  }
324 
325  /**
326  * Return tensor associated with matrix storing gradient values
327  */
328  template <int R, int D>
330  THROW_MESSAGE("Not implemented");
331  }
332 
333  /**
334  * \brief Calculate gradient and values at integration points
335  * @param side side of entity on element
336  * @param type type of entity
337  * @param data data stored on entity (dofs values, dofs indices, etc.)
338  * @return error code
339  */
343  const int nb_base_functions = data.getN().size2();
344  bool constant_diff = false;
345  if (type == MBVERTEX && data.getDiffN().size1() * data.getDiffN().size2() ==
346  DIM * nb_base_functions) {
347  constant_diff = true;
348  }
349  const int nb_dofs = data.getFieldData().size();
350  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
351  double *data_ptr, *n_ptr, *diff_n_ptr;
352  n_ptr = &data.getN()(gg, 0);
353  if (constant_diff) {
354  diff_n_ptr = &data.getDiffN()(0, 0);
355  } else {
356  diff_n_ptr = &data.getDiffN()(gg, 0);
357  }
358  data_ptr = &*data.getFieldData().data().begin();
359  for (int rr = 0; rr < RANK; rr++, data_ptr++) {
360  dataAtGaussPts(gg, rr) +=
361  cblas_ddot(nb_dofs / RANK, n_ptr, 1, data_ptr, RANK);
362  double *diff_n_ptr2 = diff_n_ptr;
363  for (unsigned int dd = 0; dd < DIM; dd++, diff_n_ptr2++) {
364  dataGradAtGaussPts(gg, DIM * rr + dd) +=
365  cblas_ddot(nb_dofs / RANK, diff_n_ptr2, DIM, data_ptr, RANK);
366  }
367  }
368  }
370  }
371 
372  MoFEMErrorCode doWork(int side, EntityType type,
374 
376 
377  if (data.getFieldData().size() == 0) {
379  }
380 
381  unsigned int nb_dofs = data.getFieldData().size();
382  if (nb_dofs == 0)
384 
385  if (nb_dofs % RANK != 0) {
386  SETERRQ4(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
387  "data inconsistency, type %d, side %d, nb_dofs %d, rank %d",
388  type, side, nb_dofs, RANK);
389  }
390  if (nb_dofs / RANK > data.getN().size2()) {
391  std::cerr << side << " " << type << " "
392  << ApproximationBaseNames[data.getBase()] << std::endl;
393  std::cerr << data.getN() << std::endl;
394  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
395  "data inconsistency nb_dofs >= data.N.size2(), i.e. %u >= %u",
396  nb_dofs, data.getN().size2());
397  }
398 
399  if (type == MBVERTEX) {
400  dataAtGaussPts.resize(data.getN().size1(), RANK, false);
401  dataGradAtGaussPts.resize(data.getN().size1(), RANK * DIM, false);
402  dataAtGaussPts.clear();
403  dataGradAtGaussPts.clear();
404  }
405 
406  CHKERR calculateValAndGrad(side, type, data);
407 
409  }
410 };
411 
412 /**
413  * \brief Specialization for field with 3 coefficients in 3 dimension
414  */
415 template <>
416 template <>
418 OpGetDataAndGradient<3, 3>::getValAtGaussPtsTensor<3>(MatrixDouble &data);
419 
420 /**
421  * \brief Specialization for field with 3 coefficients in 3 dimension
422  */
423 template <>
424 template <>
426 OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(MatrixDouble &data);
427 
428 /**
429  * \brief Specialization for field with 3 coefficients in 3 dimension
430  */
431 template <>
433  int side, EntityType type, DataForcesAndSourcesCore::EntData &data);
434 
435 /**
436  * \brief Specialization for field with for scalar field in 3 dimension
437  */
438 template <>
439 MoFEMErrorCode OpGetDataAndGradient<1, 3>::calculateValAndGrad(
440  int side, EntityType type, DataForcesAndSourcesCore::EntData &data);
441 
442 /** \brief Calculate normals at Gauss points of triangle element
443  * \ingroup mofem_forces_and_source
444  */
446 
451 
453  MatrixDouble &normals_at_gausspt,
454  MatrixDouble &tangent1_at_gausspt,
455  MatrixDouble &tangent2_at_gausspt)
456  : cOords_at_GaussPt(coords_at_gausspt),
457  nOrmals_at_GaussPt(normals_at_gausspt),
458  tAngent1_at_GaussPt(tangent1_at_gausspt),
459  tAngent2_at_GaussPt(tangent2_at_gausspt) {}
460 
462  MoFEMErrorCode doWork(int side, EntityType type,
464 
465  MoFEMErrorCode calculateNormals();
466 };
467 
468 /** \brief calculate normals at Gauss points of triangle element
469  * \ingroup mofem_forces_and_sources
470  */
472 
481 
483  MatrixDouble &coords_at_gaussptf3, MatrixDouble &normals_at_gaussptf3,
484  MatrixDouble &tangent1_at_gaussptf3, MatrixDouble &tangent2_at_gaussptf3,
485  MatrixDouble &coords_at_gaussptf4, MatrixDouble &normals_at_gaussptf4,
486  MatrixDouble &tangent1_at_gaussptf4, MatrixDouble &tangent2_at_gaussptf4)
487  : cOords_at_GaussPtF3(coords_at_gaussptf3),
488  nOrmals_at_GaussPtF3(normals_at_gaussptf3),
489  tAngent1_at_GaussPtF3(tangent1_at_gaussptf3),
490  tAngent2_at_GaussPtF3(tangent2_at_gaussptf3),
491  cOords_at_GaussPtF4(coords_at_gaussptf4),
492  nOrmals_at_GaussPtF4(normals_at_gaussptf4),
493  tAngent1_at_GaussPtF4(tangent1_at_gaussptf4),
494  tAngent2_at_GaussPtF4(tangent2_at_gaussptf4) {}
495 
497  MoFEMErrorCode doWork(int side, EntityType type,
499 
500  MoFEMErrorCode calculateNormals();
501 };
502 
503 /** \brief transform Hdiv base fluxes from reference element to physical
504  * triangle \ingroup mofem_forces_and_sources
505  */
507 
508  const VectorDouble
509  *normalRawPtr; //< Normal of the element for linear geometry
510  const MatrixDouble *normalsAtGaussPtsRawPtr; //< Normals at integration points
511  //for nonlinear geometry
512 
513  /**
514  * @brief Shift in vector for linear geometry
515  *
516  * Normal can have size larger than three, for example normal for contact
517  * prism and flat prims element have six comonents, for top and bottom
518  * triangle of the prims.
519  *
520  */
522 
524  const MatrixDouble &normals_at_pts,
525  const int normal_shift = 0)
526  : normalRawPtr(&normal), normalsAtGaussPtsRawPtr(&normals_at_pts),
527  normalShift(normal_shift) {}
528 
530  const VectorDouble *normal_raw_ptr = nullptr,
531  const MatrixDouble *normals_at_pts_ptr = nullptr,
532  const int normal_shift = 0)
533  : normalRawPtr(normal_raw_ptr),
534  normalsAtGaussPtsRawPtr(normals_at_pts_ptr), normalShift(normal_shift) {
535  }
536 
537  MoFEMErrorCode doWork(int side, EntityType type,
539 };
540 
541 /** \brief transform Hcurl base fluxes from reference element to physical
542  * triangle \ingroup mofem_forces_and_sources
543  */
545 
552 
554  const MatrixDouble &normals_at_pts,
555  const VectorDouble &tangent0,
556  const MatrixDouble &tangent0_at_pts,
557  const VectorDouble &tangent1,
558  const MatrixDouble &tangent1_at_pts)
559  : nOrmal(normal), normalsAtGaussPts(normals_at_pts), tAngent0(tangent0),
560  tangent0AtGaussPt(tangent0_at_pts), tAngent1(tangent1),
561  tangent1AtGaussPt(tangent1_at_pts) {}
562 
563  MoFEMErrorCode doWork(int side, EntityType type,
565 };
566 
567 /** \brief Calculate tangent vector on edge form HO geometry approximation
568  * \ingroup mofem_forces_and_sources
569  */
571 
573 
574  OpGetHoTangentOnEdge(MatrixDouble &tangent) : tAngent(tangent) {}
575 
576  MoFEMErrorCode doWork(int side, EntityType type,
578 };
579 
580 /** \brief transform Hcurl base fluxes from reference element to physical edge
581  * \ingroup mofem_forces_and_sources
582  */
584 
587 
589  const MatrixDouble &tangent_at_pts)
590  : tAngent(tangent), tangentAtGaussPt(tangent_at_pts) {}
591 
592  MoFEMErrorCode doWork(int side, EntityType type,
594 };
595 
596 } // namespace MoFEM
597 
598 #endif //__DATAOPERATORS_HPP
599 
600 /**
601  * \defgroup mofem_forces_and_sources Forces and sources
602  **/
bool & doTris
\deprectaed
Calculate tangent vector on edge form HO geometry approximation.
OpSetHoContravariantPiolaTransform(VectorDouble &det_jac, MatrixDouble &jac)
brief Transform local reference derivatives of shape function to global derivatives
MatrixDouble diffNinvJac
OpSetHoInvJacH1(MatrixDouble &inv_ho_jac)
FieldApproximationBase & getBase()
Get approximation base.
OpSetInvJacHdivAndHcurl(MatrixDouble3by3 &inv_jac)
base operator to do operations at Gauss Pt. level
FTensor::Index< 'j', 3 > j
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
bool & doVertices
\deprectaed If false skip vertices
transform Hcurl base fluxes from reference element to physical triangle
FTensor::Index< 'i', 3 > i
FTensor::Index< 'i', 3 > i
void unSetSymm()
unset if operator is executed for non symmetric problem
bool & doPrisms
\deprectaed
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
OpGetDataAndGradient(MatrixDouble &data_at_gauss_pt, MatrixDouble &data_grad_at_gauss_pt)
OpSetHoInvJacHdivAndHcurl(MatrixDouble &inv_ho_jac)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
transform Hcurl base fluxes from reference element to physical edge
FTensor::Tensor1< double, 2 > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, 2 > &t_disp)
Definition: ContactOps.hpp:137
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
OpSetContravariantPiolaTransform(double &volume, MatrixDouble3by3 &jac)
void setSymm()
set if operator is executed taking in account symmetry
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
FTensor::Tensor2< double *, 3, 3 > tInvJac
FTensor::Tensor2< double *, 3, 3 > tInvJac
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
OpSetCovariantPiolaTransform(MatrixDouble3by3 &inv_jac)
FTensor::Tensor2< double *, 3, 3 > tInvJac
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
virtual MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpGetCoordsAndNormalsOnFace(MatrixDouble &coords_at_gausspt, MatrixDouble &normals_at_gausspt, MatrixDouble &tangent1_at_gausspt, MatrixDouble &tangent2_at_gausspt)
FTensor::Index< 'i', 3 > i
FTensor::Tensor2< double *, R, D > getGradAtGaussPtsTensor(MatrixDouble &data)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpGetHoTangentOnEdge(MatrixDouble &tangent)
int normalShift
Shift in vector for linear geometry.
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
FTensor::Index< 'i', 3 > i
static const char *const ApproximationBaseNames[]
Definition: definitions.h:164
Get field values and gradients at Gauss points.
FTensor::Index< 'j', 3 > j
apply contravariant (Piola) transfer to Hdiv space
FTensor::Index< 'k', 3 > k
OpSetContravariantPiolaTransformOnFace(const VectorDouble &normal, const MatrixDouble &normals_at_pts, const int normal_shift=0)
bool sYmm
If true assume that matrix is symmetric structure.
bool getSymm() const
Get if operator uses symmetry of DOFs or not.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
FTensor::Tensor1< double *, R > getValAtGaussPtsTensor(MatrixDouble &data)
FTensor::Index< 'k', 3 > k
transform Hdiv base fluxes from reference element to physical triangle
bool & doEdges
\deprectaed If false skip edges
OpSetCovariantPiolaTransformOnEdge(const VectorDouble &tangent, const MatrixDouble &tangent_at_pts)
calculate normals at Gauss points of triangle element
OpSetInvJacH1(MatrixDouble3by3 &inv_jac)
OpGetCoordsAndNormalsOnPrism(MatrixDouble &coords_at_gaussptf3, MatrixDouble &normals_at_gaussptf3, MatrixDouble &tangent1_at_gaussptf3, MatrixDouble &tangent2_at_gaussptf3, MatrixDouble &coords_at_gaussptf4, MatrixDouble &normals_at_gaussptf4, MatrixDouble &tangent1_at_gaussptf4, MatrixDouble &tangent2_at_gaussptf4)
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
Apply contravariant (Piola) transfer to Hdiv space for HO geometr.
bool & doQuads
\deprectaed
transform local reference derivatives of shape function to global derivatives if higher order geometr...
#define CHKERR
Inline error check.
Definition: definitions.h:604
FTensor::Tensor2< double *, 3, 3 > tJac
OpSetContravariantPiolaTransformOnFace(const VectorDouble *normal_raw_ptr=nullptr, const MatrixDouble *normals_at_pts_ptr=nullptr, const int normal_shift=0)
Transform local reference derivatives of shape function to global derivatives.
Data on single entity (This is passed as argument to DataOperator::doWork)
OpSetHoCovariantPiolaTransform(MatrixDouble &inv_jac)
FTensor::Index< 'j', 3 > j
transform local reference derivatives of shape function to global derivatives if higher order geometr...
apply covariant transfer to Hcurl space
Calculate normals at Gauss points of triangle element.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
FTensor::Index< 'j', 3 > j
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
const VectorDouble & getFieldData() const
get dofs values
OpSetCovariantPiolaTransformOnFace(const VectorDouble &normal, const MatrixDouble &normals_at_pts, const VectorDouble &tangent0, const MatrixDouble &tangent0_at_pts, const VectorDouble &tangent1, const MatrixDouble &tangent1_at_pts)
MoFEMErrorCode calculateValAndGrad(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Calculate gradient and values at integration points.
bool & doTets
\deprectaed
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions