408 {
410
411 constexpr bool debug =
false;
412
413 constexpr int numNodes = 4;
414 constexpr int numEdges = 6;
415 constexpr int refinementLevels = 6;
416
417 auto &m_field = fe_raw_ptr->
mField;
418 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
419 auto fe_handle = fe_ptr->getFEEntityHandle();
420
421 auto set_base_quadrature = [&]() {
425 }
426 const int rule =
funRule(order_data);
427 const auto xiao_rule =
429 if (!xiao_rule) {
431 "Xiao--Gimbutas tetrahedron rule is available for polynomial "
432 "orders 0 to %d; requested %d",
434 }
435 if (xiao_rule->numBarycentricCoordinates != 4) {
437 "wrong number of tetrahedron barycentric coordinates");
438 }
439
440 const size_t nb_gauss_pts = xiao_rule->numPoints;
441 auto &gauss_pts = fe_ptr->gaussPts;
442 gauss_pts.resize(4, nb_gauss_pts, false);
443 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[1], 4, &gauss_pts(0, 0),
444 1);
445 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[2], 4, &gauss_pts(1, 0),
446 1);
447 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[3], 4, &gauss_pts(2, 0),
448 1);
449 cblas_dcopy(nb_gauss_pts, xiao_rule->weights, 1, &gauss_pts(3, 0), 1);
450 auto &data = fe_ptr->dataOnElement[
H1];
451 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
452 false);
453 double *shape_ptr =
454 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
455 cblas_dcopy(4 * nb_gauss_pts, xiao_rule->points, 1, shape_ptr, 1);
457 };
458
459 CHKERR set_base_quadrature();
460
462
463 auto get_singular_nodes = [&]() {
464 int num_nodes;
465 const EntityHandle *conn;
466 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
467 true);
468 std::bitset<numNodes> singular_nodes;
469 for (auto nn = 0; nn != numNodes; ++nn) {
471 singular_nodes.set(nn);
472 } else {
473 singular_nodes.reset(nn);
474 }
475 }
476 return singular_nodes;
477 };
478
479 auto get_singular_edges = [&]() {
480 std::bitset<numEdges> singular_edges;
481 for (int ee = 0; ee != numEdges; ee++) {
482 EntityHandle edge;
483 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
485 singular_edges.set(ee);
486 } else {
487 singular_edges.reset(ee);
488 }
489 }
490 return singular_edges;
491 };
492
493 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
495 fe_ptr->gaussPts.swap(ref_gauss_pts);
496 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
497 auto &data = fe_ptr->dataOnElement[
H1];
498 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
499 double *shape_ptr =
500 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
502 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
503 nb_gauss_pts);
505 };
506
507 auto singular_nodes = get_singular_nodes();
508 if (singular_nodes.count()) {
509 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
511 CHKERR set_gauss_pts(it_map_ref_coords->second);
513 } else {
514
515 auto refine_quadrature = [&]() {
517
518 const int max_level = refinementLevels;
519 EntityHandle tet;
520
521 moab::Core moab_ref;
522 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
523 EntityHandle nodes[4];
524 for (int nn = 0; nn != 4; nn++)
525 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
526 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
529 {
530 Range tets(tet, tet);
533 tets, 1, true, edges, moab::Interface::UNION);
536 }
537
538 Range nodes_at_front;
539 for (int nn = 0; nn != numNodes; nn++) {
540 if (singular_nodes[nn]) {
541 EntityHandle ent;
542 CHKERR moab_ref.side_element(tet, 0, nn, ent);
543 nodes_at_front.insert(ent);
544 }
545 }
546
547 auto singular_edges = get_singular_edges();
548
549 EntityHandle meshset;
550 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
551 for (int ee = 0; ee != numEdges; ee++) {
552 if (singular_edges[ee]) {
553 EntityHandle ent;
554 CHKERR moab_ref.side_element(tet, 1, ee, ent);
555 CHKERR moab_ref.add_entities(meshset, &ent, 1);
556 }
557 }
558
559
561 for (int ll = 0; ll != max_level; ll++) {
564 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
566 edges);
568 CHKERR moab_ref.get_adjacencies(
569 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
570 ref_edges = intersect(ref_edges, edges);
572 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
573 ref_edges = intersect(ref_edges, ents);
576 ->getEntitiesByTypeAndRefLevel(
578 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
582 ->updateMeshsetByEntitiesChildren(meshset,
584 meshset, MBEDGE, true);
585 }
586
587
590 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
592 tets);
593
596 }
597
599 int tt = 0;
600 for (Range::iterator tit = tets.begin(); tit != tets.end();
601 tit++, tt++) {
602 int num_nodes;
603 const EntityHandle *conn;
604 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
605 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
606 }
607
608 auto &data = fe_ptr->dataOnElement[
H1];
609 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
610 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
612 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
613 int gg = 0;
614 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
615 double *tet_coords = &ref_coords(tt, 0);
617 det *= 6;
618 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
619 for (
int dd = 0;
dd != 3;
dd++) {
620 ref_gauss_pts(dd, gg) =
621 shape_n(ggg, 0) * tet_coords[3 * 0 +
dd] +
622 shape_n(ggg, 1) * tet_coords[3 * 1 +
dd] +
623 shape_n(ggg, 2) * tet_coords[3 * 2 +
dd] +
624 shape_n(ggg, 3) * tet_coords[3 * 3 +
dd];
625 }
626 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
627 }
628 }
629
630 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
632
633
636
637 TetPolynomialBase::switchCacheBaseOff<HDIV>({fe_raw_ptr});
638 TetPolynomialBase::switchCacheBaseOn<HDIV>({fe_raw_ptr});
639
641 };
642
643 CHKERR refine_quadrature();
644 }
645 }
646 }
647
649 }
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
const Rule * getTetrahedronRule(const int order)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
static PetscBool setSingularity
static std::map< long int, MatrixDouble > mapRefCoords
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
Mesh refinement interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
