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