v0.9.1
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  private:
236  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
237 
238  };
239 
240  enum Switches {
241  NO_HO_GEOMETRY = 1 << 0,
242  NO_CONTRAVARIANT_TRANSFORM_HDIV = 1 << 1,
243  NO_COVARIANT_TRANSFORM_HCURL = 1 << 2,
244  };
245 
246  template <int SWITCH> MoFEMErrorCode OpSwitch();
247 
248 protected:
249 
250  MoFEMErrorCode getNumberOfNodes(int &num_nodes) const;
251 
252  /**
253  * \brief Calculate element area and normal of the face at integration points
254  *
255  * @return Error code
256  */
257  virtual MoFEMErrorCode calculateAreaAndNormalAtIntegrationPts();
258 
259  /**
260  * \brief Calculate element area and normal of the face
261  *
262  * Note that at that point is assumed that geometry is exclusively defined
263  * by corner nodes.
264  *
265  * @return Error code
266  */
267  virtual MoFEMErrorCode calculateAreaAndNormal();
268 
269  /**
270  * \brief Set integration points
271  * @return Error code
272  */
273  virtual MoFEMErrorCode setIntegrationPts();
274 
275  /**
276  * \brief Determine approximation space and order of base functions
277  * @return Error code
278  */
279  virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement();
280 
281  /**
282  * \brief Calculate coordinate at integration points
283  * @return Error code
284  */
285  virtual MoFEMErrorCode calculateCoordinatesAtGaussPts();
286 
287  /**
288  * \brief Calculate normal on curved elements
289  *
290  * Geometry is given by other field.
291  *
292  * @return error code
293  */
294  virtual MoFEMErrorCode calculateHoNormal();
295 
296  double aRea;
299  VectorDouble nOrmal, tangentOne, tangentTwo;
302 
310 
311  friend class UserDataOperator;
313 };
314 
315 /** \brief Face finite element switched
316  \ingroup mofem_forces_and_sources_tri_element
317 
318  */
319 template <int SWITCH>
322 
323  using FaceElementForcesAndSourcesCoreBase::
324  FaceElementForcesAndSourcesCoreBase;
325 
326  using UserDataOperator =
328 
329  MoFEMErrorCode operator()();
330 };
331 
332 /** \brief Face finite element default
333  \ingroup mofem_forces_and_sources_tri_element
334 
335  */
338 
339 template <int SWITCH>
340 MoFEMErrorCode FaceElementForcesAndSourcesCoreBase::OpSwitch() {
342 
343  const EntityType type = numeredEntFiniteElementPtr->getEntType();
344  if (type != lastEvaluatedElementEntityType) {
345  switch (type) {
346  case MBTRI:
347  getElementPolynomialBase() =
348  boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
349  break;
350  case MBQUAD:
351  getElementPolynomialBase() =
352  boost::shared_ptr<BaseFunction>(new QuadPolynomialBase());
353  break;
354  default:
356  }
357  CHKERR createDataOnElement();
358  }
359 
360  // Calculate normal and tangent vectors for face geometry
361  CHKERR calculateAreaAndNormal();
362  CHKERR getSpaceBaseAndOrderOnElement();
363 
364  CHKERR setIntegrationPts();
365  if (gaussPts.size2() == 0)
367 
368  DataForcesAndSourcesCore &data_curl = *dataOnElement[HCURL];
369  DataForcesAndSourcesCore &data_div = *dataOnElement[HDIV];
370 
371  CHKERR calculateCoordinatesAtGaussPts();
372  CHKERR calHierarchicalBaseFunctionsOnElement();
373  CHKERR calBernsteinBezierBaseFunctionsOnElement();
374 
375  switch (numeredEntFiniteElementPtr->getEntType()) {
376  case MBTRI:
377  break;
378  case MBQUAD:
379  CHKERR calculateAreaAndNormalAtIntegrationPts();
380  break;
381  default:
382  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
383  "Element type not implemented");
384  }
385 
386  if (!(NO_HO_GEOMETRY & SWITCH)) {
387  CHKERR calculateHoNormal();
388  }
389 
390  // Apply Piola transform to HDiv and HCurl spaces, uses previously
391  // calculated faces normal and tangent vectors.
392  if (!(NO_CONTRAVARIANT_TRANSFORM_HDIV & SWITCH))
393  if (dataH1.spacesOnEntities[MBTRI].test(HDIV))
394  CHKERR opContravariantTransform.opRhs(data_div);
395 
396  if (!(NO_COVARIANT_TRANSFORM_HCURL & SWITCH))
397  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL))
398  CHKERR opCovariantTransform.opRhs(data_curl);
399 
400  // Iterate over operators
401  CHKERR loopOverOperators();
402 
404 }
405 
406 template <int SWITCH>
408  return OpSwitch<SWITCH>();
409 }
410 
411 double FaceElementForcesAndSourcesCoreBase::UserDataOperator::getArea() {
412  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->aRea;
413 }
414 
415 double FaceElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() {
416  return getArea();
417 }
418 
419 VectorDouble &
420 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormal() {
421  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->nOrmal;
422 }
423 
424 VectorDouble &
425 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent1() {
426  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->tangentOne;
427 }
428 
429 VectorDouble &
430 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent2() {
431  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->tangentTwo;
432 }
433 
434 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
435  getFTensor1Normal() {
436  double *ptr = &*getNormal().data().begin();
437  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
438 }
439 
440 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
441  getFTensor1Tangent1() {
442  double *ptr = &*getTangent1().data().begin();
443  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
444 }
445 
446 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
447  getFTensor1Tangent2() {
448  double *ptr = &*getTangent2().data().begin();
449  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
450 }
451 
452 int FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNumNodes() {
453  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->num_nodes;
454 }
455 
456 const EntityHandle *
457 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getConn() {
458  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->conn;
459 }
460 
461 VectorDouble &
462 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getCoords() {
463  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->coords;
464 }
465 
466 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
467  getFTensor1Coords() {
468  double *ptr = &*getCoords().data().begin();
469  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
470  &ptr[2]);
471 }
472 
473 MatrixDouble &
474 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getCoordsAtGaussPts() {
475  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
476  ->coordsAtGaussPts;
477 }
478 
479 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
480  getFTensor1CoordsAtGaussPts() {
481  double *ptr = &*getCoordsAtGaussPts().data().begin();
482  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
483  &ptr[2]);
484 }
485 
486 MatrixDouble &
487 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getHoCoordsAtGaussPts() {
488  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
489  ->hoCoordsAtGaussPts;
490 }
491 
492 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
493  getFTensor1HoCoordsAtGaussPts() {
494  if (getHoCoordsAtGaussPts().size1() == 0 &&
495  getHoCoordsAtGaussPts().size2() != 3) {
496  return getFTensor1Coords();
497  }
498  double *ptr = &*getHoCoordsAtGaussPts().data().begin();
499  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
500  &ptr[2]);
501 }
502 
503 MatrixDouble &
504 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPts() {
505  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
506  ->normalsAtGaussPts;
507 }
508 
509 MatrixDouble &
510 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPt() {
511  return getNormalsAtGaussPts();
512 }
513 
514 ublas::matrix_row<MatrixDouble>
515 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPts(
516  const int gg) {
517  return ublas::matrix_row<MatrixDouble>(
518  static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
519  ->normalsAtGaussPts,
520  gg);
521 }
522 
523 ublas::matrix_row<MatrixDouble>
524 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPt(
525  const int gg) {
526  return getNormalsAtGaussPts(gg);
527 }
528 
529 MatrixDouble &
530 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent1AtGaussPts() {
531  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
532  ->tangentOneAtGaussPts;
533 }
534 
535 MatrixDouble &
536 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent2AtGaussPts() {
537  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
538  ->tangentTwoAtGaussPts;
539 }
540 
541 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
542  getFTensor1NormalsAtGaussPts() {
543  double *ptr = &*getNormalsAtGaussPts().data().begin();
544  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
545  &ptr[2]);
546 }
547 
548 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
549  getFTensor1Tangent1AtGaussPts() {
550  double *ptr = &*getTangent1AtGaussPts().data().begin();
551  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
552  &ptr[2]);
553 }
554 
555 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
556  getFTensor1Tangent2AtGaussPts() {
557  double *ptr = &*getTangent2AtGaussPts().data().begin();
558  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
559  &ptr[2]);
560 }
561 
563 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getFaceFE() {
564  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE);
565 }
566 
567 const FaceElementForcesAndSourcesCoreBase *FaceElementForcesAndSourcesCoreBase::
568  UserDataOperator::getFaceElementForcesAndSourcesCore() {
569  return getFaceFE();
570 }
571 
572 template <int SWITCH>
574 FaceElementForcesAndSourcesCoreBase::UserDataOperator::loopSideVolumes(
575  const string &fe_name,
577  return loopSide(fe_name, &fe_method, 3);
578 }
579 
580 } // namespace MoFEM
581 
582 #endif //__FACEELEMENTFORCESANDSOURCESCORE_HPP__
583 
584 /**
585  * \defgroup mofem_forces_and_sources_tri_element Face Element
586  * \brief Implementation of face element
587  *
588  * \ingroup mofem_forces_and_sources
589  **/
Deprecated interface functions.
field with continuous normal traction
Definition: definitions.h:178
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:482
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:513
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:177
#define CHKERR
Inline error check.
Definition: definitions.h:601
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:412
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...