v0.9.1
VolumeElementForcesAndSourcesCoreOnSide.hpp
Go to the documentation of this file.
1 /** \file VolumeElementForcesAndSourcesCoreOnSide.hpp
2  \brief Volume element.
3 
4  Those element are inherited by user to implement specific implementation of
5  particular problem.
6 
7 */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef __VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
24 #define __VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
25 
26 using namespace boost::numeric;
27 
28 namespace MoFEM {
29 
30 /**
31  * \brief Base volume element used to integrate on skeleton
32  * \ingroup mofem_forces_and_sources_volume_element
33  */
36 
38 
39  /**
40  * @brief Get the face nodes mapped on volume element
41  *
42  * \todo That this is not general, e.g., for quad number of nodes is 4.
43  *
44  * @return const std::array<int, 3>&
45  */
46  inline const std::array<int, 3> &getFaceConnMap() const;
47 
48  /**
49  * @brief Get face nodes maped on volume
50  *
51  * \todo That this is not general, e.g., for prism or hex, size of fixed array
52  * is wrong.
53  *
54  * @return const sdt::array<int, 4>&
55  */
56  inline const std::array<int, 4> &getTetConnMap() const;
57 
58  /**
59  * @brief Get node on volume opposite to volume element
60  *
61  * \todo That this is not general, e.g., for prism or hex, opoosite node is
62  * not unique.
63  *
64  * @return int
65  */
66  inline int getOppositeNode() const;
67 
68  /**
69  * @brief Sense face on volume
70  *
71  * @return int
72  */
73  inline int getFaceSense() const;
74 
75  /**
76  * @brief Face number on the volume
77  *
78  * @return int
79  */
80  inline int getFaceSideNumber() const;
81 
82  /** \brief default operator for TET element
83  * \ingroup mofem_forces_and_sources_volume_element
84  */
87 
89 
90  inline VolumeElementForcesAndSourcesCoreOnSideBase *getVolumeFE() const;
91 
92  inline FaceElementForcesAndSourcesCoreBase *getFaceFE() const;
93 
94  /**
95  * \brief get face sense in respect to volume
96  * @return error code
97  */
98  inline int getFaceSense() const;
99 
100  /**
101  * \brief get face side number in respect to volume
102  * @return error code
103  */
104  inline int getFaceSideNumber() const;
105 
106  inline bool getEdgeFace(const int ee) const;
107 
108  /**
109  * get face normal on side which is this element
110  * @return face normal
111  */
112  inline VectorDouble &getNormal();
113 
114  /** \brief get triangle tangent 1
115  */
116  inline VectorDouble &getTangent1();
117 
118  /** \brief get triangle tangent 2
119  */
120  inline VectorDouble &getTangent2();
121 
122  /** \brief get normal as tensor
123  */
124  inline auto getFTensor1Normal();
125 
126  /** \brief get tangentOne as tensor
127  */
128  inline auto getFTensor1Tangent1();
129 
130  /** \brief get tangentTwo as tensor
131  */
132  inline auto getFTensor1Tangent2();
133 
134  /** \brief if higher order geometry return normals at Gauss pts.
135 
136  Note: returned matrix has size 0 in rows and columns if no HO approximation
137  of geometry is available.
138 
139  */
140  inline MatrixDouble &getNormalsAtGaussPts();
141 
142  /** \brief if higher order geometry return normals at Gauss pts.
143  *
144  * \param gg gauss point number
145  */
146  inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPts(const int gg);
147 
148  /** \brief get normal at integration points
149 
150  Example:
151  \code
152  double nrm2;
153  FTensor::Index<'i',3> i;
154  auto t_normal = getFTensor1NormalsAtGaussPts();
155  for(int gg = gg!=data.getN().size1();gg++) {
156  nrm2 = sqrt(t_normal(i)*t_normal(i));
157  ++t_normal;
158  }
159  \endcode
160 
161  */
162  inline auto getFTensor1NormalsAtGaussPts();
163 
164  /** \brief get face coordinates at Gauss pts.
165 
166  \note Coordinates should be the same what function getCoordsAtGaussPts
167  on tets is returning. If both coordinates are different it is error, or you
168  do something very unusual.
169 
170  */
171  inline MatrixDouble &getFaceCoordsAtGaussPts();
172 
173  protected:
174 
175  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
176  };
177 
178  int getRule(int order);
179  MoFEMErrorCode setGaussPts(int order);
180 
181 private:
182  int faceSense; ///< Sense of face, could be 1 or -1
183  int faceSideNumber; ///< Face side number
184  std::array<int, 3> faceConnMap;
185  std::array<int, 4> tetConnMap;
187 };
188 
189 /**
190  * @brief Volume side finite element with switches
191  *
192  * Using SWITCH to off functions
193  *
194  * @tparam SWITCH
195  */
196 template <int SWITCH>
199 
200  using VolumeElementForcesAndSourcesCoreOnSideBase::
201  VolumeElementForcesAndSourcesCoreOnSideBase;
202 
203  using UserDataOperator =
205 
206  MoFEMErrorCode operator()();
207 };
208 
209 /** \brief Volume element used to integrate on skeleton
210  \ingroup mofem_forces_and_sources_volume_element
211 
212  */
215 
216 const std::array<int, 3> &
217 VolumeElementForcesAndSourcesCoreOnSideBase::getFaceConnMap() const {
218  return faceConnMap;
219 }
220 
221 const std::array<int, 4> &
222 VolumeElementForcesAndSourcesCoreOnSideBase::getTetConnMap() const {
223  return tetConnMap;
224 }
225 
226 int VolumeElementForcesAndSourcesCoreOnSideBase::getOppositeNode() const {
227  return oppositeNode;
228 }
229 
230 int VolumeElementForcesAndSourcesCoreOnSideBase::getFaceSense() const {
231  return faceSense;
232 }
233 
234 int VolumeElementForcesAndSourcesCoreOnSideBase::getFaceSideNumber() const {
235  return faceSideNumber;
236 }
237 
239 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getFaceFE()
240  const {
241  return static_cast<FaceElementForcesAndSourcesCoreBase *>(
242  getVolumeFE()->sidePtrFE);
243 }
244 
245 VectorDouble &
246 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getNormal() {
247  return getFaceFE()->nOrmal;
248 }
249 
250 VectorDouble &
251 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getTangent1() {
252  return getFaceFE()->tangentOne;
253 }
254 
255 VectorDouble &
256 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getTangent2() {
257  return getFaceFE()->tangentTwo;
258 }
259 
261 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getVolumeFE()
262  const {
263  return static_cast<VolumeElementForcesAndSourcesCoreOnSideBase *>(ptrFE);
264 }
265 
266 int VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
267  getFaceSense() const {
268  return getVolumeFE()->faceSense;
269 }
270 
271 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
272  getFTensor1Normal() {
273  double *ptr = &*getNormal().data().begin();
274  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
275 }
276 
277 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
278  getFTensor1Tangent1() {
279  double *ptr = &*getTangent1().data().begin();
280  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
281 }
282 
283 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
284  getFTensor1Tangent2() {
285  double *ptr = &*getTangent2().data().begin();
286  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
287 }
288 
289 inline auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
290  getFTensor1NormalsAtGaussPts() {
291  double *ptr = &*getNormalsAtGaussPts().data().begin();
292  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
293  &ptr[2]);
294 }
295 
296 int VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
297  getFaceSideNumber() const {
298  return getVolumeFE()->faceSideNumber;
299 }
300 
301 MatrixDouble &VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
302  getNormalsAtGaussPts() {
303  return getFaceFE()->normalsAtGaussPts;
304 }
305 
306 MatrixDouble &VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
307  getFaceCoordsAtGaussPts() {
308  return getFaceFE()->coordsAtGaussPts;
309 }
310 
311 bool VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getEdgeFace(
312  const int ee) const {
313  constexpr bool edges_on_faces[6][4] = {{true, false, false, true}, // e0
314  {false, true, false, true}, // e1
315  {false, false, true, true}, // e2
316  {true, false, true, false}, // e3
317  {true, true, false, false}, // e4
318  {false, true, true, false}};
319  return edges_on_faces[ee][getFaceSideNumber()];
320 }
321 
322 ublas::matrix_row<MatrixDouble> VolumeElementForcesAndSourcesCoreOnSideBase::
323  UserDataOperator::getNormalsAtGaussPts(const int gg) {
324  return ublas::matrix_row<MatrixDouble>(getNormalsAtGaussPts(), gg);
325 }
326 
327 template <int SWITCH>
330  return OpSwitch<SWITCH>();
331 }
332 
333 } // namespace MoFEM
334 
335 #endif //__VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
constexpr int order
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
ForcesAndSourcesCore::UserDataOperator UserDataOperator
structure to get information form mofem into DataForcesAndSourcesCore
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...