v0.9.1
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...
Deprecated interface functions.
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:482
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
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:601
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:412
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...