v0.13.0
HODataOperators.hpp
Go to the documentation of this file.
1 /** \file HODataOperators.hpp
2  * \brief Operators managing HO geometry
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 __HO_DATA_OPERATORS_HPP__
21 #define __HO_DATA_OPERATORS_HPP__
22 
23 namespace MoFEM {
24 
25 /**
26  * @brief Calculate jacobian on Hex or other volume which is not simplex.
27  *
28  * Using not a field data, but geometrical coordinates of element.
29  *
30  * \note Use OpCalculateVectorFieldGradient to calculate Jacobian from field
31  * data.
32  *
33  */
36 
37  OpCalculateHOJacVolume(boost::shared_ptr<MatrixDouble> jac_ptr);
38 
39  MoFEMErrorCode doWork(int side, EntityType type,
41 
42 private:
43  boost::shared_ptr<MatrixDouble> jacPtr;
44 };
45 
46 /**
47  * @brief Calculate HO coordinates at gauss points
48  *
49  */
51 
52  OpCalculateHOCoords(const std::string field_name)
54 
55  MoFEMErrorCode doWork(int side, EntityType type,
57 };
58 
59 /**
60  * \brief Set inverse jacobian to base functions
61  *
62  */
65 
67  boost::shared_ptr<MatrixDouble> inv_jac_ptr)
68  : ForcesAndSourcesCore::UserDataOperator(space), invJacPtr(inv_jac_ptr) {
69 
70  if (!inv_jac_ptr)
71  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "invJacPtr not allocated");
72 
73  if (space == L2) {
74  doVertices = false;
75  }
76  }
77 
78  MoFEMErrorCode doWork(int side, EntityType type,
80 
81 private:
82  boost::shared_ptr<MatrixDouble> invJacPtr;
84 };
85 
86 /**
87  * \brief transform local reference derivatives of shape function to global
88  derivatives if higher order geometry is given
89  *
90 
91  * \ingroup mofem_forces_and_sources
92 */
94 
96  boost::shared_ptr<MatrixDouble> inv_jac_ptr)
97  : ForcesAndSourcesCore::UserDataOperator(space), invJacPtr(inv_jac_ptr) {
98 
99  if (!invJacPtr)
101  "Pointer for invJacPtr not allocated");
102 
103  doVertices = false;
104  if (space == HDIV)
105  doEdges = false;
106  }
107 
108  MoFEMErrorCode doWork(int side, EntityType type,
110 
111 private:
112  boost::shared_ptr<MatrixDouble> invJacPtr;
115 };
116 
117 /**
118  * @brief Modify integration weights on face to take in account higher-order
119  * geometry
120  * @ingroup mofem_forces_and_sources_tri_element
121  *
122  */
127  MoFEMErrorCode doWork(int side, EntityType type,
129 };
130 
131 /**
132  * \brief Set inverse jacobian to base functions
133  *
134  */
136 
137  OpSetHOWeights(boost::shared_ptr<VectorDouble> det_ptr)
139  detPtr(det_ptr) {
140  if (!detPtr)
142  "Pointer for detPtr not allocated");
143  }
144 
145  MoFEMErrorCode doWork(int side, EntityType type,
147 
148 private:
149  boost::shared_ptr<VectorDouble> detPtr;
150 };
151 
152 /** \brief Apply contravariant (Piola) transfer to Hdiv space for HO geometr
153 
154 * \ingroup mofem_forces_and_sources
155 */
158 
160  boost::shared_ptr<VectorDouble> det_ptr,
161  boost::shared_ptr<MatrixDouble> jac_ptr)
163  jacPtr(jac_ptr) {
164  doVertices = false;
165  if (space == HDIV)
166  doEdges = false;
167  }
168 
169  MoFEMErrorCode doWork(int side, EntityType type,
171 
172 private:
173  boost::shared_ptr<VectorDouble> detPtr;
174  boost::shared_ptr<MatrixDouble> jacPtr;
175 
178 };
179 
180 /** \brief Apply covariant (Piola) transfer to Hcurl space for HO geometry
181  */
184 
186  boost::shared_ptr<MatrixDouble> jac_inv_ptr)
188  jacInvPtr(jac_inv_ptr) {
189  doVertices = false;
190  if (space == HDIV)
191  doEdges = false;
192  if (!jacInvPtr)
194  "Pointer for jacPtr not allocated");
195  }
196 
197  MoFEMErrorCode doWork(int side, EntityType type,
199 
200 private:
201  boost::shared_ptr<MatrixDouble> jacInvPtr;
202 
205 };
206 
207 /** \brief Calculate jacobian for face element
208 
209  It is assumed that face element is XY plane. Applied
210  only for 2d problems and 2d problems embedded in 3d space
211 
212  \note If you push operators for HO normal befor this operator, HO geometry is
213  taken into account when you calculate jacobian.
214 
215 */
216 template <int DIM> struct OpCalculateHOJacForFaceImpl;
217 
218 template <>
221 
222  OpCalculateHOJacForFaceImpl(boost::shared_ptr<MatrixDouble> jac_ptr);
223 
224  MoFEMErrorCode doWork(int side, EntityType type,
226 
227 protected:
228  boost::shared_ptr<MatrixDouble> jacPtr;
229 };
230 
231 template <>
233 
235 
236  MoFEMErrorCode doWork(int side, EntityType type,
238 };
239 
242 
243 /** \brief Calculate normals at Gauss points of triangle element
244  * \ingroup mofem_forces_and_source
245  */
248 
249  OpGetHONormalsOnFace(std::string field_name);
250 
251  MoFEMErrorCode doWork(int side, EntityType type,
253 };
254 
255 /** \brief transform Hdiv base fluxes from reference element to physical
256  * triangle \ingroup mofem_forces_and_sources
257  *
258  * \note Matrix which keeps normal is assumed to have three columns, and number
259  * of rows should be equal to number of integration points.
260  *
261  */
264 
266  const FieldSpace space,
267  boost::shared_ptr<MatrixDouble> normals_at_gauss_pts = nullptr)
269  normalsAtGaussPts(normals_at_gauss_pts) {}
270 
271  MoFEMErrorCode doWork(int side, EntityType type,
273 
274 private:
275  boost::shared_ptr<MatrixDouble> normalsAtGaussPts;
276 };
277 
278 /** \brief transform Hcurl base fluxes from reference element to physical edge
279  * \ingroup mofem_forces_and_sources
280  */
283 
285  const FieldSpace space = HCURL,
286  boost::shared_ptr<MatrixDouble> tangent_at_pts = nullptr)
288  tangentAtGaussPts(tangent_at_pts) {}
289 
290  MoFEMErrorCode doWork(int side, EntityType type,
292 
293 private:
294  boost::shared_ptr<MatrixDouble> tangentAtGaussPts;
295 };
296 
297 /** \brief transform Hcurl base fluxes from reference element to physical
298  * triangle \ingroup mofem_forces_and_sources
299  */
302 
304  const FieldSpace space,
305  boost::shared_ptr<MatrixDouble> normals_at_pts = nullptr,
306  boost::shared_ptr<MatrixDouble> tangent1_at_pts = nullptr,
307  boost::shared_ptr<MatrixDouble> tangent2_at_pts = nullptr)
309  normalsAtPts(normals_at_pts), tangent1AtPts(tangent1_at_pts),
310  tangent2AtPts(tangent2_at_pts) {
311  if (normals_at_pts || tangent1_at_pts || tangent2_at_pts)
312  if (normals_at_pts && tangent1_at_pts && tangent2_at_pts)
314  "All elements in constructor have to set pointer");
315  }
316 
317  MoFEMErrorCode doWork(int side, EntityType type,
319 
320 private:
321  boost::shared_ptr<MatrixDouble> normalsAtPts;
322  boost::shared_ptr<MatrixDouble> tangent1AtPts;
323  boost::shared_ptr<MatrixDouble> tangent2AtPts;
324 
327 };
328 
329 /** \brief Calculate tangent vector on edge form HO geometry approximation
330  * \ingroup mofem_forces_and_sources
331  */
334 
336  std::string field_name,
337  boost::shared_ptr<MatrixDouble> tangents_at_pts = nullptr)
339  OPROW),
340  tangentsAtPts(tangents_at_pts) {}
341 
342  MoFEMErrorCode doWork(int side, EntityType type,
344 
345 private:
346  boost::shared_ptr<MatrixDouble> tangentsAtPts;
347 };
348 
349 /**
350  * @brief Scale base functions by inverses of measure of element
351  *
352  * @tparam OP
353  */
354 template <typename OP> struct OpScaleBaseBySpaceInverseOfMeasure : public OP {
355 
357  boost::shared_ptr<VectorDouble> det_jac_ptr, const FieldSpace space = L2)
358  : OP(space), fieldSpace(space), detJacPtr(det_jac_ptr) {}
359 
360  MoFEMErrorCode doWork(int side, EntityType type,
363 
364  if (!detJacPtr) {
365  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "detJacPtr not set");
366  }
367 
368  auto scale = [&]() {
369  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
370  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
371  auto &base_fun = data.getN(base);
372  auto &diff_base_fun = data.getDiffN(base);
373  if (detJacPtr) {
374 
375  auto &det_vec = *detJacPtr;
376  const auto nb_base_fun = base_fun.size2();
377  const auto nb_int_pts = base_fun.size1();
378 
379  if (nb_int_pts != det_vec.size())
380  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
381  "Number of integration pts in detJacPtr does not mush "
382  "number of integration pts in base function");
383 
384  auto get_row = [](auto &m, auto gg) {
385  return ublas::matrix_row<decltype(m)>(m, gg);
386  };
387 
388  for (auto gg = 0; gg != nb_int_pts; ++gg)
389  get_row(base_fun) /= det_vec[gg];
390 
391  if (diff_base_fun.size1() == nb_int_pts) {
392  for (auto gg = 0; gg != nb_int_pts; ++gg)
393  get_row(diff_base_fun) /= det_vec[gg];
394  }
395  }
396  }
397  };
398 
399  if (this->getFEDim() == 3) {
400  switch (fieldSpace) {
401  case H1:
402  scale();
403  break;
404  case HCURL:
405  if (type >= MBEDGE)
406  scale();
407  break;
408  case HDIV:
409  if (type >= MBTRI)
410  scale();
411  break;
412  case L2:
413  if (type >= MBTET)
414  scale();
415  break;
416  default:
417  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "impossible case");
418  }
419  } else if (this->getFEDim() == 2) {
420  switch (fieldSpace) {
421  case H1:
422  scale();
423  break;
424  case HCURL:
425  if (type >= MBEDGE)
426  scale();
427  break;
428  case HDIV:
429  case L2:
430  if (type >= MBTRI)
431  scale();
432  break;
433  default:
434  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "impossible case");
435  }
436  } else if (this->getFEDim() == 1) {
437  switch (fieldSpace) {
438  case H1:
439  scale();
440  break;
441  case HCURL:
442  case L2:
443  if (type >= MBEDGE)
444  scale();
445  break;
446  default:
447  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "impossible case");
448  }
449  } else {
450  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "impossible case");
451  }
452 
454  }
455 
456 private:
458  boost::shared_ptr<VectorDouble> detJacPtr;
459 };
460 
461 template <typename E>
462 MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl,
463  bool hdiv, bool l2) {
465  auto material_grad_mat = boost::make_shared<MatrixDouble>();
466  auto material_det_vec = boost::make_shared<VectorDouble>();
467  auto material_inv_grad_mat = boost::make_shared<MatrixDouble>();
468  e.getOpPtrVector().push_back(
469  new OpCalculateVectorFieldGradient<3, 3>(field, material_grad_mat));
470  e.getOpPtrVector().push_back(new OpInvertMatrix<3>(
471  material_grad_mat, material_det_vec, material_inv_grad_mat));
472  e.getOpPtrVector().push_back(new OpSetHOWeights(material_det_vec));
473  if (h1)
474  e.getOpPtrVector().push_back(
475  new OpSetHOInvJacToScalarBases(H1, material_inv_grad_mat));
476  if (l2)
477  e.getOpPtrVector().push_back(
478  new OpSetHOInvJacToScalarBases(L2, material_inv_grad_mat));
479  if (hdiv) {
480  e.getOpPtrVector().push_back(new OpSetHOContravariantPiolaTransform(
481  HDIV, material_det_vec, material_grad_mat));
482  e.getOpPtrVector().push_back(
483  new OpSetHOInvJacVectorBase(HDIV, material_inv_grad_mat));
484  }
485  if (hcurl) {
486  e.getOpPtrVector().push_back(
487  new OpSetHOCovariantPiolaTransform(HCURL, material_inv_grad_mat));
488  e.getOpPtrVector().push_back(
489  new OpSetHOInvJacVectorBase(HCURL, material_inv_grad_mat));
490  }
492 }
493 
494 template <typename E>
495 MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl,
496  bool hdiv) {
498  e.meshPositionsFieldName = "none";
499  e.getOpPtrVector().push_back(new OpGetHONormalsOnFace(field));
500  e.getOpPtrVector().push_back(new OpCalculateHOCoords(field));
501  if (hcurl) {
502  e.getOpPtrVector().push_back(
504  }
505  if (hdiv) {
506  e.getOpPtrVector().push_back(
508  }
510 }
511 
512 }; // namespace MoFEM
513 
514 #endif
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ USER_BASE
user implemented approximation base
Definition: definitions.h:81
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:608
FieldSpace
approximation spaces
Definition: definitions.h:95
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ NOSPACE
Definition: definitions.h:96
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
double scale
Definition: plastic.cpp:103
FTensor::Index< 'm', 3 > m
bool & doVertices
\deprectaed If false skip vertices
bool & doEdges
\deprectaed If false skip edges
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
structure to get information form mofem into EntitiesFieldData
Calculate HO coordinates at gauss points.
OpCalculateHOCoords(const std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > jacPtr
Calculate jacobian for face element.
Calculate jacobian on Hex or other volume which is not simplex.
boost::shared_ptr< MatrixDouble > jacPtr
OpCalculateHOJacVolume(boost::shared_ptr< MatrixDouble > jac_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::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.
Calculate normals at Gauss points of triangle element.
OpGetHONormalsOnFace(std::string field_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate tangent vector on edge form HO geometry approximation.
boost::shared_ptr< MatrixDouble > tangentsAtPts
OpGetHOTangentsOnEdge(std::string field_name, boost::shared_ptr< MatrixDouble > tangents_at_pts=nullptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
transform Hcurl base fluxes from reference element to physical edge
boost::shared_ptr< MatrixDouble > tangentAtGaussPts
OpHOSetContravariantPiolaTransformOnEdge3D(const FieldSpace space=HCURL, boost::shared_ptr< MatrixDouble > tangent_at_pts=nullptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
transform Hdiv base fluxes from reference element to physical triangle
boost::shared_ptr< MatrixDouble > normalsAtGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpHOSetContravariantPiolaTransformOnFace3D(const FieldSpace space, boost::shared_ptr< MatrixDouble > normals_at_gauss_pts=nullptr)
transform Hcurl base fluxes from reference element to physical triangle
OpHOSetCovariantPiolaTransformOnFace3D(const FieldSpace space, boost::shared_ptr< MatrixDouble > normals_at_pts=nullptr, boost::shared_ptr< MatrixDouble > tangent1_at_pts=nullptr, boost::shared_ptr< MatrixDouble > tangent2_at_pts=nullptr)
boost::shared_ptr< MatrixDouble > normalsAtPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > tangent2AtPts
boost::shared_ptr< MatrixDouble > tangent1AtPts
Scale base functions by inverses of measure of element.
OpScaleBaseBySpaceInverseOfMeasure(boost::shared_ptr< VectorDouble > det_jac_ptr, const FieldSpace space=L2)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VectorDouble > detJacPtr
Apply contravariant (Piola) transfer to Hdiv space for HO geometr.
boost::shared_ptr< VectorDouble > detPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > jacPtr
OpSetHOContravariantPiolaTransform(const FieldSpace space, boost::shared_ptr< VectorDouble > det_ptr, boost::shared_ptr< MatrixDouble > jac_ptr)
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
OpSetHOCovariantPiolaTransform(const FieldSpace space, boost::shared_ptr< MatrixDouble > jac_inv_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > jacInvPtr
Set inverse jacobian to base functions.
OpSetHOInvJacToScalarBases(const FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > invJacPtr
transform local reference derivatives of shape function to global derivatives if higher order geometr...
boost::shared_ptr< MatrixDouble > invJacPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetHOInvJacVectorBase(const FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Set inverse jacobian to base functions.
OpSetHOWeights(boost::shared_ptr< VectorDouble > det_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > detPtr
Modify integration weights on face to take in account higher-order geometry.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.