mesh_data_transfer.cpp

Go to the documentation of this file.

/** \file mesh_data_transfer.cpp
 
  \brief Transfer data from source mesh to target mesh
 
  The data is transferred from tags on 3D elements of the source mesh to tags on
  3D elements in the target mesh. This code uses MBCoupler functionality of MOAB
  and is based on mbcoupler_test.cpp from MOAB.
 
  For this tool to work, source and target files should be partitioned using
  mbpart tool from MOAB, e.g.
  mbpart -z PHG -t 4 sslv116_hex.cub sslv116_hex_4mp.h5m
 
*/
 
#include <MoFEM.hpp>
#include <mbcoupler/Coupler.hpp>
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;
 
    // Create MoFEM database
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    char mesh_source_file[255] = "source.h5m";
    char mesh_target_file[255] = "target.h5m";
    char mesh_out_file[255] = "out.h5m";
    char iterp_tag_name[255] = "INTERNAL_STRESS";
 
    int interp_order = 1;
    PetscBool hybrid_interp = PETSC_TRUE;
 
    int atom_test = 0;
 
    double toler = 5.e-10;
 
    PetscOptionsBegin(m_field.get_comm(), "", "mesh data transfer tool",
                      "none");
    CHKERR PetscOptionsString("-source_file", "source mesh file name", "",
                              "source.h5m", mesh_source_file, 255,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsString("-target_file", "target mesh file name", "",
                              "target.h5m", mesh_target_file, 255,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsString("-output_file", "output mesh file name", "",
                              "out.h5m", mesh_out_file, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsString("-interp_tag", "interpolation tag name", "",
                              "INTERNAL_STRESS", iterp_tag_name, 255,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-interp_order", "interpolation order", "", 0,
                           &interp_order, PETSC_NULLPTR);
    CHKERR PetscOptionsBool("-hybrid_interp", "use hybrid interpolation", "",
                            hybrid_interp, &hybrid_interp, PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-atom_test", "atom test number", "", 0, &atom_test,
                           PETSC_NULLPTR);
 
    PetscOptionsEnd();
 
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_LOG_TAG("WORLD", "mesh_data_transfer");
 
    Coupler::Method method;
    switch (interp_order) {
    case 0:
      method = Coupler::CONSTANT;
      break;
    case 1:
      method = Coupler::LINEAR_FE;
      break;
    default:
      SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
              "Unsupported interpolation order");
    }
 
    std::vector<std::string> mesh_files(2);
    mesh_files[0] = string(mesh_source_file);
    mesh_files[1] = string(mesh_target_file);
 
    int nprocs, rank;
    ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
    CHKERRQ(ierr);
    ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    CHKERRQ(ierr);
 
    std::string read_opts, write_opts;
    read_opts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_"
                "DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS";
    if (nprocs > 1)
      read_opts += ";PARALLEL_GHOSTS=3.0.1";
    write_opts = (nprocs > 1) ? "PARALLEL=WRITE_PART" : "";
 
    std::vector<boost::shared_ptr<ParallelComm>> pcs(mesh_files.size());
    std::vector<EntityHandle> roots(mesh_files.size());
 
    for (unsigned int i = 0; i < mesh_files.size(); i++) {
      pcs[i] = boost::make_shared<ParallelComm>(&moab, MPI_COMM_WORLD);
      int index = pcs[i]->get_id();
      std::string newread_opts;
      std::ostringstream extra_opt;
      extra_opt << ";PARALLEL_COMM=" << index;
      newread_opts = read_opts + extra_opt.str();
 
      CHKERR moab.create_meshset(MESHSET_SET, roots[i]);
 
      CHKERR moab.load_file(mesh_files[i].c_str(), &roots[i],
                            newread_opts.c_str());
    }
 
    Range src_elems, targ_elems, targ_verts;
    CHKERR pcs[0]->get_part_entities(src_elems, 3);
 
    Tag interp_tag;
    CHKERR moab.tag_get_handle(iterp_tag_name, interp_tag);
 
    int interp_tag_len;
    CHKERR moab.tag_get_length(interp_tag, interp_tag_len);
 
    if (interp_tag_len != 1 && interp_tag_len != 3 && interp_tag_len != 9) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
              "Unsupported interpolation tag length: %d", interp_tag_len);
    }
 
    Coupler mbc(&moab, pcs[0].get(), src_elems, 0, true);
 
    std::vector<double> vpos; // the positions we are interested in
    int num_pts = 0;
 
    Range tmp_verts;
 
    // First get all vertices adj to partition entities in target mesh
    CHKERR pcs[1]->get_part_entities(targ_elems, 3);
    if (interp_order == 0) {
      targ_verts = targ_elems;
    } else {
      CHKERR moab.get_adjacencies(targ_elems, 0, false, targ_verts,
                                  moab::Interface::UNION);
    }
 
    // Then get non-owned verts and subtract
    CHKERR pcs[1]->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);
    targ_verts = subtract(targ_verts, tmp_verts);
 
    // get position of these entities; these are the target points
    num_pts = (int)targ_verts.size();
    vpos.resize(3 * targ_verts.size());
    CHKERR moab.get_coords(targ_verts, &vpos[0]);
 
    // Locate those points in the source mesh
    boost::shared_ptr<TupleList> tl_ptr;
    tl_ptr = boost::make_shared<TupleList>();
    CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(), false);
 
    // If some points were not located, we need to process them
    auto find_missing_points = [&](Range &targ_verts, int &num_pts,
                                   std::vector<double> &vpos,
                                   Range &missing_verts) {
      MoFEMFunctionBegin;
      int missing_pts_num = 0;
      int i = 0;
      auto vit = targ_verts.begin();
      for (; vit != targ_verts.end(); i++) {
        if (tl_ptr->vi_rd[3 * i + 1] == -1) {
          missing_verts.insert(*vit);
          vit = targ_verts.erase(vit);
          missing_pts_num++;
        } else {
          vit++;
        }
      }
 
      int missing_pts_num_global = 0;
      MPI_Allreduce(&missing_pts_num, &missing_pts_num_global, 1, MPI_INT,
                    MPI_SUM, m_field.get_comm());
      if (missing_pts_num_global) {
        MOFEM_LOG("WORLD", Sev::warning)
            << missing_pts_num_global
            << " points in target mesh were not located in source mesh. ";
      }
 
      if (missing_pts_num) {
        num_pts = (int)targ_verts.size();
        vpos.resize(3 * targ_verts.size());
        CHKERR moab.get_coords(targ_verts, &vpos[0]);
        tl_ptr->reset();
        CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(),
                                 false);
      }
      MoFEMFunctionReturn(0);
    };
 
    Range missing_verts;
    CHKERR find_missing_points(targ_verts, num_pts, vpos, missing_verts);
 
    std::vector<double> source_data(interp_tag_len * src_elems.size(), 0.0);
    std::vector<double> target_data(interp_tag_len * num_pts, 0.0);
 
    CHKERR moab.tag_get_data(interp_tag, src_elems, &source_data[0]);
 
    Tag scalar_tag, adj_count_tag;
    double def_scl = 0;
    string scalar_tag_name = string(iterp_tag_name) + "_COMP";
    CHKERR moab.tag_get_handle(scalar_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
                               scalar_tag, MB_TAG_CREAT | MB_TAG_DENSE,
                               &def_scl);
 
    string adj_count_tag_name = "ADJ_COUNT";
    double def_adj = 0;
    CHKERR moab.tag_get_handle(adj_count_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
                               adj_count_tag, MB_TAG_CREAT | MB_TAG_DENSE,
                               &def_adj);
 
    // MBCoupler functionality supports only scalar tags. For the case of
    // vector or tensor tags we need to save each component as a scalar tag
    auto create_scalar_tags = [&](const Range &src_elems,
                                  const std::vector<double> &source_data,
                                  int itag) {
      MoFEMFunctionBegin;
 
      std::vector<double> source_data_scalar(src_elems.size());
      // Populate source_data_scalar
      for (int ielem = 0; ielem < src_elems.size(); ielem++) {
        source_data_scalar[ielem] = source_data[itag + ielem * interp_tag_len];
      }
 
      // Set data on the scalar tag
      CHKERR moab.tag_set_data(scalar_tag, src_elems, &source_data_scalar[0]);
 
      if (interp_order == 1) {
        // Linear interpolation: compute average value of data on vertices
        Range src_verts;
        CHKERR moab.get_connectivity(src_elems, src_verts, true);
 
        CHKERR moab.tag_clear_data(scalar_tag, src_verts, &def_scl);
        CHKERR moab.tag_clear_data(adj_count_tag, src_verts, &def_adj);
 
        for (auto &tet : src_elems) {
          double tet_data = 0;
          CHKERR moab.tag_get_data(scalar_tag, &tet, 1, &tet_data);
 
          Range adj_verts;
          CHKERR moab.get_connectivity(&tet, 1, adj_verts, true);
 
          std::vector<double> adj_vert_data(adj_verts.size(), 0.0);
          std::vector<double> adj_vert_count(adj_verts.size(), 0.0);
 
          CHKERR moab.tag_get_data(scalar_tag, adj_verts, &adj_vert_data[0]);
          CHKERR moab.tag_get_data(adj_count_tag, adj_verts,
                                   &adj_vert_count[0]);
 
          for (int ivert = 0; ivert < adj_verts.size(); ivert++) {
            adj_vert_data[ivert] += tet_data;
            adj_vert_count[ivert] += 1;
          }
 
          CHKERR moab.tag_set_data(scalar_tag, adj_verts, &adj_vert_data[0]);
          CHKERR moab.tag_set_data(adj_count_tag, adj_verts,
                                   &adj_vert_count[0]);
        }
 
        // Reduce tags for the parallel case
        std::vector<Tag> tags = {scalar_tag, adj_count_tag};
        pcs[0]->reduce_tags(tags, tags, MPI_SUM, src_verts);
 
        std::vector<double> src_vert_data(src_verts.size(), 0.0);
        std::vector<double> src_vert_adj_count(src_verts.size(), 0.0);
 
        CHKERR moab.tag_get_data(scalar_tag, src_verts, &src_vert_data[0]);
        CHKERR moab.tag_get_data(adj_count_tag, src_verts,
                                 &src_vert_adj_count[0]);
 
        for (int ivert = 0; ivert < src_verts.size(); ivert++) {
          src_vert_data[ivert] /= src_vert_adj_count[ivert];
        }
        CHKERR moab.tag_set_data(scalar_tag, src_verts, &src_vert_data[0]);
      }
      MoFEMFunctionReturn(0);
    };
 
    for (int itag = 0; itag < interp_tag_len; itag++) {
 
      CHKERR create_scalar_tags(src_elems, source_data, itag);
 
      std::vector<double> target_data_scalar(num_pts, 0.0);
      CHKERR mbc.interpolate(method, scalar_tag_name, &target_data_scalar[0],
                             tl_ptr.get());
 
      for (int ielem = 0; ielem < num_pts; ielem++) {
        target_data[itag + ielem * interp_tag_len] = target_data_scalar[ielem];
      }
    }
 
    // Use original tag
    CHKERR moab.tag_set_data(interp_tag, targ_verts, &target_data[0]);
 
    if (missing_verts.size() && (interp_order == 1) && hybrid_interp) {
      MOFEM_LOG("WORLD", Sev::warning) << "Using hybrid interpolation for "
                                          "missing points in the target mesh.";
      Range missing_adj_elems;
      CHKERR moab.get_adjacencies(missing_verts, 3, false, missing_adj_elems,
                                  moab::Interface::UNION);
 
      int num_adj_elems = (int)missing_adj_elems.size();
      std::vector<double> vpos_adj_elems;
 
      vpos_adj_elems.resize(3 * missing_adj_elems.size());
      CHKERR moab.get_coords(missing_adj_elems, &vpos_adj_elems[0]);
 
      // Locate those points in the source mesh
      tl_ptr->reset();
      CHKERR mbc.locate_points(&vpos_adj_elems[0], num_adj_elems, 0, toler,
                               tl_ptr.get(), false);
 
      Range missing_tets;
      CHKERR find_missing_points(missing_adj_elems, num_adj_elems,
                                 vpos_adj_elems, missing_tets);
      if (missing_tets.size()) {
        MOFEM_LOG("WORLD", Sev::warning)
            << missing_tets.size()
            << "points in target mesh were not located in source mesh. ";
      }
 
      std::vector<double> target_data_adj_elems(interp_tag_len * num_adj_elems,
                                                0.0);
 
      for (int itag = 0; itag < interp_tag_len; itag++) {
        CHKERR create_scalar_tags(src_elems, source_data, itag);
 
        std::vector<double> target_data_adj_elems_scalar(num_adj_elems, 0.0);
        CHKERR mbc.interpolate(method, scalar_tag_name,
                               &target_data_adj_elems_scalar[0], tl_ptr.get());
 
        for (int ielem = 0; ielem < num_adj_elems; ielem++) {
          target_data_adj_elems[itag + ielem * interp_tag_len] =
              target_data_adj_elems_scalar[ielem];
        }
      }
 
      CHKERR moab.tag_set_data(interp_tag, missing_adj_elems,
                               &target_data_adj_elems[0]);
 
      // FIXME: add implementation for parallel case
      for (auto &vert : missing_verts) {
        Range adj_elems;
        CHKERR moab.get_adjacencies(&vert, 1, 3, false, adj_elems,
                                    moab::Interface::UNION);
 
        std::vector<double> adj_elems_data(adj_elems.size() * interp_tag_len,
                                           0.0);
        CHKERR moab.tag_get_data(interp_tag, adj_elems, &adj_elems_data[0]);
 
        std::vector<double> vert_data(interp_tag_len, 0.0);
        for (int itag = 0; itag < interp_tag_len; itag++) {
          for (int i = 0; i < adj_elems.size(); i++) {
            vert_data[itag] += adj_elems_data[i * interp_tag_len + itag];
          }
          vert_data[itag] /= adj_elems.size();
        }
        CHKERR moab.tag_set_data(interp_tag, &vert, 1, &vert_data[0]);
      }
    }
 
    CHKERR moab.tag_delete(scalar_tag);
    CHKERR moab.tag_delete(adj_count_tag);
 
    if (atom_test == 1 || atom_test == 2) {
      // Check the integrated stress for the SSLV116 test
 
      if (interp_order == 1) {
        // If the data was interpolated to target mesh vertices, we need to
        // average the vertex values for each 3D element and save on its tag
 
        // We need to reduce tags for the parallel case
        std::vector<Tag> tags;
        tags.push_back(interp_tag);
        pcs[1]->reduce_tags(tags, tags, MPI_SUM, targ_verts);
 
        for (auto &tet : targ_elems) {
          Range adj_verts;
          CHKERR moab.get_connectivity(&tet, 1, adj_verts, true);
 
          std::vector<double> adj_vert_data(adj_verts.size() * interp_tag_len,
                                            0.0);
          CHKERR moab.tag_get_data(interp_tag, adj_verts, &adj_vert_data[0]);
 
          std::vector<double> tet_data(interp_tag_len, 0.0);
          for (int itag = 0; itag < interp_tag_len; itag++) {
            for (int i = 0; i < adj_verts.size(); i++) {
              tet_data[itag] += adj_vert_data[i * interp_tag_len + itag];
            }
            tet_data[itag] /= adj_verts.size();
          }
          CHKERR moab.tag_set_data(interp_tag, &tet, 1, &tet_data[0]);
        }
      }
 
      std::vector<double> data_integ(interp_tag_len, 0.0);
      for (auto &tet : targ_elems) {
 
        std::vector<double> tet_data(interp_tag_len, 0.0);
        CHKERR moab.tag_get_data(interp_tag, &tet, 1, &tet_data[0]);
 
        const EntityHandle *vert_conn;
        int vert_num;
        CHKERR moab.get_connectivity(tet, vert_conn, vert_num, true);
        std::vector<double> vpos(3 * vert_num);
        CHKERR moab.get_coords(vert_conn, vert_num, &vpos[0]);
        double vol = Tools::tetVolume(vpos.data());
 
        for (int itag = 0; itag < interp_tag_len; itag++) {
          data_integ[itag] += tet_data[itag] * vol;
        }
      }
 
      std::vector<double> global_data_integ(interp_tag_len, 0.0);
      MPI_Allreduce(&data_integ[0], &global_data_integ[0], interp_tag_len,
                    MPI_DOUBLE, MPI_SUM, m_field.get_comm());
 
      std::string temp = "Integrated stress for sslv116 test: ";
      for (int i = 0; i < 9; ++i) {
        std::stringstream ss;
        ss << std::scientific << std::setprecision(3) << global_data_integ[i];
        temp += ss.str() + " ";
      }
      MOFEM_LOG("WORLD", Sev::inform) << temp;
 
      double non_zero_val = 1.655e12;
      double non_zero_tol = 5.1e-2;
      if (interp_order == 1) {
        non_zero_tol = 1e-3;
      }
      std::vector<int> non_zero_inds(3);
      std::vector<int> zero_inds(6);
      if (atom_test == 1) {
        non_zero_inds = {0, 4, 8};
        zero_inds = {1, 2, 3, 5, 6, 7};
      } else {
        non_zero_inds = {0, 1, 2};
        zero_inds = {3, 4, 5, 6, 7, 8};
      }
      bool non_zero_check =
          all_of(begin(non_zero_inds), end(non_zero_inds), [&](int i) {
            return abs(global_data_integ[i] - non_zero_val) / non_zero_val <
                   non_zero_tol;
          });
      bool zero_check = all_of(begin(zero_inds), end(zero_inds), [&](int i) {
        return abs(global_data_integ[i]) < 1e-12;
      });
      if (!non_zero_check || !zero_check) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                "Wrong value of the integrated stress");
      }
    } else if (atom_test) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
              "Unknown test number: %d", atom_test);
    }
 
    Range part_sets;
    // Only save the target mesh
    part_sets.insert(roots[1]);
    std::string new_write_opts;
    std::ostringstream extra_opt;
    if (nprocs > 1) {
      int index = pcs[1]->get_id();
      extra_opt << ";PARALLEL_COMM=" << index;
    }
    new_write_opts = write_opts + extra_opt.str();
    CHKERR moab.write_file(mesh_out_file, NULL, new_write_opts.c_str(),
                           part_sets);
    if (0 == rank) {
      MOFEM_LOG("WORLD", Sev::inform) << "Wrote file " << mesh_out_file;
    }
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}