v0.9.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 
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 
130  enum Switches {
131  NO_HO_GEOMETRY = 1 << 0,
132  NO_COVARIANT_TRANSFORM_HCURL = 1 << 2
133  };
134 
135  template <int SWITCH> MoFEMErrorCode OpSwitch();
136 
137 protected:
141 
142  double lEngth;
143 
144  int numNodes;
149 
150  MoFEMErrorCode calculateEdgeDirection();
151  MoFEMErrorCode setIntegrationPts();
152  MoFEMErrorCode calculateCoordsAtIntegrationPts();
153  MoFEMErrorCode calculateHoCoordsAtIntegrationPts();
154 
156 };
157 
158 /** \brief Edge finite element
159  * \ingroup mofem_forces_and_sources_edge_element
160  *
161  * User is implementing own operator at Gauss points level, by own object
162  * derived from EdgeElementForcesAndSourcesCoreL::UserDataOperator. Arbitrary
163  * number of operator added pushing objects to rowOpPtrVector and
164  * rowColOpPtrVector.
165  *
166  */
167 template <int SWITCH>
170 
171  using EdgeElementForcesAndSourcesCoreBase::
172  EdgeElementForcesAndSourcesCoreBase;
173 
174  MoFEMErrorCode operator()();
175 };
176 
177 /** \brief Edge finite element default
178  \ingroup mofem_forces_and_sources_edge_element
179 
180  */
183 
184 template <int SWITCH>
185 MoFEMErrorCode EdgeElementForcesAndSourcesCoreBase::OpSwitch() {
187 
188  if (numeredEntFiniteElementPtr->getEntType() != MBEDGE)
190 
191  CHKERR createDataOnElement();
192 
193  DataForcesAndSourcesCore &data_curl = *dataOnElement[HCURL];
194 
195  CHKERR calculateEdgeDirection();
196  CHKERR getSpacesAndBaseOnEntities(dataH1);
197  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
198  dataH1.dataOnEntities[MBEDGE][0].getSense() =
199  1; // set sense to 1, this is this entity
200 
201  // Hcurl
202  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
203  CHKERR getEntityDataOrder<MBEDGE>(data_curl, HCURL);
204  data_curl.dataOnEntities[MBEDGE][0].getSense() =
205  1; // set sense to 1, this is this entity
206  data_curl.spacesOnEntities[MBEDGE].set(HCURL);
207  }
208 
209  CHKERR setIntegrationPts();
210  CHKERR calculateCoordsAtIntegrationPts();
211  CHKERR calHierarchicalBaseFunctionsOnElement();
212  CHKERR calBernsteinBezierBaseFunctionsOnElement();
213  if (!(SWITCH & NO_HO_GEOMETRY))
214  CHKERR calculateHoCoordsAtIntegrationPts();
215 
216  if (!(SWITCH & NO_COVARIANT_TRANSFORM_HCURL))
217  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL))
218  CHKERR opCovariantTransform.opRhs(data_curl);
219 
220  // Iterate over operators
221  CHKERR loopOverOperators();
222 
224 }
225 
226 template <int SWITCH>
228  return OpSwitch<SWITCH>();
229 }
230 
231 const EntityHandle *
232 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getConn() {
233  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->cOnn;
234 }
235 
236 double EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getLength() {
237  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->lEngth;
238 }
239 
240 double EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure() {
241  return getLength();
242 }
243 
244 VectorDouble &
245 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getDirection() {
246  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->dIrection;
247 }
248 
249 VectorDouble &
250 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getCoords() {
251  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)->cOords;
252 }
253 
254 MatrixDouble &
255 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getCoordsAtGaussPts() {
256  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)
257  ->coordsAtGaussPts;
258 }
259 
260 auto EdgeElementForcesAndSourcesCoreBase::UserDataOperator::
261  getFTensor1CoordsAtGaussPts() {
262  double *ptr = &*getCoordsAtGaussPts().data().begin();
263  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
264  &ptr[2]);
265 }
266 
267 MatrixDouble &
268 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getTangetAtGaussPts() {
269  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE)
270  ->tangentAtGaussPts;
271 }
272 
274 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getEdgeFE() {
275  return static_cast<EdgeElementForcesAndSourcesCoreBase *>(ptrFE);
276 }
277 
279 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getFTensor1Direction() {
280  return FTensor::Tensor1<double, 3>(getDirection()[0], getDirection()[1],
281  getDirection()[2]);
282 }
283 
285 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::getTensor1Coords() {
286  double *ptr = &*getCoords().data().begin();
287  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
288  &ptr[2]);
289 }
290 
292 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::
293  getFTensor1TangentAtGaussPts() {
294  double *ptr = &*getTangetAtGaussPts().data().begin();
295  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(ptr, &ptr[1],
296  &ptr[2]);
297 }
298 
299 template <int SWITCH>
301 EdgeElementForcesAndSourcesCoreBase::UserDataOperator::loopSideFaces(
302  const string &fe_name,
304  return loopSide(fe_name, &fe_side, 2);
305 }
306 
307 } // namespace MoFEM
308 
309 #endif //__EDGEELEMENTFORCESANDSURCESCORE_HPP__
310 
311 /**
312  * \defgroup mofem_forces_and_sources_edge_element Edge Element
313  *
314  * \brief Implementation of edge element.
315  *
316  * \ingroup mofem_forces_and_sources
317  */
Calculate tangent vector on edge form HO geometry approximation.
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:477
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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:172
#define CHKERR
Inline error check.
Definition: definitions.h:596
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:407
continuous field
Definition: definitions.h:171
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...