v0.9.0
FaceElementForcesAndSourcesCore.hpp
Go to the documentation of this file.
1 /** \file FaceElementForcesAndSourcesCore.hpp
2  \brief Implementation of face element.
3 
4 */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 #ifndef __FACEELEMENTFORCESANDSOURCESCORE_HPP__
21 #define __FACEELEMENTFORCESANDSOURCESCORE_HPP__
22 
23 using namespace boost::numeric;
24 
25 namespace MoFEM {
26 
27 template <int SWITCH> struct VolumeElementForcesAndSourcesCoreOnSideSwitch;
28 
29 /** \brief Face finite element
30  \ingroup mofem_forces_and_sources_tri_element
31 
32  User is implementing own operator at Gauss point level, by own object
33  derived from FaceElementForcesAndSourcesCoreL::UserDataOperator. Arbitrary
34  number of operator added pushing objects to OpPtrVector
35 
36  */
38 
39  std::string meshPositionsFieldName; ///< Name of the field with geometry
41 
42  /** \brief default operator for TRI element
43  * \ingroup mofem_forces_and_sources_tri_element
44  */
46 
48 
49  /**
50  * \brief get area of face
51  * @return area of face
52  */
53  inline double getArea();
54 
55  /**
56  * \brief get measure of element
57  * @return area of face
58  */
59  inline double getMeasure();
60 
61  /** \brief get triangle normal
62  */
63  inline VectorDouble &getNormal();
64 
65  /** \brief get triangle tangent 1
66  */
67  inline VectorDouble &getTangent1();
68 
69  /** \brief get triangle tangent 2
70  */
71  inline VectorDouble &getTangent2();
72 
73  /** \brief get normal as tensor
74  */
75  inline auto getFTensor1Normal();
76 
77  /** \brief get tangentOne as tensor
78  */
79  inline auto getFTensor1Tangent1();
80 
81  /** \brief get tangentTwo as tensor
82  */
83  inline auto getFTensor1Tangent2();
84 
85  /** \brief get element number of nodes
86  */
87  inline int getNumNodes();
88 
89  /** \brief get element connectivity
90  */
91  inline const EntityHandle *getConn();
92 
93  /** \brief get triangle coordinates
94  */
95  inline VectorDouble &getCoords();
96 
97  /**
98  * \brief get get coords at gauss points
99 
100  \code
101  FTensor::Index<'i',3> i;
102  FTensor::Tensor1<double,3> t_center;
103  auto t_coords = getFTensor1Coords();
104  t_center(i) = 0;
105  for(int nn = 0;nn!=3;nn++) {
106  t_center(i) += t_coords(i);
107  ++t_coords;
108  }
109  t_center(i) /= 3;
110  \endcode
111 
112  */
113  inline auto getFTensor1Coords();
114 
115  /** \brief Gauss points and weight, matrix (nb. of points x 3)
116 
117  Column 0-2 integration points coordinate x and y, respectively. At rows are
118  integration points.
119 
120  */
121  inline MatrixDouble &getCoordsAtGaussPts();
122 
123  /** \brief get coordinates at Gauss pts.
124  */
125  inline auto getFTensor1CoordsAtGaussPts();
126 
127  /** \brief coordinate at Gauss points (if hierarchical approximation of
128  element geometry)
129 
130  Note: returned matrix has size 0 in rows and columns if no HO approximation
131  of geometry is available.
132 
133  */
134  inline MatrixDouble &getHoCoordsAtGaussPts();
135 
136  /** \brief get coordinates at Gauss pts (takes in account ho approx. of
137  * geometry)
138  */
139  inline auto getFTensor1HoCoordsAtGaussPts();
140 
141  /** \brief if higher order geometry return normals at Gauss pts.
142 
143  Note: returned matrix has size 0 in rows and columns if no HO approximation
144  of geometry is available.
145 
146  */
147  inline MatrixDouble &getNormalsAtGaussPts();
148 
149  /**
150  * @deprecated Use getNormalsAtGaussPts
151  */
152  DEPRECATED inline MatrixDouble &getNormalsAtGaussPt();
153 
154  /** \brief if higher order geometry return normals at Gauss pts.
155  *
156  * \param gg gauss point number
157  */
158  inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPts(const int gg);
159 
160  /**
161  * @deprecated Cotrect name is getNormalsAtGaussPts
162  */
163  DEPRECATED inline ublas::matrix_row<MatrixDouble>
164  getNormalsAtGaussPt(const int gg);
165 
166  /** \brief if higher order geometry return tangent vector to triangle at
167  Gauss pts.
168 
169  Note: returned matrix has size 0 in rows and columns if no HO approximation
170  of geometry is avaliable.
171 
172  */
173  inline MatrixDouble &getTangent1AtGaussPts();
174 
175  /** \brief if higher order geometry return tangent vector to triangle at
176  Gauss pts.
177 
178  Note: returned matrix has size 0 in rows and columns if no HO approximation
179  of geometry is avaliable.
180 
181  */
182  inline MatrixDouble &getTangent2AtGaussPts();
183 
184  /** \brief get normal at integration points
185 
186  Example:
187  \code
188  double nrm2;
189  FTensor::Index<'i',3> i;
190  auto t_normal = getFTensor1NormalsAtGaussPts();
191  for(int gg = gg!=data.getN().size1();gg++) {
192  nrm2 = sqrt(t_normal(i)*t_normal(i));
193  ++t_normal;
194  }
195  \endcode
196 
197  */
198  inline auto getFTensor1NormalsAtGaussPts();
199 
200  /** \brief get tangent 1 at integration points
201 
202  */
203  inline auto getFTensor1Tangent1AtGaussPts();
204 
205  /** \brief get tangent 2 at integration points
206 
207  */
208  inline auto getFTensor1Tangent2AtGaussPts();
209 
210  /** \brief return pointer to Generic Triangle Finite Element object
211  */
212  inline const FaceElementForcesAndSourcesCoreBase *getFaceFE();
213 
214  /**
215  * @deprecated Use getFaceFE
216  */
218  getFaceElementForcesAndSourcesCore();
219 
220  /**
221  *
222  * User call this function to loop over elements on the side of face. This
223  * function calls MoFEM::VolumeElementForcesAndSourcesCoreOnSide with is
224  * operator to do calculations.
225  *
226  * @param fe_name Name of the element
227  * @param method Finite element object
228  * @return error code
229  */
230  template <int SWITCH>
231  MoFEMErrorCode loopSideVolumes(
232  const string &fe_name,
234  };
235 
236  enum Switches {
237  NO_HO_GEOMETRY = 1 << 0,
238  NO_CONTRAVARIANT_TRANSFORM_HDIV = 1 << 1,
239  NO_COVARIANT_TRANSFORM_HCURL = 1 << 2,
240  };
241 
242  template <int SWITCH> MoFEMErrorCode OpSwitch();
243 
244 protected:
245  MoFEMErrorCode getNumberOfNodes(int &num_nodes) const;
246 
247  /**
248  * \brief Calculate element area and normal of the face at integration points
249  *
250  * @return Error code
251  */
252  virtual MoFEMErrorCode calculateAreaAndNormalAtIntegrationPts();
253 
254  /**
255  * \brief Calculate element area and normal of the face
256  *
257  * Note that at that point is assumed that geometry is exclusively defined
258  * by corner nodes.
259  *
260  * @return Error code
261  */
262  virtual MoFEMErrorCode calculateAreaAndNormal();
263 
264  /**
265  * \brief Set integration points
266  * @return Error code
267  */
268  virtual MoFEMErrorCode setIntegrationPts();
269 
270  /**
271  * \brief Determine approximation space and order of base functions
272  * @return Error code
273  */
274  virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement();
275 
276  /**
277  * \brief Calculate coordinate at integration points
278  * @return Error code
279  */
280  virtual MoFEMErrorCode calculateCoordinatesAtGaussPts();
281 
282  /**
283  * \brief Calculate normal on curved elements
284  *
285  * Geometry is given by other field.
286  *
287  * @return error code
288  */
289  virtual MoFEMErrorCode calculateHoNormal();
290 
291  double aRea;
294  VectorDouble nOrmal, tangentOne, tangentTwo;
297 
305 
306  friend class UserDataOperator;
308 };
309 
310 /** \brief Face finite element switched
311  \ingroup mofem_forces_and_sources_tri_element
312 
313  */
314 template <int SWITCH>
317 
318  using FaceElementForcesAndSourcesCoreBase::
319  FaceElementForcesAndSourcesCoreBase;
320 
321  using UserDataOperator =
323 
324  MoFEMErrorCode operator()();
325 };
326 
327 /** \brief Face finite element default
328  \ingroup mofem_forces_and_sources_tri_element
329 
330  */
333 
334 template <int SWITCH>
335 MoFEMErrorCode FaceElementForcesAndSourcesCoreBase::OpSwitch() {
337 
338  const EntityType type = numeredEntFiniteElementPtr->getEntType();
339  if (type != lastEvaluatedElementEntityType) {
340  switch (type) {
341  case MBTRI:
342  getElementPolynomialBase() =
343  boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
344  break;
345  case MBQUAD:
346  getElementPolynomialBase() =
347  boost::shared_ptr<BaseFunction>(new QuadPolynomialBase());
348  break;
349  default:
351  }
352  CHKERR createDataOnElement();
353  }
354 
355  // Calculate normal and tangent vectors for face geometry
356  CHKERR calculateAreaAndNormal();
357  CHKERR getSpaceBaseAndOrderOnElement();
358 
359  CHKERR setIntegrationPts();
360  if (gaussPts.size2() == 0)
362 
363  DataForcesAndSourcesCore &data_curl = *dataOnElement[HCURL];
364  DataForcesAndSourcesCore &data_div = *dataOnElement[HDIV];
365 
366  CHKERR calculateCoordinatesAtGaussPts();
367  CHKERR calHierarchicalBaseFunctionsOnElement();
368  CHKERR calBernsteinBezierBaseFunctionsOnElement();
369 
370  switch (numeredEntFiniteElementPtr->getEntType()) {
371  case MBTRI:
372  break;
373  case MBQUAD:
374  CHKERR calculateAreaAndNormalAtIntegrationPts();
375  break;
376  default:
377  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
378  "Element type not implemented");
379  }
380 
381  if (!(NO_HO_GEOMETRY & SWITCH)) {
382  CHKERR calculateHoNormal();
383  }
384 
385  // Apply Piola transform to HDiv and HCurl spaces, uses previously
386  // calculated faces normal and tangent vectors.
387  if (!(NO_CONTRAVARIANT_TRANSFORM_HDIV & SWITCH))
388  if (dataH1.spacesOnEntities[MBTRI].test(HDIV))
389  CHKERR opContravariantTransform.opRhs(data_div);
390 
391  if (!(NO_COVARIANT_TRANSFORM_HCURL & SWITCH))
392  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL))
393  CHKERR opCovariantTransform.opRhs(data_curl);
394 
395  // Iterate over operators
396  CHKERR loopOverOperators();
397 
399 }
400 
401 template <int SWITCH>
403  return OpSwitch<SWITCH>();
404 }
405 
406 double FaceElementForcesAndSourcesCoreBase::UserDataOperator::getArea() {
407  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->aRea;
408 }
409 
410 double FaceElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() {
411  return getArea();
412 }
413 
414 VectorDouble &
415 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormal() {
416  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->nOrmal;
417 }
418 
419 VectorDouble &
420 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent1() {
421  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->tangentOne;
422 }
423 
424 VectorDouble &
425 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent2() {
426  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->tangentTwo;
427 }
428 
429 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
430  getFTensor1Normal() {
431  double *ptr = &*getNormal().data().begin();
432  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
433 }
434 
435 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
436  getFTensor1Tangent1() {
437  double *ptr = &*getTangent1().data().begin();
438  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
439 }
440 
441 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
442  getFTensor1Tangent2() {
443  double *ptr = &*getTangent2().data().begin();
444  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
445 }
446 
447 int FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNumNodes() {
448  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->num_nodes;
449 }
450 
451 const EntityHandle *
452 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getConn() {
453  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->conn;
454 }
455 
456 VectorDouble &
457 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getCoords() {
458  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->coords;
459 }
460 
461 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
462  getFTensor1Coords() {
463  double *ptr = &*getCoords().data().begin();
464  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
465  &ptr[2]);
466 }
467 
468 MatrixDouble &
469 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getCoordsAtGaussPts() {
470  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
471  ->coordsAtGaussPts;
472 }
473 
474 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
475  getFTensor1CoordsAtGaussPts() {
476  double *ptr = &*getCoordsAtGaussPts().data().begin();
477  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
478  &ptr[2]);
479 }
480 
481 MatrixDouble &
482 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getHoCoordsAtGaussPts() {
483  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
484  ->hoCoordsAtGaussPts;
485 }
486 
487 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
488  getFTensor1HoCoordsAtGaussPts() {
489  if (getHoCoordsAtGaussPts().size1() == 0 &&
490  getHoCoordsAtGaussPts().size2() != 3) {
491  return getFTensor1Coords();
492  }
493  double *ptr = &*getHoCoordsAtGaussPts().data().begin();
494  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
495  &ptr[2]);
496 }
497 
498 MatrixDouble &
499 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPts() {
500  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
501  ->normalsAtGaussPts;
502 }
503 
504 MatrixDouble &
505 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPt() {
506  return getNormalsAtGaussPts();
507 }
508 
509 ublas::matrix_row<MatrixDouble>
510 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPts(
511  const int gg) {
512  return ublas::matrix_row<MatrixDouble>(
513  static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
514  ->normalsAtGaussPts,
515  gg);
516 }
517 
518 ublas::matrix_row<MatrixDouble>
519 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPt(
520  const int gg) {
521  return getNormalsAtGaussPts(gg);
522 }
523 
524 MatrixDouble &
525 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent1AtGaussPts() {
526  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
527  ->tangentOneAtGaussPts;
528 }
529 
530 MatrixDouble &
531 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent2AtGaussPts() {
532  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
533  ->tangentTwoAtGaussPts;
534 }
535 
536 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
537  getFTensor1NormalsAtGaussPts() {
538  double *ptr = &*getNormalsAtGaussPts().data().begin();
539  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
540  &ptr[2]);
541 }
542 
543 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
544  getFTensor1Tangent1AtGaussPts() {
545  double *ptr = &*getTangent1AtGaussPts().data().begin();
546  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
547  &ptr[2]);
548 }
549 
550 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
551  getFTensor1Tangent2AtGaussPts() {
552  double *ptr = &*getTangent2AtGaussPts().data().begin();
553  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
554  &ptr[2]);
555 }
556 
558 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getFaceFE() {
559  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE);
560 }
561 
562 const FaceElementForcesAndSourcesCoreBase *FaceElementForcesAndSourcesCoreBase::
563  UserDataOperator::getFaceElementForcesAndSourcesCore() {
564  return getFaceFE();
565 }
566 
567 template <int SWITCH>
569 FaceElementForcesAndSourcesCoreBase::UserDataOperator::loopSideVolumes(
570  const string &fe_name,
572  return loopSide(fe_name, &fe_method, 3);
573 }
574 
575 } // namespace MoFEM
576 
577 #endif //__FACEELEMENTFORCESANDSOURCESCORE_HPP__
578 
579 /**
580  * \defgroup mofem_forces_and_sources_tri_element Face Element
581  * \brief Implementation of face element
582  *
583  * \ingroup mofem_forces_and_sources
584  **/
field with continuous normal traction
Definition: definitions.h:173
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
transform Hcurl base fluxes from reference element to physical triangle
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
std::string meshPositionsFieldName
Name of the field with geometry.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
T data[Tensor_Dim]
OpSetContravariantPiolaTransformOnFace opContravariantTransform
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Calculate base functions on triangle.
transform Hdiv base fluxes from reference element to physical triangle
field with continuous tangents
Definition: definitions.h:172
#define CHKERR
Inline error check.
Definition: definitions.h:596
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Calculate base functions on triangle.
#define DEPRECATED
Definition: definitions.h:29
structure to get information form mofem into DataForcesAndSourcesCore
Calculate normals at Gauss points of triangle element.
#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
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...