389 {
391
392 constexpr bool debug =
false;
393
394 constexpr int numNodes = 4;
395 constexpr int numEdges = 6;
396 constexpr int refinementLevels = 6;
397
398 auto &m_field = fe_raw_ptr->
mField;
399 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
400 auto fe_handle = fe_ptr->getFEEntityHandle();
401
402 auto set_base_quadrature = [&]() {
404 int rule =
funRule(order_data);
408 "wrong dimension");
409 }
413 }
415 auto &gauss_pts = fe_ptr->gaussPts;
416 gauss_pts.resize(4, nb_gauss_pts, false);
417 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
418 &gauss_pts(0, 0), 1);
419 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
420 &gauss_pts(1, 0), 1);
421 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
422 &gauss_pts(2, 0), 1);
424 &gauss_pts(3, 0), 1);
425 auto &data = fe_ptr->dataOnElement[
H1];
426 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
427 false);
428 double *shape_ptr =
429 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
430 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
431 1);
432 } else {
435 }
437 };
438
439 CHKERR set_base_quadrature();
440
442
443 auto get_singular_nodes = [&]() {
444 int num_nodes;
446 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
447 true);
448 std::bitset<numNodes> singular_nodes;
449 for (auto nn = 0; nn != numNodes; ++nn) {
451 singular_nodes.set(nn);
452 } else {
453 singular_nodes.reset(nn);
454 }
455 }
456 return singular_nodes;
457 };
458
459 auto get_singular_edges = [&]() {
460 std::bitset<numEdges> singular_edges;
461 for (int ee = 0; ee != numEdges; ee++) {
463 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
465 singular_edges.set(ee);
466 } else {
467 singular_edges.reset(ee);
468 }
469 }
470 return singular_edges;
471 };
472
473 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
475 fe_ptr->gaussPts.swap(ref_gauss_pts);
476 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
477 auto &data = fe_ptr->dataOnElement[
H1];
478 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
479 double *shape_ptr =
480 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
482 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
483 nb_gauss_pts);
485 };
486
487 auto singular_nodes = get_singular_nodes();
488 if (singular_nodes.count()) {
489 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
491 CHKERR set_gauss_pts(it_map_ref_coords->second);
493 } else {
494
495 auto refine_quadrature = [&]() {
497
498 const int max_level = refinementLevels;
500
501 moab::Core moab_ref;
502 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
504 for (int nn = 0; nn != 4; nn++)
505 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
506 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
509 {
510 Range tets(tet, tet);
513 tets, 1, true, edges, moab::Interface::UNION);
516 }
517
518 Range nodes_at_front;
519 for (int nn = 0; nn != numNodes; nn++) {
520 if (singular_nodes[nn]) {
522 CHKERR moab_ref.side_element(tet, 0, nn, ent);
523 nodes_at_front.insert(ent);
524 }
525 }
526
527 auto singular_edges = get_singular_edges();
528
530 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
531 for (int ee = 0; ee != numEdges; ee++) {
532 if (singular_edges[ee]) {
534 CHKERR moab_ref.side_element(tet, 1, ee, ent);
535 CHKERR moab_ref.add_entities(meshset, &ent, 1);
536 }
537 }
538
539
541 for (int ll = 0; ll != max_level; ll++) {
544 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
546 edges);
548 CHKERR moab_ref.get_adjacencies(
549 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
550 ref_edges = intersect(ref_edges, edges);
552 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
553 ref_edges = intersect(ref_edges, ents);
556 ->getEntitiesByTypeAndRefLevel(
558 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
562 ->updateMeshsetByEntitiesChildren(meshset,
564 meshset, MBEDGE, true);
565 }
566
567
570 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
572 tets);
573
576 }
577
579 int tt = 0;
580 for (Range::iterator tit = tets.begin(); tit != tets.end();
581 tit++, tt++) {
582 int num_nodes;
584 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
585 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
586 }
587
588 auto &data = fe_ptr->dataOnElement[
H1];
589 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
590 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
592 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
593 int gg = 0;
594 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
595 double *tet_coords = &ref_coords(tt, 0);
597 det *= 6;
598 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
599 for (
int dd = 0;
dd != 3;
dd++) {
600 ref_gauss_pts(dd, gg) =
601 shape_n(ggg, 0) * tet_coords[3 * 0 +
dd] +
602 shape_n(ggg, 1) * tet_coords[3 * 1 +
dd] +
603 shape_n(ggg, 2) * tet_coords[3 * 2 +
dd] +
604 shape_n(ggg, 3) * tet_coords[3 * 3 +
dd];
605 }
606 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
607 }
608 }
609
610 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
612
613
616
617 TetPolynomialBase::switchCacheBaseOff<HDIV>({fe_raw_ptr});
618 TetPolynomialBase::switchCacheBaseOn<HDIV>({fe_raw_ptr});
619
621 };
622
623 CHKERR refine_quadrature();
624 }
625 }
626 }
627
629 }
#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.
