v0.15.0
Loading...
Searching...
No Matches
VolumeElementForcesAndSourcesCoreOnSide.hpp
Go to the documentation of this file.
1/** \file VolumeElementForcesAndSourcesCoreOnSide.hpp
2 \brief Side volume element for skeleton integration.
3
4 This header defines VolumeElementForcesAndSourcesCoreOnSide, which extends
5 VolumeElementForcesAndSourcesCore to evaluate volume finite elements on a
6 selected face (the "skeleton"). It provides connectivity maps, face
7 orientation information, and access to face geometry (normals, tangents)
8 for implementing mortar methods, flux terms, and interface conditions.
9
10 The nested UserDataOperator adapts the standard volume operator interface
11 for side evaluations. Users inherit from UserDataOperator and override
12 methods like doWork() to assemble residuals and tangent contributions using
13 face normals at Gauss points, face coordinates, and the skeleton sense.
14
15 Typical usage:
16 1. Instantiate VolumeElementForcesAndSourcesCoreOnSide for each side
17 2. Attach a derived UserDataOperator to handle side-specific assembly
18 3. Loop over face elements on the skeleton
19 4. For each face, find adjacent volume elements (sides) and loop over both
20 5. Access getFaceSideNumber(), getSkeletonSense(), getNormal(), etc. in your
21 operator to retrieve face geometry and orientation for each side
22
23 \author Anonymous author contributing to MoFEM under MiT License
24 \date June 2024
25*/
26
27#ifndef __VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
28#define __VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
29
30using namespace boost::numeric;
31
32namespace MoFEM {
33
34/**
35 * \brief Base volume element used to integrate on skeleton
36 * \ingroup mofem_forces_and_sources_volume_element
37 */
40
42
43 /**
44 * @brief Get the face nodes mapped on the volume element
45 *
46 * \todo This is not general; e.g., for quads the number of nodes is 4.
47 *
48 * @return const std::array<int, 4>&
49 */
50 inline const std::array<int, 4> &getFaceConnMap() const;
51
52 /**
53 * @brief Get face nodes mapped on volume
54 *
55 * \todo This is not general; e.g., for prisms or hexes the size of the fixed
56 * array is incorrect.
57 *
58 * @return const std::array<int, 8>&
59 */
60 inline const std::array<int, 8> &getTetConnMap() const;
61
62 /**
63 * @brief Get node on volume opposite to volume element
64 *
65 * \todo This is not general; for prisms or hexes the opposite node is
66 * not unique.
67 *
68 * @return int
69 */
70 inline int getOppositeNode() const;
71
72 /**
73 * @brief Sense face on volume
74 * @deprecated use getSkeletonSense()
75 *
76 * @return int
77 */
78 DEPRECATED inline int getFaceSense() const;
79
80 /**
81 * @brief Sense face on volume
82 *
83 * @return int
84 */
85 inline int getSkeletonSense() const;
86
87 /**
88 * @brief Face number on the volume
89 *
90 * @return int
91 */
92 inline int getFaceSideNumber() const;
93
94 /** \brief default operator for TET element
95 * \ingroup mofem_forces_and_sources_volume_element
96 */
97 struct UserDataOperator;
98
99 /**
100 * \brief Compute quadrature rule index for a given polynomial order
101 *
102 * Note: This function for volume side by default return -1, since the integartion
103 * rule is set by
104 *
105 *
106 * @param order Requested polynomial order
107 * @return Integration rule index compatible with the requested order
108 */
109 int getRule(int order);
110
111 /**
112 * \brief Configure Gauss points for the element side
113 *
114 * @param order Requested polynomial order for integration
115 * @return Error code (0 on success)
116 */
118
119 /**
120 * \brief Execute user operator
121 * @return Error code (0 on success)
122 */
124
125private:
126
127 /**
128 * \brief Implementation detail for operator()
129 *
130 * @return Error code (0 on success)
131 */
133
134 int faceSense; ///< Sense of face, could be 1 or -1
135 int faceSideNumber; ///< Face side number
136 int nbNodesOnFace; ///< Number of nodes on face
137 std::array<int, 4> faceConnMap;
138 std::array<int, 8> tetConnMap;
140 bool operatorImplCalled = false; ///< flag to check if operatorImpl was called
141};
142
143/** \brief default operator for TET element
144 * \ingroup mofem_forces_and_sources_volume_element
145 */
148
149 using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
150
151 /**
152 * \brief Access the owning volume finite element
153 *
154 * Volume element is side element, to face element.
155 *
156 * @return Pointer to the parent volume finite element
157 */
159
160 /**
161 * \brief Access the face finite element
162 * Face element, execute side element, that is volume, on which the
163 * volume side operators are executed.
164 *
165 * @return Pointer to the face finite element attached to this side
166 */
168
169 /**
170 * \brief Get face sense with respect to the volume
171 * \deprecated Use getSkeletonSense()
172 * @return Face sense (+1 or -1)
173 */
174 DEPRECATED inline int getFaceSense() const;
175
176 /**
177 * @brief Get the skeleton sense
178 *
179 * Calls getFaceSense().
180 *
181 * @return Skeleton sense (+1 or -1)
182 */
183 inline int getSkeletonSense() const;
184
185 /**
186 * \brief Get face side number with respect to the volume
187 * @return Face side number
188 */
189 inline int getFaceSideNumber() const;
190
191 /**
192 * \brief Query if an edge belongs to the current face
193 * @param ee Local edge index (0-5 for tetrahedra)
194 * @return True when the edge lies on the active face
195 */
196 inline bool getEdgeFace(const int ee) const;
197
198 /**
199 * \brief Get the face normal for this side element
200 * @return Face normal vector
201 */
202 inline VectorDouble &getNormal();
203
204 /** \brief Get first triangle tangent
205 */
206 inline VectorDouble &getTangent1();
207
208 /** \brief Get second triangle tangent
209 */
210 inline VectorDouble &getTangent2();
211
212 /** \brief Get normal as tensor
213 */
214 inline auto getFTensor1Normal();
215
216 /** \brief Get tangent one as tensor
217 */
218 inline auto getFTensor1Tangent1();
219
220 /** \brief Get tangent two as tensor
221 */
222 inline auto getFTensor1Tangent2();
223
224 /** \brief If higher-order geometry is available, return normals at Gauss points.
225 *
226 * \note The returned matrix has zero rows and columns when no higher-order
227 * approximation of geometry is available.
228 */
230
231 /** \brief If higher-order geometry is available, return normals at Gauss points.
232 *
233 * \param gg Gauss point number
234 */
235 inline ublas::matrix_row<MatrixDouble> getNormalsAtGaussPts(const int gg);
236
237 /** \brief Get normal at integration points
238 *
239 * Example:
240 * \code
241 * double nrm2;
242 * FTensor::Index<'i',3> i;
243 * auto t_normal = getFTensor1NormalsAtGaussPts();
244 * for(int gg = gg!=data.getN().size1();gg++) {
245 * nrm2 = sqrt(t_normal(i)*t_normal(i));
246 * ++t_normal;
247 * }
248 * \endcode
249 */
250 inline auto getFTensor1NormalsAtGaussPts();
251
252 /** \brief Get face coordinates at Gauss points.
253 *
254 * \note Coordinates should match the values returned by getCoordsAtGaussPts
255 * on tetrahedra. If they differ, it indicates an error or an unusual setup.
256 */
258
259protected:
260 /**
261 * \brief Assign base class pointer for the operator
262 * @param ptr Pointer to the forces and sources core instance
263 * @return Error code (0 on success)
264 */
266};
267
268const std::array<int, 4> &
272
273const std::array<int, 8> &
277
281
285
289
293
295VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getFaceFE() const {
296 return static_cast<FaceElementForcesAndSourcesCore *>(
298}
299
301VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getNormal() {
302 return getFaceFE()->nOrmal;
303}
304
306VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getTangent1() {
307 return getFaceFE()->tangentOne;
308}
309
311VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getTangent2() {
312 return getFaceFE()->tangentTwo;
313}
314
316VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getVolumeFE() const {
317 return static_cast<VolumeElementForcesAndSourcesCoreOnSide *>(ptrFE);
318}
319
320int VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getFaceSense()
321 const {
322 return getSkeletonSense();
323}
324
325int VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
326 getSkeletonSense() const {
327 return getVolumeFE()->faceSense;
328}
329
330auto VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
331 getFTensor1Normal() {
332 double *ptr = &*getNormal().data().begin();
333 return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
334}
335
336auto VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
337 getFTensor1Tangent1() {
338 double *ptr = &*getTangent1().data().begin();
339 return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
340}
341
342auto VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
343 getFTensor1Tangent2() {
344 double *ptr = &*getTangent2().data().begin();
345 return FTensor::Tensor1<double *, 3>(ptr, &ptr[1], &ptr[2]);
346}
347
348inline auto VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
349 getFTensor1NormalsAtGaussPts() {
350 double *ptr = &*getNormalsAtGaussPts().data().begin();
352 &ptr[2]);
353}
354
355int VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
356 getFaceSideNumber() const {
357 return getVolumeFE()->faceSideNumber;
358}
359
360MatrixDouble &VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
361 getNormalsAtGaussPts() {
362 return getFaceFE()->normalsAtGaussPts;
363}
364
365MatrixDouble &VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::
366 getFaceCoordsAtGaussPts() {
367 return getFaceFE()->coordsAtGaussPts;
368}
369
370bool VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getEdgeFace(
371 const int ee) const {
372 constexpr bool edges_on_faces[6][4] = {{true, false, false, true}, // e0
373 {false, true, false, true}, // e1
374 {false, false, true, true}, // e2
375 {true, false, true, false}, // e3
376 {true, true, false, false}, // e4
377 {false, true, true, false}};
378 return edges_on_faces[ee][getFaceSideNumber()];
379}
380
381ublas::matrix_row<MatrixDouble>
382VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator::getNormalsAtGaussPts(
383 const int gg) {
384 return ublas::matrix_row<MatrixDouble>(getNormalsAtGaussPts(), gg);
385}
386
387/**
388 * @deprecated Use VolumeElementForcesAndSourcesCore.
389 */
392
393/**
394 * @deprecated Do not use this template; it is obsolete.
395 */
396template <int SWITCH>
404
405} // namespace MoFEM
406
407#endif //__VOLUMEELEMENTFORCESANDSOURCESCORE_ONSIDE_HPP__
#define DEPRECATED
Definition definitions.h:17
constexpr int order
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
DEPRECATED typedef VolumeElementForcesAndSourcesCoreOnSide VolumeElementForcesAndSourcesCoreOnSideBase
structure to get information from mofem into EntitiesFieldData
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
bool getEdgeFace(const int ee) const
Query if an edge belongs to the current face.
FaceElementForcesAndSourcesCore * getFaceFE() const
Access the face finite element Face element, execute side element, that is volume,...
int getFaceSideNumber() const
Get face side number with respect to the volume.
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
Assign base class pointer for the operator.
DEPRECATED int getFaceSense() const
Get face sense with respect to the volume.
MatrixDouble & getNormalsAtGaussPts()
If higher-order geometry is available, return normals at Gauss points.
VolumeElementForcesAndSourcesCoreOnSide * getVolumeFE() const
Access the owning volume finite element.
MoFEMErrorCode operatorImpl()
Implementation detail for operator()
const std::array< int, 8 > & getTetConnMap() const
Get face nodes mapped on volume.
const std::array< int, 4 > & getFaceConnMap() const
Get the face nodes mapped on the volume element.
int getRule(int order)
Compute quadrature rule index for a given polynomial order.
int getOppositeNode() const
Get node on volume opposite to volume element.
MoFEMErrorCode setGaussPts(int order)
Configure Gauss points for the element side.