v0.13.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, 4>&
45  */
46  inline const std::array<int, 4> &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, 8> &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  */
85  struct UserDataOperator;
86 
87  int getRule(int order);
88  MoFEMErrorCode setGaussPts(int order);
89 
90 private:
91  int faceSense; ///< Sense of face, could be 1 or -1
92  int faceSideNumber; ///< Face side number
93  std::array<int, 4> faceConnMap;
94  std::array<int, 8> tetConnMap;
96 };
97 
98 /** \brief default operator for TET element
99  * \ingroup mofem_forces_and_sources_volume_element
100  */
103 
105 
106  inline VolumeElementForcesAndSourcesCoreOnSideBase *getVolumeFE() const;
107 
108  inline FaceElementForcesAndSourcesCoreBase *getFaceFE() const;
109 
110  /**
111  * \brief get face sense in respect to volume
112  * @return error code
113  */
114  inline int getFaceSense() const;
115 
116  /**
117  * \brief get face side number in respect to volume
118  * @return error code
119  */
120  inline int getFaceSideNumber() const;
121 
122  inline bool getEdgeFace(const int ee) const;
123 
124  /**
125  * get face normal on side which is this element
126  * @return face normal
127  */
128  inline VectorDouble &getNormal();
129 
130  /** \brief get triangle tangent 1
131  */
132  inline VectorDouble &getTangent1();
133 
134  /** \brief get triangle tangent 2
135  */
136  inline VectorDouble &getTangent2();
137 
138  /** \brief get normal as tensor
139  */
140  inline auto getFTensor1Normal();
141 
142  /** \brief get tangentOne as tensor
143  */
144  inline auto getFTensor1Tangent1();
145 
146  /** \brief get tangentTwo as tensor
147  */
148  inline auto getFTensor1Tangent2();
149 
150  /** \brief if higher order geometry return normals at Gauss pts.
151 
152  Note: returned matrix has size 0 in rows and columns if no HO approximation
153  of geometry is available.
154 
155  */
156  inline MatrixDouble &getNormalsAtGaussPts();
157 
158  /** \brief if higher order geometry return normals at Gauss pts.
159  *
160  * \param gg gauss point number
161  */
162  inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPts(const int gg);
163 
164  /** \brief get normal at integration points
165 
166  Example:
167  \code
168  double nrm2;
169  FTensor::Index<'i',3> i;
170  auto t_normal = getFTensor1NormalsAtGaussPts();
171  for(int gg = gg!=data.getN().size1();gg++) {
172  nrm2 = sqrt(t_normal(i)*t_normal(i));
173  ++t_normal;
174  }
175  \endcode
176 
177  */
178  inline auto getFTensor1NormalsAtGaussPts();
179 
180  /** \brief get face coordinates at Gauss pts.
181 
182  \note Coordinates should be the same what function getCoordsAtGaussPts
183  on tets is returning. If both coordinates are different it is error, or you
184  do something very unusual.
185 
186  */
187  inline MatrixDouble &getFaceCoordsAtGaussPts();
188 
189 protected:
190  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
191 };
192 
193 /**
194  * @brief Volume side finite element with switches
195  *
196  * Using SWITCH to off functions
197  *
198  * @tparam SWITCH
199  */
200 template <int SWITCH>
203 
204  using VolumeElementForcesAndSourcesCoreOnSideBase::
205  VolumeElementForcesAndSourcesCoreOnSideBase;
206 
209 
210  MoFEMErrorCode operator()();
211 };
212 
213 /** \brief Volume element used to integrate on skeleton
214  \ingroup mofem_forces_and_sources_volume_element
215 
216  */
219 
220 const std::array<int, 4> &
221 VolumeElementForcesAndSourcesCoreOnSideBase::getFaceConnMap() const {
222  return faceConnMap;
223 }
224 
225 const std::array<int, 8> &
226 VolumeElementForcesAndSourcesCoreOnSideBase::getTetConnMap() const {
227  return tetConnMap;
228 }
229 
230 int VolumeElementForcesAndSourcesCoreOnSideBase::getOppositeNode() const {
231  return oppositeNode;
232 }
233 
234 int VolumeElementForcesAndSourcesCoreOnSideBase::getFaceSense() const {
235  return faceSense;
236 }
237 
238 int VolumeElementForcesAndSourcesCoreOnSideBase::getFaceSideNumber() const {
239  return faceSideNumber;
240 }
241 
243 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getFaceFE()
244  const {
245  return static_cast<FaceElementForcesAndSourcesCoreBase *>(
246  getVolumeFE()->sidePtrFE);
247 }
248 
249 VectorDouble &
250 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getNormal() {
251  return getFaceFE()->nOrmal;
252 }
253 
254 VectorDouble &
255 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getTangent1() {
256  return getFaceFE()->tangentOne;
257 }
258 
259 VectorDouble &
260 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getTangent2() {
261  return getFaceFE()->tangentTwo;
262 }
263 
265 VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getVolumeFE()
266  const {
267  return static_cast<VolumeElementForcesAndSourcesCoreOnSideBase *>(ptrFE);
268 }
269 
270 int VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
271  getFaceSense() const {
272  return getVolumeFE()->faceSense;
273 }
274 
275 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
276  getFTensor1Normal() {
277  double *ptr = &*getNormal().data().begin();
278  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
279 }
280 
281 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
282  getFTensor1Tangent1() {
283  double *ptr = &*getTangent1().data().begin();
284  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
285 }
286 
287 auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
288  getFTensor1Tangent2() {
289  double *ptr = &*getTangent2().data().begin();
290  return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
291 }
292 
293 inline auto VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
294  getFTensor1NormalsAtGaussPts() {
295  double *ptr = &*getNormalsAtGaussPts().data().begin();
296  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
297  &ptr[2]);
298 }
299 
300 int VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
301  getFaceSideNumber() const {
302  return getVolumeFE()->faceSideNumber;
303 }
304 
305 MatrixDouble &VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
306  getNormalsAtGaussPts() {
307  return getFaceFE()->normalsAtGaussPts;
308 }
309 
310 MatrixDouble &VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::
311  getFaceCoordsAtGaussPts() {
312  return getFaceFE()->coordsAtGaussPts;
313 }
314 
315 bool VolumeElementForcesAndSourcesCoreOnSideBase::UserDataOperator::getEdgeFace(
316  const int ee) const {
317  constexpr bool edges_on_faces[6][4] = {{true, false, false, true}, // e0
318  {false, true, false, true}, // e1
319  {false, false, true, true}, // e2
320  {true, false, true, false}, // e3
321  {true, true, false, false}, // e4
322  {false, true, true, false}};
323  return edges_on_faces[ee][getFaceSideNumber()];
324 }
325 
326 ublas::matrix_row<MatrixDouble> VolumeElementForcesAndSourcesCoreOnSideBase::
327  UserDataOperator::getNormalsAtGaussPts(const int gg) {
328  return ublas::matrix_row<MatrixDouble>(getNormalsAtGaussPts(), gg);
329 }
330 
331 template <int SWITCH>
334  return opSwitch<SWITCH>();
335 }
336 
337 } // namespace MoFEM
338 
339 #endif //__VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume 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
structure to get information form mofem into EntitiesFieldData