v0.15.0
Loading...
Searching...
No Matches
TriPolynomialBase.cpp
Go to the documentation of this file.
1/** \file TriPolynomialBase.cpp
2\brief Implementation of bases on triangle
3
4*/
5
6using namespace MoFEM;
7
8MoFEMErrorCode
9TriPolynomialBase::query_interface(boost::typeindex::type_index type_index,
10 UnknownInterface **iface) const {
11 MoFEMFunctionBegin;
12 *iface = const_cast<TriPolynomialBase *>(this);
13 MoFEMFunctionReturn(0);
14}
15
31
32MoFEMErrorCode
33TriPolynomialBase::getValueH1BernsteinBezierBase(MatrixDouble &pts) {
34 MoFEMFunctionBegin;
35 EntitiesFieldData &data = cTx->dAta;
36 const std::string &field_name = cTx->fieldName;
37 int nb_gauss_pts = pts.size2();
38
39 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
40 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
41 if (!ptr)
42 ptr.reset(new MatrixInt());
43 return *ptr;
44 };
45
46 auto get_base = [field_name](auto &data) -> MatrixDouble & {
47 auto &ptr = data.getBBNSharedPtr(field_name);
48 if (!ptr)
49 ptr.reset(new MatrixDouble());
50 return *ptr;
51 };
52
53 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
54 auto &ptr = data.getBBDiffNSharedPtr(field_name);
55 if (!ptr)
56 ptr.reset(new MatrixDouble());
57 return *ptr;
58 };
59
60 auto get_alpha_by_name_ptr =
61 [](auto &data,
62 const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
63 return data.getBBAlphaIndicesSharedPtr(field_name);
64 };
65
66 auto get_base_by_name_ptr =
67 [](auto &data,
68 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
69 return data.getBBNSharedPtr(field_name);
70 };
71
72 auto get_diff_base_by_name_ptr =
73 [](auto &data,
74 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
75 return data.getBBDiffNSharedPtr(field_name);
76 };
77
78 auto get_alpha_by_order_ptr =
79 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
80 return data.getBBAlphaIndicesByOrderSharedPtr(o);
81 };
82
83 auto get_base_by_order_ptr =
84 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
85 return data.getBBNByOrderSharedPtr(o);
86 };
87
88 auto get_diff_base_by_order_ptr =
89 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
90 return data.getBBDiffNByOrderSharedPtr(o);
91 };
92
93 auto &vert_get_n = get_base(data.dataOnEntities[MBVERTEX][0]);
94 auto &vert_get_diff_n = get_diff_base(data.dataOnEntities[MBVERTEX][0]);
95 vert_get_n.resize(nb_gauss_pts, 3, false);
96 vert_get_diff_n.resize(nb_gauss_pts, 6, false);
97
98 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
99 (unsigned int)nb_gauss_pts)
100 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
101 "Base functions or nodes has wrong number of integration points "
102 "for base %s",
103 ApproximationBaseNames[NOBASE]);
104 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
105
106 auto &vertex_alpha = get_alpha(data.dataOnEntities[MBVERTEX][0]);
107 vertex_alpha.resize(3, 3, false);
108 vertex_alpha.clear();
109 for (int n = 0; n != 3; ++n)
110 vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
111
112 CHKERR BernsteinBezier::baseFunctionsTri(
113 1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
114 &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &vert_get_n(0, 0),
115 &vert_get_diff_n(0, 0));
116 for (int n = 0; n != 3; ++n) {
117 const int f = boost::math::factorial<double>(
118 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
119 for (int g = 0; g != nb_gauss_pts; ++g) {
120 vert_get_n(g, n) *= f;
121 for (int d = 0; d != 2; ++d)
122 vert_get_diff_n(g, 2 * n + d) *= f;
123 }
124 }
125
126 // edges
127 if (data.spacesOnEntities[MBEDGE].test(H1)) {
128 if (data.dataOnEntities[MBEDGE].size() != 3)
129 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
130 "Wrong size of ent data");
131
132 constexpr int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
133 for (int ee = 0; ee != 3; ++ee) {
134 auto &ent_data = data.dataOnEntities[MBEDGE][ee];
135
136 if (ent_data.getSense() == 0)
137 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
138 "Sense of the edge unknown");
139
140 const int sense = ent_data.getSense();
141 const int order = ent_data.getOrder();
142 const int nb_dofs = NBEDGE_H1(ent_data.getOrder());
143
144 if (nb_dofs) {
145 if (get_alpha_by_order_ptr(ent_data, order)) {
146 get_alpha_by_name_ptr(ent_data, field_name) =
147 get_alpha_by_order_ptr(ent_data, order);
148 get_base_by_name_ptr(ent_data, field_name) =
149 get_base_by_order_ptr(ent_data, order);
150 get_diff_base_by_name_ptr(ent_data, field_name) =
151 get_diff_base_by_order_ptr(ent_data, order);
152 } else {
153 auto &get_n = get_base(ent_data);
154 auto &get_diff_n = get_diff_base(ent_data);
155 get_n.resize(nb_gauss_pts, nb_dofs, false);
156 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
157
158 auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
159 edge_alpha.resize(nb_dofs, 3, false);
160 CHKERR BernsteinBezier::generateIndicesEdgeTri(ee, order,
161 &edge_alpha(0, 0));
162 if (sense == -1) {
163 for (int i = 0; i != edge_alpha.size1(); ++i) {
164 int a = edge_alpha(i, edges_nodes[ee][0]);
165 edge_alpha(i, edges_nodes[ee][0]) =
166 edge_alpha(i, edges_nodes[ee][1]);
167 edge_alpha(i, edges_nodes[ee][1]) = a;
168 }
169 }
170 CHKERR BernsteinBezier::baseFunctionsTri(
171 order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
172 &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
173 &get_diff_n(0, 0));
174
175 get_alpha_by_order_ptr(ent_data, order) =
176 get_alpha_by_name_ptr(ent_data, field_name);
177 get_base_by_order_ptr(ent_data, order) =
178 get_base_by_name_ptr(ent_data, field_name);
179 get_diff_base_by_order_ptr(ent_data, order) =
180 get_diff_base_by_name_ptr(ent_data, field_name);
181 }
182 }
183 }
184 } else {
185 for (int ee = 0; ee != 3; ++ee) {
186 auto &ent_data = data.dataOnEntities[MBEDGE][ee];
187 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
188 auto &get_n = get_base(ent_data);
189 auto &get_diff_n = get_diff_base(ent_data);
190 get_n.resize(nb_gauss_pts, 0, false);
191 get_diff_n.resize(nb_gauss_pts, 0, false);
192 }
193 }
194
195 if (data.spacesOnEntities[MBTRI].test(H1)) {
196 if (data.dataOnEntities[MBTRI].size() != 1)
197 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
198 "Wrong size ent of ent data");
199
200 auto &ent_data = data.dataOnEntities[MBTRI][0];
201 const int order = ent_data.getOrder();
202 const int nb_dofs = NBFACETRI_H1(order);
203 if (get_alpha_by_order_ptr(ent_data, order)) {
204 get_alpha_by_name_ptr(ent_data, field_name) =
205 get_alpha_by_order_ptr(ent_data, order);
206 get_base_by_name_ptr(ent_data, field_name) =
207 get_base_by_order_ptr(ent_data, order);
208 get_diff_base_by_name_ptr(ent_data, field_name) =
209 get_diff_base_by_order_ptr(ent_data, order);
210 } else {
211
212 auto &get_n = get_base(ent_data);
213 auto &get_diff_n = get_diff_base(ent_data);
214 get_n.resize(nb_gauss_pts, nb_dofs, false);
215 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
216 if (nb_dofs) {
217 auto &tri_alpha = get_alpha(ent_data);
218 tri_alpha.resize(nb_dofs, 3, false);
219
220 CHKERR BernsteinBezier::generateIndicesTriTri(order, &tri_alpha(0, 0));
221 CHKERR BernsteinBezier::baseFunctionsTri(
222 order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
223 &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
224 &get_diff_n(0, 0));
225
226 get_alpha_by_order_ptr(ent_data, order) =
227 get_alpha_by_name_ptr(ent_data, field_name);
228 get_base_by_order_ptr(ent_data, order) =
229 get_base_by_name_ptr(ent_data, field_name);
230 get_diff_base_by_order_ptr(ent_data, order) =
231 get_diff_base_by_name_ptr(ent_data, field_name);
232 }
233 }
234 } else {
235 auto &ent_data = data.dataOnEntities[MBTRI][0];
236 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
237 auto &get_n = get_base(ent_data);
238 auto &get_diff_n = get_diff_base(ent_data);
239 get_n.resize(nb_gauss_pts, 0, false);
240 get_diff_n.resize(nb_gauss_pts, 0, false);
241 }
242
243 MoFEMFunctionReturn(0);
244}
245
246MoFEMErrorCode TriPolynomialBase::getValueH1AinsworthBase(MatrixDouble &pts) {
247 MoFEMFunctionBegin;
248
249 EntitiesFieldData &data = cTx->dAta;
250 const FieldApproximationBase base = cTx->bAse;
251
252 if (cTx->basePolynomialsType0 == NULL)
253 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
254 "Polynomial type not set");
255
256 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
257 double *diffL, const int dim) =
258 cTx->basePolynomialsType0;
259
260 int nb_gauss_pts = pts.size2();
261
262 if (data.spacesOnEntities[MBEDGE].test(H1)) {
263 // edges
264 if (data.dataOnEntities[MBEDGE].size() != 3)
265 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
266 "expected size of data.dataOnEntities[MBEDGE] is 3");
267
268 int sense[3], order[3];
269 double *H1edgeN[3], *diffH1edgeN[3];
270 for (int ee = 0; ee < 3; ee++) {
271
272 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
273 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
274 "sense of the edge unknown");
275
276 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
277 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
278 int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
279 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
280 false);
281 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
282 2 * nb_dofs, false);
283 H1edgeN[ee] = &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
284 diffH1edgeN[ee] =
285 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
286 }
287 CHKERR H1_EdgeShapeFunctions_MBTRI(
288 sense, order,
289 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
290 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
291 H1edgeN, diffH1edgeN, nb_gauss_pts, base_polynomials);
292 }
293
294 if (data.spacesOnEntities[MBTRI].test(H1)) {
295 // face
296 if (data.dataOnEntities[MBTRI].size() != 1) {
297 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
298 "expected that size data.dataOnEntities[MBTRI] is one");
299 }
300
301 int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][0].getOrder());
302 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, nb_dofs,
303 false);
304 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
305 2 * nb_dofs, false);
306 const int face_nodes[] = {0, 1, 2};
307 CHKERR H1_FaceShapeFunctions_MBTRI(
308 face_nodes, data.dataOnEntities[MBTRI][0].getOrder(),
309 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
310 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
311 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
312 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
313 nb_gauss_pts, base_polynomials);
314 }
315
316 MoFEMFunctionReturn(0);
317}
318
319MoFEMErrorCode TriPolynomialBase::getValueL2(MatrixDouble &pts) {
320 MoFEMFunctionBegin;
321
322 switch (cTx->bAse) {
323 case AINSWORTH_LEGENDRE_BASE:
324 case AINSWORTH_LOBATTO_BASE:
325 case DEMKOWICZ_JACOBI_BASE:
326 CHKERR getValueL2AinsworthBase(pts);
327 break;
328 case AINSWORTH_BERNSTEIN_BEZIER_BASE:
329 CHKERR getValueL2BernsteinBezierBase(pts);
330 break;
331 default:
332 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
333 }
334
335 MoFEMFunctionReturn(0);
336}
337
338MoFEMErrorCode TriPolynomialBase::getValueL2AinsworthBase(MatrixDouble &pts) {
339 MoFEMFunctionBegin;
340
341 EntitiesFieldData &data = cTx->dAta;
342 const FieldApproximationBase base = cTx->bAse;
343 if (cTx->basePolynomialsType0 == NULL)
344 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
345 "Polynomial type not set");
346 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
347 double *diffL, const int dim) =
348 cTx->basePolynomialsType0;
349
350 int nb_gauss_pts = pts.size2();
351
352 data.dataOnEntities[MBTRI][0].getN(base).resize(
353 nb_gauss_pts, NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()),
354 false);
355 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(
356 nb_gauss_pts, 2 * NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()),
357 false);
358
359 CHKERR L2_Ainsworth_ShapeFunctions_MBTRI(
360 data.dataOnEntities[MBTRI][0].getOrder(),
361 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
362 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
363 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
364 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
365 nb_gauss_pts, base_polynomials);
366
367 MoFEMFunctionReturn(0);
368}
369
370MoFEMErrorCode
371TriPolynomialBase::getValueL2BernsteinBezierBase(MatrixDouble &pts) {
372 MoFEMFunctionBegin;
373
374 EntitiesFieldData &data = cTx->dAta;
375 const std::string &field_name = cTx->fieldName;
376 int nb_gauss_pts = pts.size2();
377
378 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
379 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
380 if (!ptr)
381 ptr.reset(new MatrixInt());
382 return *ptr;
383 };
384
385 auto get_base = [field_name](auto &data) -> MatrixDouble & {
386 auto &ptr = data.getBBNSharedPtr(field_name);
387 if (!ptr)
388 ptr.reset(new MatrixDouble());
389 return *ptr;
390 };
391
392 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
393 auto &ptr = data.getBBDiffNSharedPtr(field_name);
394 if (!ptr)
395 ptr.reset(new MatrixDouble());
396 return *ptr;
397 };
398
399 auto get_alpha_by_name_ptr =
400 [](auto &data,
401 const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
402 return data.getBBAlphaIndicesSharedPtr(field_name);
403 };
404
405 auto get_base_by_name_ptr =
406 [](auto &data,
407 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
408 return data.getBBNSharedPtr(field_name);
409 };
410
411 auto get_diff_base_by_name_ptr =
412 [](auto &data,
413 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
414 return data.getBBDiffNSharedPtr(field_name);
415 };
416
417 auto get_alpha_by_order_ptr =
418 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
419 return data.getBBAlphaIndicesByOrderSharedPtr(o);
420 };
421
422 auto get_base_by_order_ptr =
423 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
424 return data.getBBNByOrderSharedPtr(o);
425 };
426
427 auto get_diff_base_by_order_ptr =
428 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
429 return data.getBBDiffNByOrderSharedPtr(o);
430 };
431
432 if (data.spacesOnEntities[MBTRI].test(L2)) {
433 if (data.dataOnEntities[MBTRI].size() != 1)
434 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
435 "Wrong size ent of ent data");
436
437 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
438 (unsigned int)nb_gauss_pts)
439 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
440 "Base functions or nodes has wrong number of integration points "
441 "for base %s",
442 ApproximationBaseNames[NOBASE]);
443 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
444
445 auto &ent_data = data.dataOnEntities[MBTRI][0];
446 const int order = ent_data.getOrder();
447 const int nb_dofs = NBFACETRI_L2(order);
448
449 if (get_alpha_by_order_ptr(ent_data, order)) {
450 get_alpha_by_name_ptr(ent_data, field_name) =
451 get_alpha_by_order_ptr(ent_data, order);
452 get_base_by_name_ptr(ent_data, field_name) =
453 get_base_by_order_ptr(ent_data, order);
454 get_diff_base_by_name_ptr(ent_data, field_name) =
455 get_diff_base_by_order_ptr(ent_data, order);
456 } else {
457
458 auto &get_n = get_base(ent_data);
459 auto &get_diff_n = get_diff_base(ent_data);
460 get_n.resize(nb_gauss_pts, nb_dofs, false);
461 get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
462 if (nb_dofs) {
463
464 if (order == 0) {
465
466 if (nb_dofs != 1)
467 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
468 "Inconsistent number of DOFs");
469
470 auto &tri_alpha = get_alpha(ent_data);
471 tri_alpha.clear();
472 get_n(0, 0) = 1;
473 get_diff_n.clear();
474
475 } else {
476
477 if (nb_dofs != 3 + 3 * NBEDGE_H1(order) + NBFACETRI_H1(order))
478 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
479 "Inconsistent number of DOFs");
480
481 auto &tri_alpha = get_alpha(ent_data);
482 tri_alpha.resize(nb_dofs, 3, false);
483
484 CHKERR BernsteinBezier::generateIndicesVertexTri(order,
485 &tri_alpha(0, 0));
486
487 if (order > 1) {
488 std::array<int, 3> face_n{order, order, order};
489 std::array<int *, 3> face_edge_ptr{
490 &tri_alpha(3, 0), &tri_alpha(3 + NBEDGE_H1(order), 0),
491 &tri_alpha(3 + 2 * NBEDGE_H1(order), 0)};
492 CHKERR BernsteinBezier::generateIndicesEdgeTri(
493 face_n.data(), face_edge_ptr.data());
494 if (order > 2)
495 CHKERR BernsteinBezier::generateIndicesTriTri(
496 order, &tri_alpha(3 + 3 * NBEDGE_H1(order), 0));
497 }
498 CHKERR BernsteinBezier::baseFunctionsTri(
499 order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
500 &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
501 &get_diff_n(0, 0));
502 }
503
504 get_alpha_by_order_ptr(ent_data, order) =
505 get_alpha_by_name_ptr(ent_data, field_name);
506 get_base_by_order_ptr(ent_data, order) =
507 get_base_by_name_ptr(ent_data, field_name);
508 get_diff_base_by_order_ptr(ent_data, order) =
509 get_diff_base_by_name_ptr(ent_data, field_name);
510 }
511 }
512 } else {
513 auto &ent_data = data.dataOnEntities[MBTRI][0];
514 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
515 auto &get_n = get_base(ent_data);
516 auto &get_diff_n = get_diff_base(ent_data);
517 get_n.resize(nb_gauss_pts, 0, false);
518 get_diff_n.resize(nb_gauss_pts, 0, false);
519 }
520
521 MoFEMFunctionReturn(0);
522}
523
524MoFEMErrorCode TriPolynomialBase::getValueHdivAinsworthBase(MatrixDouble &pts) {
525 MoFEMFunctionBegin;
526
527 EntitiesFieldData &data = cTx->dAta;
528 const FieldApproximationBase base = cTx->bAse;
529 if (cTx->basePolynomialsType0 == NULL)
530 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
531 "Polynomial type not set");
532 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
533 double *diffL, const int dim) =
534 cTx->basePolynomialsType0;
535
536 int nb_gauss_pts = pts.size2();
537
538 double *phi_f_e[3];
539 double *phi_f;
540
541 N_face_edge.resize(1, 3, false);
542 N_face_bubble.resize(1, false);
543 int face_order = data.dataOnEntities[MBTRI][0].getOrder();
544 auto nb_dofs_on_face =
545 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(
546 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order)) +
547 NBFACETRI_AINSWORTH_FACE_HDIV(
548 AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
549 if (nb_dofs_on_face) {
550 auto face_edge_dofs = NBFACETRI_AINSWORTH_EDGE_HDIV(
551 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order));
552 // three edges on face
553 for (int ee = 0; ee < 3; ee++) {
554 N_face_edge(0, ee).resize(nb_gauss_pts, 3 * face_edge_dofs, false);
555 phi_f_e[ee] = &((N_face_edge(0, ee))(0, 0));
556 }
557 auto face_bubble_dofs = NBFACETRI_AINSWORTH_FACE_HDIV(
558 AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
559 N_face_bubble[0].resize(nb_gauss_pts, 3 * face_bubble_dofs, false);
560 phi_f = &*(N_face_bubble[0].data().begin());
561
562 int face_nodes[3] = {0, 1, 2};
563 if (face_edge_dofs)
564 CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
565 face_nodes,
566 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order),
567 &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, phi_f_e,
568 NULL, nb_gauss_pts, 3, base_polynomials);
569 if (face_bubble_dofs)
570 CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
571 face_nodes,
572 AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order),
573 &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, phi_f, NULL,
574 nb_gauss_pts, 3, base_polynomials);
575
576 // set shape functions into data structure
577 if (data.dataOnEntities[MBTRI].size() != 1) {
578 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
579 }
580 data.dataOnEntities[MBTRI][0].getN(base).resize(
581 nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(face_order), false);
582
583 auto max_face_order =
584 std::max(face_order,
585 AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order));
586
587 int col = 0;
588 for (int oo = 0; oo < max_face_order; oo++) {
589 if (oo < AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order))
590 for (int ee = 0; ee < 3; ee++) {
591 for (int dd = 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
592 dd < 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++, col++) {
593 for (int gg = 0; gg < nb_gauss_pts; gg++) {
594 data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
595 N_face_edge(0, ee)(gg, dd);
596 }
597 }
598 }
599
600 if (oo < AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order))
601 for (int dd = 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo);
602 dd < 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++, col++) {
603 for (int gg = 0; gg < nb_gauss_pts; gg++) {
604 data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
605 N_face_bubble[0](gg, dd);
606 }
607 }
608 }
609 }
610
611 MoFEMFunctionReturn(0);
612}
613
614MoFEMErrorCode TriPolynomialBase::getValueHdivDemkowiczBase(MatrixDouble &pts) {
615 MoFEMFunctionBegin;
616
617 EntitiesFieldData &data = cTx->dAta;
618 // set shape functions into data structure
619 if (data.dataOnEntities[MBTRI].size() != 1) {
620 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
621 }
622 const FieldApproximationBase base = cTx->bAse;
623 if (base != DEMKOWICZ_JACOBI_BASE) {
624 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
625 "This should be used only with DEMKOWICZ_JACOBI_BASE "
626 "but base is %s",
627 ApproximationBaseNames[base]);
628 }
629 int order = data.dataOnEntities[MBTRI][0].getOrder();
630 int nb_gauss_pts = pts.size2();
631 data.dataOnEntities[MBTRI][0].getN(base).resize(
632 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
633 double *phi_f = &*data.dataOnEntities[MBTRI][0].getN(base).data().begin();
634 if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
635 MoFEMFunctionReturnHot(0);
636 int face_nodes[3] = {0, 1, 2};
637 CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
638 face_nodes, order,
639 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
640 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
641 NULL, nb_gauss_pts, 3);
642
643 MoFEMFunctionReturn(0);
644}
645
661
662MoFEMErrorCode
663TriPolynomialBase::getValueHcurlAinsworthBase(MatrixDouble &pts) {
664 MoFEMFunctionBegin;
665
666 EntitiesFieldData &data = cTx->dAta;
667 const FieldApproximationBase base = cTx->bAse;
668 if (data.dataOnEntities[MBTRI].size() != 1)
669 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
670
671 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
672 double *diffL, const int dim) =
673 cTx->basePolynomialsType0;
674
675 int nb_gauss_pts = pts.size2();
676
677 // Calculation H-curl on triangle faces
678 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
679 if (data.dataOnEntities[MBEDGE].size() != 3) {
680 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
681 }
682 int sense[3], order[3];
683 double *hcurl_edge_n[3];
684 double *diff_hcurl_edge_n[3];
685 for (int ee = 0; ee < 3; ee++) {
686 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
687 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
688 "data inconsistency");
689 }
690 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
691 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
692 int nb_dofs =
693 NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
694 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
695 3 * nb_dofs, false);
696 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
697 nb_gauss_pts, 2 * 3 * nb_dofs, false);
698 hcurl_edge_n[ee] =
699 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
700 diff_hcurl_edge_n[ee] =
701 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
702 }
703 CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(
704 sense, order,
705 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
706 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
707 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
708 } else {
709 for (int ee = 0; ee < 3; ee++) {
710 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
711 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
712 false);
713 }
714 }
715
716 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
717
718 // cerr << data.dataOnEntities[MBVERTEX][0].getN(base) << endl;
719 // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
720 //
721 // face
722 if (data.dataOnEntities[MBTRI].size() != 1) {
723 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
724 }
725 int order = data.dataOnEntities[MBTRI][0].getOrder();
726 int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order);
727 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
728 false);
729 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
730 3 * 2 * nb_dofs, false);
731 // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
732 int face_nodes[] = {0, 1, 2};
733 CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(
734 face_nodes, order,
735 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
736 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
737 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
738 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
739 nb_gauss_pts, base_polynomials);
740 // cerr << data.dataOnEntities[MBTRI][0].getN(base) << endl;
741 } else {
742 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
743 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
744 }
745
746 MoFEMFunctionReturn(0);
747}
748
749MoFEMErrorCode
750TriPolynomialBase::getValueHcurlAinsworthBrokenBase(MatrixDouble &pts) {
751 MoFEMFunctionBegin;
752
753 EntitiesFieldData &data = cTx->dAta;
754 const FieldApproximationBase base = cTx->bAse;
755 if (data.dataOnEntities[MBTRI].size() != 1)
756 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
757
758 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
759 double *diffL, const int dim) =
760 cTx->basePolynomialsType0;
761
762 int nb_gauss_pts = pts.size2();
763
764 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
765
766 // face
767 if (data.dataOnEntities[MBTRI].size() != 1) {
768 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
769 }
770
771 int order = data.dataOnEntities[MBTRI][0].getOrder();
772 auto nb_edge_dofs = NBEDGE_AINSWORTH_HCURL(order);
773 int nb_dofs_face = NBFACETRI_AINSWORTH_HCURL(order);
774
775 // three faces on triangle
776 int nb_dofs = 3 * nb_edge_dofs + nb_dofs_face;
777 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
778 false);
779 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
780 3 * 2 * nb_dofs, false);
781
782 MatrixDouble edge_bases[] = {
783 MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
784 MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
785 MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false)};
786 MatrixDouble diff_edge_bases[] = {
787 MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
788 MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
789 MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false)};
790
791 std::array<double *, 3> phi{&*edge_bases[0].data().begin(),
792 &*edge_bases[1].data().begin(),
793 &*edge_bases[2].data().begin()};
794 std::array<double *, 3> diff_phi{&*diff_edge_bases[0].data().begin(),
795 &*diff_edge_bases[1].data().begin(),
796 &*diff_edge_bases[2].data().begin()};
797
798 std::array<int, 3> edge_order{order, order, order};
799 std::array<int, 3> edge_sense{1, 1, 1};
800 CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(
801 edge_sense.data(), edge_order.data(),
802 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
803 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
804 phi.data(), diff_phi.data(), nb_gauss_pts, base_polynomials);
805
806 MatrixDouble face_bases(nb_gauss_pts, 3 * nb_dofs_face, false);
807 MatrixDouble diff_face_bases(nb_gauss_pts, 3 * 2 * nb_dofs_face, false);
808 int face_nodes[] = {0, 1, 2};
809 CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(
810 face_nodes, order,
811 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
812 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
813 &*face_bases.data().begin(), &*diff_face_bases.data().begin(),
814 nb_gauss_pts, base_polynomials);
815
816 // edges
817 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_edge_base[] = {
818 getFTensor1FromPtr<3>(&*edge_bases[0].data().begin()),
819 getFTensor1FromPtr<3>(&*edge_bases[1].data().begin()),
820 getFTensor1FromPtr<3>(&*edge_bases[2].data().begin())};
821 FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_edge_diff_base[] = {
822 getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[0].data().begin()),
823 getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[1].data().begin()),
824 getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[2].data().begin())};
825
826 // face
827 auto t_vol_base = getFTensor1FromPtr<3>(&*face_bases.data().begin());
828 auto t_vol_diff_base =
829 getFTensor2HVecFromPtr<3, 2>(&*diff_face_bases.data().begin());
830
831 auto t_base = getFTensor1FromPtr<3>(
832 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin());
833 auto t_diff_base = getFTensor2HVecFromPtr<3, 2>(
834 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin());
835
836 FTENSOR_INDEX(3, i);
837 FTENSOR_INDEX(2, j);
838
839 for (int gg = 0; gg != nb_gauss_pts; gg++) {
840 for (int oo = 0; oo < order; oo++) {
841 // edges
842 for (int dd = NBEDGE_AINSWORTH_HCURL(oo);
843 dd != NBEDGE_AINSWORTH_HCURL(oo + 1); ++dd) {
844 for (int ee = 0; ee != 3; ++ee) {
845 t_base(i) = t_edge_base[ee](i);
846 ++t_base;
847 ++t_edge_base[ee];
848 t_diff_base(i, j) = t_edge_diff_base[ee](i, j);
849 ++t_diff_base;
850 ++t_edge_diff_base[ee];
851 }
852 }
853 // face
854 for (int dd = NBFACETRI_AINSWORTH_HCURL(oo);
855 dd != NBFACETRI_AINSWORTH_HCURL(oo + 1); ++dd) {
856 t_base(i) = t_vol_base(i);
857 ++t_base;
858 ++t_vol_base;
859 t_diff_base(i, j) = t_vol_diff_base(i, j);
860 ++t_diff_base;
861 ++t_vol_diff_base;
862 }
863 }
864 }
865
866 } else {
867 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
868 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
869 }
870
871 MoFEMFunctionReturn(0);
872}
873
874MoFEMErrorCode
875TriPolynomialBase::getValueHcurlDemkowiczBase(MatrixDouble &pts) {
876 MoFEMFunctionBegin;
877
878 EntitiesFieldData &data = cTx->dAta;
879 const FieldApproximationBase base = cTx->bAse;
880
881 int nb_gauss_pts = pts.size2();
882
883 // Calculation H-curl on triangle faces
884 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
885
886 if (data.dataOnEntities[MBEDGE].size() != 3)
887 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
888 "wrong number of data structures on edges, should be three but "
889 "is %ld",
890 data.dataOnEntities[MBEDGE].size());
891
892 int sense[3], order[3];
893 double *hcurl_edge_n[3];
894 double *diff_hcurl_edge_n[3];
895
896 for (int ee = 0; ee != 3; ++ee) {
897
898 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
899 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
900 "orientation (sense) of edge is not set");
901
902 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
903 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
904 int nb_dofs =
905 NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
906 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
907 3 * nb_dofs, false);
908 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
909 nb_gauss_pts, 2 * 3 * nb_dofs, false);
910
911 hcurl_edge_n[ee] =
912 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
913 diff_hcurl_edge_n[ee] =
914 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
915 }
916
917 CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(
918 sense, order,
919 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
920 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
921 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
922
923 } else {
924
925 // No DOFs on faces, resize base function matrices, indicating that no
926 // dofs on them.
927 for (int ee = 0; ee != 3; ++ee) {
928 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
929 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
930 false);
931 }
932 }
933
934 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
935
936 // face
937 if (data.dataOnEntities[MBTRI].size() != 1)
938 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
939 "No data struture to keep base functions on face");
940
941 int order = data.dataOnEntities[MBTRI][0].getOrder();
942 int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
943 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
944 false);
945 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
946 3 * 2 * nb_dofs, false);
947
948 int face_nodes[] = {0, 1, 2};
949 CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(
950 face_nodes, order,
951 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
952 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
953 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
954 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
955 nb_gauss_pts);
956
957 } else {
958
959 // No DOFs on faces, resize base function matrices, indicating that no
960 // dofs on them.
961 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
962 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
963 }
964
965 MoFEMFunctionReturn(0);
966}
967
968MoFEMErrorCode
969TriPolynomialBase::getValueHcurlDemkowiczBrokenBase(MatrixDouble &pts) {
970 MoFEMFunctionBegin;
971
972 EntitiesFieldData &data = cTx->dAta;
973 const FieldApproximationBase base = cTx->bAse;
974
975 int nb_gauss_pts = pts.size2();
976
977 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
978
979 // face
980 if (data.dataOnEntities[MBTRI].size() != 1)
981 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
982 "No data struture to keep base functions on face");
983
984 int order = data.dataOnEntities[MBTRI][0].getOrder();
985 int nb_edge_dofs = NBEDGE_DEMKOWICZ_HCURL(order);
986 int nb_volume_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
987 int nb_dofs = 3 * nb_edge_dofs + nb_volume_dofs;
988 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
989 false);
990 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
991 3 * 2 * nb_dofs, false);
992
993 MatrixDouble edge_bases[] = {
994 MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
995 MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
996 MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false)};
997 MatrixDouble diff_edge_bases[] = {
998 MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
999 MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
1000 MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false)};
1001
1002 std::array<int, 3> edge_order{order, order, order};
1003 std::array<int, 3> edge_sense{1, 1, 1};
1004 std::array<double *, 3> phi{&*edge_bases[0].data().begin(),
1005 &*edge_bases[1].data().begin(),
1006 &*edge_bases[2].data().begin()};
1007 std::array<double *, 3> diff_phi{&*diff_edge_bases[0].data().begin(),
1008 &*diff_edge_bases[1].data().begin(),
1009 &*diff_edge_bases[2].data().begin()};
1010 CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(
1011 edge_sense.data(), edge_order.data(),
1012 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1013 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1014 phi.data(), diff_phi.data(), nb_gauss_pts);
1015
1016 MatrixDouble face_bases(nb_gauss_pts, 3 * nb_volume_dofs, false);
1017 MatrixDouble diff_face_bases(nb_gauss_pts, 3 * 2 * nb_volume_dofs, false);
1018 int face_nodes[] = {0, 1, 2};
1019 CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(
1020 face_nodes, order,
1021 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1022 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1023 &*face_bases.data().begin(), &*diff_face_bases.data().begin(),
1024 nb_gauss_pts);
1025
1026 // edges
1027 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_edge_base[] = {
1028 getFTensor1FromPtr<3>(&*edge_bases[0].data().begin()),
1029 getFTensor1FromPtr<3>(&*edge_bases[1].data().begin()),
1030 getFTensor1FromPtr<3>(&*edge_bases[2].data().begin())};
1031 FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_edge_diff_base[] = {
1032 getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[0].data().begin()),
1033 getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[1].data().begin()),
1034 getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[2].data().begin())};
1035 // face
1036 auto t_vol_base = getFTensor1FromPtr<3>(&*face_bases.data().begin());
1037 auto t_vol_diff_base =
1038 getFTensor2HVecFromPtr<3, 2>(&*diff_face_bases.data().begin());
1039
1040 auto t_base = getFTensor1FromPtr<3>(
1041 &*data.dataOnEntities[MBTRI][0].getN(base).data().begin());
1042 auto t_diff_base = getFTensor2HVecFromPtr<3, 2>(
1043 &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin());
1044
1045 FTENSOR_INDEX(3, i);
1046 FTENSOR_INDEX(2, j);
1047
1048 // Fill the base functions and their derivatives
1049 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1050 for (int oo = 0; oo < order; oo++) {
1051 // edges
1052 for (int dd = NBEDGE_DEMKOWICZ_HCURL(oo);
1053 dd != NBEDGE_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1054 for (int ee = 0; ee != 3; ++ee) {
1055 t_base(i) = t_edge_base[ee](i);
1056 ++t_base;
1057 ++t_edge_base[ee];
1058 t_diff_base(i, j) = t_edge_diff_base[ee](i, j);
1059 ++t_diff_base;
1060 ++t_edge_diff_base[ee];
1061 }
1062 }
1063 // faces
1064 for (int dd = NBFACETRI_DEMKOWICZ_HCURL(oo);
1065 dd != NBFACETRI_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1066 t_base(i) = t_vol_base(i);
1067 ++t_base;
1068 ++t_vol_base;
1069 t_diff_base(i, j) = t_vol_diff_base(i, j);
1070 ++t_diff_base;
1071 ++t_vol_diff_base;
1072 }
1073 }
1074 }
1075
1076 } else {
1077
1078 // No DOFs on faces, resize base function matrices, indicating that no
1079 // dofs on them.
1080 data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
1081 data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1082 }
1083
1084 MoFEMFunctionReturn(0);
1085}
1086
1087MoFEMErrorCode TriPolynomialBase::getValueHcurl(MatrixDouble &pts) {
1088 MoFEMFunctionBegin;
1089 switch (cTx->spaceContinuity) {
1090 case CONTINUOUS:
1091
1092 switch (cTx->bAse) {
1093 case AINSWORTH_LEGENDRE_BASE:
1094 case AINSWORTH_LOBATTO_BASE:
1095 MoFEMFunctionReturnHot(getValueHcurlAinsworthBase(pts));
1096 case DEMKOWICZ_JACOBI_BASE:
1097 MoFEMFunctionReturnHot(getValueHcurlDemkowiczBase(pts));
1098 default:
1099 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1100 }
1101
1102 break;
1103
1104 case DISCONTINUOUS:
1105 switch (cTx->bAse) {
1106 case AINSWORTH_LEGENDRE_BASE:
1107 case AINSWORTH_LOBATTO_BASE:
1108 MoFEMFunctionReturnHot(getValueHcurlAinsworthBrokenBase(pts));
1109 case DEMKOWICZ_JACOBI_BASE:
1110 MoFEMFunctionReturnHot(getValueHcurlDemkowiczBrokenBase(pts));
1111 default:
1112 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1113 }
1114
1115 default:
1116 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1117 }
1118
1119 MoFEMFunctionReturn(0);
1120}
1121
1122MoFEMErrorCode
1123TriPolynomialBase::getValue(MatrixDouble &pts,
1124 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
1125 MoFEMFunctionBegin;
1126
1127 cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
1128
1129 int nb_gauss_pts = pts.size2();
1130 if (!nb_gauss_pts) {
1131 MoFEMFunctionReturnHot(0);
1132 }
1133
1134 if (pts.size1() < 2)
1135 SETERRQ(
1136 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1137 "Wrong dimension of pts, should be at least 3 rows with coordinates");
1138
1139 const FieldApproximationBase base = cTx->bAse;
1140 EntitiesFieldData &data = cTx->dAta;
1141
1142 if (base != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1143 if (cTx->copyNodeBase == LASTBASE) {
1144 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 3,
1145 false);
1146 CHKERR ShapeMBTRI(
1147 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1148 &pts(0, 0), &pts(1, 0), nb_gauss_pts);
1149 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(3, 2, false);
1150 CHKERR ShapeDiffMBTRI(
1151 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1152 } else {
1153 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1154 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1155 data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(base) =
1156 data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(cTx->copyNodeBase);
1157 }
1158 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1159 (unsigned int)nb_gauss_pts) {
1160 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1161 "Base functions or nodes has wrong number of integration points "
1162 "for base %s",
1163 ApproximationBaseNames[base]);
1164 }
1165 }
1166
1167 auto set_node_derivative_for_all_gauss_pts = [&]() {
1168 MoFEMFunctionBegin;
1169 // In linear geometry derivatives are constant,
1170 // this in expense of efficiency makes implementation
1171 // consistent between vertices and other types of entities
1172 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 6,
1173 false);
1174 for (int gg = 0; gg != nb_gauss_pts; ++gg)
1175 std::copy(Tools::diffShapeFunMBTRI.begin(),
1176 Tools::diffShapeFunMBTRI.end(),
1177 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
1178 MoFEMFunctionReturn(0);
1179 };
1180
1181 CHKERR set_node_derivative_for_all_gauss_pts();
1182
1183 switch (cTx->spaceContinuity) {
1184 case CONTINUOUS:
1185
1186 switch (cTx->sPace) {
1187 case H1:
1188 MoFEMFunctionReturnHot(getValueH1(pts));
1189 case HDIV:
1190 MoFEMFunctionReturnHot(getValueHdiv(pts));
1191 case HCURL:
1192 MoFEMFunctionReturnHot(getValueHcurl(pts));
1193 case L2:
1194 MoFEMFunctionReturnHot(getValueL2(pts));
1195 default:
1196 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1197 }
1198
1199 case DISCONTINUOUS:
1200 switch (cTx->sPace) {
1201 case HCURL:
1202 MoFEMFunctionReturnHot(getValueHcurl(pts));
1203 default:
1204 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1205 }
1206
1207 default:
1208 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1209 }
1210
1211 MoFEMFunctionReturn(0);
1212}
1213
1214MoFEMErrorCode TriPolynomialBase::setDofsSideMap(
1215
1216 const FieldSpace space, const FieldContinuity continuity,
1217 const FieldApproximationBase base, DofsSideMap &dofs_side_map
1218
1219) {
1220 MoFEMFunctionBeginHot;
1221
1222 switch (continuity) {
1223 case DISCONTINUOUS:
1224
1225 switch (space) {
1226 case HCURL:
1227 MoFEMFunctionReturnHot(
1228 setDofsSideMapHcurl(space, continuity, base, dofs_side_map));
1229 default:
1230 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
1231 FieldSpaceNames[space]);
1232 }
1233 break;
1234
1235 default:
1236 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1237 "Unknown (or not implemented) continuity");
1238 }
1239
1240 MoFEMFunctionReturnHot(0);
1241}
1242
1243MoFEMErrorCode TriPolynomialBase::setDofsSideMapHcurl(
1244 const FieldSpace space, const FieldContinuity continuity,
1245 const FieldApproximationBase base, DofsSideMap &dofs_side_map) {
1246 MoFEMFunctionBeginHot;
1247
1248 // That has to be consistent with implementation of getValueHdiv for
1249 // particular base functions.
1250
1251 auto set_ainsworth = [&dofs_side_map]() {
1252 MoFEMFunctionBegin;
1253
1254 dofs_side_map.clear();
1255
1256 int dof = 0;
1257 for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
1258
1259 // edges
1260 for (int dd = NBEDGE_AINSWORTH_HCURL(oo);
1261 dd != NBEDGE_AINSWORTH_HCURL(oo + 1); ++dd) {
1262 for (int ee = 0; ee != 3; ++ee) {
1263 dofs_side_map.insert(DofsSideMapData{MBEDGE, ee, dof});
1264 ++dof;
1265 }
1266 }
1267 // face
1268 for (int dd = NBFACETRI_AINSWORTH_HCURL(oo);
1269 dd != NBFACETRI_AINSWORTH_HCURL(oo + 1); ++dd) {
1270 dofs_side_map.insert(DofsSideMapData{MBTRI, 0, dof});
1271 ++dof;
1272 }
1273 }
1274
1275 MoFEMFunctionReturn(0);
1276 };
1277
1278 auto set_demkowicz = [&dofs_side_map]() {
1279 MoFEMFunctionBegin;
1280
1281 dofs_side_map.clear();
1282
1283 int dof = 0;
1284 for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
1285
1286 // edges
1287 for (int dd = NBEDGE_DEMKOWICZ_HCURL(oo);
1288 dd != NBEDGE_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1289 for (int ee = 0; ee != 3; ++ee) {
1290 dofs_side_map.insert(DofsSideMapData{MBEDGE, ee, dof});
1291 ++dof;
1292 }
1293 }
1294 // faces
1295 for (int dd = NBFACETRI_DEMKOWICZ_HCURL(oo);
1296 dd != NBFACETRI_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1297 dofs_side_map.insert(DofsSideMapData{MBTRI, 0, dof});
1298 ++dof;
1299 }
1300 }
1301
1302 MoFEMFunctionReturn(0);
1303 };
1304
1305 switch (continuity) {
1306 case DISCONTINUOUS:
1307 switch (base) {
1308 case AINSWORTH_LEGENDRE_BASE:
1309 case AINSWORTH_LOBATTO_BASE:
1310 MoFEMFunctionReturnHot(set_ainsworth());
1311 case DEMKOWICZ_JACOBI_BASE:
1312 MoFEMFunctionReturnHot(set_demkowicz());
1313 default:
1314 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1315 }
1316 break;
1317
1318 default:
1319 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1320 "Unknown (or not implemented) continuity");
1321 }
1322
1323 MoFEMFunctionReturnHot(0);
1324}
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2049
a
constexpr double a
Definition approx_sphere.cpp:30
FieldApproximationBase
FieldApproximationBase
approximation base
Definition definitions.h:58
LASTBASE
@ LASTBASE
Definition definitions.h:69
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
NOBASE
@ NOBASE
Definition definitions.h:59
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
FieldSpace
FieldSpace
approximation spaces
Definition definitions.h:82
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
H1
@ H1
continuous field
Definition definitions.h:85
HCURL
@ HCURL
field with continuous tangents
Definition definitions.h:86
HDIV
@ HDIV
field with continuous normal traction
Definition definitions.h:87
FieldContinuity
FieldContinuity
Field continuity.
Definition definitions.h:99
CONTINUOUS
@ CONTINUOUS
Regular field.
Definition definitions.h:100
DISCONTINUOUS
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
Definition definitions.h:101
FieldSpaceNames
static const char *const FieldSpaceNames[]
Definition definitions.h:92
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
ApproximationBaseNames
static const char *const ApproximationBaseNames[]
Definition definitions.h:72
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:18
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition fem_tools.c:194
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition fem_tools.c:182
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:108
NBFACETRI_AINSWORTH_HCURL
#define NBFACETRI_AINSWORTH_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:100
NBFACETRI_L2
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
Definition h1_hdiv_hcurl_l2.h:42
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition h1_hdiv_hcurl_l2.h:55
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:139
L2_Ainsworth_ShapeFunctions_MBTRI
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTRI(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on triangle for L2 space.
Definition l2.c:19
NBFACETRI_DEMKOWICZ_HCURL
#define NBFACETRI_DEMKOWICZ_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:109
H1_EdgeShapeFunctions_MBTRI
PetscErrorCode H1_EdgeShapeFunctions_MBTRI(int *sense, int *p, double *N, double *diffN, double *edgeN[3], double *diff_edgeN[3], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H1_EdgeShapeFunctions_MBTRI.
Definition h1.c:17
NBFACETRI_AINSWORTH_FACE_HDIV
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:131
H1_FaceShapeFunctions_MBTRI
PetscErrorCode H1_FaceShapeFunctions_MBTRI(const int *face_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition h1.c:191
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:130
NBEDGE_AINSWORTH_HCURL
#define NBEDGE_AINSWORTH_HCURL(P)
Definition h1_hdiv_hcurl_l2.h:97
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition h1_hdiv_hcurl_l2.h:60
NBFACETRI_AINSWORTH_HDIV
#define NBFACETRI_AINSWORTH_HDIV(P)
Definition h1_hdiv_hcurl_l2.h:132
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
lambda
static double lambda
Definition incompressible_elasticity.cpp:181
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::MatrixInt
UBlasMatrix< int > MatrixInt
Definition Types.hpp:76
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on teriangle.
Definition Hcurl.cpp:2105
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition Hdiv.cpp:634
MoFEM::Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition Hcurl.cpp:1249
MoFEM::Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on face.
Definition Hcurl.cpp:237
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:47
MoFEM::getFTensor2HVecFromPtr< 3, 2 >
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition Templates.hpp:957
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth nme:nme847.
Definition Hdiv.cpp:174
MoFEM::Hcurl_Demkowicz_FaceBaseFunctions_MBTRI
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(int *faces_nodes, int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Face base interior function.
Definition Hcurl.cpp:2433
MoFEM::getFTensor1FromPtr
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition Templates.hpp:838
phi
static double phi
Definition poisson_2d_dis_galerkin.cpp:30
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
g
constexpr double g
Definition shallow_wave.cpp:63
MoFEM::AinsworthOrderHooks::broken_nbfacetri_face_hdiv
static boost::function< int(int)> broken_nbfacetri_face_hdiv
Definition Hdiv.hpp:26
MoFEM::AinsworthOrderHooks::broken_nbfacetri_edge_hdiv
static boost::function< int(int)> broken_nbfacetri_edge_hdiv
Definition Hdiv.hpp:25
MoFEM::BaseFunction::DofsSideMap
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Definition BaseFunction.hpp:55
MoFEM::BernsteinBezier::baseFunctionsTri
static MoFEMErrorCode baseFunctionsTri(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
Definition BernsteinBezier.cpp:429
MoFEM::BernsteinBezier::generateIndicesVertexTri
static MoFEMErrorCode generateIndicesVertexTri(const int N, int *alpha)
Definition BernsteinBezier.cpp:90
MoFEM::BernsteinBezier::generateIndicesTriTri
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
Definition BernsteinBezier.cpp:124
MoFEM::BernsteinBezier::generateIndicesEdgeTri
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])
Definition BernsteinBezier.cpp:99
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition EntPolynomialBaseCtx.hpp:22
MoFEM::EntPolynomialBaseCtx::copyNodeBase
const FieldApproximationBase copyNodeBase
Definition EntPolynomialBaseCtx.hpp:41
MoFEM::EntPolynomialBaseCtx::basePolynomialsType0
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Definition EntPolynomialBaseCtx.hpp:27
MoFEM::EntPolynomialBaseCtx::spaceContinuity
const FieldContinuity spaceContinuity
Definition EntPolynomialBaseCtx.hpp:38
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition EntPolynomialBaseCtx.hpp:39
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition EntitiesFieldData.hpp:40
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition EntitiesFieldData.hpp:49
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition EntitiesFieldData.hpp:57
MoFEM::Field::maxBrokenDofsOrder
static constexpr int maxBrokenDofsOrder
max number of broken dofs
Definition FieldMultiIndices.hpp:282
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition Tools.hpp:104
MoFEM::TriPolynomialBase
Calculate base functions on triangle.
Definition TriPolynomialBase.hpp:16
MoFEM::TriPolynomialBase::getValueH1BernsteinBezierBase
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:33
MoFEM::TriPolynomialBase::getValueL2BernsteinBezierBase
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:371
MoFEM::TriPolynomialBase::getValueH1
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:16
MoFEM::TriPolynomialBase::getValueH1AinsworthBase
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:246
MoFEM::TriPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:319
MoFEM::TriPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition TriPolynomialBase.cpp:1123
MoFEM::TriPolynomialBase::N_face_edge
ublas::matrix< MatrixDouble > N_face_edge
Definition TriPolynomialBase.hpp:57
MoFEM::TriPolynomialBase::getValueHdivAinsworthBase
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:524
MoFEM::TriPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:614
MoFEM::TriPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition TriPolynomialBase.cpp:9
MoFEM::TriPolynomialBase::N_face_bubble
ublas::vector< MatrixDouble > N_face_bubble
Definition TriPolynomialBase.hpp:58
MoFEM::TriPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:875
MoFEM::TriPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:1087
MoFEM::TriPolynomialBase::setDofsSideMapHcurl
static MoFEMErrorCode setDofsSideMapHcurl(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &dofs_side_map)
Definition TriPolynomialBase.cpp:1243
MoFEM::TriPolynomialBase::setDofsSideMap
static MoFEMErrorCode setDofsSideMap(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &)
Set map of dof to side number.
Definition TriPolynomialBase.cpp:1214
MoFEM::TriPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition TriPolynomialBase.hpp:47
MoFEM::TriPolynomialBase::getValueHcurlAinsworthBase
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:663
MoFEM::TriPolynomialBase::getValueHcurlAinsworthBrokenBase
MoFEMErrorCode getValueHcurlAinsworthBrokenBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:750
MoFEM::TriPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:646
MoFEM::TriPolynomialBase::getValueL2AinsworthBase
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:338
MoFEM::TriPolynomialBase::getValueHcurlDemkowiczBrokenBase
MoFEMErrorCode getValueHcurlDemkowiczBrokenBase(MatrixDouble &pts)
Definition TriPolynomialBase.cpp:969
MoFEM::UnknownInterface
base class for all interface classes
Definition UnknownInterface.hpp:34
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
Generated by Doxygen 1.12.0 and hosted at University of Glasgow