v0.12.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
40 
41  /** \brief default operator for TRI element
42  * \ingroup mofem_forces_and_sources_tri_element
43  */
44  struct UserDataOperator;
45 
46  enum Switches {
47  };
48 
49  template <int SWITCH> MoFEMErrorCode opSwitch();
50 
51 protected:
53 
54  /**
55  * \brief Calculate element area and normal of the face at integration points
56  *
57  * @return Error code
58  */
59  virtual MoFEMErrorCode calculateAreaAndNormalAtIntegrationPts();
60 
61  /**
62  * \brief Calculate element area and normal of the face
63  *
64  * Note that at that point is assumed that geometry is exclusively defined
65  * by corner nodes.
66  *
67  * @return Error code
68  */
69  virtual MoFEMErrorCode calculateAreaAndNormal();
70 
71  /**
72  * \brief Set integration points
73  * @return Error code
74  */
75  virtual MoFEMErrorCode setIntegrationPts();
76 
77  /**
78  * \brief Determine approximation space and order of base functions
79  * @return Error code
80  */
81  virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement();
82 
83  /**
84  * \brief Calculate coordinate at integration points
85  * @return Error code
86  */
87  virtual MoFEMErrorCode calculateCoordinatesAtGaussPts();
88 
89  double aRea;
90  int num_nodes;
92  VectorDouble nOrmal, tangentOne, tangentTwo;
94 
98 
99  friend class UserDataOperator;
101 };
102 
103 /** \brief default operator for TRI element
104  * \ingroup mofem_forces_and_sources_tri_element
105  */
108 
110 
111  /**
112  * \brief get area of face
113  * @return area of face
114  */
115  inline double getArea();
116 
117  /**
118  * \brief get measure of element
119  * @return area of face
120  */
121  inline double getMeasure();
122 
123  /** \brief get triangle normal
124  */
125  inline VectorDouble &getNormal();
126 
127  /** \brief get triangle tangent 1
128  */
129  inline VectorDouble &getTangent1();
130 
131  /** \brief get triangle tangent 2
132  */
133  inline VectorDouble &getTangent2();
134 
135  /** \brief get normal as tensor
136  */
137  inline auto getFTensor1Normal();
138 
139  /** \brief get tangentOne as tensor
140  */
141  inline auto getFTensor1Tangent1();
142 
143  /** \brief get tangentTwo as tensor
144  */
145  inline auto getFTensor1Tangent2();
146 
147  /** \brief get element number of nodes
148  */
149  inline int getNumNodes();
150 
151  /** \brief get element connectivity
152  */
153  inline const EntityHandle *getConn();
154 
155  /** \brief get triangle coordinates
156  */
157  inline VectorDouble &getCoords();
158 
159  /**
160  * \brief get get coords at gauss points
161 
162  \code
163  FTensor::Index<'i',3> i;
164  FTensor::Tensor1<double,3> t_center;
165  auto t_coords = getFTensor1Coords();
166  t_center(i) = 0;
167  for(int nn = 0;nn!=3;nn++) {
168  t_center(i) += t_coords(i);
169  ++t_coords;
170  }
171  t_center(i) /= 3;
172  \endcode
173 
174  */
175  inline auto getFTensor1Coords();
176 
177  /** \brief if higher order geometry return normals at Gauss pts.
178 
179  Note: returned matrix has size 0 in rows and columns if no HO approximation
180  of geometry is available.
181 
182  */
183  inline MatrixDouble &getNormalsAtGaussPts();
184 
185  /** \brief if higher order geometry return normals at Gauss pts.
186  *
187  * \param gg gauss point number
188  */
189  inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPts(const int gg);
190 
191  /** \brief if higher order geometry return tangent vector to triangle at
192  Gauss pts.
193 
194  Note: returned matrix has size 0 in rows and columns if no HO approximation
195  of geometry is avaliable.
196 
197  */
198  inline MatrixDouble &getTangent1AtGaussPts();
199 
200  /** \brief if higher order geometry return tangent vector to triangle at
201  Gauss pts.
202 
203  Note: returned matrix has size 0 in rows and columns if no HO approximation
204  of geometry is avaliable.
205 
206  */
207  inline MatrixDouble &getTangent2AtGaussPts();
208 
209  /** \brief get normal at integration points
210 
211  Example:
212  \code
213  double nrm2;
214  FTensor::Index<'i',3> i;
215  auto t_normal = getFTensor1NormalsAtGaussPts();
216  for(int gg = gg!=data.getN().size1();gg++) {
217  nrm2 = sqrt(t_normal(i)*t_normal(i));
218  ++t_normal;
219  }
220  \endcode
221 
222  */
223  inline auto getFTensor1NormalsAtGaussPts();
224 
225  /** \brief get tangent 1 at integration points
226 
227  */
228  inline auto getFTensor1Tangent1AtGaussPts();
229 
230  /** \brief get tangent 2 at integration points
231 
232  */
233  inline auto getFTensor1Tangent2AtGaussPts();
234 
235  /** \brief return pointer to Generic Triangle Finite Element object
236  */
237  inline FaceElementForcesAndSourcesCoreBase *getFaceFE();
238 
239 
240  /**
241  *
242  * User call this function to loop over elements on the side of face. This
243  * function calls MoFEM::VolumeElementForcesAndSourcesCoreOnSide with is
244  * operator to do calculations.
245  *
246  * @param fe_name Name of the element
247  * @param method Finite element object
248  * @return error code
249  */
250  template <int SWITCH>
251  MoFEMErrorCode loopSideVolumes(
252  const string &fe_name,
254 
255 private:
256  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
257 };
258 
259 /** \brief Face finite element switched
260  \ingroup mofem_forces_and_sources_tri_element
261 
262  */
263 template <int SWITCH>
266 
269 
272 
273  MoFEMErrorCode operator()();
274 };
275 
276 /** \brief Face finite element default
277  \ingroup mofem_forces_and_sources_tri_element
278 
279  */
282 
283 template <int SWITCH>
284 MoFEMErrorCode FaceElementForcesAndSourcesCoreBase::opSwitch() {
286 
287  const auto type = numeredEntFiniteElementPtr->getEntType();
288  if (type != lastEvaluatedElementEntityType) {
289  switch (type) {
290  case MBTRI:
291  getElementPolynomialBase() =
292  boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
293  break;
294  case MBQUAD:
295  getElementPolynomialBase() =
296  boost::shared_ptr<BaseFunction>(new QuadPolynomialBase());
297  break;
298  default:
300  }
301  CHKERR createDataOnElement();
302  }
303 
304  // Calculate normal and tangent vectors for face geometry
305  CHKERR calculateAreaAndNormal();
306  CHKERR getSpaceBaseAndOrderOnElement();
307 
308  CHKERR setIntegrationPts();
309  if (gaussPts.size2() == 0)
311 
312  CHKERR calculateCoordinatesAtGaussPts();
313  CHKERR calHierarchicalBaseFunctionsOnElement();
314  CHKERR calBernsteinBezierBaseFunctionsOnElement();
315  CHKERR calculateAreaAndNormalAtIntegrationPts();
316 
317  // Iterate over operators
318  CHKERR loopOverOperators();
319 
321 }
322 
323 template <int SWITCH>
325  return opSwitch<SWITCH>();
326 }
327 
328 double FaceElementForcesAndSourcesCoreBase::UserDataOperator::getArea() {
329  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->aRea;
330 }
331 
332 double FaceElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() {
333  return getArea();
334 }
335 
336 VectorDouble &
337 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormal() {
338  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->nOrmal;
339 }
340 
341 VectorDouble &
342 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent1() {
343  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->tangentOne;
344 }
345 
346 VectorDouble &
347 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent2() {
348  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->tangentTwo;
349 }
350 
351 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
352  getFTensor1Normal() {
353  double *ptr = &*getNormal().data().begin();
354  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
355 }
356 
357 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
358  getFTensor1Tangent1() {
359  double *ptr = &*getTangent1().data().begin();
360  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
361 }
362 
363 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
364  getFTensor1Tangent2() {
365  double *ptr = &*getTangent2().data().begin();
366  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
367 }
368 
369 int FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNumNodes() {
370  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->num_nodes;
371 }
372 
373 const EntityHandle *
374 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getConn() {
375  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->conn;
376 }
377 
378 VectorDouble &
379 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getCoords() {
380  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)->coords;
381 }
382 
383 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
384  getFTensor1Coords() {
385  double *ptr = &*getCoords().data().begin();
386  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
387  &ptr[2]);
388 }
389 
390 MatrixDouble &
391 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPts() {
392  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
393  ->normalsAtGaussPts;
394 }
395 
396 ublas::matrix_row<MatrixDouble>
397 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPts(
398  const int gg) {
399  return ublas::matrix_row<MatrixDouble>(
400  static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
401  ->normalsAtGaussPts,
402  gg);
403 }
404 
405 MatrixDouble &
406 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent1AtGaussPts() {
407  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
408  ->tangentOneAtGaussPts;
409 }
410 
411 MatrixDouble &
412 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getTangent2AtGaussPts() {
413  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE)
414  ->tangentTwoAtGaussPts;
415 }
416 
417 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
418  getFTensor1NormalsAtGaussPts() {
419  double *ptr = &*getNormalsAtGaussPts().data().begin();
420  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
421  &ptr[2]);
422 }
423 
424 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
425  getFTensor1Tangent1AtGaussPts() {
426  double *ptr = &*getTangent1AtGaussPts().data().begin();
427  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
428  &ptr[2]);
429 }
430 
431 auto FaceElementForcesAndSourcesCoreBase::UserDataOperator::
432  getFTensor1Tangent2AtGaussPts() {
433  double *ptr = &*getTangent2AtGaussPts().data().begin();
434  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
435  &ptr[2]);
436 }
437 
439 FaceElementForcesAndSourcesCoreBase::UserDataOperator::getFaceFE() {
440  return static_cast<FaceElementForcesAndSourcesCoreBase *>(ptrFE);
441 }
442 
443 
444 
445 template <int SWITCH>
447 FaceElementForcesAndSourcesCoreBase::UserDataOperator::loopSideVolumes(
448  const string &fe_name,
450  return loopSide(fe_name, &fe_method, 3);
451 }
452 
453 } // namespace MoFEM
454 
455 #endif //__FACEELEMENTFORCESANDSOURCESCORE_HPP__
456 
457 /**
458  * \defgroup mofem_forces_and_sources_tri_element Face Element
459  * \brief Implementation of face element
460  *
461  * \ingroup mofem_forces_and_sources
462  **/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
T data[Tensor_Dim]
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Deprecated interface functions.
std::string meshPositionsFieldName
Name of the field with geometry.
structure to get information form mofem into DataForcesAndSourcesCore
Calculate base functions on triangle.
Calculate base functions on triangle.