v0.14.0
Loading...
Searching...
No Matches
FaceElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1/** \file FaceElementForcesAndSourcesCore.cpp
2
3\brief Implementation of face element
4
5*/
6
7namespace MoFEM {
8
10 Interface &m_field)
11 : ForcesAndSourcesCore(m_field),
12 meshPositionsFieldName("MESH_NODE_POSITIONS"), aRea(elementMeasure) {}
13
17
18 auto type = numeredEntFiniteElementPtr->getEntType();
19
23
24 auto get_ftensor_from_vec_3d = [](VectorDouble &v) {
26 &v[2]);
27 };
28
29 auto get_ftensor_n_diff = [&]() {
30 const auto &m = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
32 &m(0, 1));
33 };
34
35 auto get_ftensor_from_mat_3d = [](MatrixDouble &m) {
37 &m(0, 0), &m(0, 1), &m(0, 2));
38 };
39
40 if (type == MBTRI) {
41
42 const size_t nb_gauss_pts = gaussPts.size2();
43 normalsAtGaussPts.resize(nb_gauss_pts, 3);
44 tangentOneAtGaussPts.resize(nb_gauss_pts, 3);
45 tangentTwoAtGaussPts.resize(nb_gauss_pts, 3);
46
47 auto t_tan1 = get_ftensor_from_mat_3d(tangentOneAtGaussPts);
48 auto t_tan2 = get_ftensor_from_mat_3d(tangentTwoAtGaussPts);
49 auto t_normal = get_ftensor_from_mat_3d(normalsAtGaussPts);
50
51 auto t_n =
54 &tangentOne[2]);
56 &tangentTwo[2]);
57
58 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
59 t_normal(i) = t_n(i);
60 t_tan1(i) = t_t1(i);
61 t_tan2(i) = t_t2(i);
62 ++t_tan1;
63 ++t_tan2;
64 ++t_normal;
65 }
66
67 } else if (type == MBQUAD) {
68
70 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
71 coords.resize(num_nodes * 3, false);
72 CHKERR mField.get_moab().get_coords(conn, num_nodes,
73 &*coords.data().begin());
74
75 const size_t nb_gauss_pts = gaussPts.size2();
76 normalsAtGaussPts.resize(nb_gauss_pts, 3);
77 tangentOneAtGaussPts.resize(nb_gauss_pts, 3);
78 tangentTwoAtGaussPts.resize(nb_gauss_pts, 3);
79 normalsAtGaussPts.clear();
82
83 auto t_t1 = get_ftensor_from_mat_3d(tangentOneAtGaussPts);
84 auto t_t2 = get_ftensor_from_mat_3d(tangentTwoAtGaussPts);
85 auto t_normal = get_ftensor_from_mat_3d(normalsAtGaussPts);
86
89
90 auto t_diff = get_ftensor_n_diff();
91 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
92 auto t_coords = get_ftensor_from_vec_3d(coords);
93 for (int nn = 0; nn != num_nodes; ++nn) {
94 t_t1(i) += t_coords(i) * t_diff(N0);
95 t_t2(i) += t_coords(i) * t_diff(N1);
96 ++t_diff;
97 ++t_coords;
98 }
99 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
100
101 ++t_t1;
102 ++t_t2;
103 ++t_normal;
104 }
105 } else {
106 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
107 "Element type not implemented");
108 }
109
111}
112
115
117 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
118 coords.resize(num_nodes * 3, false);
119 CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
120 nOrmal.resize(3, false);
121 tangentOne.resize(3, false);
122 tangentTwo.resize(3, false);
123
124 auto calc_normal = [&](const double *diff_ptr) {
127 &coords[0], &coords[1], &coords[2]);
129 &nOrmal[0], &nOrmal[1], &nOrmal[2]);
131 &tangentOne[0], &tangentOne[1], &tangentOne[2]);
133 &tangentTwo[0], &tangentTwo[1], &tangentTwo[2]);
135 diff_ptr, &diff_ptr[1]);
136
140
143 t_t1(i) = 0;
144 t_t2(i) = 0;
145
146 for (int nn = 0; nn != num_nodes; ++nn) {
147 t_t1(i) += t_coords(i) * t_diff(N0);
148 t_t2(i) += t_coords(i) * t_diff(N1);
149 ++t_coords;
150 ++t_diff;
151 }
152 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
153 aRea = sqrt(t_normal(i) * t_normal(i));
155 };
156
157 const double *diff_ptr;
158 switch (numeredEntFiniteElementPtr->getEntType()) {
159 case MBTRI:
160 diff_ptr = Tools::diffShapeFunMBTRI.data();
161 CHKERR calc_normal(diff_ptr);
162 // FIXME: Normal should be divided not the area for triangle!!
163 aRea /= 2;
164 break;
165 case MBQUAD:
166 diff_ptr = Tools::diffShapeFunMBQUADAtCenter.data();
167 CHKERR calc_normal(diff_ptr);
168 break;
169 default:
170 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
171 "Element type not implemented");
172 }
173
175}
176
179 // Set integration points
180 int order_data = getMaxDataOrder();
181 int order_row = getMaxRowOrder();
182 int order_col = getMaxColOrder();
183
184 const auto type = numeredEntFiniteElementPtr->getEntType();
185
186 auto get_rule_by_type = [&]() {
187 switch (type) {
188 case MBQUAD:
189 return getRule(order_row + 1, order_col + 1, order_data + 1);
190 default:
191 return getRule(order_row, order_col, order_data);
192 }
193 };
194
195 const int rule = get_rule_by_type();
196
197 auto set_integration_pts_for_tri = [&]() {
199 if (rule < QUAD_2D_TABLE_SIZE) {
200 if (QUAD_2D_TABLE[rule]->dim != 2) {
201 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
202 }
203 if (QUAD_2D_TABLE[rule]->order < rule) {
204 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
205 "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
206 }
207 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
208 gaussPts.resize(3, nb_gauss_pts, false);
209 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
210 &gaussPts(0, 0), 1);
211 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
212 &gaussPts(1, 0), 1);
213 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
214 &gaussPts(2, 0), 1);
215 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
216 false);
217 double *shape_ptr =
218 &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
219 cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
220 1);
221 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
222 std::copy(
224 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
225 } else {
226 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
227 "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
228 }
230 };
231
232 auto calc_base_for_tri = [&]() {
234 const size_t nb_gauss_pts = gaussPts.size2();
235 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
236 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
237 base.resize(nb_gauss_pts, 3, false);
238 diff_base.resize(3, 2, false);
239 CHKERR ShapeMBTRI(&*base.data().begin(), &gaussPts(0, 0), &gaussPts(1, 0),
240 nb_gauss_pts);
241 std::copy(
243 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
245 };
246
247 auto calc_base_for_quad = [&]() {
249 const size_t nb_gauss_pts = gaussPts.size2();
250 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
251 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
252 base.resize(nb_gauss_pts, 4, false);
253 diff_base.resize(nb_gauss_pts, 8, false);
254 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
255 const double ksi = gaussPts(0, gg);
256 const double zeta = gaussPts(1, gg);
257 base(gg, 0) = N_MBQUAD0(ksi, zeta);
258 base(gg, 1) = N_MBQUAD1(ksi, zeta);
259 base(gg, 2) = N_MBQUAD2(ksi, zeta);
260 base(gg, 3) = N_MBQUAD3(ksi, zeta);
261 diff_base(gg, 0) = diffN_MBQUAD0x(zeta);
262 diff_base(gg, 1) = diffN_MBQUAD0y(ksi);
263 diff_base(gg, 2) = diffN_MBQUAD1x(zeta);
264 diff_base(gg, 3) = diffN_MBQUAD1y(ksi);
265 diff_base(gg, 4) = diffN_MBQUAD2x(zeta);
266 diff_base(gg, 5) = diffN_MBQUAD2y(ksi);
267 diff_base(gg, 6) = diffN_MBQUAD3x(zeta);
268 diff_base(gg, 7) = diffN_MBQUAD3y(ksi);
269 }
271 };
272
273 if (rule >= 0) {
274 switch (type) {
275 case MBTRI:
276 CHKERR set_integration_pts_for_tri();
277 break;
278 case MBQUAD:
280 rule);
281 CHKERR calc_base_for_quad();
282 break;
283 default:
284 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
285 "Element type not implemented: %d", type);
286 }
287
288 } else {
289 // If rule is negative, set user defined integration points
290 CHKERR setGaussPts(order_row, order_col, order_data);
291 const size_t nb_gauss_pts = gaussPts.size2();
292 if (nb_gauss_pts) {
293 switch (type) {
294 case MBTRI:
295 CHKERR calc_base_for_tri();
296 break;
297 case MBQUAD:
298 CHKERR calc_base_for_quad();
299 break;
300 default:
301 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
302 "Element type not implemented: %d", type);
303 }
304 }
305 }
307}
308
312 // Get spaces order/base and sense of entities.
313
315
316 // H1
317 if (dataH1.spacesOnEntities[MBEDGE].test(H1)) {
318 CHKERR getEntitySense<MBEDGE>(dataH1);
319 CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
320 }
321 if (dataH1.spacesOnEntities[MBTRI].test(H1)) {
322 CHKERR getEntitySense<MBTRI>(dataH1);
323 CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
324 }
325 if (dataH1.spacesOnEntities[MBQUAD].test(H1)) {
326 CHKERR getEntitySense<MBQUAD>(dataH1);
327 CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
328 }
329
330 // Hcurl
331 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
332 CHKERR getEntitySense<MBEDGE>(dataHcurl);
333 CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
334 dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
335 }
336 if (dataH1.spacesOnEntities[MBTRI].test(HCURL)) {
337 CHKERR getEntitySense<MBTRI>(dataHcurl);
338 CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
339 dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
340 }
341 if (dataH1.spacesOnEntities[MBQUAD].test(HCURL)) {
342 CHKERR getEntitySense<MBQUAD>(dataHcurl);
343 CHKERR getEntityDataOrder<MBQUAD>(dataHcurl, HCURL);
344 dataHcurl.spacesOnEntities[MBQUAD].set(HCURL);
345 }
346
347 // Hdiv
348 if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
349 CHKERR getEntitySense<MBTRI>(dataHdiv);
350 CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
351 dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
352 }
353 if (dataH1.spacesOnEntities[MBQUAD].test(HDIV)) {
354 CHKERR getEntitySense<MBQUAD>(dataHdiv);
355 CHKERR getEntityDataOrder<MBQUAD>(dataHdiv, HDIV);
356 dataHdiv.spacesOnEntities[MBQUAD].set(HDIV);
357 }
358
359 // L2
360 if (dataH1.spacesOnEntities[MBTRI].test(L2)) {
361 CHKERR getEntitySense<MBTRI>(dataL2);
362 CHKERR getEntityDataOrder<MBTRI>(dataL2, L2);
363 dataL2.spacesOnEntities[MBTRI].set(L2);
364 }
365 if (dataH1.spacesOnEntities[MBQUAD].test(L2)) {
366 CHKERR getEntitySense<MBQUAD>(dataL2);
367 CHKERR getEntityDataOrder<MBQUAD>(dataL2, L2);
368 dataL2.spacesOnEntities[MBQUAD].set(L2);
369 }
370
372}
373
377
378 const size_t nb_nodes =
379 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).size2();
380 double *shape_functions =
381 &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
382 const size_t nb_gauss_pts = gaussPts.size2();
383 coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
384 for (int gg = 0; gg != nb_gauss_pts; ++gg)
385 for (int dd = 0; dd != 3; ++dd)
386 coordsAtGaussPts(gg, dd) = cblas_ddot(
387 nb_nodes, &shape_functions[nb_nodes * gg], 1, &coords[dd], 3);
388
390}
391
392MoFEMErrorCode FaceElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
395 if (!(ptrFE = dynamic_cast<FaceElementForcesAndSourcesCore *>(ptr)))
396 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
397 "User operator and finite element do not work together");
399}
400
403
404 const auto type = numeredEntFiniteElementPtr->getEntType();
405 if (type != lastEvaluatedElementEntityType) {
406 switch (type) {
407 case MBTRI:
409 boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
410 break;
411 case MBQUAD:
413 boost::shared_ptr<BaseFunction>(new QuadPolynomialBase());
414 break;
415 default:
417 }
420 }
421
422 // Calculate normal and tangent vectors for face geometry
425
427 if (gaussPts.size2() == 0)
429
434
435 // Iterate over operators
437
439}
440
442FaceElementForcesAndSourcesCore::UserDataOperator::loopSideVolumes(
443 const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method) {
444 return loopSide(fe_name, &fe_method, 3);
445}
446
447} // namespace MoFEM
static Number< 1 > N1
static Number< 0 > N0
@ NOBASE
Definition: definitions.h:59
#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
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#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
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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 diffN_MBQUAD2y(x)
Definition: fem_tools.h:66
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:63
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:61
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:64
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:68
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:62
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:57
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:67
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:65
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:59
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:58
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
double zeta
Viscous hardening.
Definition: plastic.cpp:170
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
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
virtual MoFEMErrorCode calculateAreaAndNormal()
Calculate element area and normal of the face.
MoFEMErrorCode operator()()
function is run for every finite element
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
virtual MoFEMErrorCode calculateAreaAndNormalAtIntegrationPts()
Calculate element area and normal of the face at integration points.
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
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
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.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
Calculate base functions on triangle.
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:577
static constexpr std::array< double, 8 > diffShapeFunMBQUADAtCenter
Definition: Tools.hpp:212
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
Calculate base functions on triangle.
Base volume element used to integrate on skeleton.
int npoints
Definition: quad.h:29