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 element or nodal tags on the source mesh to
  element or nodal tags on 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";
    char output_tag_name[255] = "";
 
    int interp_order = 1;
    PetscBool hybrid_interp = PETSC_TRUE;
    PetscBool src_tag_nodes_to_elem_avg = PETSC_FALSE;
    PetscBool src_tag_nodes_to_elem_set = PETSC_FALSE;
    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 PetscOptionsString("-output_tag", "output tag name",
                              "", "", output_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 PetscOptionsBool(
        "-src_nodes_to_elem", "0: use source nodes directly; 1: average nodes to elements; omitted: use source elements",
        "", src_tag_nodes_to_elem_avg, &src_tag_nodes_to_elem_avg,
        &src_tag_nodes_to_elem_set);
    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");
    }
 
    const bool src_tag_nodes =
        src_tag_nodes_to_elem_set && !src_tag_nodes_to_elem_avg;
    if (src_tag_nodes && interp_order != 1) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
              "-src_nodes_to_elem 0 requires -interp_order 1; use "
              "-src_nodes_to_elem 1 to average nodal data for elemental output");
    }
 
    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);
    }
 
    Tag output_tag = interp_tag;
    if (output_tag_name[0]) {
      std::vector<double> def_val(interp_tag_len, 0.0);
      auto rval = moab.tag_get_handle(output_tag_name, interp_tag_len,
                                      MB_TYPE_DOUBLE, output_tag,
                                      MB_TAG_CREAT | MB_TAG_DENSE,
                                      def_val.data());
      if (rval == MB_ALREADY_ALLOCATED) {
        SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
                "Output tag '%s' has a conflicting default; use another name",
                output_tag_name);
      }
      CHKERR rval;
    }
 
    if (src_tag_nodes_to_elem_avg) {
      for (auto &elem : src_elems) {
        Range adj_verts;
        CHKERR moab.get_connectivity(&elem, 1, adj_verts, true);
 
        std::vector<double> adj_vert_data(adj_verts.size() * interp_tag_len, 0.0);
        std::vector<double> elem_data(interp_tag_len, 0.0);
        CHKERR moab.tag_get_data(interp_tag, adj_verts, adj_vert_data.data());
 
        // Average vertex tag values and store the result on the source element.
        for (int itag = 0; itag < interp_tag_len; ++itag) {
          for (size_t ivert = 0; ivert != adj_verts.size(); ++ivert)
            elem_data[itag] += adj_vert_data[ivert * interp_tag_len + itag];
          elem_data[itag] /= adj_verts.size();
        }
        CHKERR moab.tag_set_data(interp_tag, &elem, 1, elem_data.data());
      }
    }
 
    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);
 
    Range src_entities;
    if (src_tag_nodes) {
      CHKERR moab.get_connectivity(src_elems, src_entities, true);
    } else {
      src_entities = src_elems;
    }
    std::vector<double> source_data(interp_tag_len * src_entities.size(), 0.0);
    std::vector<double> target_data(interp_tag_len * num_pts, 0.0);
 
    CHKERR moab.tag_get_data(interp_tag, src_entities, source_data.data());
 
    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_entities.size());
      for (size_t i = 0; i != src_entities.size(); ++i)
        source_data_scalar[i] = source_data[i * interp_tag_len + itag];
 
      CHKERR moab.tag_set_data(scalar_tag, src_entities,
                               source_data_scalar.data());
 
      if (!src_tag_nodes && 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];
      }
    }
 
    // Store the interpolated values on the target tag
    CHKERR moab.tag_set_data(output_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(output_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(output_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(output_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(output_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(output_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(output_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(output_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;
}