v0.9.0
FatPrismElementForcesAndSourcesCore.hpp
Go to the documentation of this file.
1 /** \file ForcesAndSourcesCore.hpp
2 
3  \brief Implementation of elements on entities.
4 
5  Those element are inherited by user to implement specific implementation of
6  particular problem.
7 
8 */
9 
10 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #ifndef __FATPRISMELEMENTFORCESANDSURCESCORE_HPP__
25 #define __FATPRISMELEMENTFORCESANDSURCESCORE_HPP__
26 
27 using namespace boost::numeric;
28 
29 namespace MoFEM {
30 
31 /** \brief FatPrism finite element
32  \ingroup mofem_forces_and_sources_prism_element
33 
34  User is implementing own operator at Gauss points level, by own object
35  derived from FatPrismElementForcesAndSourcesCoreL::UserDataOperator. Arbitrary
36  number of operator added pushing objects to rowOpPtrVector and
37  rowColOpPtrVector.
38 
39  \todo Need to implement operators that will make this element work as Volume
40  element
41 
42  */
45 
47 
48  virtual int getRuleTrianglesOnly(int order) { return 2 * order; };
49  virtual int getRuleThroughThickness(int order) { return 2 * order; };
50 
51  virtual MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only) {
53  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
55  }
56 
57  virtual MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness) {
59  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
61  }
62 
63  /** \brief default operator for Flat Prism element
64  * \ingroup mofem_forces_and_sources_prism_element
65  */
68 
70 
71  /** \brief get face aRea
72  \param dd if dd == 0 it is for face F3 if dd == 1 is for face F4
73  */
74  inline double getArea(const int dd);
75 
76  inline double getAreaF3();
77  inline double getAreaF4();
78 
79  /** \brief get triangle normal
80 
81  Normal has 6 elements, first 3 are for face F3 another three for face F4
82 
83  */
84  inline VectorDouble &getNormal();
85 
86  inline VectorAdaptor getNormalF3();
87 
88  inline VectorAdaptor getNormalF4();
89 
90  /** \brief get Gauss pts. in the prism
91  */
92  inline MatrixDouble &getGaussPts();
93 
94  /** \brief get Gauss pts. on triangles
95  */
96  inline MatrixDouble &getGaussPtsTrianglesOnly();
97 
98  /** \brief get Gauss pts. through thickness
99  */
100  inline MatrixDouble &getGaussPtsThroughThickness();
101 
102  /** \brief get coordinates at Gauss pts.
103 
104  Matrix has size (nb integration points)x(coordinates on F3 and F4 = 6),
105  i.e. coordinates on face F3 and F4
106 
107  */
108  inline MatrixDouble &getCoordsAtGaussPts();
109 
110  /** \brief get coordinates at Gauss pts.
111 
112  Matrix has size (nb integration points)x(coordinates on F3 and F4 = 6),
113  i.e. coordinates on face F3 and F4
114 
115  */
116  inline MatrixDouble &getCoordsAtGaussPtsTrianglesOnly();
117 
118  /** \brief coordinate at Gauss points on face 3 (if hierarchical
119  * approximation of element geometry)
120  */
121  inline MatrixDouble &getHoCoordsAtGaussPtsF3();
122 
123  /** \brief coordinate at Gauss points on face 4 (if hierarchical
124  * approximation of element geometry)
125  */
126  inline MatrixDouble &getHoCoordsAtGaussPtsF4();
127 
128  /** \brief if higher order geometry return normals at face F3 at Gauss pts.
129  *
130  * Face 3 is top face in canonical triangle numeration, see \cite
131  * tautges2010canonical
132  *
133  */
134  inline MatrixDouble &getNormalsAtGaussPtF3();
135 
136  /** \brief if higher order geometry return normals at face F4 at Gauss pts.
137  *
138  * Face 4 is top face in canonical triangle numeration, see \cite
139  * tautges2010canonical
140  *
141  */
142  inline MatrixDouble &getNormalsAtGaussPtF4();
143 
144  /** \brief if higher order geometry return normals at Gauss pts.
145  *
146  * Face 3 is top face in canonical triangle numeration, see \cite
147  * tautges2010canonical
148  *
149  * \param gg gauss point number
150  */
151  inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPtF3(const int gg);
152 
153  /** \brief if higher order geometry return normals at Gauss pts.
154  *
155  * Face 3 is top face in canonical triangle numeration, see \cite
156  * tautges2010canonical
157  *
158  * \param gg gauss point number
159  */
160  inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPtF4(const int gg);
161 
162  /** \brief if higher order geometry return tangent vector to triangle at
163  * Gauss pts.
164  */
165  inline MatrixDouble &getTangent1AtGaussPtF3();
166 
167  /** \brief if higher order geometry return tangent vector to triangle at
168  * Gauss pts.
169  */
170  inline MatrixDouble &getTangent2AtGaussPtF3();
171 
172  /** \brief if higher order geometry return tangent vector to triangle at
173  * Gauss pts.
174  */
175  inline MatrixDouble &getTangent1AtGaussPtF4();
176 
177  /** \brief if higher order geometry return tangent vector to triangle at
178  * Gauss pts.
179  */
180  inline MatrixDouble &getTangent2AtGaussPtF4();
181 
182  inline DataForcesAndSourcesCore &getTrianglesOnlyDataStructure();
183 
184  inline DataForcesAndSourcesCore &getTroughThicknessDataStructure();
185 
186  /** \brief return pointer to fat prism finite element
187  */
188  inline const FatPrismElementForcesAndSourcesCore *getPrismFE();
189  };
190 
191  MoFEMErrorCode operator()();
192 
193 protected:
194  double aRea[2];
196 
200 
203 
213 
214  friend class UserDataOperator;
215 };
216 
217 inline double
218 FatPrismElementForcesAndSourcesCore::UserDataOperator::getArea(const int dd) {
219  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)->aRea[0];
220 }
221 
222 inline double
223 FatPrismElementForcesAndSourcesCore::UserDataOperator::getAreaF3() {
224  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)->aRea[0];
225 }
226 inline double
227 FatPrismElementForcesAndSourcesCore::UserDataOperator::getAreaF4() {
228  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)->aRea[1];
229 }
230 
231 inline VectorDouble &
232 FatPrismElementForcesAndSourcesCore::UserDataOperator::getNormal() {
233  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)->normal;
234 }
235 
236 inline VectorAdaptor
237 FatPrismElementForcesAndSourcesCore::UserDataOperator::getNormalF3() {
238  double *data =
239  &(static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)->normal[0]);
240  return VectorAdaptor(3, ublas::shallow_array_adaptor<double>(3, data));
241 }
242 
243 inline VectorAdaptor
244 FatPrismElementForcesAndSourcesCore::UserDataOperator::getNormalF4() {
245  double *data =
246  &(static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)->normal[3]);
247  return VectorAdaptor(3, ublas::shallow_array_adaptor<double>(3, data));
248 }
249 
250 inline MatrixDouble &
251 FatPrismElementForcesAndSourcesCore::UserDataOperator::getGaussPts() {
252  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)->gaussPts;
253 }
254 
255 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
256  getGaussPtsTrianglesOnly() {
257  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
258  ->gaussPtsTrianglesOnly;
259 }
260 
261 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
262  getGaussPtsThroughThickness() {
263  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
264  ->gaussPtsThroughThickness;
265 }
266 
267 inline MatrixDouble &
268 FatPrismElementForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts() {
269  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
270  ->coordsAtGaussPts;
271 }
272 
273 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
274  getCoordsAtGaussPtsTrianglesOnly() {
275  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
276  ->coordsAtGaussPtsTrianglesOnly;
277 }
278 
279 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
280  getHoCoordsAtGaussPtsF3() {
281  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
282  ->hoCoordsAtGaussPtsF3;
283 }
284 
285 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
286  getHoCoordsAtGaussPtsF4() {
287  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
288  ->hoCoordsAtGaussPtsF4;
289 }
290 
291 inline MatrixDouble &
292 FatPrismElementForcesAndSourcesCore::UserDataOperator::getNormalsAtGaussPtF3() {
293  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
294  ->nOrmals_at_GaussPtF3;
295 }
296 
297 inline MatrixDouble &
298 FatPrismElementForcesAndSourcesCore::UserDataOperator::getNormalsAtGaussPtF4() {
299  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
300  ->nOrmals_at_GaussPtF4;
301 }
302 
303 inline ublas::matrix_row<MatrixDouble>
304 FatPrismElementForcesAndSourcesCore::UserDataOperator::getNormalsAtGaussPtF3(
305  const int gg) {
306  return ublas::matrix_row<MatrixDouble>(
307  static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
308  ->nOrmals_at_GaussPtF3,
309  gg);
310 }
311 
312 inline ublas::matrix_row<MatrixDouble>
313 FatPrismElementForcesAndSourcesCore::UserDataOperator::getNormalsAtGaussPtF4(
314  const int gg) {
315  return ublas::matrix_row<MatrixDouble>(
316  static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
317  ->nOrmals_at_GaussPtF4,
318  gg);
319 }
320 
321 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
322  getTangent1AtGaussPtF3() {
323  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
324  ->tAngent1_at_GaussPtF3;
325 }
326 
327 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
328  getTangent2AtGaussPtF3() {
329  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
330  ->tAngent2_at_GaussPtF3;
331 }
332 
333 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
334  getTangent1AtGaussPtF4() {
335  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
336  ->tAngent1_at_GaussPtF4;
337 }
338 
339 inline MatrixDouble &FatPrismElementForcesAndSourcesCore::UserDataOperator::
340  getTangent2AtGaussPtF4() {
341  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
342  ->tAngent2_at_GaussPtF4;
343 }
344 
345 inline DataForcesAndSourcesCore &FatPrismElementForcesAndSourcesCore::
346  UserDataOperator::getTrianglesOnlyDataStructure() {
347  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
348  ->dataH1TrianglesOnly;
349 }
350 
351 inline DataForcesAndSourcesCore &FatPrismElementForcesAndSourcesCore::
352  UserDataOperator::getTroughThicknessDataStructure() {
353  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE)
354  ->dataH1TroughThickness;
355 }
356 
358 FatPrismElementForcesAndSourcesCore::UserDataOperator::getPrismFE() {
359  return static_cast<FatPrismElementForcesAndSourcesCore *>(ptrFE);
360 }
361 
362 } // namespace MoFEM
363 
364 #endif //__FATPRISMELEMENTFORCESANDSURCESCORE_HPP__
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#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
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:101
calculate normals at Gauss points of triangle element
virtual MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
virtual MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
FatPrism finite elementUser is implementing own operator at Gauss points level, by own object derived...
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72