406 {
408
409 constexpr bool debug =
false;
410
411 constexpr int numNodes = 4;
412 constexpr int numEdges = 6;
413 constexpr int refinementLevels = 6;
414
415 auto &m_field = fe_raw_ptr->
mField;
416 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
417 auto fe_handle = fe_ptr->getFEEntityHandle();
418
419 auto set_base_quadrature = [&]() {
421 int rule =
funRule(order_data);
425 "wrong dimension");
426 }
430 }
432 auto &gauss_pts = fe_ptr->gaussPts;
433 gauss_pts.resize(4, nb_gauss_pts, false);
434 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
435 &gauss_pts(0, 0), 1);
436 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
437 &gauss_pts(1, 0), 1);
438 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
439 &gauss_pts(2, 0), 1);
441 &gauss_pts(3, 0), 1);
442 auto &data = fe_ptr->dataOnElement[
H1];
443 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
444 false);
445 double *shape_ptr =
446 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
447 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
448 1);
449 } else {
452 }
454 };
455
456 CHKERR set_base_quadrature();
457
459
460 auto get_singular_nodes = [&]() {
461 int num_nodes;
463 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
464 true);
465 std::bitset<numNodes> singular_nodes;
466 for (auto nn = 0; nn != numNodes; ++nn) {
468 singular_nodes.set(nn);
469 } else {
470 singular_nodes.reset(nn);
471 }
472 }
473 return singular_nodes;
474 };
475
476 auto get_singular_edges = [&]() {
477 std::bitset<numEdges> singular_edges;
478 for (int ee = 0; ee != numEdges; ee++) {
480 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
482 singular_edges.set(ee);
483 } else {
484 singular_edges.reset(ee);
485 }
486 }
487 return singular_edges;
488 };
489
490 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
492 fe_ptr->gaussPts.swap(ref_gauss_pts);
493 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
494 auto &data = fe_ptr->dataOnElement[
H1];
495 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
496 double *shape_ptr =
497 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
499 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
500 nb_gauss_pts);
502 };
503
504 auto singular_nodes = get_singular_nodes();
505 if (singular_nodes.count()) {
506 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
508 CHKERR set_gauss_pts(it_map_ref_coords->second);
510 } else {
511
512 auto refine_quadrature = [&]() {
514
515 const int max_level = refinementLevels;
517
518 moab::Core moab_ref;
519 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
521 for (int nn = 0; nn != 4; nn++)
522 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
523 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
526 {
527 Range tets(tet, tet);
530 tets, 1, true, edges, moab::Interface::UNION);
533 }
534
535 Range nodes_at_front;
536 for (int nn = 0; nn != numNodes; nn++) {
537 if (singular_nodes[nn]) {
539 CHKERR moab_ref.side_element(tet, 0, nn, ent);
540 nodes_at_front.insert(ent);
541 }
542 }
543
544 auto singular_edges = get_singular_edges();
545
547 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
548 for (int ee = 0; ee != numEdges; ee++) {
549 if (singular_edges[ee]) {
551 CHKERR moab_ref.side_element(tet, 1, ee, ent);
552 CHKERR moab_ref.add_entities(meshset, &ent, 1);
553 }
554 }
555
556
558 for (int ll = 0; ll != max_level; ll++) {
561 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
563 edges);
565 CHKERR moab_ref.get_adjacencies(
566 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
567 ref_edges = intersect(ref_edges, edges);
569 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
570 ref_edges = intersect(ref_edges, ents);
573 ->getEntitiesByTypeAndRefLevel(
575 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
579 ->updateMeshsetByEntitiesChildren(meshset,
581 meshset, MBEDGE, true);
582 }
583
584
587 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
589 tets);
590
593 }
594
596 int tt = 0;
597 for (Range::iterator tit = tets.begin(); tit != tets.end();
598 tit++, tt++) {
599 int num_nodes;
601 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
602 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
603 }
604
605 auto &data = fe_ptr->dataOnElement[
H1];
606 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
607 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
609 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
610 int gg = 0;
611 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
612 double *tet_coords = &ref_coords(tt, 0);
614 det *= 6;
615 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
616 for (
int dd = 0;
dd != 3;
dd++) {
617 ref_gauss_pts(dd, gg) =
618 shape_n(ggg, 0) * tet_coords[3 * 0 +
dd] +
619 shape_n(ggg, 1) * tet_coords[3 * 1 +
dd] +
620 shape_n(ggg, 2) * tet_coords[3 * 2 +
dd] +
621 shape_n(ggg, 3) * tet_coords[3 * 3 +
dd];
622 }
623 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
624 }
625 }
626
627 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
629
630
633
634 TetPolynomialBase::switchCacheBaseOff<HDIV>({fe_raw_ptr});
635 TetPolynomialBase::switchCacheBaseOn<HDIV>({fe_raw_ptr});
636
638 };
639
640 CHKERR refine_quadrature();
641 }
642 }
643 }
644
646 }
#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)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
#define QUAD_3D_TABLE_SIZE
static QUAD *const QUAD_3D_TABLE[]
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.
