Inheritance diagram for MyPostProc:

Collaboration diagram for MyPostProc:

Public Member Functions
MoFEMErrorCode	generateReferenceElementMesh ()

MoFEMErrorCode	setGaussPts (int order)

MoFEMErrorCode	preProcess ()

MoFEMErrorCode	postProcess ()

Protected Attributes
ublas::matrix< int >	refEleMap

MatrixDouble	shapeFunctions

Detailed Description

Examples: plot_base.cpp.

Definition at line 34 of file plot_base.cpp.

Member Function Documentation

◆ generateReferenceElementMesh()

MoFEMErrorCode MyPostProc::generateReferenceElementMesh ( )

Examples: plot_base.cpp.

Definition at line 397 of file plot_base.cpp.

                                                        {
  MoFEMFunctionBegin;
  moab::Core core_ref;
  moab::Interface &moab_ref = core_ref;
 
  char ref_mesh_file_name[255];
 
  if (SPACE_DIM == 2)
    strcpy(ref_mesh_file_name, "ref_mesh2d.h5m");
  else if (SPACE_DIM == 3)
    strcpy(ref_mesh_file_name, "ref_mesh3d.h5m");
  else
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Dimension not implemented");
 
  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-ref_file", ref_mesh_file_name,
                               255, PETSC_NULL);
  CHKERR moab_ref.load_file(ref_mesh_file_name, 0, "");
 
  // Get elements
  Range elems;
  CHKERR moab_ref.get_entities_by_dimension(0, SPACE_DIM, elems);
 
  // Add mid-nodes on edges
  EntityHandle meshset;
  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
  CHKERR moab_ref.add_entities(meshset, elems);
  CHKERR moab_ref.convert_entities(meshset, true, false, false);
  CHKERR moab_ref.delete_entities(&meshset, 1);
 
  // Get nodes on the mesh
  Range elem_nodes;
  CHKERR moab_ref.get_connectivity(elems, elem_nodes, false);
 
  // Map node entity and Gauss pint number
  std::map<EntityHandle, int> nodes_pts_map;
 
  // Set gauss points coordinates from the reference mesh
  gaussPts.resize(SPACE_DIM + 1, elem_nodes.size(), false);
  gaussPts.clear();
  Range::iterator nit = elem_nodes.begin();
  for (int gg = 0; nit != elem_nodes.end(); nit++, gg++) {
    double coords[3];
    CHKERR moab_ref.get_coords(&*nit, 1, coords);
    for (auto d : {0, 1, 2})
      gaussPts(d, gg) = coords[d];
    nodes_pts_map[*nit] = gg;
  }
 
  if (SPACE_DIM == 2) {
    // Set size of adjacency matrix (note ho order nodes 3 nodes and 3 nodes on
    // edges)
    refEleMap.resize(elems.size(), 3 + 3);
  } else if (SPACE_DIM == 3) {
    refEleMap.resize(elems.size(), 4 + 6);
  }
 
  // Set adjacency matrix
  Range::iterator tit = elems.begin();
  for (int tt = 0; tit != elems.end(); ++tit, ++tt) {
    const EntityHandle *conn;
    int num_nodes;
    CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
    for (int nn = 0; nn != num_nodes; ++nn) {
      refEleMap(tt, nn) = nodes_pts_map[conn[nn]];
    }
  }
 
  MoFEMFunctionReturn(0);
}

◆ postProcess()

MoFEMErrorCode MyPostProc::postProcess ( )

Examples: plot_base.cpp.

