390 {
392
393 constexpr bool debug =
false;
394
395 constexpr int numNodes = 4;
396 constexpr int numEdges = 6;
397 constexpr int refinementLevels = 6;
398
399 auto &m_field = fe_raw_ptr->
mField;
400 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
401 auto fe_handle = fe_ptr->getFEEntityHandle();
402
403 auto set_base_quadrature = [&]() {
407 }
408 const int rule =
funRule(order_data);
409 const auto xiao_rule = IntRules::XiaoGimbutas::getTetrahedronRule(rule);
410 if (!xiao_rule) {
412 "Xiao--Gimbutas tetrahedron rule is available for polynomial "
413 "orders 0 to %d; requested %d",
414 IntRules::XiaoGimbutas::tetrahedronRuleCount, rule);
415 }
416 if (xiao_rule->numBarycentricCoordinates != 4) {
418 "wrong number of tetrahedron barycentric coordinates");
419 }
420
421 const size_t nb_gauss_pts = xiao_rule->numPoints;
422 auto &gauss_pts = fe_ptr->gaussPts;
423 gauss_pts.resize(4, nb_gauss_pts, false);
424 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[1], 4, &gauss_pts(0, 0), 1);
425 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[2], 4, &gauss_pts(1, 0), 1);
426 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[3], 4, &gauss_pts(2, 0), 1);
427 cblas_dcopy(nb_gauss_pts, xiao_rule->weights, 1, &gauss_pts(3, 0), 1);
428 auto &data = fe_ptr->dataOnElement[
H1];
429 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
430 false);
431 double *shape_ptr =
432 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
433 cblas_dcopy(4 * nb_gauss_pts, xiao_rule->points, 1, shape_ptr, 1);
435 };
436
437 CHKERR set_base_quadrature();
438
440
441 auto get_singular_nodes = [&]() {
442 int num_nodes;
443 const EntityHandle *conn;
444 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
445 true);
446 std::bitset<numNodes> singular_nodes;
447 for (auto nn = 0; nn != numNodes; ++nn) {
449 singular_nodes.set(nn);
450 } else {
451 singular_nodes.reset(nn);
452 }
453 }
454 return singular_nodes;
455 };
456
457 auto get_singular_edges = [&]() {
458 std::bitset<numEdges> singular_edges;
459 for (int ee = 0; ee != numEdges; ee++) {
460 EntityHandle edge;
461 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
463 singular_edges.set(ee);
464 } else {
465 singular_edges.reset(ee);
466 }
467 }
468 return singular_edges;
469 };
470
471 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
473 fe_ptr->gaussPts.swap(ref_gauss_pts);
474 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
475 auto &data = fe_ptr->dataOnElement[
H1];
476 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
477 double *shape_ptr =
478 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
480 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
481 nb_gauss_pts);
483 };
484
485 auto singular_nodes = get_singular_nodes();
486 if (singular_nodes.count()) {
487 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
489 CHKERR set_gauss_pts(it_map_ref_coords->second);
491 } else {
492
493 auto refine_quadrature = [&]() {
495
496 const int max_level = refinementLevels;
497 EntityHandle tet;
498
499 moab::Core moab_ref;
500 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
501 EntityHandle nodes[4];
502 for (int nn = 0; nn != 4; nn++)
503 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
504 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
507 {
508 Range tets(tet, tet);
511 tets, 1, true, edges, moab::Interface::UNION);
514 }
515
516 Range nodes_at_front;
517 for (int nn = 0; nn != numNodes; nn++) {
518 if (singular_nodes[nn]) {
519 EntityHandle ent;
520 CHKERR moab_ref.side_element(tet, 0, nn, ent);
521 nodes_at_front.insert(ent);
522 }
523 }
524
525 auto singular_edges = get_singular_edges();
526
527 EntityHandle meshset;
528 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
529 for (int ee = 0; ee != numEdges; ee++) {
530 if (singular_edges[ee]) {
531 EntityHandle ent;
532 CHKERR moab_ref.side_element(tet, 1, ee, ent);
533 CHKERR moab_ref.add_entities(meshset, &ent, 1);
534 }
535 }
536
537
539 for (int ll = 0; ll != max_level; ll++) {
542 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
544 edges);
546 CHKERR moab_ref.get_adjacencies(
547 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
548 ref_edges = intersect(ref_edges, edges);
550 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
551 ref_edges = intersect(ref_edges, ents);
554 ->getEntitiesByTypeAndRefLevel(
556 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
560 ->updateMeshsetByEntitiesChildren(meshset,
562 meshset, MBEDGE, true);
563 }
564
565
568 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
570 tets);
571
574 }
575
577 int tt = 0;
578 for (Range::iterator tit = tets.begin(); tit != tets.end();
579 tit++, tt++) {
580 int num_nodes;
581 const EntityHandle *conn;
582 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
583 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
584 }
585
586 auto &data = fe_ptr->dataOnElement[
H1];
587 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
588 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
590 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
591 int gg = 0;
592 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
593 double *tet_coords = &ref_coords(tt, 0);
595 det *= 6;
596 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
597 for (
int dd = 0;
dd != 3;
dd++) {
598 ref_gauss_pts(dd, gg) =
599 shape_n(ggg, 0) * tet_coords[3 * 0 +
dd] +
600 shape_n(ggg, 1) * tet_coords[3 * 1 +
dd] +
601 shape_n(ggg, 2) * tet_coords[3 * 2 +
dd] +
602 shape_n(ggg, 3) * tet_coords[3 * 3 +
dd];
603 }
604 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
605 }
606 }
607
608 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
610
611
614
615 TetPolynomialBase::switchCacheBaseOff<HDIV>({fe_raw_ptr});
616 TetPolynomialBase::switchCacheBaseOn<HDIV>({fe_raw_ptr});
617
619 };
620
621 CHKERR refine_quadrature();
622 }
623 }
624 }
625
627 }
#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.
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.
