243 {
245
248 const int nb_gauss_pts = pts.size2();
249
251 (unsigned int)nb_gauss_pts)
253 "Base functions or nodes has wrong number of integration points "
254 "for base %s",
257
259 auto &ptr = data.getBBAlphaIndicesSharedPtr(
field_name);
260 if (!ptr)
262 return *ptr;
263 };
264
267 if (!ptr)
269 return *ptr;
270 };
271
273 auto &ptr = data.getBBDiffNSharedPtr(
field_name);
274 if (!ptr)
276 return *ptr;
277 };
278
279 auto get_alpha_by_name_ptr =
280 [](auto &data,
281 const std::string &
field_name) -> boost::shared_ptr<MatrixInt> & {
282 return data.getBBAlphaIndicesSharedPtr(
field_name);
283 };
284
285 auto get_base_by_name_ptr =
286 [](auto &data,
287 const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
289 };
290
291 auto get_diff_base_by_name_ptr =
292 [](auto &data,
293 const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
295 };
296
297 auto get_alpha_by_order_ptr =
298 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
299 return data.getBBAlphaIndicesByOrderSharedPtr(o);
300 };
301
302 auto get_base_by_order_ptr =
303 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
304 return data.getBBNByOrderSharedPtr(o);
305 };
306
307 auto get_diff_base_by_order_ptr =
308 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
309 return data.getBBDiffNByOrderSharedPtr(o);
310 };
311
312
313
314
315 auto make_own = [
field_name](
auto &ent_data) {
316 auto &n_ptr = ent_data.getBBNSharedPtr(
field_name);
317 if (!n_ptr)
318 return;
319 for (auto &ptr : ent_data.getBBNByOrderArray()) {
320 if (ptr && ptr.get() == n_ptr.get()) {
321 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
322 ent_data.getBBDiffNSharedPtr(
field_name).reset();
323 n_ptr.reset();
324 return;
325 }
326 }
327 };
328
329
330
331 auto set_empty = [&](auto &ent_data) {
332 make_own(ent_data);
333 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
334 get_base(ent_data).resize(nb_gauss_pts, 0, false);
335 get_diff_base(ent_data).resize(nb_gauss_pts, 0, false);
336 };
337
338
339
340
341
342
343
344 auto check_uniform_order = [&]() {
346 auto &node_order = data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder();
347 if (node_order.size() != 4)
349 "Wrong number of vertex orders for BB base: %d",
350 static_cast<int>(node_order.size()));
352 for (
int n = 0;
n != 4; ++
n)
354 auto max_order_on = [&](
const EntityType
type) {
356 for (
auto &ent_data : data.dataOnEntities[
type])
358 };
359 max_order_on(MBEDGE);
360 max_order_on(MBTRI);
361 max_order_on(MBTET);
364 "Bernstein-Bezier H1 base requires approximation order >= 1");
365 if (
static_cast<size_t>(
order) >=
368 "Bernstein-Bezier order %d exceeds the supported maximum %d",
370 static_cast<int>(
372 1);
373 for (
int n = 0;
n != 4; ++
n)
374 if (node_order[
n] < 1)
376 "Bernstein-Bezier H1 base requires vertex order between 1 "
377 "and the element order %d; vertex %d has order %d",
379 auto check_type = [&](
const EntityType
type,
const int nb_dofs_at_order) {
383 int side = 0;
384 for (
auto &ent_data : data.dataOnEntities[
type]) {
385 if (ent_data.getOrder() !=
order)
387 "Bernstein-Bezier base in H1 space requires uniform "
388 "approximation order on the element: %s %d has order %d, "
389 "element order is %d. BB functions are not hierarchical, "
390 "so mixed degrees do not span a complete polynomial space. "
391 "Set one order for the whole H1 field (L2 fields are not "
392 "restricted).",
393 moab::CN::EntityTypeName(
type), side, ent_data.getOrder(),
395 ++side;
396 }
398 };
403 };
404 CHKERR check_uniform_order();
405
407 auto &vertex_alpha = get_alpha(vert_ent_data);
408 vertex_alpha.resize(4, 4, false);
409 vertex_alpha.clear();
410 for (
int n = 0;
n != 4; ++
n)
412
413 auto &vert_get_n = get_base(vert_ent_data);
414 auto &vert_get_diff_n = get_diff_base(vert_ent_data);
415 vert_get_n.resize(nb_gauss_pts, 4, false);
416 vert_get_diff_n.resize(nb_gauss_pts, 12, false);
418 1,
lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
420 &vert_get_diff_n(0, 0));
421 for (
int n = 0;
n != 4; ++
n) {
422 const double f = boost::math::factorial<double>(
424 for (
int g = 0;
g != nb_gauss_pts; ++
g) {
425 vert_get_n(
g,
n) *=
f;
426 for (
int d = 0;
d != 3; ++
d)
427 vert_get_diff_n(
g, 3 *
n + d) *=
f;
428 }
429 }
430
431
435 "Wrong size of ent data");
436
437 constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
438 {0, 3}, {1, 3}, {2, 3}};
439 for (int ee = 0; ee != 6; ++ee) {
441 const int sense = ent_data.getSense();
442 if (sense == 0)
444 "Sense of the edge unknown");
445 const int order = ent_data.getOrder();
447
448 if (nb_dofs) {
449 if (get_alpha_by_order_ptr(ent_data,
order)) {
451 get_alpha_by_order_ptr(ent_data,
order);
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 make_own(ent_data);
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, 3 * nb_dofs, false);
462
463 auto &edge_alpha = get_alpha(ent_data);
464 edge_alpha.resize(nb_dofs, 4, false);
466 &edge_alpha(0, 0));
467 if (sense == -1) {
468 for (
int i = 0;
i != edge_alpha.size1(); ++
i) {
469 int a = edge_alpha(
i, edges_nodes[ee][0]);
470 edge_alpha(
i, edges_nodes[ee][0]) =
471 edge_alpha(
i, edges_nodes[ee][1]);
472 edge_alpha(
i, edges_nodes[ee][1]) =
a;
473 }
474 }
476 order,
lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
478 &get_diff_n(0, 0));
479
480 get_alpha_by_order_ptr(ent_data,
order) =
482 get_base_by_order_ptr(ent_data,
order) =
484 get_diff_base_by_order_ptr(ent_data,
order) =
485 get_diff_base_by_name_ptr(ent_data,
field_name);
486 }
487 } else {
488 set_empty(ent_data);
489 }
490 }
491 } else {
492 for (int ee = 0; ee != 6; ++ee)
494 }
495
496
500 "Wrong size of ent data");
505
506 for (int ff = 0; ff != 4; ++ff) {
508 const int order = ent_data.getOrder();
510
511 if (nb_dofs) {
512 if (get_alpha_by_order_ptr(ent_data,
order)) {
514 get_alpha_by_order_ptr(ent_data,
order);
516 get_base_by_order_ptr(ent_data,
order);
517 get_diff_base_by_name_ptr(ent_data,
field_name) =
518 get_diff_base_by_order_ptr(ent_data,
order);
519 } else {
520
521 make_own(ent_data);
522 auto &get_n = get_base(ent_data);
523 auto &get_diff_n = get_diff_base(ent_data);
524 get_n.resize(nb_gauss_pts, nb_dofs, false);
525 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
526
527 auto &face_alpha = get_alpha(ent_data);
528 face_alpha.resize(nb_dofs, 4, false);
529
531 &face_alpha(0, 0));
532 senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(),
false);
534 constexpr int tri_nodes[4][3] = {
535 {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
536 for (
int d = 0;
d != nb_dofs; ++
d)
537 for (
int n = 0;
n != 3; ++
n)
539 face_alpha(d, tri_nodes[ff][
n]);
542 order,
lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
544 &get_diff_n(0, 0));
545
546 get_alpha_by_order_ptr(ent_data,
order) =
548 get_base_by_order_ptr(ent_data,
order) =
550 get_diff_base_by_order_ptr(ent_data,
order) =
551 get_diff_base_by_name_ptr(ent_data,
field_name);
552 }
553 } else {
554 set_empty(ent_data);
555 }
556 }
557 } else {
558 for (int ff = 0; ff != 4; ++ff)
560 }
561
565 "Wrong size ent of ent data");
566
568 const int order = ent_data.getOrder();
570 if (get_alpha_by_order_ptr(ent_data,
order)) {
572 get_alpha_by_order_ptr(ent_data,
order);
574 get_base_by_order_ptr(ent_data,
order);
575 get_diff_base_by_name_ptr(ent_data,
field_name) =
576 get_diff_base_by_order_ptr(ent_data,
order);
577 } else {
578
579 make_own(ent_data);
580 auto &get_n = get_base(ent_data);
581 auto &get_diff_n = get_diff_base(ent_data);
582 get_n.resize(nb_gauss_pts, nb_dofs, false);
583 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
584 if (nb_dofs) {
585 auto &tet_alpha = get_alpha(ent_data);
586 tet_alpha.resize(nb_dofs, 4, false);
587
590 order,
lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
592 &get_diff_n(0, 0));
593
594 get_alpha_by_order_ptr(ent_data,
order) =
596 get_base_by_order_ptr(ent_data,
order) =
598 get_diff_base_by_order_ptr(ent_data,
order) =
599 get_diff_base_by_name_ptr(ent_data,
field_name);
600 } else {
601 ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
602 }
603 }
604 } else {
606 }
607
609}
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
UBlasMatrix< double > MatrixDouble
UBlasMatrix< int > MatrixInt
constexpr auto field_name
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
static MoFEMErrorCode baseFunctionsTet(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)
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
const std::string fieldName
static constexpr size_t MaxBernsteinBezierOrder
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Public Member Functions inherited from MoFEM::BaseFunction