EshelbianPlasticity::SetIntegrationAtFrontVolume Struct Reference

Collaboration diagram for EshelbianPlasticity::SetIntegrationAtFrontVolume:

Classes
struct	VolFe

Public Member Functions
	SetIntegrationAtFrontVolume ()=default
MoFEMErrorCode	operator() (ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)

Static Public Attributes
static constexpr int	numNodes = 4
static constexpr int	numEdges = 6
static constexpr int	refinementLevels = 3
static boost::shared_ptr< Range >	frontNodes
static boost::shared_ptr< Range >	frontEdges
static std::map< long int, MatrixDouble >	mapRefCoords

Detailed Description

Definition at line 103 of file EshelbianPlasticity.cpp.

Constructor & Destructor Documentation

SetIntegrationAtFrontVolume()

EshelbianPlasticity::SetIntegrationAtFrontVolume::SetIntegrationAtFrontVolume ( )

default

Member Function Documentation

operator()()

MoFEMErrorCode EshelbianPlasticity::SetIntegrationAtFrontVolume::operator() ( ForcesAndSourcesCore *  fe_raw_ptr, int  order_row, int  order_col, int  order_data  )

inline

Definition at line 107 of file EshelbianPlasticity.cpp.

                                                           {
    MoFEMFunctionBegin;
 
    auto &m_field = fe_raw_ptr->mField;
    auto vol_fe_ptr = dynamic_cast<VolFe *>(fe_raw_ptr);
    auto fe_handle = vol_fe_ptr->getFEEntityHandle();
 
    auto set_base_quadrature = [&]() {
      MoFEMFunctionBegin;
      int rule = 2 * order_data + 1;
      if (rule < QUAD_3D_TABLE_SIZE) {
        if (QUAD_3D_TABLE[rule]->dim != 3) {
          SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                  "wrong dimension");
        }
        if (QUAD_3D_TABLE[rule]->order < rule) {
          SETERRQ2(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
                   "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
        }
        const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
        auto &gauss_pts = vol_fe_ptr->gaussPts;
        gauss_pts.resize(4, nb_gauss_pts, false);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
                    &gauss_pts(0, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
                    &gauss_pts(1, 0), 1);
        cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
                    &gauss_pts(2, 0), 1);
        cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
                    &gauss_pts(3, 0), 1);
        auto &data = vol_fe_ptr->dataOnElement[H1];
        data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
                                                              false);
        double *shape_ptr =
            &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1,
                    shape_ptr, 1);
      } else {
        SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                 "rule > quadrature order %d < %d", rule,
                 QUAD_3D_TABLE_SIZE);
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_base_quadrature();
 
    if (frontNodes->size() && EshelbianCore::setSingularity) {
 
      auto get_singular_nodes = [&]() {
        int num_nodes;
        const EntityHandle *conn;
        CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
                                                   true);
        std::bitset<numNodes> singular_nodes;
        for (auto nn = 0; nn != numNodes; ++nn) {
          if (frontNodes->find(conn[nn]) != frontNodes->end()) {
            singular_nodes.set(nn);
          } else {
            singular_nodes.reset(nn);
          }
        }
        return singular_nodes;
      };
 
      auto get_singular_edges = [&]() {
        std::bitset<numEdges> singular_edges;
        for (int ee = 0; ee != numEdges; ee++) {
          EntityHandle edge;
          CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
          if (frontEdges->find(edge) != frontEdges->end()) {
            singular_edges.set(ee);
          } else {
            singular_edges.reset(ee);
          }
        }
        return singular_edges;
      };
 
      auto set_gauss_pts = [&](auto &ref_gauss_pts) {
        MoFEMFunctionBegin;
        vol_fe_ptr->gaussPts.swap(ref_gauss_pts);
        const size_t nb_gauss_pts = vol_fe_ptr->gaussPts.size2();
        // auto &dataH1 = vol_fe_ptr->dataOnElement[H1];
        // data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
        // double *shape_ptr =
        //     &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
        // CHKERR ShapeMBTET(shape_ptr, &vol_fe_ptr->gaussPts(0, 0),
        //                   &vol_fe_ptr->gaussPts(1, 0), &gaussPts(2, 0),
        //                   nb_gauss_pts);
        MoFEMFunctionReturn(0);
      };
 
      auto singular_nodes = get_singular_nodes();
      if (singular_nodes.count()) {
        auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
        if (it_map_ref_coords != mapRefCoords.end()) {
          CHKERR set_gauss_pts(it_map_ref_coords->second);
          MoFEMFunctionReturnHot(0);
        } else {
 
          auto refine_quadrature = [&]() {
            MoFEMFunctionBegin;
 
            const int max_level = refinementLevels;
            EntityHandle tet;
 
            moab::Core moab_ref;
            double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
            EntityHandle nodes[4];
            for (int nn = 0; nn != 4; nn++)
              CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
            CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
            MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
            MoFEM::Interface &m_field_ref = mofem_ref_core;
            CHKERR m_field_ref.getInterface<BitRefManager>()
                ->setBitRefLevelByDim(0, 3, BitRefLevel().set(0));
 
            Range nodes_at_front;
            for (int nn = 0; nn != numNodes; nn++) {
              if (singular_nodes[nn]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tet, 0, nn, ent);
                nodes_at_front.insert(ent);
              }
            }
 
            auto singular_edges = get_singular_edges();
 
            EntityHandle meshset;
            CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
            for (int ee = 0; ee != 6; ee++) {
              if (singular_edges[ee]) {
                EntityHandle ent;
                CHKERR moab_ref.side_element(tet, 1, ee, ent);
                CHKERR moab_ref.add_entities(meshset, &ent, 1);
              }
            }
 
            // refine mesh
            auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
            for (int ll = 0; ll != max_level; ll++) {
              Range edges;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
                                                 BitRefLevel().set(), MBEDGE,
                                                 edges);
              Range ref_edges;
              CHKERR moab_ref.get_adjacencies(
                  nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
              ref_edges = intersect(ref_edges, edges);
              Range ents;
              CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
              ref_edges = intersect(ref_edges, ents);
              Range tets;
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->getEntitiesByTypeAndRefLevel(
                      BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
              CHKERR m_ref->addVerticesInTheMiddleOfEdges(
                  ref_edges, BitRefLevel().set(ll + 1));
              CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
              CHKERR m_field_ref.getInterface<BitRefManager>()
                  ->updateMeshsetByEntitiesChildren(meshset,
                                                    BitRefLevel().set(ll + 1),
                                                    meshset, MBEDGE, true);
            }
 
            // get ref coords
            Range tets;
            CHKERR m_field_ref.getInterface<BitRefManager>()
                ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
                                               BitRefLevel().set(), MBTET,
                                               tets);
            MatrixDouble ref_coords(tets.size(), 12, false);
            int tt = 0;
            for (Range::iterator tit = tets.begin(); tit != tets.end();
                 tit++, tt++) {
              int num_nodes;
              const EntityHandle *conn;
              CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
              CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
            }
 
            auto &data = vol_fe_ptr->dataOnElement[H1];
            const size_t nb_gauss_pts = vol_fe_ptr->gaussPts.size2();
            MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
            MatrixDouble &shape_n =
                data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
            int gg = 0;
            for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
              double *tet_coords = &ref_coords(tt, 0);
              double det = Tools::tetVolume(tet_coords);
              det *= 6;
              for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
                for (int dd = 0; dd != 3; dd++) {
                  ref_gauss_pts(dd, gg) =
                      shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
                      shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
                      shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
                      shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
                }
                ref_gauss_pts(3, gg) = vol_fe_ptr->gaussPts(3, ggg) * det;
              }
            }
 
            mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
            CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
 
            MoFEMFunctionReturn(0);
          };
        }
      }
    }
 
    MoFEMFunctionReturn(0);
  }