Definition at line 622 of file plot_base.cpp.

                                       {
  MoFEMFunctionBeginHot;
 
  auto resolve_shared_ents = [&]() {
    MoFEMFunctionBegin;
 
    ParallelComm *pcomm_post_proc_mesh =
        ParallelComm::get_pcomm(&(postProcMesh), MYPCOMM_INDEX);
    if (pcomm_post_proc_mesh == NULL) {
      // wrapRefMeshComm =
      // boost::make_shared<WrapMPIComm>(T::mField.get_comm(), false);
      pcomm_post_proc_mesh = new ParallelComm(
          &(postProcMesh),
          PETSC_COMM_WORLD /*(T::wrapRefMeshComm)->get_comm()*/);
    }
 
    CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR resolve_shared_ents();
 
  MoFEMFunctionReturnHot(0);
}

◆ preProcess()

MoFEMErrorCode MyPostProc::preProcess ( )

Examples: plot_base.cpp.

Definition at line 612 of file plot_base.cpp.

                                      {
  MoFEMFunctionBegin;
  moab::Interface &moab = coreMesh;
  ParallelComm *pcomm_post_proc_mesh =
      ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
  if (pcomm_post_proc_mesh != NULL)
    delete pcomm_post_proc_mesh;
  MoFEMFunctionReturn(0);
};

◆ setGaussPts()

MoFEMErrorCode MyPostProc::setGaussPts ( int order )

Examples: plot_base.cpp.

Definition at line 468 of file plot_base.cpp.

                                                {
  MoFEMFunctionBegin;
 
  const int num_nodes = gaussPts.size2();
 
  // Calculate shape functions
 
  switch (numeredEntFiniteElementPtr->getEntType()) {
  case MBTRI:
    shapeFunctions.resize(num_nodes, 3);
    CHKERR Tools::shapeFunMBTRI(&*shapeFunctions.data().begin(),
                                &gaussPts(0, 0), &gaussPts(1, 0), num_nodes);
    break;
  case MBQUAD: {
    shapeFunctions.resize(num_nodes, 4);
    for (int gg = 0; gg != num_nodes; gg++) {
      double ksi = gaussPts(0, gg);
      double eta = gaussPts(1, gg);
      shapeFunctions(gg, 0) = N_MBQUAD0(ksi, eta);
      shapeFunctions(gg, 1) = N_MBQUAD1(ksi, eta);
      shapeFunctions(gg, 2) = N_MBQUAD2(ksi, eta);
      shapeFunctions(gg, 3) = N_MBQUAD3(ksi, eta);
    }
  } break;
  case MBTET: {
    shapeFunctions.resize(num_nodes, 8);
    CHKERR Tools::shapeFunMBTET(&*shapeFunctions.data().begin(),
                                &gaussPts(0, 0), &gaussPts(1, 0),
                                &gaussPts(2, 0), num_nodes);
  } break;
  case MBHEX: {
    shapeFunctions.resize(num_nodes, 8);
    for (int gg = 0; gg != num_nodes; gg++) {
      double ksi = gaussPts(0, gg);
      double eta = gaussPts(1, gg);
      double zeta = gaussPts(2, gg);
      shapeFunctions(gg, 0) = N_MBHEX0(ksi, eta, zeta);
      shapeFunctions(gg, 1) = N_MBHEX1(ksi, eta, zeta);
      shapeFunctions(gg, 2) = N_MBHEX2(ksi, eta, zeta);
      shapeFunctions(gg, 3) = N_MBHEX3(ksi, eta, zeta);
      shapeFunctions(gg, 4) = N_MBHEX4(ksi, eta, zeta);
      shapeFunctions(gg, 5) = N_MBHEX5(ksi, eta, zeta);
      shapeFunctions(gg, 6) = N_MBHEX6(ksi, eta, zeta);
      shapeFunctions(gg, 7) = N_MBHEX7(ksi, eta, zeta);
    }
  } break;
  default:
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Not implemented element type");
  }
 
  // Create physical nodes
 
  // MoAB interface allowing for creating nodes and elements in the bulk
  ReadUtilIface *iface;
  CHKERR postProcMesh.query_interface(iface);
 
  std::vector<double *> arrays; /// pointers to memory allocated by MoAB for
                                /// storing X, Y, and Z coordinates
  EntityHandle startv;          // Starting handle for vertex
  // Allocate memory for num_nodes, and return starting handle, and access to
  // memort.
  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
 
  mapGaussPts.resize(gaussPts.size2());
  for (int gg = 0; gg != num_nodes; ++gg)
    mapGaussPts[gg] = startv + gg;
 
  Tag th;
  int def_in_the_loop = -1;
  CHKERR postProcMesh.tag_get_handle("NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER, th,
                                     MB_TAG_CREAT | MB_TAG_SPARSE,
                                     &def_in_the_loop);
 
  // Create physical elements
 
  const int num_el = refEleMap.size1();
  const int num_nodes_on_ele = refEleMap.size2();
 
  EntityHandle starte; // Starting handle to first created element
  EntityHandle *conn;  // Access to MOAB memory with connectivity of elements
 
  // Create tris/tets in the bulk in MoAB database
  if (SPACE_DIM == 2)
    CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTRI, 0,
                                      starte, conn);
  else if (SPACE_DIM == 3)
    CHKERR iface->get_element_connect(num_el, num_nodes_on_ele, MBTET, 0,
                                      starte, conn);
  else
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Dimension not implemented");
 
  // At this point elements (memory for elements) is allocated, at code bellow
  // actual connectivity of elements is set.
  for (unsigned int tt = 0; tt != refEleMap.size1(); ++tt) {
    for (int nn = 0; nn != num_nodes_on_ele; ++nn)
      conn[num_nodes_on_ele * tt + nn] = mapGaussPts[refEleMap(tt, nn)];
  }
 
  // Finalise elements creation. At that point MOAB updates adjacency tables,
  // and elements are ready to use.
  CHKERR iface->update_adjacencies(starte, num_el, num_nodes_on_ele, conn);
 
  auto physical_elements = Range(starte, starte + num_el - 1);
  CHKERR postProcMesh.tag_clear_data(th, physical_elements, &(nInTheLoop));
 
  EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
  int fe_num_nodes;
  {
    const EntityHandle *conn;
    mField.get_moab().get_connectivity(fe_ent, conn, fe_num_nodes, true);
    coords.resize(3 * fe_num_nodes, false);
    CHKERR mField.get_moab().get_coords(conn, fe_num_nodes, &coords[0]);
  }
 
  // Set physical coordinates to physical nodes
  FTensor::Index<'i', 3> i;
  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_n(
      &*shapeFunctions.data().begin());
 
  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_coords(
      arrays[0], arrays[1], arrays[2]);
  const double *t_coords_ele_x = &coords[0];
  const double *t_coords_ele_y = &coords[1];
  const double *t_coords_ele_z = &coords[2];
  for (int gg = 0; gg != num_nodes; ++gg) {
    FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_ele_coords(
        t_coords_ele_x, t_coords_ele_y, t_coords_ele_z);
    t_coords(i) = 0;
    for (int nn = 0; nn != fe_num_nodes; ++nn) {
      t_coords(i) += t_n * t_ele_coords(i);
      for (auto ii : {0, 1, 2})
        if (std::abs(t_coords(ii)) < std::numeric_limits<float>::epsilon())
          t_coords(ii) = 0;
      ++t_ele_coords;
      ++t_n;
    }
    ++t_coords;
  }
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

◆ refEleMap

ublas::matrix<int> MyPostProc::refEleMap

protected

Examples: plot_base.cpp.

Definition at line 44 of file plot_base.cpp.

◆ shapeFunctions

MatrixDouble MyPostProc::shapeFunctions