v0.13.0
EdgeElementForcesAndSourcesCore.hpp
Go to the documentation of this file.
1 /** \file EdgeElementForcesAndSourcesCore.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 __EDGEELEMENTFORCESANDSURCESCORE_HPP__
25 #define __EDGEELEMENTFORCESANDSURCESCORE_HPP__
26 
27 using namespace boost::numeric;
28 
29 namespace MoFEM {
30 
31 template <int SWITCH> struct FaceElementForcesAndSourcesCoreOnSideSwitch;
32 
33 struct FaceElementForcesAndSourcesCoreOnSideBase;
34 
35 /** \brief Edge finite element
36  * \ingroup mofem_forces_and_sources_edge_element
37  *
38  * User is implementing own operator at Gauss points level, by own object
39  * derived from EdgeElementForcesAndSourcesCoreL::UserDataOperator. Arbitrary
40  * number of operator added pushing objects to rowOpPtrVector and
41  * rowColOpPtrVector.
42  *
43  */
45 
47 
48  /** \brief default operator for EDGE element
49  \ingroup mofem_forces_and_sources_edge_element
50  */
51  struct UserDataOperator;
52 
53  enum Switches {};
54 
55  template <int SWITCH> MoFEMErrorCode opSwitch();
56 
57 protected:
59 
61 
62  double lEngth;
63 
64  int numNodes;
68 
69  MoFEMErrorCode calculateEdgeDirection();
70  MoFEMErrorCode setIntegrationPts();
71  MoFEMErrorCode calculateCoordsAtIntegrationPts();
72 
74 };
75 
76 /** \brief default operator for EDGE element
77  \ingroup mofem_forces_and_sources_edge_element
78  */
81 
83 
84  /** \brief get element connectivity
85  */
86  inline const EntityHandle *getConn();
87 
88  /**
89  * \brief get edge length
90  */
91  inline double getLength();
92 
93  /**
94  * \brief get measure of element
95  * @return length of face
96  */
97  inline double getMeasure();
98 
99  /**
100  * \brief get edge direction
101  */
102  inline VectorDouble &getDirection();
103 
104  /**
105  * \brief get edge normal
106  * NOTE: it should be used only in 2D analysis
107  */
108  inline auto getFTensor1Normal();
109 
110  /**
111  * @brief get ftensor1 edge normal
112  *
113  * @param vec vector in third direction
114  * @return auto
115  */
116  inline auto getFTensor1Normal(const FTensor::Tensor1<double, 3> &vec);
117 
118  /**
119  * \brief get edge node coordinates
120  */
121  inline VectorDouble &getCoords();
122 
123  /**
124  * \brief get tangent vector to edge curve at integration points
125  */
126  inline MatrixDouble &getTangetAtGaussPts();
127 
128  /**
129  * \brief get pointer to this finite element
130  */
131  inline const EdgeElementForcesAndSourcesCoreBase *getEdgeFE();
132 
133  inline FTensor::Tensor1<double, 3> getFTensor1Direction();
134 
135  /**
136  * \brief get get coords at gauss points
137 
138  \code
139  FTensor::Index<'i',3> i;
140  auto t_center;
141  auto t_coords = getTensor1Coords();
142  t_center(i) = 0;
143  for(int nn = 0;nn!=2;nn++) {
144  t_center(i) += t_coords(i);
145  ++t_coords;
146  }
147  t_center(i) /= 2;
148  \endcode
149 
150  */
151  inline FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> getTensor1Coords();
152 
154  getFTensor1TangentAtGaussPts();
155 
156  template <int SWITCH>
158  loopSideFaces(const string &fe_name,
160 
161 protected:
162  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
163 };
164 
165 /** \brief Edge finite element
166  * \ingroup mofem_forces_and_sources_edge_element
167  *
168  * User is implementing own operator at Gauss points level, by own object
169  * derived from EdgeElementForcesAndSourcesCoreL::UserDataOperator. Arbitrary
170  * number of operator added pushing objects to rowOpPtrVector and
171  * rowColOpPtrVector.
172  *
173  */
174 template <int SWITCH>
177 
180 
181  MoFEMErrorCode operator()();
182 };
183 
184 /** \brief Edge finite element default
185  \ingroup mofem_forces_and_sources_edge_element
186 
187  */
190 
191 template <int SWITCH>
192 MoFEMErrorCode EdgeElementForcesAndSourcesCoreBase::opSwitch() {
194 
195  if (numeredEntFiniteElementPtr->getEntType() != MBEDGE)
197 
198  CHKERR calculateEdgeDirection();
199  CHKERR getSpacesAndBaseOnEntities(dataH1);
200  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
201  dataH1.dataOnEntities[MBEDGE][0].getSense() =
202  1; // set sense to 1, this is this entity
203 
204  // L2
205  if (dataH1.spacesOnEntities[MBEDGE].test(L2)) {
206  auto &data_l2 = *dataOnElement[L2];
207  CHKERR getEntityDataOrder<MBEDGE>(data_l2, L2);
208  data_l2.dataOnEntities[MBEDGE][0].getSense() =
209  1; // set sense to 1, this is this entity
210  data_l2.spacesOnEntities[MBEDGE].set(L2);
211  }
212 
213  // Hcurl
214  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
215  auto &data_curl = *dataOnElement[HCURL];
216  CHKERR getEntityDataOrder<MBEDGE>(data_curl, HCURL);
217  data_curl.dataOnEntities[MBEDGE][0].getSense() =
218  1; // set sense to 1, this is this entity
219  data_curl.spacesOnEntities[MBEDGE].set(HCURL);
220  }
221 
222  CHKERR setIntegrationPts();
223  CHKERR calculateCoordsAtIntegrationPts();
224  CHKERR calHierarchicalBaseFunctionsOnElement();
225  CHKERR calBernsteinBezierBaseFunctionsOnElement();
226 
227  // Iterate over operators
228  CHKERR loopOverOperators();
229 
231 }
232 
233 template <int SWITCH>
235  return opSwitch<SWITCH>();
236 }
237 
238 const EntityHandle *
239 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getConn() {
240  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->cOnn;
241 }
242 
243 double EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getLength() {
244  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->lEngth;
245 }
246 
247 double EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() {
248  return getLength();
249 }
250 
251 VectorDouble &
252 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getDirection() {
253  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->dIrection;
254 }
255 
256 VectorDouble &
257 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getCoords() {
258  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->cOords;
259 }
260 
261 MatrixDouble &
262 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getTangetAtGaussPts() {
263  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)
264  ->tangentAtGaussPts;
265 }
266 
268 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getEdgeFE() {
269  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE);
270 }
271 
273 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getFTensor1Direction() {
274  return FTensor::Tensor1<double, 3>(getDirection()[0], getDirection()[1],
275  getDirection()[2]);
276 }
277 
279 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getTensor1Coords() {
280  double *ptr = &*getCoords().data().begin();
281  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
282  &ptr[2]);
283 }
284 
286 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::
287  getFTensor1TangentAtGaussPts() {
288  double *ptr = &*getTangetAtGaussPts().data().begin();
289  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
290  &ptr[2]);
291 }
292 
293 auto EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getFTensor1Normal(
294  const FTensor::Tensor1<double, 3> &vec) {
299  auto t_dir = getFTensor1Direction();
300  t_normal(i) = FTensor::levi_civita(i, j, k) * t_dir(j) * vec(k);
301  return t_normal;
302 }
303 
304 auto EdgeElementForcesAndSourcesCoreBase::UserDataOperator::
305  getFTensor1Normal() {
306  FTensor::Tensor1<double, 3> t_normal{0., 0., 1.};
307  return getFTensor1Normal(t_normal);
308 }
309 
310 template <int SWITCH>
312 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::loopSideFaces(
313  const string &fe_name,
315  return loopSide(fe_name, &fe_side, 2);
316 }
317 
318 } // namespace MoFEM
319 
320 #endif //__EDGEELEMENTFORCESANDSURCESCORE_HPP__
321 
322 /**
323  * \defgroup mofem_forces_and_sources_edge_element Edge Element
324  *
325  * \brief Implementation of edge element.
326  *
327  * \ingroup mofem_forces_and_sources
328  */
ForcesAndSourcesCore::UserDataOperator UserDataOperator
T data[Tensor_Dim]
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
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
Deprecated interface functions.
structure to get information form mofem into EntitiesFieldData