v0.13.2
Loading...
Searching...
No Matches
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 0., 0., 1.};
24
26 Interface &m_field)
27 : ForcesAndSourcesCore(m_field),
28 meshPositionsFieldName("MESH_NODE_POSITIONS") {
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));
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
135MoFEMErrorCode EdgeElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
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
145EdgeElementForcesAndSourcesCore::UserDataOperator::loopSideFaces(
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
@ 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: Common.hpp:10
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.
static FTensor::Tensor1< double, 3 > tFaceOrientation
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