v0.9.0
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 
174  int getRule(int order);
175  MoFEMErrorCode setGaussPts(int order);
176 
177 private:
178  int faceSense; ///< Sense of face, could be 1 or -1
179  int faceSideNumber; ///< Face side number
180  std::array<int, 3> faceConnMap;
181  std::array<int, 4> tetConnMap;
183 };
184 
185 /**
186  * @brief Volume side finite element with switches
187  *
188  * Using SWITCH to off functions
189  *
190  * @tparam SWITCH
191  */
192 template <int SWITCH>
195 
196  using VolumeElementForcesAndSourcesCoreOnSideBase::
197  VolumeElementForcesAndSourcesCoreOnSideBase;
198 
199  using UserDataOperator =
201 
202  MoFEMErrorCode operator()();
203 };
204 
205 /** \brief Volume element used to integrate on skeleton
206  \ingroup mofem_forces_and_sources_volume_element
207 
208  */
211 
212 const std::array<int, 3> &
213 VolumeElementForcesAndSourcesCoreOnSideBase::getFaceConnMap() const {
214  return faceConnMap;
215 }
216 
217 const std::array<int, 4> &
218 VolumeElementForcesAndSourcesCoreOnSideBase::getTetConnMap() const {
219  return tetConnMap;
220 }
221 
222 int VolumeElementForcesAndSourcesCoreOnSideBase::getOppositeNode() const {
223  return oppositeNode;
224 }
225 
226 int VolumeElementForcesAndSourcesCoreOnSideBase::getFaceSense() const {
227  return faceSense;
228 }
229 
230 int VolumeElementForcesAndSourcesCoreOnSideBase::getFaceSideNumber() const {
231  return faceSideNumber;
232 }
233 
235 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getFaceFE()
236  const {
237  return static_cast<FaceElementForcesAndSourcesCoreBase *>(
238  getVolumeFE()->sidePtrFE);
239 }
240 
241 VectorDouble &
242 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getNormal() {
243  return getFaceFE()->nOrmal;
244 }
245 
246 VectorDouble &
247 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getTangent1() {
248  return getFaceFE()->tangentOne;
249 }
250 
251 VectorDouble &
252 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getTangent2() {
253  return getFaceFE()->tangentTwo;
254 }
255 
257 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getVolumeFE()
258  const {
259  return static_cast<VolumeElementForcesAndSourcesCoreOnSideBase *>(ptrFE);
260 }
261 
262 int VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
263  getFaceSense() const {
264  return getVolumeFE()->faceSense;
265 }
266 
267 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
268  getFTensor1Normal() {
269  double *ptr = &*getNormal().data().begin();
270  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
271 }
272 
273 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
274  getFTensor1Tangent1() {
275  double *ptr = &*getTangent1().data().begin();
276  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
277 }
278 
279 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
280  getFTensor1Tangent2() {
281  double *ptr = &*getTangent2().data().begin();
282  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
283 }
284 
285 inline auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
286  getFTensor1NormalsAtGaussPts() {
287  double *ptr = &*getNormalsAtGaussPts().data().begin();
288  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
289  &ptr[2]);
290 }
291 
292 int VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
293  getFaceSideNumber() const {
294  return getVolumeFE()->faceSideNumber;
295 }
296 
297 MatrixDouble &VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
298  getNormalsAtGaussPts() {
299  return getFaceFE()->normalsAtGaussPts;
300 }
301 
302 MatrixDouble &VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
303  getFaceCoordsAtGaussPts() {
304  return getFaceFE()->coordsAtGaussPts;
305 }
306 
307 bool VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getEdgeFace(
308  const int ee) const {
309  constexpr bool edges_on_faces[6][4] = {{true, false, false, true}, // e0
310  {false, true, false, true}, // e1
311  {false, false, true, true}, // e2
312  {true, false, true, false}, // e3
313  {true, true, false, false}, // e4
314  {false, true, true, false}};
315  return edges_on_faces[ee][getFaceSideNumber()];
316 }
317 
318 ublas::matrix_row<MatrixDouble> VolumeElementForcesAndSourcesCoreOnSideBase::
319  UserDataOperator::getNormalsAtGaussPts(const int gg) {
320  return ublas::matrix_row<MatrixDouble>(getNormalsAtGaussPts(), gg);
321 }
322 
323 template <int SWITCH>
326  return OpSwitch<SWITCH>();
327 }
328 
329 } // namespace MoFEM
330 
331 #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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
ForcesAndSourcesCore::UserDataOperator UserDataOperator
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...