v0.10.0
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 base
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 
41 
42  /** \brief default operator for TET element
43  * \ingroup mofem_forces_and_sources_volume_element
44  */
46 
48 
49  /** \brief get element number of nodes
50  */
51  inline int getNumNodes();
52 
53  /** \brief get element connectivity
54  */
55  inline const EntityHandle *getConn();
56 
57  /** \brief element volume (linear geometry)
58  */
59  inline double getVolume() const;
60 
61  /** \brief element volume (linear geometry)
62  */
63  inline double &getVolume();
64 
65  /**
66  * \brief get element Jacobian
67  */
68  inline FTensor::Tensor2<double *, 3, 3> &getJac();
69 
70  /**
71  * \brief get element inverse Jacobian
72  */
73  inline FTensor::Tensor2<double *, 3, 3> &getInvJac();
74 
75  /**
76  * \brief get measure of element
77  * @return volume
78  */
79  inline double getMeasure() const;
80 
81  /**
82  * \brief get measure of element
83  * @return volume
84  */
85  inline double &getMeasure();
86 
87  /** \brief nodal coordinates
88  */
89  inline VectorDouble &getCoords();
90 
91  /** \brief Gauss points and weight, matrix (nb. of points x 3)
92 
93  Column 0-2 integration points coordinate x, y and z, respectively. At rows
94  are integration points.
95 
96  */
97  inline MatrixDouble &getCoordsAtGaussPts();
98 
99  /** \brief coordinate at Gauss points (if hierarchical approximation of
100  * element geometry)
101  */
102  inline MatrixDouble &getHoCoordsAtGaussPts();
103 
104  inline MatrixDouble &getHoGaussPtsJac();
105 
106  inline MatrixDouble &getHoGaussPtsInvJac();
107 
108  inline VectorDouble &getHoGaussPtsDetJac();
109 
110  inline auto getFTenosr0HoMeasure();
111 
112  /**
113  * \brief Get coordinates at integration points assuming linear geometry
114  *
115  * \code
116  * auto t_coords = getFTensor1CoordsAtGaussPts();
117  * for(int gg = 0;gg!=nb_int_ptrs;gg++) {
118  * // do something
119  * ++t_coords;
120  * }
121  * \endcode
122  *
123  */
124  inline auto getFTensor1CoordsAtGaussPts();
125 
126  /**
127  * \brief Get coordinates at integration points taking geometry from field
128  *
129  * This is HO geometry given by arbitrary order polynomial
130  * \code
131  * auto t_coords = getFTensor1HoCoordsAtGaussPts();
132  * for(int gg = 0;gg!=nb_int_ptrs;gg++) {
133  * // do something
134  * ++t_coords;
135  * }
136  * \endcode
137  *
138  */
139  inline auto getFTensor1HoCoordsAtGaussPts();
140 
141  inline auto getFTensor2HoGaussPtsJac();
142 
143  inline auto getFTensor2HoGaussPtsInvJac();
144 
145  /** \brief return pointer to Generic Volume Finite Element object
146  */
147  inline VolumeElementForcesAndSourcesCoreBase *getVolumeFE() const;
148 
149  /**
150  * \brief Get divergence of base functions at integration point
151  *
152  * Works only for H-div space.
153  *
154  * How to use it:
155  * \code
156  * VectorDouble div_vec(data.getFieldData().size());
157  * for(int gg = 0;gg<nb_gauss_pts;gg++) {
158  * CHKERR getDivergenceOfHDivBaseFunctions(side,type,data,gg,div_vec);
159  * // do something with vec_div
160  * }
161  * \endcode
162  * where vec_div has size of nb. of dofs, at each element represents
163  * divergence of base functions.
164  *
165  * @param side side (local) number of entity on element
166  * @param type type of entity
167  * @param data data structure
168  * @param gg gauss pts
169  * @param div divergence vector, size of vector is equal to number of base
170  * functions
171  * @return error code
172  */
174  getDivergenceOfHDivBaseFunctions(int side, EntityType type,
176  int gg, VectorDouble &div);
177 
178  /**
179  * \brief Get curl of base functions at integration point
180  *
181  * \f[
182  * \nabla \times \mathbf{\phi}
183  * \f]
184  *
185  * Works only for H-curl and H-div space.
186  *
187  * How to use it:
188  * \code
189  * MatrixDouble curl_mat(data.getFieldData().size(),3);
190  * for(int gg = 0;gg<nb_gauss_pts;gg++) {
191  * CHKERR getCurlOfHCurlBaseFunctions(side,type,data,gg,curl_mat);
192  * FTensor::Tensor1<FTensor::PackPtr<double*, 3>, 3> t_curl(
193  * &curl_mat(0,0), &curl_mat(0,1), &curl_mat(0,2)
194  * );
195  * // do something with curl
196  * }
197  * \endcode
198  * where curl_mat is matrix which number of rows is equal to nb. of base
199  * functions. Number of columns is 3, since we work in 3d here. Rows
200  * represents curl of base functions.
201  *
202  * @param side side (local) number of entity on element
203  * @param type type of entity
204  * @param data data structure
205  * @param gg gauss pts
206  * @param curl curl matrix, nb. of rows is equal to number of base
207  * functions, columns are curl of base vector
208  * @return error code
209  */
211  getCurlOfHCurlBaseFunctions(int side, EntityType type,
212  DataForcesAndSourcesCore::EntData &data, int gg,
213  MatrixDouble &curl);
214 
215  protected:
216 
217  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
218 
219  };
220 
221  enum Switches {
222  NO_HO_GEOMETRY = 1 << 0 | 1 << 2,
223  NO_TRANSFORM = 1 << 1 | 1 << 2,
224  NO_HO_TRANSFORM = 1 << 2
225  };
226 
227  template <int SWITCH> MoFEMErrorCode OpSwitch();
228 
229 protected:
230 
232  const EntityType type = MBTET);
233 
234  // Note that functions below could be overloaded by user to change default
235  // behavior of the element.
236 
237  /**
238  * \brief Set integration points
239  * @return Error code
240  */
241  virtual MoFEMErrorCode setIntegrationPts();
242 
243  /**
244  * \brief Calculate element volume and Jacobian
245  *
246  * Note that at that point is assumed that geometry is exclusively defined by
247  * corner nodes.
248  *
249  * @return Error code
250  */
251  virtual MoFEMErrorCode calculateVolumeAndJacobian();
252 
253  /**
254  * \brief Calculate coordinate at integration points
255  * @return Error code
256  */
257  virtual MoFEMErrorCode calculateCoordinatesAtGaussPts();
258 
259  /**
260  * \brief Determine approximation space and order of base functions
261  * @return Error code
262  */
263  virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement();
264 
265  /**
266  * \brief Transform base functions based on geometric element Jacobian.
267  *
268  * This function apply transformation to base functions and its derivatives.
269  * For example when base functions for H-div are present the
270  * Piola-Transformarion is applied to base functions and their derivatives.
271  *
272  * @return Error code
273  */
274  virtual MoFEMErrorCode transformBaseFunctions();
275 
276  /** \brief Calculate Jacobian for HO geometry
277  *
278  * MoFEM use hierarchical approximate base to describe geometry of the body.
279  * This function transform derivatives of base functions when HO geometry is
280  * set and calculate Jacobian, inverse of Jacobian and determinant of
281  * transformation.
282  *
283  */
284  virtual MoFEMErrorCode calculateHoJacobian();
285 
286  /**
287  * \brief Transform base functions based on ho-geometry element Jacobian.
288  *
289  * This function apply transformation to base functions and its derivatives.
290  * For example when base functions for H-div are present the
291  * Piola-Transformarion is applied to base functions and their derivatives.
292  *
293  * @return Error code
294  */
295  virtual MoFEMErrorCode transformHoBaseFunctions();
296 
300 
305 
310 
312  opHOatGaussPoints; ///< higher order geometry data at Gauss pts
317 
318  double vOlume;
319 
321  const EntityHandle *conn;
324 
326 
327  friend class UserDataOperator;
328 };
329 
330 /**
331  * @brief Volume finite element with switches
332  *
333  * Using SWITCH to off functions
334  *
335  * @tparam SWITCH
336  */
337 template <int SWITCH>
340 
342  const EntityType type = MBTET): VolumeElementForcesAndSourcesCoreBase(m_field, MBTET) {}
345 
346  MoFEMErrorCode operator()();
347 };
348 
349 /** \brief Volume finite element default
350  \ingroup mofem_forces_and_sources_volume_element
351 
352  */
355 
356 template <int SWITCH>
357 MoFEMErrorCode VolumeElementForcesAndSourcesCoreBase::OpSwitch() {
359 
360  if (numeredEntFiniteElementPtr->getEntType() != MBTET)
362  CHKERR createDataOnElement();
363 
364  CHKERR calculateVolumeAndJacobian();
365  CHKERR getSpaceBaseAndOrderOnElement();
366  CHKERR setIntegrationPts();
367  if (gaussPts.size2() == 0)
369  CHKERR calculateCoordinatesAtGaussPts();
370  CHKERR calHierarchicalBaseFunctionsOnElement();
371  CHKERR calBernsteinBezierBaseFunctionsOnElement();
372 
373  if (!(NO_TRANSFORM & SWITCH))
374  CHKERR transformBaseFunctions();
375 
376  if (!(NO_HO_GEOMETRY & SWITCH))
377  CHKERR calculateHoJacobian();
378 
379  if (!(NO_HO_TRANSFORM & SWITCH))
380  CHKERR transformHoBaseFunctions();
381 
382  // Iterate over operators
383  CHKERR loopOverOperators();
384 
386 }
387 
388 template <int SWITCH>
390  return OpSwitch<SWITCH>();
391 }
392 
393 int VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getNumNodes() {
394  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->num_nodes;
395 }
396 
397 const EntityHandle *
398 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getConn() {
399  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->conn;
400 }
401 
402 double
403 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getVolume() const {
404  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->vOlume;
405 }
406 
407 double &VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getVolume() {
408  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->vOlume;
409 }
410 
412 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getJac() {
413  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->tJac;
414 }
415 
417 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getInvJac() {
418  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->tInvJac;
419 }
420 
421 double
422 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() const {
423  return getVolume();
424 }
425 
426 double &VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() {
427  return getVolume();
428 }
429 
430 VectorDouble &
431 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getCoords() {
432  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->coords;
433 }
434 
435 MatrixDouble &
436 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getCoordsAtGaussPts() {
437  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
438  ->coordsAtGaussPts;
439 }
440 
441 MatrixDouble &VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
442  getHoCoordsAtGaussPts() {
443  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
444  ->hoCoordsAtGaussPts;
445 }
446 
447 MatrixDouble &
448 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getHoGaussPtsJac() {
449  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
450  ->hoGaussPtsJac;
451 }
452 
453 MatrixDouble &
454 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getHoGaussPtsInvJac() {
455  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
456  ->hoGaussPtsInvJac;
457 }
458 
459 VectorDouble &
460 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getHoGaussPtsDetJac() {
461  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
462  ->hoGaussPtsDetJac;
463 }
464 
465 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
466  getFTenosr0HoMeasure() {
468  &*getHoGaussPtsDetJac().data().begin());
469 }
470 
471 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
472  getFTensor1CoordsAtGaussPts() {
474  &getCoordsAtGaussPts()(0, 0), &getCoordsAtGaussPts()(0, 1),
475  &getCoordsAtGaussPts()(0, 2));
476 }
477 
478 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
479  getFTensor1HoCoordsAtGaussPts() {
480  double *ptr = &*getHoCoordsAtGaussPts().data().begin();
481  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, ptr + 1,
482  ptr + 2);
483 }
484 
485 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
486  getFTensor2HoGaussPtsJac() {
487  double *ptr = &*getHoGaussPtsJac().data().begin();
489  ptr, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
490  ptr + 8);
491 }
492 
493 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
494  getFTensor2HoGaussPtsInvJac() {
495  double *ptr = &*getHoGaussPtsInvJac().data().begin();
497  ptr, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
498  ptr + 8);
499 }
500 
502 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getVolumeFE() const {
503  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE);
504 }
505 
506 } // namespace MoFEM
507 
508 #endif //__VOLUMEELEMENTFORCESANDSOURCESCORE_HPP__
509 
510 /**
511  * \defgroup mofem_forces_and_sources_volume_element Volume Element
512  * \brief Implementation of general volume element.
513  *
514  * \ingroup mofem_forces_and_sources
515  **/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
T data[Tensor_Dim]
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:102
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Data on single entity (This is passed as argument to DataOperator::doWork)
Deprecated interface functions.
Data operator to do calculations at integration points.
structure to get information form mofem into DataForcesAndSourcesCore
apply contravariant (Piola) transfer to Hdiv space
apply covariant transfer to Hcurl space
Apply contravariant (Piola) transfer to Hdiv space for HO geometr.
Apply covariant (Piola) transfer to Hcurl space for HO geometry.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Transform local reference derivatives of shape function to global derivatives.
brief Transform local reference derivatives of shape function to global derivatives
OpGetDataAndGradient< 3, 3 > opHOatGaussPoints
higher order geometry data at Gauss pts
VolumeElementForcesAndSourcesCoreSwitch(Interface &m_field, const EntityType type=MBTET)