149                                                                            {
  151 
  152#ifndef NDEBUG
  155        << "Refinement for hexes is not implemented";
  156#endif
  157 
  158  moab::Core core_ref;
  159  moab::Interface &moab_ref = core_ref;
  160 
  161  auto create_reference_element = [&moab_ref]() {
  163    constexpr double base_coords[] = {0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0,
  164                                      0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1};
  166    for (int nn = 0; nn < 8; nn++)
  167      CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
 
  169    CHKERR moab_ref.create_element(MBHEX, nodes, 8, hex);
 
  171  };
  172 
  173  auto add_ho_nodes = [&]() {
  176    CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, 
true);
 
  178    CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
 
  179    CHKERR moab_ref.add_entities(meshset, hexes);
 
  180    CHKERR moab_ref.convert_entities(meshset, 
true, 
true, 
true);
 
  181    CHKERR moab_ref.delete_entities(&meshset, 1);
 
  183  };
  184 
  185  auto set_gauss_pts = [&](std::map<EntityHandle, int> &little_map) {
  188    CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, 
true);
 
  190    CHKERR moab_ref.get_connectivity(hexes, hexes_nodes, 
false);
 
  192    gauss_pts.resize(hexes_nodes.size(), 4, false);
  193    size_t gg = 0;
  194    for (auto node : hexes_nodes) {
  195      CHKERR moab_ref.get_coords(&node, 1, &gauss_pts(gg, 0));
 
  196      little_map[node] = gg;
  197      ++gg;
  198    }
  199    gauss_pts = trans(gauss_pts);
  201  };
  202 
  203  auto set_ref_hexes = [&](std::map<EntityHandle, int> &little_map) {
  206    CHKERR moab_ref.get_entities_by_type(0, MBHEX, hexes, 
true);
 
  207    size_t hh = 0;
  209    for (auto hex : hexes) {
  211      int num_nodes;
  212      CHKERR moab_ref.get_connectivity(hex, conn, num_nodes, 
false);
 
  213      if (ref_hexes.size2() != num_nodes) {
  214        ref_hexes.resize(hexes.size(), num_nodes);
  215      }
  216      for (int nn = 0; nn != num_nodes; ++nn) {
  217        ref_hexes(hh, nn) = little_map[conn[nn]];
  218      }
  219      ++hh;
  220    }
  222  };
  223 
  224  auto set_shape_functions = [&]() {
  228    const auto nb_gauss_pts = gauss_pts.size2();
  229    shape_functions.resize(nb_gauss_pts, 8);
  230    for (int gg = 0; gg != nb_gauss_pts; ++gg) {
  231      const double ksi = gauss_pts(0, gg);
  232      const double zeta = gauss_pts(1, gg);
 
  233      const double eta = gauss_pts(2, gg);
 
  242    }
  244  };
  245 
  249 
  250  CHKERR create_reference_element();
 
  253  std::map<EntityHandle, int> little_map;
  254  CHKERR set_gauss_pts(little_map);
 
  255  CHKERR set_ref_hexes(little_map);
 
  256  CHKERR set_shape_functions();
 
  257 
  259}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
double zeta
Viscous hardening.
std::vector< ublas::matrix< int > > levelRef
std::vector< MatrixDouble > levelShapeFunctions
std::vector< MatrixDouble > levelGaussPtsOnRefMesh