679 {
681
682 constexpr bool debug =
false;
683
684 constexpr int numNodes = 3;
685 constexpr int numEdges = 3;
686 constexpr int refinementLevels = 6;
687
688 auto &m_field = fe_raw_ptr->
mField;
689 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
690 auto fe_handle = fe_ptr->getFEEntityHandle();
691
692 auto set_base_quadrature = [&]() {
696 if (!xiao_rule) {
698 "Xiao--Gimbutas triangle rule is available for polynomial "
699 "orders 0 to %d; requested %d",
701 }
702 if (xiao_rule->numBarycentricCoordinates != 3) {
704 "wrong number of triangle barycentric coordinates");
705 }
706
707 const size_t nb_gauss_pts = xiao_rule->numPoints;
708 auto &gauss_pts = fe_ptr->gaussPts;
709 gauss_pts.resize(3, nb_gauss_pts, false);
710 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[1], 3, &gauss_pts(0, 0),
711 1);
712 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[2], 3, &gauss_pts(1, 0),
713 1);
714 cblas_dcopy(nb_gauss_pts, xiao_rule->weights, 1, &gauss_pts(2, 0), 1);
716 };
717
718 CHKERR set_base_quadrature();
719
721
722 auto get_singular_nodes = [&]() {
723 int num_nodes;
724 const EntityHandle *conn;
725 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
726 true);
727 std::bitset<numNodes> singular_nodes;
728 for (auto nn = 0; nn != numNodes; ++nn) {
730 singular_nodes.set(nn);
731 } else {
732 singular_nodes.reset(nn);
733 }
734 }
735 return singular_nodes;
736 };
737
738 auto get_singular_edges = [&]() {
739 std::bitset<numEdges> singular_edges;
740 for (int ee = 0; ee != numEdges; ee++) {
741 EntityHandle edge;
742 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
744 singular_edges.set(ee);
745 } else {
746 singular_edges.reset(ee);
747 }
748 }
749 return singular_edges;
750 };
751
752 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
754 fe_ptr->gaussPts.swap(ref_gauss_pts);
756 };
757
758 auto singular_nodes = get_singular_nodes();
759 if (singular_nodes.count()) {
760 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
762 CHKERR set_gauss_pts(it_map_ref_coords->second);
764 } else {
765
766 auto refine_quadrature = [&]() {
768
769 const int max_level = refinementLevels;
770
771 moab::Core moab_ref;
772 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
773 EntityHandle nodes[numNodes];
774 for (int nn = 0; nn != numNodes; nn++)
775 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
776 EntityHandle tri;
777 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
780 {
781 Range tris(tri, tri);
784 tris, 1, true, edges, moab::Interface::UNION);
787 }
788
789 Range nodes_at_front;
790 for (int nn = 0; nn != numNodes; nn++) {
791 if (singular_nodes[nn]) {
792 EntityHandle ent;
793 CHKERR moab_ref.side_element(tri, 0, nn, ent);
794 nodes_at_front.insert(ent);
795 }
796 }
797
798 auto singular_edges = get_singular_edges();
799
800 EntityHandle meshset;
801 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
802 for (int ee = 0; ee != numEdges; ee++) {
803 if (singular_edges[ee]) {
804 EntityHandle ent;
805 CHKERR moab_ref.side_element(tri, 1, ee, ent);
806 CHKERR moab_ref.add_entities(meshset, &ent, 1);
807 }
808 }
809
810
812 for (int ll = 0; ll != max_level; ll++) {
815 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
817 edges);
819 CHKERR moab_ref.get_adjacencies(
820 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
821 ref_edges = intersect(ref_edges, edges);
823 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
824 ref_edges = intersect(ref_edges, ents);
827 ->getEntitiesByTypeAndRefLevel(
829 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
833 ->updateMeshsetByEntitiesChildren(meshset,
835 meshset, MBEDGE, true);
836 }
837
838
841 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
843 tris);
844
848 }
849
851 int tt = 0;
852 for (Range::iterator tit = tris.begin(); tit != tris.end();
853 tit++, tt++) {
854 int num_nodes;
855 const EntityHandle *conn;
856 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
857 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
858 }
859
860 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
861 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
864 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
865 int gg = 0;
866 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
867 double *tri_coords = &ref_coords(tt, 0);
870 auto det = t_normal.
l2();
871 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
872 for (
int dd = 0;
dd != 2;
dd++) {
873 ref_gauss_pts(dd, gg) =
874 shape_n(ggg, 0) * tri_coords[3 * 0 +
dd] +
875 shape_n(ggg, 1) * tri_coords[3 * 1 +
dd] +
876 shape_n(ggg, 2) * tri_coords[3 * 2 +
dd];
877 }
878 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
879 }
880 }
881
882 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
884
886 };
887
888 CHKERR refine_quadrature();
889 }
890 }
891 }
892
894 }
#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 * getTriangleRule(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.
