656 {
658
659 constexpr bool debug =
false;
660
661 constexpr int numNodes = 3;
662 constexpr int numEdges = 3;
663 constexpr int refinementLevels = 6;
664
665 auto &m_field = fe_raw_ptr->
mField;
666 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
667 auto fe_handle = fe_ptr->getFEEntityHandle();
668
669 auto set_base_quadrature = [&]() {
672 const auto xiao_rule = IntRules::XiaoGimbutas::getTriangleRule(rule);
673 if (!xiao_rule) {
675 "Xiao--Gimbutas triangle rule is available for polynomial "
676 "orders 0 to %d; requested %d",
677 IntRules::XiaoGimbutas::triangleRuleCount, rule);
678 }
679 if (xiao_rule->numBarycentricCoordinates != 3) {
681 "wrong number of triangle barycentric coordinates");
682 }
683
684 const size_t nb_gauss_pts = xiao_rule->numPoints;
685 auto &gauss_pts = fe_ptr->gaussPts;
686 gauss_pts.resize(3, nb_gauss_pts, false);
687 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[1], 3, &gauss_pts(0, 0), 1);
688 cblas_dcopy(nb_gauss_pts, &xiao_rule->points[2], 3, &gauss_pts(1, 0), 1);
689 cblas_dcopy(nb_gauss_pts, xiao_rule->weights, 1, &gauss_pts(2, 0), 1);
691 };
692
693 CHKERR set_base_quadrature();
694
696
697 auto get_singular_nodes = [&]() {
698 int num_nodes;
699 const EntityHandle *conn;
700 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
701 true);
702 std::bitset<numNodes> singular_nodes;
703 for (auto nn = 0; nn != numNodes; ++nn) {
705 singular_nodes.set(nn);
706 } else {
707 singular_nodes.reset(nn);
708 }
709 }
710 return singular_nodes;
711 };
712
713 auto get_singular_edges = [&]() {
714 std::bitset<numEdges> singular_edges;
715 for (int ee = 0; ee != numEdges; ee++) {
716 EntityHandle edge;
717 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
719 singular_edges.set(ee);
720 } else {
721 singular_edges.reset(ee);
722 }
723 }
724 return singular_edges;
725 };
726
727 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
729 fe_ptr->gaussPts.swap(ref_gauss_pts);
731 };
732
733 auto singular_nodes = get_singular_nodes();
734 if (singular_nodes.count()) {
735 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
737 CHKERR set_gauss_pts(it_map_ref_coords->second);
739 } else {
740
741 auto refine_quadrature = [&]() {
743
744 const int max_level = refinementLevels;
745
746 moab::Core moab_ref;
747 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
748 EntityHandle nodes[numNodes];
749 for (int nn = 0; nn != numNodes; nn++)
750 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
751 EntityHandle tri;
752 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
755 {
756 Range tris(tri, tri);
759 tris, 1, true, edges, moab::Interface::UNION);
762 }
763
764 Range nodes_at_front;
765 for (int nn = 0; nn != numNodes; nn++) {
766 if (singular_nodes[nn]) {
767 EntityHandle ent;
768 CHKERR moab_ref.side_element(tri, 0, nn, ent);
769 nodes_at_front.insert(ent);
770 }
771 }
772
773 auto singular_edges = get_singular_edges();
774
775 EntityHandle meshset;
776 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
777 for (int ee = 0; ee != numEdges; ee++) {
778 if (singular_edges[ee]) {
779 EntityHandle ent;
780 CHKERR moab_ref.side_element(tri, 1, ee, ent);
781 CHKERR moab_ref.add_entities(meshset, &ent, 1);
782 }
783 }
784
785
787 for (int ll = 0; ll != max_level; ll++) {
790 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
792 edges);
794 CHKERR moab_ref.get_adjacencies(
795 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
796 ref_edges = intersect(ref_edges, edges);
798 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
799 ref_edges = intersect(ref_edges, ents);
802 ->getEntitiesByTypeAndRefLevel(
804 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
808 ->updateMeshsetByEntitiesChildren(meshset,
810 meshset, MBEDGE, true);
811 }
812
813
816 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
818 tris);
819
823 }
824
826 int tt = 0;
827 for (Range::iterator tit = tris.begin(); tit != tris.end();
828 tit++, tt++) {
829 int num_nodes;
830 const EntityHandle *conn;
831 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
832 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
833 }
834
835 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
836 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
839 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
840 int gg = 0;
841 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
842 double *tri_coords = &ref_coords(tt, 0);
845 auto det = t_normal.
l2();
846 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
847 for (
int dd = 0;
dd != 2;
dd++) {
848 ref_gauss_pts(dd, gg) =
849 shape_n(ggg, 0) * tri_coords[3 * 0 +
dd] +
850 shape_n(ggg, 1) * tri_coords[3 * 1 +
dd] +
851 shape_n(ggg, 2) * tri_coords[3 * 2 +
dd];
852 }
853 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
854 }
855 }
856
857 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
859
861 };
862
863 CHKERR refine_quadrature();
864 }
865 }
866 }
867
869 }
#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.
