v0.14.0
EdgeElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file EdgeElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of edge element
4 
5 */
6 
7 
8 
9 #ifdef __cplusplus
10 extern "C" {
11 #endif
12 #include <cblas.h>
13 #include <lapack_wrap.h>
14 // #include <gm_rule.h>
15 #include <quad.h>
16 #ifdef __cplusplus
17 }
18 #endif
19 
20 namespace MoFEM {
21 
23  0., 0., 1.};
24 
26  Interface &m_field)
27  : ForcesAndSourcesCore(m_field),
28  meshPositionsFieldName("MESH_NODE_POSITIONS"), lEngth(elementMeasure) {
30  boost::shared_ptr<BaseFunction>(new EdgePolynomialBase());
32  "Problem with creation data on element");
33 }
34 
38  CHKERR mField.get_moab().get_connectivity(ent, cOnn, numNodes, true);
39  cOords.resize(numNodes * 3, false);
40  CHKERR mField.get_moab().get_coords(cOnn, numNodes, &*cOords.data().begin());
41  dIrection.resize(3, false);
42  cblas_dcopy(3, &cOords[3], 1, &*dIrection.data().begin(), 1);
43  cblas_daxpy(3, -1., &cOords[0], 1, &*dIrection.data().begin(), 1);
44  lEngth = cblas_dnrm2(3, &*dIrection.data().begin(), 1);
46 }
47 
50  int order_data = getMaxDataOrder();
51  int order_row = getMaxRowOrder();
52  int order_col = getMaxColOrder();
53  int rule = getRule(order_row, order_col, order_data);
54 
55  int nb_gauss_pts;
56  if (rule >= 0) {
57  if (rule < QUAD_1D_TABLE_SIZE) {
58  if (QUAD_1D_TABLE[rule]->dim != 1) {
59  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
60  }
61  if (QUAD_1D_TABLE[rule]->order < rule) {
62  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
63  "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
64  }
65  nb_gauss_pts = QUAD_1D_TABLE[rule]->npoints;
66  gaussPts.resize(2, nb_gauss_pts, false);
67  cblas_dcopy(nb_gauss_pts, &QUAD_1D_TABLE[rule]->points[1], 2,
68  &gaussPts(0, 0), 1);
69  cblas_dcopy(nb_gauss_pts, QUAD_1D_TABLE[rule]->weights, 1,
70  &gaussPts(1, 0), 1);
71  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 2,
72  false);
73  double *shape_ptr =
74  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
75  cblas_dcopy(2 * nb_gauss_pts, QUAD_1D_TABLE[rule]->points, 1, shape_ptr,
76  1);
77  } else {
78  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
79  "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
80  nb_gauss_pts = 0;
81  }
82  } else {
83  // If rule is negative, set user defined integration points
84  CHKERR setGaussPts(order_row, order_col, order_data);
85  nb_gauss_pts = gaussPts.size2();
86  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 2,
87  false);
88  if (nb_gauss_pts) {
89  for (int gg = 0; gg != nb_gauss_pts; gg++) {
90  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0) =
91  N_MBEDGE0(gaussPts(0, gg));
92  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 1) =
93  N_MBEDGE1(gaussPts(0, gg));
94  }
95  }
96  }
98 }
99 
104  const auto nb_gauss_pts = gaussPts.size2();
105  coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
106  tangentAtGaussPts.resize(nb_gauss_pts, 3, false);
107 
108  if (nb_gauss_pts) {
110  &*gaussPts.data().begin());
112  &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
113  &coordsAtGaussPts(0, 2));
115  &tangentAtGaussPts(0, 0), &tangentAtGaussPts(0, 1),
116  &tangentAtGaussPts(0, 2));
117 
118  FTensor::Tensor1<double, 3> t_coords_node0(cOords[0], cOords[1], cOords[2]);
119  FTensor::Tensor1<double, 3> t_coords_node1(cOords[3], cOords[4], cOords[5]);
121  t_dir(i) = t_coords_node1(i) - t_coords_node0(i);
122 
123  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
124  t_coords(i) = N_MBEDGE0(t_ksi) * t_coords_node0(i) +
125  N_MBEDGE1(t_ksi) * t_coords_node1(i);
126  t_tangent(i) = t_dir(i);
127  ++t_tangent;
128  ++t_coords;
129  ++t_ksi;
130  }
131  }
133 }
134 
136  ForcesAndSourcesCore *ptr) {
138  if (!(ptrFE = dynamic_cast<EdgeElementForcesAndSourcesCore *>(ptr)))
139  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
140  "User operator and finite element do not work together");
142 }
143 
146  const string fe_name, FaceElementForcesAndSourcesCoreOnSide &fe_side) {
147  return loopSide(fe_name, &fe_side, 2);
148 }
149 
152 
153  if (numeredEntFiniteElementPtr->getEntType() != MBEDGE)
155 
158  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
159  dataH1.dataOnEntities[MBEDGE][0].getSense() =
160  1; // set sense to 1, this is this entity
161 
162  // L2
163  if (dataH1.spacesOnEntities[MBEDGE].test(L2)) {
164  auto &data_l2 = *dataOnElement[L2];
165  CHKERR getEntityDataOrder<MBEDGE>(data_l2, L2);
166  data_l2.dataOnEntities[MBEDGE][0].getSense() =
167  1; // set sense to 1, this is this entity
168  data_l2.spacesOnEntities[MBEDGE].set(L2);
169  }
170 
171  // Hcurl
172  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
173  auto &data_curl = *dataOnElement[HCURL];
174  CHKERR getEntityDataOrder<MBEDGE>(data_curl, HCURL);
175  data_curl.dataOnEntities[MBEDGE][0].getSense() =
176  1; // set sense to 1, this is this entity
177  data_curl.spacesOnEntities[MBEDGE].set(HCURL);
178  }
179 
184 
185  // Iterate over operators
187 
189 }
190 
191 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::ForcesAndSourcesCore::calBernsteinBezierBaseFunctionsOnElement
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
Definition: ForcesAndSourcesCore.cpp:1132
MoFEM::EdgeElementForcesAndSourcesCore::setIntegrationPts
MoFEMErrorCode setIntegrationPts()
Definition: EdgeElementForcesAndSourcesCore.cpp:48
MoFEM::ForcesAndSourcesCore::getMaxRowOrder
int getMaxRowOrder() const
Get max order of approximation for field in rows.
Definition: ForcesAndSourcesCore.cpp:131
FTensor::Tensor1< double, 3 >
EntityHandle
MoFEM::EdgeElementForcesAndSourcesCore::EdgeElementForcesAndSourcesCore
EdgeElementForcesAndSourcesCore(Interface &m_field)
Definition: EdgeElementForcesAndSourcesCore.cpp:25
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::EdgeElementForcesAndSourcesCore::tangentAtGaussPts
MatrixDouble tangentAtGaussPts
Definition: EdgeElementForcesAndSourcesCore.hpp:48
MoFEM::ForcesAndSourcesCore::setGaussPts
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
Definition: ForcesAndSourcesCore.cpp:1961
MoFEM::ForcesAndSourcesCore::getMaxDataOrder
int getMaxDataOrder() const
Get max order of approximation for data fields.
Definition: ForcesAndSourcesCore.cpp:120
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
MoFEM::EdgeElementForcesAndSourcesCore::cOnn
const EntityHandle * cOnn
Definition: EdgeElementForcesAndSourcesCore.hpp:54
MoFEM::EdgeElementForcesAndSourcesCore::calculateEdgeDirection
MoFEMErrorCode calculateEdgeDirection()
Definition: EdgeElementForcesAndSourcesCore.cpp:35
MoFEM::ForcesAndSourcesCore::getElementPolynomialBase
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
Definition: ForcesAndSourcesCore.hpp:90
MoFEM::ForcesAndSourcesCore::getMaxColOrder
int getMaxColOrder() const
Get max order of approximation for field in columns.
Definition: ForcesAndSourcesCore.cpp:135
MoFEM::ForcesAndSourcesCore::calHierarchicalBaseFunctionsOnElement
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
Definition: ForcesAndSourcesCore.cpp:1112
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:49
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::EdgeElementForcesAndSourcesCore::lEngth
double & lEngth
Definition: EdgeElementForcesAndSourcesCore.hpp:51
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::EdgeElementForcesAndSourcesCore::numNodes
int numNodes
Definition: EdgeElementForcesAndSourcesCore.hpp:53
QUAD_1D_TABLE
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ForcesAndSourcesCore::coordsAtGaussPts
MatrixDouble coordsAtGaussPts
coordinated at gauss points
Definition: ForcesAndSourcesCore.hpp:544
MoFEM::ForcesAndSourcesCore::dataOnElement
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
Definition: ForcesAndSourcesCore.hpp:461
MoFEM::ForcesAndSourcesCore::createDataOnElement
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
Definition: ForcesAndSourcesCore.cpp:1337
lapack_wrap.h
MoFEM::EdgeElementForcesAndSourcesCore::cOords
VectorDouble cOords
Definition: EdgeElementForcesAndSourcesCore.hpp:56
QUAD_::npoints
int npoints
Definition: quad.h:29
MoFEM::ForcesAndSourcesCore::getRule
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
Definition: ForcesAndSourcesCore.cpp:1955
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator::loopSideFaces
MoFEMErrorCode loopSideFaces(const string fe_name, FaceElementForcesAndSourcesCoreOnSide &fe_side)
Definition: EdgeElementForcesAndSourcesCore.cpp:145
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::ForcesAndSourcesCore::getSpacesAndBaseOnEntities
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
Definition: ForcesAndSourcesCore.cpp:990
FTensor::Index< 'i', 3 >
MoFEM::EdgeElementForcesAndSourcesCore::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: EdgeElementForcesAndSourcesCore.cpp:150
MoFEM::EdgeElementForcesAndSourcesCore::tFaceOrientation
static FTensor::Tensor1< double, 3 > tFaceOrientation
Definition: EdgeElementForcesAndSourcesCore.hpp:44
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::ForcesAndSourcesCore::dataH1
EntitiesFieldData & dataH1
Definition: ForcesAndSourcesCore.hpp:471
MoFEM::EdgePolynomialBase
Calculate base functions on tetrahedral.
Definition: EdgePolynomialBase.hpp:19
QUAD_1D_TABLE_SIZE
#define QUAD_1D_TABLE_SIZE
Definition: quad.h:163
FTensor::Tensor0
Definition: Tensor0.hpp:16
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:109
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:57
MoFEM::ForcesAndSourcesCore::loopOverOperators
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
Definition: ForcesAndSourcesCore.cpp:1379
MoFEM::EdgeElementForcesAndSourcesCore::calculateCoordsAtIntegrationPts
MoFEMErrorCode calculateCoordsAtIntegrationPts()
Definition: EdgeElementForcesAndSourcesCore.cpp:101
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:383
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::EdgeElementForcesAndSourcesCore::dIrection
VectorDouble dIrection
Definition: EdgeElementForcesAndSourcesCore.hpp:55
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::ForcesAndSourcesCore::UserDataOperator::ptrFE
ForcesAndSourcesCore * ptrFE
Definition: ForcesAndSourcesCore.hpp:985
MoFEM::FaceElementForcesAndSourcesCoreOnSide
Base face element used to integrate on skeleton.
Definition: FaceElementForcesAndSourcesCoreOnSide.hpp:19
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator::setPtrFE
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
Definition: EdgeElementForcesAndSourcesCore.cpp:135
quad.h