662 {
664
665 constexpr bool debug =
false;
666
667 constexpr int numNodes = 3;
668 constexpr int numEdges = 3;
669 constexpr int refinementLevels = 6;
670
671 auto &m_field = fe_raw_ptr->
mField;
672 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
673 auto fe_handle = fe_ptr->getFEEntityHandle();
674
675 auto set_base_quadrature = [&]() {
677 int rule =
funRule(order_data);
681 }
685 }
687 fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
688 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[1], 3,
689 &fe_ptr->gaussPts(0, 0), 1);
690 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[rule]->points[2], 3,
691 &fe_ptr->gaussPts(1, 0), 1);
693 &fe_ptr->gaussPts(2, 0), 1);
694 auto &data = fe_ptr->dataOnElement[
H1];
695 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3,
696 false);
697 double *shape_ptr =
698 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
699 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
700 1);
701 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(3, 2,
false);
702 std::copy(
704 data->dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).data().begin());
705
706 } else {
709 }
711 };
712
713 CHKERR set_base_quadrature();
714
716
717 auto get_singular_nodes = [&]() {
718 int num_nodes;
720 CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
721 true);
722 std::bitset<numNodes> singular_nodes;
723 for (auto nn = 0; nn != numNodes; ++nn) {
725 singular_nodes.set(nn);
726 } else {
727 singular_nodes.reset(nn);
728 }
729 }
730 return singular_nodes;
731 };
732
733 auto get_singular_edges = [&]() {
734 std::bitset<numEdges> singular_edges;
735 for (int ee = 0; ee != numEdges; ee++) {
737 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
739 singular_edges.set(ee);
740 } else {
741 singular_edges.reset(ee);
742 }
743 }
744 return singular_edges;
745 };
746
747 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
749 fe_ptr->gaussPts.swap(ref_gauss_pts);
750 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
751 auto &data = fe_ptr->dataOnElement[
H1];
752 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
753 double *shape_ptr =
754 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
756 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
758 };
759
760 auto singular_nodes = get_singular_nodes();
761 if (singular_nodes.count()) {
762 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
764 CHKERR set_gauss_pts(it_map_ref_coords->second);
766 } else {
767
768 auto refine_quadrature = [&]() {
770
771 const int max_level = refinementLevels;
772
773 moab::Core moab_ref;
774 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
776 for (int nn = 0; nn != numNodes; nn++)
777 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
779 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
782 {
783 Range tris(tri, tri);
786 tris, 1, true, edges, moab::Interface::UNION);
789 }
790
791 Range nodes_at_front;
792 for (int nn = 0; nn != numNodes; nn++) {
793 if (singular_nodes[nn]) {
795 CHKERR moab_ref.side_element(tri, 0, nn, ent);
796 nodes_at_front.insert(ent);
797 }
798 }
799
800 auto singular_edges = get_singular_edges();
801
803 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
804 for (int ee = 0; ee != numEdges; ee++) {
805 if (singular_edges[ee]) {
807 CHKERR moab_ref.side_element(tri, 1, ee, ent);
808 CHKERR moab_ref.add_entities(meshset, &ent, 1);
809 }
810 }
811
812
814 for (int ll = 0; ll != max_level; ll++) {
817 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
819 edges);
821 CHKERR moab_ref.get_adjacencies(
822 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
823 ref_edges = intersect(ref_edges, edges);
825 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
826 ref_edges = intersect(ref_edges, ents);
829 ->getEntitiesByTypeAndRefLevel(
831 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
835 ->updateMeshsetByEntitiesChildren(meshset,
837 meshset, MBEDGE, true);
838 }
839
840
843 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
845 tris);
846
850 }
851
853 int tt = 0;
854 for (Range::iterator tit = tris.begin(); tit != tris.end();
855 tit++, tt++) {
856 int num_nodes;
858 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
859 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
860 }
861
862 auto &data = fe_ptr->dataOnElement[
H1];
863 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
864 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
866 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
867 int gg = 0;
868 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
869 double *tri_coords = &ref_coords(tt, 0);
872 auto det = t_normal.
l2();
873 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
874 for (
int dd = 0;
dd != 2;
dd++) {
875 ref_gauss_pts(dd, gg) =
876 shape_n(ggg, 0) * tri_coords[3 * 0 +
dd] +
877 shape_n(ggg, 1) * tri_coords[3 * 1 +
dd] +
878 shape_n(ggg, 2) * tri_coords[3 * 2 +
dd];
879 }
880 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
881 }
882 }
883
884 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
886
888 };
889
890 CHKERR refine_quadrature();
891 }
892 }
893 }
894
896 }
#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_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.
