project_roughness.cpp

Go to the documentation of this file.

/** \file project_roughness.cpp
  \example project_roughness.cpp
  \brief project rough surface heights on mesh
*/
 
/* This file is part of MoFEM tool.
 * Mandatory to execute this tool:
 *   - A valid mesh file (-file_name)
 *   - A surface data file (-surface_data_file) with exactly nb_side_nodes × nb_side_nodes entries
 *
 * Example command usages:
 *   project_roughness -file_name mesh.cub -surface_data_file surface.data \
 *                     -mesh_side_nodes_num 55 -mesh_side_length 1.0 \
 *                     -mesh_scale 1 -mesh_height 2 -tolerance 1e-5
*/
// #include <BasicFiniteElements.hpp>
 
#include <MoFEM.hpp>
 
namespace bio = boost::iostreams;
using bio::stream;
using bio::tee_device;
 
using namespace MoFEM;
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    int rank;
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 
    // Read parameters from line command
    PetscBool flg_mesh_file = PETSC_TRUE;
    PetscBool flg_out_file = PETSC_TRUE;
    PetscBool flg_data_file = PETSC_TRUE;
 
    char mesh_file_name[255];
    char out_file_name[255] = "out.h5m";
    char data_file_name[255];
 
    int nb_side_nodes = 1;
    double side_length = 1.0;
    double mesh_height = 1.0;
    double mesh_scale = 1.0;
    double tol = 1e-5;
    int atom_test = 0;
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "Project Roughness", nullptr);
    CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
                              mesh_file_name, 255, &flg_mesh_file);
    CHKERR PetscOptionsString("-output_file", "out file name", "", "out.h5m",
                              out_file_name, 255, &flg_out_file);
    CHKERR PetscOptionsString("-surface_data_file", "data file name", "",
                              "data.txt", data_file_name, 255, &flg_data_file);
 
    CHKERR PetscOptionsInt("-mesh_side_nodes_num",
                           "number of nodes on each side", "", nb_side_nodes,
                           &nb_side_nodes, PETSC_NULLPTR);
    CHKERR PetscOptionsReal("-mesh_side_length", "side length", "", side_length,
                            &side_length, PETSC_NULLPTR);
    CHKERR PetscOptionsReal("-mesh_scale", "mesh scale", "", mesh_scale,
                            &mesh_scale, PETSC_NULLPTR);
    CHKERR PetscOptionsReal("-mesh_height", "vertical height of the mesh", "",
                            mesh_height, &mesh_height, PETSC_NULLPTR);
    CHKERR PetscOptionsReal("-tolerance", "Tolerance for coordinate matching",
                            "", tol, &tol, PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-atom_test", "atom test number", "", 0, &atom_test,
                           PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    MOFEM_LOG_CHANNEL("SELF");
 
    MOFEM_LOG_C("SELF", Sev::inform, "Mesh side length: %.2f", side_length);
    MOFEM_LOG_C("SELF", Sev::inform, "Mesh height: %.2f", mesh_height);
    MOFEM_LOG_C("SELF", Sev::inform, "Number of nodes on each side: %d",
                nb_side_nodes);
 
    std::ifstream infile(data_file_name, std::ios::in);
    std::vector<double> surface_data;
    std::vector<bool> surface_check;
 
    if (!infile.is_open()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
              "Could not open surface data file");
    }
 
    std::string line;
    while (std::getline(infile, line)) {
      std::istringstream iss(line);
      double num = 0.0;
      if (iss >> num)
        surface_data.push_back(num);
      else
        continue;
    }
    int expected_surface_nodes = nb_side_nodes * nb_side_nodes;
 
    if ((int)surface_data.size() != expected_surface_nodes) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Surface data size mismatch: expected %d nodes (%d x %d), got %d",
              expected_surface_nodes, nb_side_nodes, nb_side_nodes,
              (int)surface_data.size());
    }
    double max_surf_height =
        *std::max_element(std::begin(surface_data), std::end(surface_data));
    for (auto &node_data : surface_data) {
      node_data -= max_surf_height;
    }
 
    surface_check.resize(surface_data.size());
    std::fill(surface_check.begin(), surface_check.end(), false);
 
    if (flg_mesh_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
    }
 
    // Read mesh to MOAB
    const char *option;
    option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
    CHKERR moab.load_file(mesh_file_name, 0, option);
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
 
    MoFEM::Core core(moab);
    MOFEM_LOG_C("SELF", Sev::inform, "%s",
                "Projecting surface data on the mesh");
 
    Range elastic_elems, elastic_nodes;
    CHKERR moab.get_entities_by_dimension(
        0, 3, elastic_elems); // Todo: currently for 3D meshes only
    CHKERR moab.get_connectivity(elastic_elems, elastic_nodes);
 
    const double edge_length = side_length / (double)(nb_side_nodes - 1);
 
    MOFEM_LOG_C("SELF", Sev::inform, "Actual edge length per element: %.8f",
                edge_length);
 
    auto flat_ind = [&](int i, int j) { return j * nb_side_nodes + i; };
 
    for (Range::iterator nit = elastic_nodes.begin();
         nit != elastic_nodes.end(); nit++) {
      double coords[3];
      CHKERR moab.get_coords(&*nit, 1, coords);
      double x = coords[0];
      double y = coords[1];
      double z = coords[2];
      double height_coef = (mesh_height + z) / mesh_height;
 
      double ksi, eta;
      double pos_x = x / edge_length;
      double pos_y = y / edge_length;
 
      double pos_x_round = std::round(pos_x);
      double pos_y_round = std::round(pos_y);
 
      if (std::abs(pos_x_round - pos_x) < tol) {
        pos_x = pos_x_round;
      }
      if (std::abs(pos_y_round - pos_y) < tol) {
        pos_y = pos_y_round;
      }
 
      int i, j, i_next, j_next;
      i = i_next = (int)std::floor(pos_x);
      j = j_next = (int)std::floor(pos_y);
 
      if (std::abs(pos_x - std::floor(pos_x)) < tol) {
        ksi = -1.0;
      } else {
        ksi = (pos_x - std::floor(pos_x)) * 2.0 - 1.0;
        if (ksi < -1. - tol || ksi > 1. + tol) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Wrong ksi for i = %d, j = %d", i, j);
        }
        i_next++;
      }
      if (std::abs(pos_y - std::floor(pos_y)) < tol) {
        eta = -1.0;
      } else {
        eta = (pos_y - std::floor(pos_y)) * 2.0 - 1.0;
        if (eta < -1. - tol || eta > 1. + tol) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Wrong eta for i = %d, j = %d", i, j);
        }
        j_next++;
      }
      if (std::abs(z) < tol) {
        if (std::abs(pos_x - std::floor(pos_x)) > tol ||
            std::abs(pos_y - std::floor(pos_y)) > tol) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Cannot find data for surface node i = %d, j = %d", i, j);
        } else {
          if (surface_check[flat_ind(i, j)]) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                    "Data used twice for surface node i = %d, j = %d", i, j);
          } else {
            surface_check[flat_ind(i, j)] = true;
          }
        }
      }
      std::array<double, 4> shape_quad;
      shape_quad[0] = 0.25 * (1.0 - ksi) * (1.0 - eta);
      shape_quad[1] = 0.25 * (1.0 + ksi) * (1.0 - eta);
      shape_quad[2] = 0.25 * (1.0 + ksi) * (1.0 + eta);
      shape_quad[3] = 0.25 * (1.0 - ksi) * (1.0 + eta);
 
      std::array<double, 4> gap_quad;
      gap_quad[0] = surface_data[flat_ind(i, j)];
      gap_quad[1] = surface_data[flat_ind(i_next, j)];
      gap_quad[2] = surface_data[flat_ind(i_next, j_next)];
      gap_quad[3] = surface_data[flat_ind(i, j_next)];
 
      for (int k : {0, 1, 2, 3}) {
        coords[2] += shape_quad[k] * gap_quad[k] * height_coef;
      }
      CHKERR moab.set_coords(&*nit, 1, coords);
    }
 
    for (int i = 0; i < nb_side_nodes; i++) {
      for (int j = 0; j < nb_side_nodes; j++) {
        if (!surface_check[flat_ind(i, j)]) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Data not used for surface node i = %d, j = %d", i, j);
        }
      }
    }
    MOFEM_LOG_C("SELF", Sev::inform, "%s", "Scaling the mesh");
 
    for (Range::iterator nit = elastic_nodes.begin();
         nit != elastic_nodes.end(); nit++) {
      double coords[3];
      CHKERR moab.get_coords(&*nit, 1, coords);
      for (int i : {0, 1, 2}) {
        coords[i] *= mesh_scale;
      }
      CHKERR moab.set_coords(&*nit, 1, coords);
    }
 
    Range tets;
    CHKERR moab.get_entities_by_type(0, MBTET, tets);
 
    if (tets.empty()) {
    MOFEM_LOG_C("SELF", Sev::warning, "%s",
                 "No tetrahedral elements found, skipping quality and volume checks");
    } else {
      // variables for quality check and volume computation
      double total_volume = 0.0;
      double min_quality = std::numeric_limits<double>::max();
 
      for (auto tet : tets) {
        const EntityHandle *conn = nullptr;
        int num_nodes = 0;
        double coords[12];
 
        CHKERR moab.get_connectivity(tet, conn, num_nodes, true);
        CHKERR moab.get_coords(conn, num_nodes, coords);
 
        double q = Tools::volumeLengthQuality(coords);
        if (!std::isnormal(q) || q <= 0) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "invalid tet quality %.6e", q);
        }
        min_quality = std::min(min_quality, q);
 
        double vol = Tools::tetVolume(coords);
        total_volume += vol;
      }
      MOFEM_LOG_C("SELF", Sev::inform, "Min quality of the final mesh = %.6e",
                  min_quality);
      MOFEM_LOG_C("SELF", Sev::inform, "Total volume of the final mesh: %.4f",
                  total_volume);
 
      if (atom_test == 1) {
        double expected_volume = 1.98;
        double diff = std::abs(total_volume - expected_volume);
 
        if (diff > tol) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Atom test %d failed: Unexpected total tet volume computed: "
                  "%.6e (Expected %.6e)",
                  atom_test, total_volume, expected_volume);
        }
      }
    }
    MOFEM_LOG_C("SELF", Sev::inform, "Writing output file: %s", out_file_name);
    EntityHandle rootset = moab.get_root_set();
    CHKERR moab.write_file(out_file_name, "MOAB", "", &rootset, 1);
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
  return 0;
}