v0.9.1
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 
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 
48 
49  /** \brief default operator for EDGE element
50  \ingroup mofem_forces_and_sources_edge_element
51  */
53 
55 
56  /** \brief get element connectivity
57  */
58  inline const EntityHandle *getConn();
59 
60  /**
61  * \brief get edge length
62  */
63  inline double getLength();
64 
65  /**
66  * \brief get measure of element
67  * @return length of face
68  */
69  inline double getMeasure();
70 
71  /**
72  * \brief get edge direction
73  */
74  inline VectorDouble &getDirection();
75 
76  /**
77  * \brief get edge node coordinates
78  */
79  inline VectorDouble &getCoords();
80 
81  /**
82  * \brief get coordinate at integration point
83  */
84  inline MatrixDouble &getCoordsAtGaussPts();
85 
86  /** \brief get coordinates at Gauss pts.
87  */
88  inline auto getFTensor1CoordsAtGaussPts();
89 
90  /**
91  * \brief get tangent vector to edge curve at integration points
92  */
93  inline MatrixDouble &getTangetAtGaussPts();
94 
95  /**
96  * \brief get pointer to this finite element
97  */
98  inline const EdgeElementForcesAndSourcesCoreBase *getEdgeFE();
99 
100  inline FTensor::Tensor1<double, 3> getFTensor1Direction();
101 
102  /**
103  * \brief get get coords at gauss points
104 
105  \code
106  FTensor::Index<'i',3> i;
107  auto t_center;
108  auto t_coords = getTensor1Coords();
109  t_center(i) = 0;
110  for(int nn = 0;nn!=2;nn++) {
111  t_center(i) += t_coords(i);
112  ++t_coords;
113  }
114  t_center(i) /= 2;
115  \endcode
116 
117  */
119  getTensor1Coords();
120 
122  getFTensor1TangentAtGaussPts();
123 
124  template <int SWITCH>
126  loopSideFaces(const string &fe_name,
128 
129  protected:
130  MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr);
131 
132  };
133 
134  enum Switches {
135  NO_HO_GEOMETRY = 1 << 0,
136  NO_COVARIANT_TRANSFORM_HCURL = 1 << 2
137  };
138 
139  template <int SWITCH> MoFEMErrorCode OpSwitch();
140 
141 protected:
145 
146  double lEngth;
147 
148  int numNodes;
153 
154  MoFEMErrorCode calculateEdgeDirection();
155  MoFEMErrorCode setIntegrationPts();
156  MoFEMErrorCode calculateCoordsAtIntegrationPts();
157  MoFEMErrorCode calculateHoCoordsAtIntegrationPts();
158 
160 };
161 
162 /** \brief Edge finite element
163  * \ingroup mofem_forces_and_sources_edge_element
164  *
165  * User is implementing own operator at Gauss points level, by own object
166  * derived from EdgeElementForcesAndSourcesCoreL::UserDataOperator. Arbitrary
167  * number of operator added pushing objects to rowOpPtrVector and
168  * rowColOpPtrVector.
169  *
170  */
171 template <int SWITCH>
174 
175  using EdgeElementForcesAndSourcesCoreBase::
176  EdgeElementForcesAndSourcesCoreBase;
177 
178  MoFEMErrorCode operator()();
179 };
180 
181 /** \brief Edge finite element default
182  \ingroup mofem_forces_and_sources_edge_element
183 
184  */
187 
188 template <int SWITCH>
189 MoFEMErrorCode EdgeElementForcesAndSourcesCoreBase::OpSwitch() {
191 
192  if (numeredEntFiniteElementPtr->getEntType() != MBEDGE)
194 
195  CHKERR createDataOnElement();
196 
197  DataForcesAndSourcesCore &data_curl = *dataOnElement[HCURL];
198 
199  CHKERR calculateEdgeDirection();
200  CHKERR getSpacesAndBaseOnEntities(dataH1);
201  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
202  dataH1.dataOnEntities[MBEDGE][0].getSense() =
203  1; // set sense to 1, this is this entity
204 
205  // Hcurl
206  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
207  CHKERR getEntityDataOrder<MBEDGE>(data_curl, HCURL);
208  data_curl.dataOnEntities[MBEDGE][0].getSense() =
209  1; // set sense to 1, this is this entity
210  data_curl.spacesOnEntities[MBEDGE].set(HCURL);
211  }
212 
213  CHKERR setIntegrationPts();
214  CHKERR calculateCoordsAtIntegrationPts();
215  CHKERR calHierarchicalBaseFunctionsOnElement();
216  CHKERR calBernsteinBezierBaseFunctionsOnElement();
217  if (!(SWITCH & NO_HO_GEOMETRY))
218  CHKERR calculateHoCoordsAtIntegrationPts();
219 
220  if (!(SWITCH & NO_COVARIANT_TRANSFORM_HCURL))
221  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL))
222  CHKERR opCovariantTransform.opRhs(data_curl);
223 
224  // Iterate over operators
225  CHKERR loopOverOperators();
226 
228 }
229 
230 template <int SWITCH>
232  return OpSwitch<SWITCH>();
233 }
234 
235 const EntityHandle *
236 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getConn() {
237  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->cOnn;
238 }
239 
240 double EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getLength() {
241  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->lEngth;
242 }
243 
244 double EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() {
245  return getLength();
246 }
247 
248 VectorDouble &
249 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getDirection() {
250  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->dIrection;
251 }
252 
253 VectorDouble &
254 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getCoords() {
255  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->cOords;
256 }
257 
258 MatrixDouble &
259 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getCoordsAtGaussPts() {
260  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)
261  ->coordsAtGaussPts;
262 }
263 
264 auto EdgeElementForcesAndSourcesCoreBase::UserDataOperator::
265  getFTensor1CoordsAtGaussPts() {
266  double *ptr = &*getCoordsAtGaussPts().data().begin();
267  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
268  &ptr[2]);
269 }
270 
271 MatrixDouble &
272 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getTangetAtGaussPts() {
273  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)
274  ->tangentAtGaussPts;
275 }
276 
278 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getEdgeFE() {
279  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE);
280 }
281 
283 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getFTensor1Direction() {
284  return FTensor::Tensor1<double, 3>(getDirection()[0], getDirection()[1],
285  getDirection()[2]);
286 }
287 
289 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getTensor1Coords() {
290  double *ptr = &*getCoords().data().begin();
291  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
292  &ptr[2]);
293 }
294 
296 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::
297  getFTensor1TangentAtGaussPts() {
298  double *ptr = &*getTangetAtGaussPts().data().begin();
299  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
300  &ptr[2]);
301 }
302 
303 template <int SWITCH>
305 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::loopSideFaces(
306  const string &fe_name,
308  return loopSide(fe_name, &fe_side, 2);
309 }
310 
311 } // namespace MoFEM
312 
313 #endif //__EDGEELEMENTFORCESANDSURCESCORE_HPP__
314 
315 /**
316  * \defgroup mofem_forces_and_sources_edge_element Edge Element
317  *
318  * \brief Implementation of edge element.
319  *
320  * \ingroup mofem_forces_and_sources
321  */
Calculate tangent vector on edge form HO geometry approximation.
Deprecated interface functions.
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
transform Hcurl base fluxes from reference element to physical edge
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:513
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
T data[Tensor_Dim]
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
field with continuous tangents
Definition: definitions.h:177
#define CHKERR
Inline error check.
Definition: definitions.h:601
ForcesAndSourcesCore::UserDataOperator UserDataOperator
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
structure to get information form mofem into DataForcesAndSourcesCore
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
continuous field
Definition: definitions.h:176
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...