Member Data Documentation

frontEdges

boost::shared_ptr<Range> EshelbianPlasticity::SetIntegrationAtFrontVolume::frontEdges

inlinestatic

Definition at line 335 of file EshelbianPlasticity.cpp.

frontNodes

boost::shared_ptr<Range> EshelbianPlasticity::SetIntegrationAtFrontVolume::frontNodes

inlinestatic

Definition at line 334 of file EshelbianPlasticity.cpp.

mapRefCoords

std::map<long int, MatrixDouble> EshelbianPlasticity::SetIntegrationAtFrontVolume::mapRefCoords

inlinestatic

Definition at line 336 of file EshelbianPlasticity.cpp.

numEdges

constexpr int EshelbianPlasticity::SetIntegrationAtFrontVolume::numEdges = 6

inlinestaticconstexpr

Definition at line 331 of file EshelbianPlasticity.cpp.

numNodes

constexpr int EshelbianPlasticity::SetIntegrationAtFrontVolume::numNodes = 4

inlinestaticconstexpr

Definition at line 330 of file EshelbianPlasticity.cpp.

refinementLevels

constexpr int EshelbianPlasticity::SetIntegrationAtFrontVolume::refinementLevels = 3

inlinestaticconstexpr

Definition at line 332 of file EshelbianPlasticity.cpp.

---

The documentation for this struct was generated from the following file:

• users_modules/eshelbian_plasticit/src/impl/EshelbianPlasticity.cpp