v0.13.1
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
10extern "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
20namespace MoFEM {
21
23 Interface &m_field)
24 : ForcesAndSourcesCore(m_field),
25 meshPositionsFieldName("MESH_NODE_POSITIONS") {
27 boost::shared_ptr<BaseFunction>(new EdgePolynomialBase());
29 "Problem with creation data on element");
30}
34 EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
35 CHKERR mField.get_moab().get_connectivity(ent, cOnn, numNodes, true);
36 cOords.resize(numNodes * 3, false);
37 CHKERR mField.get_moab().get_coords(cOnn, numNodes, &*cOords.data().begin());
38 dIrection.resize(3, false);
39 cblas_dcopy(3, &cOords[3], 1, &*dIrection.data().begin(), 1);
40 cblas_daxpy(3, -1., &cOords[0], 1, &*dIrection.data().begin(), 1);
41 lEngth = cblas_dnrm2(3, &*dIrection.data().begin(), 1);
43}
44
47 int order_data = getMaxDataOrder();
48 int order_row = getMaxRowOrder();
49 int order_col = getMaxColOrder();
50 int rule = getRule(order_row, order_col, order_data);
51
52 int nb_gauss_pts;
53 if (rule >= 0) {
54 if (rule < QUAD_1D_TABLE_SIZE) {
55 if (QUAD_1D_TABLE[rule]->dim != 1) {
56 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
57 }
58 if (QUAD_1D_TABLE[rule]->order < rule) {
59 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
60 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
61 }
62 nb_gauss_pts = QUAD_1D_TABLE[rule]->npoints;
63 gaussPts.resize(2, nb_gauss_pts, false);
64 cblas_dcopy(nb_gauss_pts, &QUAD_1D_TABLE[rule]->points[1], 2,
65 &gaussPts(0, 0), 1);
66 cblas_dcopy(nb_gauss_pts, QUAD_1D_TABLE[rule]->weights, 1,
67 &gaussPts(1, 0), 1);
68 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 2,
69 false);
70 double *shape_ptr =
71 &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
72 cblas_dcopy(2 * nb_gauss_pts, QUAD_1D_TABLE[rule]->points, 1, shape_ptr,
73 1);
74 } else {
75 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
76 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
77 nb_gauss_pts = 0;
78 }
79 } else {
80 // If rule is negative, set user defined integration points
81 CHKERR setGaussPts(order_row, order_col, order_data);
82 nb_gauss_pts = gaussPts.size2();
83 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 2,
84 false);
85 if (nb_gauss_pts) {
86 for (int gg = 0; gg != nb_gauss_pts; gg++) {
87 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 0) =
88 N_MBEDGE0(gaussPts(0, gg));
89 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE)(gg, 1) =
90 N_MBEDGE1(gaussPts(0, gg));
91 }
92 }
93 }
95}
96
101 const auto nb_gauss_pts = gaussPts.size2();
102 coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
103 tangentAtGaussPts.resize(nb_gauss_pts, 3, false);
104
105 if (nb_gauss_pts) {
107 &*gaussPts.data().begin());
109 &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
110 &coordsAtGaussPts(0, 2));
113 &tangentAtGaussPts(0, 2));
114
115 FTensor::Tensor1<double, 3> t_coords_node0(cOords[0], cOords[1], cOords[2]);
116 FTensor::Tensor1<double, 3> t_coords_node1(cOords[3], cOords[4], cOords[5]);
118 t_dir(i) = t_coords_node1(i) - t_coords_node0(i);
119
120 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
121 t_coords(i) = N_MBEDGE0(t_ksi) * t_coords_node0(i) +
122 N_MBEDGE1(t_ksi) * t_coords_node1(i);
123 t_tangent(i) = t_dir(i);
124 ++t_tangent;
125 ++t_coords;
126 ++t_ksi;
127 }
128 }
130}
131
135 if (!(ptrFE = dynamic_cast<EdgeElementForcesAndSourcesCore *>(ptr)))
136 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
137 "User operator and finite element do not work together");
139}
140
143 const string fe_name, FaceElementForcesAndSourcesCoreOnSide &fe_side) {
144 return loopSide(fe_name, &fe_side, 2);
145}
146
149
150 if (numeredEntFiniteElementPtr->getEntType() != MBEDGE)
152
155 CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
156 dataH1.dataOnEntities[MBEDGE][0].getSense() =
157 1; // set sense to 1, this is this entity
158
159 // L2
160 if (dataH1.spacesOnEntities[MBEDGE].test(L2)) {
161 auto &data_l2 = *dataOnElement[L2];
162 CHKERR getEntityDataOrder<MBEDGE>(data_l2, L2);
163 data_l2.dataOnEntities[MBEDGE][0].getSense() =
164 1; // set sense to 1, this is this entity
165 data_l2.spacesOnEntities[MBEDGE].set(L2);
166 }
167
168 // Hcurl
169 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
170 auto &data_curl = *dataOnElement[HCURL];
171 CHKERR getEntityDataOrder<MBEDGE>(data_curl, HCURL);
172 data_curl.dataOnEntities[MBEDGE][0].getSense() =
173 1; // set sense to 1, this is this entity
174 data_curl.spacesOnEntities[MBEDGE].set(HCURL);
175 }
176
181
182 // Iterate over operators
184
186}
187
188} // namespace MoFEM
@ NOBASE
Definition: definitions.h:59
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:105
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:106
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
#define QUAD_1D_TABLE_SIZE
Definition: quad.h:163
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
MoFEMErrorCode loopSideFaces(const string fe_name, FaceElementForcesAndSourcesCoreOnSide &fe_side)
MoFEMErrorCode operator()()
function is run for every finite element
Calculate base functions on tetrahedral.
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Base face element used to integrate on skeleton.
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
int getMaxRowOrder() const
Get max order of approximation for field in rows.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
MatrixDouble coordsAtGaussPts
coordinated at gauss points
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
int getMaxColOrder() const
Get max order of approximation for field in columns.
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
int npoints
Definition: quad.h:29