v0.8.23
VolumeElementForcesAndSourcesCore.hpp
Go to the documentation of this file.
1 /** \file VolumeElementForcesAndSourcesCore.hpp
2  \brief Volume element.
3 
4  Those element are inherited by user to implement specific implementation of
5  particular problem.
6 
7 */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef __VOLUMEELEMENTFORCESANDSOURCESCORE_HPP__
24 #define __VOLUMEELEMENTFORCESANDSOURCESCORE_HPP__
25 
26 using namespace boost::numeric;
27 
28 namespace MoFEM {
29 
30 /** \brief Volume finite element
31  \ingroup mofem_forces_and_sources_volume_element
32 
33  User is implementing own operator at Gauss point level, by class
34  derived from VolumeElementForcesAndSourcesCore::UserDataOperator. Arbitrary
35  number of operator can be added by pushing objects to OpPtrVector
36 
37  */
39 
43 
48 
54 
56  opHOatGaussPoints; ///< higher order geometry data at Gauss pts
61 
63  const EntityType type = MBTET);
65 
66  double vOlume;
67 
68  int num_nodes;
72 
74 
75  /** \brief default operator for TET element
76  * \ingroup mofem_forces_and_sources_volume_element
77  */
79 
82 
83  UserDataOperator(const std::string &field_name, const char type)
84  : ForcesAndSourcesCore::UserDataOperator(field_name, type) {}
85 
86  UserDataOperator(const std::string &row_field_name,
87  const std::string &col_field_name, const char type,
88  const bool symm = true)
89  : ForcesAndSourcesCore::UserDataOperator(row_field_name, col_field_name,
90  type, symm) {}
91 
92  /** \brief get element number of nodes
93  */
94  inline int getNumNodes() {
95  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)->num_nodes;
96  }
97 
98  /** \brief get element connectivity
99  */
100  inline const EntityHandle *getConn() {
101  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)->conn;
102  }
103 
104  /** \brief element volume (linear geometry)
105  */
106  inline double getVolume() {
107  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)->vOlume;
108  }
109 
110  /**
111  * \brief get element Jacobian
112  */
114  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)->tJac;
115  }
116 
117  /**
118  * \brief get element inverse Jacobian
119  */
121  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)->tInvJac;
122  }
123 
124  /**
125  * \brief get measure of element
126  * @return area of face
127  */
128  inline double getMeasure() { return getVolume(); }
129 
130  /** \brief nodal coordinates
131  */
133  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)->coords;
134  }
135 
136  /** \brief Gauss points and weight, matrix (nb. of points x 3)
137 
138  Column 0-2 integration points coordinate x, y and z, respectively. At rows
139  are integration points.
140 
141  */
143  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)
144  ->coordsAtGaussPts;
145  }
146 
147  /** \brief coordinate at Gauss points (if hierarchical approximation of
148  * element geometry)
149  */
151  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)
152  ->hoCoordsAtGaussPts;
153  }
154 
156  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)
157  ->hoGaussPtsJac;
158  }
159 
161  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)
162  ->hoGaussPtsInvJac;
163  }
164 
166  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE)
167  ->hoGaussPtsDetJac;
168  }
169 
170  inline auto getFTenosr0HoMeasure() {
172  &*getHoGaussPtsDetJac().data().begin());
173  }
174 
175  /**
176  * \brief Get coordinates at integration points assuming linear geometry
177  *
178  * \code
179  * auto t_coords = getFTensor1CoordsAtGaussPts();
180  * for(int gg = 0;gg!=nb_int_ptrs;gg++) {
181  * // do something
182  * ++t_coords;
183  * }
184  * \endcode
185  *
186  */
189  &getCoordsAtGaussPts()(0, 0), &getCoordsAtGaussPts()(0, 1),
190  &getCoordsAtGaussPts()(0, 2));
191  }
192 
193  /**
194  * \brief Get coordinates at integration points taking geometry from field
195  *
196  * This is HO geometry given by arbitrary order polynomial
197  * \code
198  * auto t_coords = getFTensor1HoCoordsAtGaussPts();
199  * for(int gg = 0;gg!=nb_int_ptrs;gg++) {
200  * // do something
201  * ++t_coords;
202  * }
203  * \endcode
204  *
205  */
207  double *ptr = &*getHoCoordsAtGaussPts().data().begin();
208  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, ptr + 1,
209  ptr + 2);
210  }
211 
212  inline auto getFTensor2HoGaussPtsJac() {
213  double *ptr = &*getHoGaussPtsJac().data().begin();
215  ptr, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
216  ptr + 8);
217  }
218 
220  double *ptr = &*getHoGaussPtsInvJac().data().begin();
222  ptr, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
223  ptr + 8);
224  }
225 
226  /** \brief return pointer to Generic Volume Finite Element object
227  */
229  return static_cast<VolumeElementForcesAndSourcesCore *>(ptrFE);
230  }
231 
232  /**
233  * \brief Get divergence of base functions at integration point
234  *
235  * Works only for H-div space.
236  *
237  * How to use it:
238  * \code
239  * VectorDouble div_vec(data.getFieldData().size());
240  * for(int gg = 0;gg<nb_gauss_pts;gg++) {
241  * CHKERR getDivergenceOfHDivBaseFunctions(side,type,data,gg,div_vec);
242  * // do something with vec_div
243  * }
244  * \endcode
245  * where vec_div has size of nb. of dofs, at each element represents
246  * divergence of base functions.
247  *
248  * @param side side (local) number of entity on element
249  * @param type type of entity
250  * @param data data structure
251  * @param gg gauss pts
252  * @param div divergence vector, size of vector is equal to number of base
253  * functions
254  * @return error code
255  */
257  getDivergenceOfHDivBaseFunctions(int side, EntityType type,
259  int gg, VectorDouble &div);
260 
261  /**
262  * \brief Get curl of base functions at integration point
263  *
264  * \f[
265  * \nabla \times \mathbf{\phi}
266  * \f]
267  *
268  * Works only for H-curl and H-div space.
269  *
270  * How to use it:
271  * \code
272  * MatrixDouble curl_mat(data.getFieldData().size(),3);
273  * for(int gg = 0;gg<nb_gauss_pts;gg++) {
274  * CHKERR getCurlOfHCurlBaseFunctions(side,type,data,gg,curl_mat);
275  * FTensor::Tensor1<FTensor::PackPtr<double*, 3>, 3> t_curl(
276  * &curl_mat(0,0), &curl_mat(0,1), &curl_mat(0,2)
277  * );
278  * // do something with curl
279  * }
280  * \endcode
281  * where curl_mat is matrix which number of rows is equal to nb. of base
282  * functions. Number of columns is 3, since we work in 3d here. Rows
283  * represents curl of base functions.
284  *
285  * @param side side (local) number of entity on element
286  * @param type type of entity
287  * @param data data structure
288  * @param gg gauss pts
289  * @param curl curl matrix, nb. of rows is equal to number of base
290  * functions, columns are curl of base vector
291  * @return error code
292  */
294  getCurlOfHCurlBaseFunctions(int side, EntityType type,
295  DataForcesAndSourcesCore::EntData &data, int gg,
296  MatrixDouble &curl);
297 
298  /// \deprecated use getFTensor1CoordsAtGaussPts
300  return getFTensor1CoordsAtGaussPts();
301  }
302 
303  /// \deprecated use getFTensor1HoCoordsAtGaussPts
305  return getFTensor1HoCoordsAtGaussPts();
306  }
307 
308  };
309 
310  unsigned int nbGaussPts; ///< Number of integration points
311 
312  // Note that functions below could be overloaded by user to change default
313  // behavior of the element.
314 
315  /**
316  * \brief Set integration points
317  * @return Error code
318  */
319  virtual MoFEMErrorCode setIntegrationPts();
320 
321  /// \deprecated function with spelling mistake, use setIntegrationPts
323  return setIntegrationPts();
324  }
325 
326  /**
327  * \brief Calculate element volume and Jacobian
328  *
329  * Note that at that point is assumed that geometry is exclusively defined by
330  * corner nodes.
331  *
332  * @return Error code
333  */
334  virtual MoFEMErrorCode calculateVolumeAndJacobian();
335 
336  /**
337  * \brief Calculate coordinate at integration points
338  * @return Error code
339  */
340  virtual MoFEMErrorCode calculateCoordinatesAtGaussPts();
341 
342  /**
343  * \brief Determine approximation space and order of base functions
344  * @return Error code
345  */
346  virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement();
347 
348  /**
349  * \brief Transform base functions based on geometric element Jacobian.
350  *
351  * This function apply transformation to base functions and its derivatives.
352  * For example when base functions for H-div are present the
353  * Piola-Transformarion is applied to base functions and their derivatives.
354  *
355  * @return Error code
356  */
357  virtual MoFEMErrorCode transformBaseFunctions();
358 
359  /** \brief Calculate Jacobian for HO geometry
360  *
361  * MoFEM use hierarchical approximate base to describe geometry of the body.
362  * This function transform derivatives of base functions when HO geometry is
363  * set and calculate Jacobian, inverse of Jacobian and determinant of
364  * transformation.
365  *
366  */
367  virtual MoFEMErrorCode calculateHoJacobian();
368 
369  /**
370  * \brief Transform base functions based on ho-geometry element Jacobian.
371  *
372  * This function apply transformation to base functions and its derivatives.
373  * For example when base functions for H-div are present the
374  * Piola-Transformarion is applied to base functions and their derivatives.
375  *
376  * @return Error code
377  */
378  virtual MoFEMErrorCode transformHoBaseFunctions();
379 
380  MoFEMErrorCode operator()();
381 };
382 
384 
385 /**
386  * \brief Volume element used to integrate on skeleton
387  * \ingroup mofem_forces_and_sources_volume_element
388  */
391 
394  const EntityType type = MBTET)
395  : VolumeElementForcesAndSourcesCore(m_field, type), faceFEPtr(NULL) {}
397 
398  inline MoFEMErrorCode
401  faceFEPtr = const_cast<FaceElementForcesAndSourcesCore *>(face_fe_ptr);
403  }
404 
405  int getRule(int order) { return -1; };
406 
407  int faceSense; ///< Sense of face, could be 1 or -1
408  int faceSideNumber; ///< Face side number
409  int faceConnMap[3];
410  int tetConnMap[4];
412 
413  MoFEMErrorCode setGaussPts(int order);
414 
415  /** \brief default operator for TET element
416  * \ingroup mofem_forces_and_sources_volume_element
417  */
420 
421  UserDataOperator(const std::string &field_name, const char type)
423  type) {}
424 
425  UserDataOperator(const std::string &row_field_name,
426  const std::string &col_field_name, const char type)
428  row_field_name, col_field_name, type) {}
429 
430  /** \brief return pointer to Generic Volume Finite Element object
431  */
433  return static_cast<VolumeElementForcesAndSourcesCoreOnSide *>(ptrFE);
434  }
435 
437  return getVolumeFE()->faceFEPtr;
438  }
439 
440  /**
441  * \brief get face sense in respect to volume
442  * @return error code
443  */
444  inline int getFaceSense() const { return getVolumeFE()->faceSense; }
445 
446  /**
447  * \brief get face side number in respect to volume
448  * @return error code
449  */
450  inline int getFaceSideNumber() const {
451  return getVolumeFE()->faceSideNumber;
452  }
453 
454  inline bool getEdgeFace(const int ee) const {
455  const bool edges_on_faces[6][4] = {{true, false, false, true}, // e0
456  {false, true, false, true}, // e1
457  {false, false, true, true}, // e2
458  {true, false, true, false}, // e3
459  {true, true, false, false}, // e4
460  {false, true, true, false}};
461  return edges_on_faces[ee][getFaceSideNumber()];
462  }
463 
464  /**
465  * get face normal on side which is this element
466  * @return face normal
467  */
468  VectorDouble &getNormal();
469 
470  /** \brief get normal as tensor
471  */
472  inline auto getFTensor1Normal() {
473  double *ptr = &*getNormal().data().begin();
474  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
475  }
476 
477  /** \brief if higher order geometry return normals at Gauss pts.
478 
479  Note: returned matrix has size 0 in rows and columns if no HO approximation
480  of geometry is available.
481 
482  */
483  MatrixDouble &getNormalsAtGaussPts();
484 
485  /** \brief if higher order geometry return normals at Gauss pts.
486  *
487  * \param gg gauss point number
488  */
489  ublas::matrix_row<MatrixDouble> getNormalsAtGaussPts(const int gg);
490 
491  /// \deprecated use getNormalsAtGaussPts
492  DEPRECATED inline auto getNormalsAtGaussPt(const int gg) {
493  return getNormalsAtGaussPts(gg);
494  }
495 
496  /** \brief get normal at integration points
497 
498  Example:
499  \code
500  double nrm2;
501  FTensor::Index<'i',3> i;
502  auto t_normal = getFTensor1NormalsAtGaussPts();
503  for(int gg = gg!=data.getN().size1();gg++) {
504  nrm2 = sqrt(t_normal(i)*t_normal(i));
505  ++t_normal;
506  }
507  \endcode
508 
509  */
511  double *ptr = &*getNormalsAtGaussPts().data().begin();
512  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
513  &ptr[2]);
514  }
515 
516  /// \deprecated use getFTensor1NormalsAtGaussPts
518  return getFTensor1NormalsAtGaussPts();
519  }
520 
521  /** \brief get face coordinates at Gauss pts.
522 
523  \note Coordinates should be the same what function getCoordsAtGaussPts
524  on tets is returning. If both coordinates are different it is error, or you
525  do something very unusual.
526 
527  */
528  MatrixDouble &getFaceCoordsAtGaussPts();
529  };
530 
531 };
532 
533 /// \deprecated Use VolumeElementForcesAndSourcesCore
534 DEPRECATED typedef VolumeElementForcesAndSourcesCore
536 
537 } // namespace MoFEM
538 
539 #endif //__VOLUMEELEMENTFORCESANDSOURCESCORE_HPP__
540 
541 /**
542  * \defgroup mofem_forces_and_sources_volume_element Volume Element
543  * \brief Implementation of general volume element.
544  *
545  * \ingroup mofem_forces_and_sources
546  **/
brief Transform local reference derivatives of shape function to global derivatives
OpGetDataAndGradient< 3, 3 > opHOatGaussPoints
higher order geometry data at Gauss pts
auto getFTensor1HoCoordsAtGaussPts()
Get coordinates at integration points taking geometry from field.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:500
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
UserDataOperator(const std::string &row_field_name, const std::string &col_field_name, const char type, const bool symm=true)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
Volume finite elementUser is implementing own operator at Gauss point level, by class derived from Vo...
VolumeElementForcesAndSourcesCoreOnSide(Interface &m_field, const EntityType type=MBTET)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpSetHoContravariantPiolaTransform opHoContravariantTransform
const VolumeElementForcesAndSourcesCoreOnSide * getVolumeFE() const
return pointer to Generic Volume Finite Element object
apply contravariant (Piola) transfer to Hdiv space
unsigned int nbGaussPts
Number of integration points.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
DEPRECATED typedef VolumeElementForcesAndSourcesCore VolumeElementForcesAndSurcesCore
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:98
Apply contravariant (Piola) transfer to Hdiv space for HO geometr.
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
FieldSpace
approximation spaces
Definition: definitions.h:167
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Transform local reference derivatives of shape function to global derivatives.
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
UserDataOperator(const std::string &row_field_name, const std::string &col_field_name, const char type)
FTensor::Tensor2< double *, 3, 3 > & getInvJac()
get element inverse Jacobian
const VolumeElementForcesAndSourcesCore * getVolumeFE()
return pointer to Generic Volume Finite Element object
#define DEPRECATED
Definition: definitions.h:29
MatrixDouble & getHoCoordsAtGaussPts()
coordinate at Gauss points (if hierarchical approximation of element geometry)
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
Data on single entity (This is passed as argument to DataOperator::doWork)
OpSetContravariantPiolaTransform opContravariantPiolaTransform
structure to get information form mofem into DataForcesAndSourcesCore
transform local reference derivatives of shape function to global derivatives if higher order geometr...
apply covariant transfer to Hcurl space
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:76
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
const int order
approximation order
MoFEMErrorCode setFaceFEPtr(const FaceElementForcesAndSourcesCore *face_fe_ptr)