v0.8.23
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 /** \brief Face finite element
28  \ingroup mofem_forces_and_sources_tri_element
29 
30  User is implementing own operator at Gauss point level, by own object
31  derived from FaceElementForcesAndSourcesCoreL::UserDataOperator. Arbitrary
32  number of operator added pushing objects to OpPtrVector
33 
34  */
36 
37  double aRea;
38  ;
39  int num_nodes;
41  VectorDouble nOrmal, tangentOne, tangentTwo;
44 
45  std::string meshPositionsFieldName; ///< Name of the field with geometry
46 
54 
56 
57  /** \brief default operator for TRI element
58  * \ingroup mofem_forces_and_sources_tri_element
59  */
61 
64 
65  UserDataOperator(const std::string &field_name, const char type)
66  : ForcesAndSourcesCore::UserDataOperator(field_name, type) {}
67 
68  UserDataOperator(const std::string &row_field_name,
69  const std::string &col_field_name, const char type,
70  const bool symm = true)
71  : ForcesAndSourcesCore::UserDataOperator(row_field_name, col_field_name,
72  type, symm){};
73 
74  /**
75  * \brief get area of face
76  * @return area of face
77  */
78  inline double getArea() {
79  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)->aRea;
80  }
81 
82  /**
83  * \brief get measure of element
84  * @return area of face
85  */
86  inline double getMeasure() { return getArea(); }
87 
88  /** \brief get triangle normal
89  */
90  inline VectorDouble &getNormal() {
91  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)->nOrmal;
92  }
93 
94  /** \brief get triangle tangent 1
95  */
97  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)->tangentOne;
98  }
99 
100  /** \brief get triangle tangent 2
101  */
103  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)->tangentTwo;
104  }
105 
106  /** \brief get normal as tensor
107  */
108  inline auto getFTensor1Normal() {
109  double *ptr = &*getNormal().data().begin();
110  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
111  }
112 
113  /// \deprecated use getFTensor1Normal()
114  DEPRECATED inline auto getTensor1Normal() { return getFTensor1Normal(); }
115 
116  /** \brief get tangentOne as tensor
117  */
118  inline auto getFTensor1Tangent1() {
119  double *ptr = &*getTangent1().data().begin();
120  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
121  }
122 
123  /// \deprecated use getFTensor1Tangent1
125  return getFTensor1Tangent1();
126  }
127 
128  /** \brief get tangentTwo as tensor
129  */
130  inline auto getFTensor2Tangent1() {
131  double *ptr = &*getTangent2().data().begin();
132  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
133  }
134 
135  /// \deprecated use getFTensor2Tangent1
137  return getFTensor2Tangent1();
138  }
139 
140  /** \brief get element number of nodes
141  */
142  inline int getNumNodes() {
143  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)->num_nodes;
144  }
145 
146  /** \brief get element connectivity
147  */
148  inline const EntityHandle *getConn() {
149  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)->conn;
150  }
151 
152  /** \brief get triangle coordinates
153  */
155  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)->coords;
156  }
157 
158  /**
159  * \brief get get coords at gauss points
160 
161  \code
162  FTensor::Index<'i',3> i;
163  FTensor::Tensor1<double,3> t_center;
164  auto t_coords = getFTensor1Coords();
165  t_center(i) = 0;
166  for(int nn = 0;nn!=3;nn++) {
167  t_center(i) += t_coords(i);
168  ++t_coords;
169  }
170  t_center(i) /= 3;
171  \endcode
172 
173  */
174  inline auto getFTensor1Coords() {
175  double *ptr = &*getCoords().data().begin();
176  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
177  &ptr[2]);
178  }
179 
180  /// \deprecated use getFTensor1Coords
181  DEPRECATED inline auto getTensor1Coords() { return getFTensor1Coords(); }
182 
183  /** \brief Gauss points and weight, matrix (nb. of points x 3)
184 
185  Column 0-2 integration points coordinate x and y, respectively. At rows are
186  integration points.
187 
188  */
190  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)
191  ->coordsAtGaussPts;
192  }
193 
194  /** \brief get coordinates at Gauss pts.
195  */
197  double *ptr = &*getCoordsAtGaussPts().data().begin();
198  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
199  &ptr[2]);
200  }
201 
202  /// \deprecated use getFTensor1CoordsAtGaussPts
204  return getFTensor1CoordsAtGaussPts();
205  }
206 
207  /** \brief coordinate at Gauss points (if hierarchical approximation of
208  element geometry)
209 
210  Note: returned matrix has size 0 in rows and columns if no HO approximation
211  of geometry is available.
212 
213  */
215  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)
216  ->hoCoordsAtGaussPts;
217  }
218 
219  /** \brief get coordinates at Gauss pts (takes in account ho approx. of
220  * geometry)
221  */
223  if (getHoCoordsAtGaussPts().size1() == 0 &&
224  getHoCoordsAtGaussPts().size2() != 3) {
225  return getFTensor1Coords();
226  }
227  double *ptr = &*getHoCoordsAtGaussPts().data().begin();
228  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
229  &ptr[2]);
230  }
231 
232  /// \deprecated use getTensor1HoCoordsAtGaussPts
234  return getFTensor1HoCoordsAtGaussPts();
235  }
236 
237  /** \brief if higher order geometry return normals at Gauss pts.
238 
239  Note: returned matrix has size 0 in rows and columns if no HO approximation
240  of geometry is available.
241 
242  */
244  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)
245  ->normalsAtGaussPts;
246  }
247 
248  /// \deprecated use getNormalsAtGaussPts
250  return getNormalsAtGaussPts();
251  }
252 
253  /** \brief if higher order geometry return normals at Gauss pts.
254  *
255  * \param gg gauss point number
256  */
257  inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPts(const int gg) {
258  return ublas::matrix_row<MatrixDouble>(
259  static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)
260  ->normalsAtGaussPts,
261  gg);
262  }
263 
264  /// \deprecated use getNormalsAtGaussPts
265  DEPRECATED inline auto getNormalsAtGaussPt(const int gg) {
266  return getNormalsAtGaussPts(gg);
267  }
268 
269  /** \brief if higher order geometry return tangent vector to triangle at
270  Gauss pts.
271 
272  Note: returned matrix has size 0 in rows and columns if no HO approximation
273  of geometry is avaliable.
274 
275  */
277  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)
278  ->tangentOneAtGaussPts;
279  }
280 
281  /// \deprecated use getTangent1AtGaussPts
283  return getTangent1AtGaussPts();
284  }
285 
286  /** \brief if higher order geometry return tangent vector to triangle at
287  Gauss pts.
288 
289  Note: returned matrix has size 0 in rows and columns if no HO approximation
290  of geometry is avaliable.
291 
292  */
294  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE)
295  ->tangentTwoAtGaussPts;
296  }
297 
298  /// \deprecated use getTangent2AtGaussPts
300  return getTangent2AtGaussPts();
301  }
302 
303  /** \brief get normal at integration points
304 
305  Example:
306  \code
307  double nrm2;
308  FTensor::Index<'i',3> i;
309  auto t_normal = getFTensor1NormalsAtGaussPts();
310  for(int gg = gg!=data.getN().size1();gg++) {
311  nrm2 = sqrt(t_normal(i)*t_normal(i));
312  ++t_normal;
313  }
314  \endcode
315 
316  */
318  double *ptr = &*getNormalsAtGaussPts().data().begin();
319  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
320  &ptr[2]);
321  }
322 
323  /// \deprecated use getFTensor1NormalsAtGaussPt
325  return getFTensor1NormalsAtGaussPts();
326  }
327 
328  /** \brief get tangent 1 at integration points
329 
330  */
332  double *ptr = &*getTangent1AtGaussPts().data().begin();
333  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
334  &ptr[2]);
335  }
336 
337  /// \deprecated use getFTensor1Tangent1AtGaussPt
339  return getFTensor1Tangent1AtGaussPts();
340  }
341 
342  /** \brief get tangent 2 at integration points
343 
344  */
346  double *ptr = &*getTangent2AtGaussPts().data().begin();
347  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
348  &ptr[2]);
349  }
350 
351  /// \deprecated use getFTensor1Tangent2AtGaussPt
353  return getFTensor1Tangent2AtGaussPts();
354  }
355 
356  /** \deprecated use getFaceFE instead
357  */
360  return getFaceFE();
361  }
362 
363  /** \deprecated use getFaceFE instead
364  */
366  return getFaceFE();
367  }
368 
369  /** \brief return pointer to Generic Triangle Finite Element object
370  */
372  return static_cast<FaceElementForcesAndSourcesCore *>(ptrFE);
373  }
374 
375  /**
376  *
377  * User call this function to loop over elements on the side of face. This
378  * function calls MoFEM::VolumeElementForcesAndSourcesCoreOnSide with is
379  * operator to do calculations.
380  *
381  * @param fe_name Name of the element
382  * @param method Finite element object
383  * @return error code
384  */
386  loopSideVolumes(const string &fe_name,
388  };
389 
390  /**
391  * \brief Calculate element area and normal of the face
392  *
393  * Note that at that point is assumed that geometry is exclusively defined by
394  * corner nodes.
395  *
396  * @return Error code
397  */
398  virtual MoFEMErrorCode calculateAreaAndNormal();
399 
400  int nbGaussPts; ///< Number of integration points
401 
402  /**
403  * \brief Set integration points
404  * @return Error code
405  */
406  virtual MoFEMErrorCode setIntegrationPts();
407 
408  /// \deprecated method with spelling mistake, use setIntegrationPts
409  DEPRECATED MoFEMErrorCode setIntegartionPts() { return setIntegrationPts(); }
410 
411  /**
412  * \brief Determine approximation space and order of base functions
413  * @return Error code
414  */
415  virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement();
416 
417  /**
418  * \brief Calculate coordinate at integration points
419  * @return Error code
420  */
421  virtual MoFEMErrorCode calculateCoordinatesAtGaussPts();
422 
423  /**
424  * \brief Calculate normal on curved elements
425  *
426  * Geometry is given by other field.
427  *
428  * @return error code
429  */
430  virtual MoFEMErrorCode calculateHoNormal();
431 
432  MoFEMErrorCode operator()();
433 };
434 
435 /** \brief Calculate inverse of jacobian for face element
436 
437  It is assumed that face element is XY plane. Applied
438  only for 2d problems.
439 
440  \todo Generalize function for arbitrary face orientation in 3d space
441 
442  \ingroup mofem_forces_and_sources_tri_element
443 
444 */
448 
449  // /**
450  // * \deprecated Field name do not needed to construct class, change v0.5.17.
451  // */
452  // DEPRECATED OpCalculateInvJacForFace(const std::string
453  // &field_name,MatrixDouble &inv_jac):
454  // FaceElementForcesAndSourcesCore::UserDataOperator(H1),
455  // invJac(inv_jac) {}
456 
459  }
460 
461  MoFEMErrorCode doWork(int side, EntityType type,
463 };
464 
465 /** \brief Transform local reference derivatives of shape functions to global
466 derivatives
467 
468 \ingroup mofem_forces_and_sources_tri_element
469 
470 */
473 
476  }
477 
478  MoFEMErrorCode doWork(int side, EntityType type,
480 
481 private:
484 };
485 
486 /**
487  * \brief brief Transform local reference derivatives of shape function to
488  global derivatives for face
489 
490  * \ingroup mofem_forces_and_sources_tri_element
491  */
494 
497  invJac(inv_jac) {}
498 
499  MoFEMErrorCode doWork(int side, EntityType type,
501 
502 private:
505 };
506 
507 } // namespace MoFEM
508 
509 #endif //__FACEELEMENTFORCESANDSOURCESCORE_HPP__
510 
511 /**
512  * \defgroup mofem_forces_and_sources_tri_element Face Element
513  * \brief Implementation of face element
514  *
515  * \ingroup mofem_forces_and_sources
516  **/
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
DEPRECATED const FaceElementForcesAndSourcesCore * getFaceElementForcesAndSourcesCore()
transform Hdiv base fluxes from reference element to physical triangle
MatrixDouble & getHoCoordsAtGaussPts()
coordinate at Gauss points (if hierarchical approximation of element geometry)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MatrixDouble & getTangent2AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
ublas::matrix_row< MatrixDouble > getNormalsAtGaussPts(const int gg)
if higher order geometry return normals at Gauss pts.
MatrixDouble & getTangent1AtGaussPts()
if higher order geometry return tangent vector to triangle at Gauss pts.
UserDataOperator(const std::string &row_field_name, const std::string &col_field_name, const char type, const bool symm=true)
OpSetCovariantPiolaTransformOnTriangle opCovariantTransform
Calculate inverse of jacobian for face element.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
std::string meshPositionsFieldName
Name of the field with geometry.
field with continuous tangents
Definition: definitions.h:171
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
FieldSpace
approximation spaces
Definition: definitions.h:167
DEPRECATED const FaceElementForcesAndSourcesCore * getTriFE()
Transform local reference derivatives of shape functions to global derivatives.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define DEPRECATED
Definition: definitions.h:29
brief Transform local reference derivatives of shape function to global derivatives for face
OpSetContravariantPiolaTransformOnTriangle opContravariantTransform
auto getFTensor1HoCoordsAtGaussPts()
get coordinates at Gauss pts (takes in account ho approx. of geometry)
Data on single entity (This is passed as argument to DataOperator::doWork)
structure to get information form mofem into DataForcesAndSourcesCore
Calculate normals at Gauss points of triangle element.
UserDataOperator(const std::string &field_name, const char type)
transform Hcurl base fluxes from reference element to physical triangle
continuous field
Definition: definitions.h:170
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:76
const FaceElementForcesAndSourcesCore * getFaceFE()
return pointer to Generic Triangle Finite Element object
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...