645 {
647
648 constexpr bool debug =
false;
649
650 constexpr int numNodes = 3;
651 constexpr int numEdges = 3;
652 constexpr int refinementLevels = 6;
653
654 auto &m_field = fe_raw_ptr->
mField;
655 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
656 auto fe_handle = fe_ptr->getFEEntityHandle();
657
658 auto set_base_quadrature = [&]() {
660 int rule =
funRule(order_data);
664 }
668 }
670 fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
671 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[1], 3,
672 &fe_ptr->gaussPts(0, 0), 1);
673 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[2], 3,
674 &fe_ptr->gaussPts(1, 0), 1);
676 &fe_ptr->gaussPts(2, 0), 1);
677 auto &data = fe_ptr->dataOnElement[
H1];
678 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3,
679 false);
680 double *shape_ptr =
681 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
682 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
683 1);
684 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(3, 2,
false);
685 std::copy(
687 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).data().begin());
688
689 } else {
692 }
694 };
695
696 CHKERR set_base_quadrature();
697
699
700 auto get_singular_nodes = [&]() {
701 int num_nodes;
703 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
704 true);
705 std::bitset<numNodes> singular_nodes;
706 for (auto nn = 0; nn != numNodes; ++nn) {
708 singular_nodes.set(nn);
709 } else {
710 singular_nodes.reset(nn);
711 }
712 }
713 return singular_nodes;
714 };
715
716 auto get_singular_edges = [&]() {
717 std::bitset<numEdges> singular_edges;
718 for (int ee = 0; ee != numEdges; ee++) {
720 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
722 singular_edges.set(ee);
723 } else {
724 singular_edges.reset(ee);
725 }
726 }
727 return singular_edges;
728 };
729
730 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
732 fe_ptr->gaussPts.swap(ref_gauss_pts);
733 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
734 auto &data = fe_ptr->dataOnElement[
H1];
735 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
736 double *shape_ptr =
737 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
739 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
741 };
742
743 auto singular_nodes = get_singular_nodes();
744 if (singular_nodes.count()) {
745 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
747 CHKERR set_gauss_pts(it_map_ref_coords->second);
749 } else {
750
751 auto refine_quadrature = [&]() {
753
754 const int max_level = refinementLevels;
755
756 moab::Core moab_ref;
757 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
759 for (int nn = 0; nn != numNodes; nn++)
760 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
762 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
765 {
766 Range tris(tri, tri);
769 tris, 1, true, edges, moab::Interface::UNION);
772 }
773
774 Range nodes_at_front;
775 for (int nn = 0; nn != numNodes; nn++) {
776 if (singular_nodes[nn]) {
778 CHKERR moab_ref.side_element(tri, 0, nn, ent);
779 nodes_at_front.insert(ent);
780 }
781 }
782
783 auto singular_edges = get_singular_edges();
784
786 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
787 for (int ee = 0; ee != numEdges; ee++) {
788 if (singular_edges[ee]) {
790 CHKERR moab_ref.side_element(tri, 1, ee, ent);
791 CHKERR moab_ref.add_entities(meshset, &ent, 1);
792 }
793 }
794
795
797 for (int ll = 0; ll != max_level; ll++) {
800 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
802 edges);
804 CHKERR moab_ref.get_adjacencies(
805 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
806 ref_edges = intersect(ref_edges, edges);
808 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
809 ref_edges = intersect(ref_edges, ents);
812 ->getEntitiesByTypeAndRefLevel(
814 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
818 ->updateMeshsetByEntitiesChildren(meshset,
820 meshset, MBEDGE, true);
821 }
822
823
826 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
828 tris);
829
833 }
834
836 int tt = 0;
837 for (Range::iterator tit = tris.begin(); tit != tris.end();
838 tit++, tt++) {
839 int num_nodes;
841 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
842 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
843 }
844
845 auto &data = fe_ptr->dataOnElement[
H1];
846 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
847 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
849 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
850 int gg = 0;
851 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
852 double *tri_coords = &ref_coords(tt, 0);
855 auto det = t_normal.
l2();
856 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
857 for (
int dd = 0;
dd != 2;
dd++) {
858 ref_gauss_pts(dd, gg) =
859 shape_n(ggg, 0) * tri_coords[3 * 0 +
dd] +
860 shape_n(ggg, 1) * tri_coords[3 * 1 +
dd] +
861 shape_n(ggg, 2) * tri_coords[3 * 2 +
dd];
862 }
863 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
864 }
865 }
866
867 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
869
871 };
872
873 CHKERR refine_quadrature();
874 }
875 }
876 }
877
879 }
#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.
constexpr int order
Order displacement.
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_2D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
static QUAD *const QUAD_2D_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.
