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  protected:
219 
220  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
221 
222  };
223 
224  enum Switches {
225  NO_HO_GEOMETRY = 1 << 0 | 1 << 2,
226  NO_TRANSFORM = 1 << 1 | 1 << 2,
227  NO_HO_TRANSFORM = 1 << 2
228  };
229 
230  template <int SWITCH> MoFEMErrorCode OpSwitch();
231 
232 protected:
233  // Note that functions below could be overloaded by user to change default
234  // behavior of the element.
235 
236  /**
237  * \brief Set integration points
238  * @return Error code
239  */
240  virtual MoFEMErrorCode setIntegrationPts();
241 
242  /**
243  * \brief Calculate element volume and Jacobian
244  *
245  * Note that at that point is assumed that geometry is exclusively defined by
246  * corner nodes.
247  *
248  * @return Error code
249  */
250  virtual MoFEMErrorCode calculateVolumeAndJacobian();
251 
252  /**
253  * \brief Calculate coordinate at integration points
254  * @return Error code
255  */
256  virtual MoFEMErrorCode calculateCoordinatesAtGaussPts();
257 
258  /**
259  * \brief Determine approximation space and order of base functions
260  * @return Error code
261  */
262  virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement();
263 
264  /**
265  * \brief Transform base functions based on geometric element Jacobian.
266  *
267  * This function apply transformation to base functions and its derivatives.
268  * For example when base functions for H-div are present the
269  * Piola-Transformarion is applied to base functions and their derivatives.
270  *
271  * @return Error code
272  */
273  virtual MoFEMErrorCode transformBaseFunctions();
274 
275  /** \brief Calculate Jacobian for HO geometry
276  *
277  * MoFEM use hierarchical approximate base to describe geometry of the body.
278  * This function transform derivatives of base functions when HO geometry is
279  * set and calculate Jacobian, inverse of Jacobian and determinant of
280  * transformation.
281  *
282  */
283  virtual MoFEMErrorCode calculateHoJacobian();
284 
285  /**
286  * \brief Transform base functions based on ho-geometry element Jacobian.
287  *
288  * This function apply transformation to base functions and its derivatives.
289  * For example when base functions for H-div are present the
290  * Piola-Transformarion is applied to base functions and their derivatives.
291  *
292  * @return Error code
293  */
294  virtual MoFEMErrorCode transformHoBaseFunctions();
295 
299 
304 
309 
311  opHOatGaussPoints; ///< higher order geometry data at Gauss pts
316 
317  double vOlume;
318 
323 
325 
326  friend class UserDataOperator;
327 };
328 
329 /**
330  * @brief Volume finite element with switches
331  *
332  * Using SWITCH to off functions
333  *
334  * @tparam SWITCH
335  */
336 template <int SWITCH>
339 
340  using VolumeElementForcesAndSourcesCoreBase::
341  VolumeElementForcesAndSourcesCoreBase;
342  using UserDataOperator =
344 
345  MoFEMErrorCode operator()();
346 };
347 
348 /** \brief Volume finite element default
349  \ingroup mofem_forces_and_sources_volume_element
350 
351  */
354 
355 template <int SWITCH>
356 MoFEMErrorCode VolumeElementForcesAndSourcesCoreBase::OpSwitch() {
358 
359  if (numeredEntFiniteElementPtr->getEntType() != MBTET)
361  CHKERR createDataOnElement();
362 
363  CHKERR calculateVolumeAndJacobian();
364  CHKERR getSpaceBaseAndOrderOnElement();
365  CHKERR setIntegrationPts();
366  if (gaussPts.size2() == 0)
368  CHKERR calculateCoordinatesAtGaussPts();
369  CHKERR calHierarchicalBaseFunctionsOnElement();
370  CHKERR calBernsteinBezierBaseFunctionsOnElement();
371 
372  if (!(NO_TRANSFORM & SWITCH))
373  CHKERR transformBaseFunctions();
374 
375  if (!(NO_HO_GEOMETRY & SWITCH))
376  CHKERR calculateHoJacobian();
377 
378  if (!(NO_HO_TRANSFORM & SWITCH))
379  CHKERR transformHoBaseFunctions();
380 
381  // Iterate over operators
382  CHKERR loopOverOperators();
383 
385 }
386 
387 template <int SWITCH>
389  return OpSwitch<SWITCH>();
390 }
391 
392 int VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getNumNodes() {
393  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->num_nodes;
394 }
395 
396 const EntityHandle *
397 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getConn() {
398  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->conn;
399 }
400 
401 double
402 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getVolume() const {
403  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->vOlume;
404 }
405 
406 double &VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getVolume() {
407  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->vOlume;
408 }
409 
411 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getJac() {
412  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->tJac;
413 }
414 
416 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getInvJac() {
417  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->tInvJac;
418 }
419 
420 double
421 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() const {
422  return getVolume();
423 }
424 
425 double &VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() {
426  return getVolume();
427 }
428 
429 VectorDouble &
430 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getCoords() {
431  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)->coords;
432 }
433 
434 MatrixDouble &
435 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getCoordsAtGaussPts() {
436  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
437  ->coordsAtGaussPts;
438 }
439 
440 MatrixDouble &VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
441  getHoCoordsAtGaussPts() {
442  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
443  ->hoCoordsAtGaussPts;
444 }
445 
446 MatrixDouble &
447 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getHoGaussPtsJac() {
448  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
449  ->hoGaussPtsJac;
450 }
451 
452 MatrixDouble &
453 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getHoGaussPtsInvJac() {
454  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
455  ->hoGaussPtsInvJac;
456 }
457 
458 VectorDouble &
459 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getHoGaussPtsDetJac() {
460  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE)
461  ->hoGaussPtsDetJac;
462 }
463 
464 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
465  getFTenosr0HoMeasure() {
467  &*getHoGaussPtsDetJac().data().begin());
468 }
469 
470 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
471  getFTensor1CoordsAtGaussPts() {
473  &getCoordsAtGaussPts()(0, 0), &getCoordsAtGaussPts()(0, 1),
474  &getCoordsAtGaussPts()(0, 2));
475 }
476 
477 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
478  getFTensor1HoCoordsAtGaussPts() {
479  double *ptr = &*getHoCoordsAtGaussPts().data().begin();
480  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, ptr + 1,
481  ptr + 2);
482 }
483 
484 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
485  getFTensor2HoGaussPtsJac() {
486  double *ptr = &*getHoGaussPtsJac().data().begin();
488  ptr, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
489  ptr + 8);
490 }
491 
492 auto VolumeElementForcesAndSourcesCoreBase::UserDataOperator::
493  getFTensor2HoGaussPtsInvJac() {
494  double *ptr = &*getHoGaussPtsInvJac().data().begin();
496  ptr, ptr + 1, ptr + 2, ptr + 3, ptr + 4, ptr + 5, ptr + 6, ptr + 7,
497  ptr + 8);
498 }
499 
501 VolumeElementForcesAndSourcesCoreBase::UserDataOperator::getVolumeFE() const {
502  return static_cast<VolumeElementForcesAndSourcesCoreBase *>(ptrFE);
503 }
504 
505 } // namespace MoFEM
506 
507 #endif //__VOLUMEELEMENTFORCESANDSOURCESCORE_HPP__
508 
509 /**
510  * \defgroup mofem_forces_and_sources_volume_element Volume Element
511  * \brief Implementation of general volume element.
512  *
513  * \ingroup mofem_forces_and_sources
514  **/
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:481
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:512
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:97
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:600
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:411
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...