21int main(
int argc, 
char *argv[]) {
 
   26    moab::Core mb_instance;
 
   27    moab::Interface &moab = mb_instance;
 
   33    char mesh_source_file[255] = 
"source.h5m";
 
   34    char mesh_target_file[255] = 
"target.h5m";
 
   35    char mesh_out_file[255] = 
"out.h5m";
 
   36    char iterp_tag_name[255] = 
"INTERNAL_STRESS";
 
   39    PetscBool hybrid_interp = PETSC_TRUE;
 
   43    double toler = 5.e-10;
 
   45    PetscOptionsBegin(m_field.
get_comm(), 
"", 
"mesh data transfer tool",
 
   47    CHKERR PetscOptionsString(
"-source_file", 
"source mesh file name", 
"",
 
   48                              "source.h5m", mesh_source_file, 255,
 
   50    CHKERR PetscOptionsString(
"-target_file", 
"target mesh file name", 
"",
 
   51                              "target.h5m", mesh_target_file, 255,
 
   53    CHKERR PetscOptionsString(
"-output_file", 
"output mesh file name", 
"",
 
   54                              "out.h5m", mesh_out_file, 255, PETSC_NULLPTR);
 
   55    CHKERR PetscOptionsString(
"-interp_tag", 
"interpolation tag name", 
"",
 
   56                              "INTERNAL_STRESS", iterp_tag_name, 255,
 
   58    CHKERR PetscOptionsInt(
"-interp_order", 
"interpolation order", 
"", 0,
 
   59                           &interp_order, PETSC_NULLPTR);
 
   60    CHKERR PetscOptionsBool(
"-hybrid_interp", 
"use hybrid interpolation", 
"",
 
   61                            hybrid_interp, &hybrid_interp, PETSC_NULLPTR);
 
   62    CHKERR PetscOptionsInt(
"-atom_test", 
"atom test number", 
"", 0, &
atom_test,
 
   70    Coupler::Method method;
 
   71    switch (interp_order) {
 
   73      method = Coupler::CONSTANT;
 
   76      method = Coupler::LINEAR_FE;
 
   80              "Unsupported interpolation order");
 
   83    std::vector<std::string> mesh_files(2);
 
   84    mesh_files[0] = string(mesh_source_file);
 
   85    mesh_files[1] = string(mesh_target_file);
 
   88    ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
 
   90    ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
   93    std::string read_opts, write_opts;
 
   94    read_opts = 
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_" 
   95                "DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS";
 
   97      read_opts += 
";PARALLEL_GHOSTS=3.0.1";
 
   98    write_opts = (nprocs > 1) ? 
"PARALLEL=WRITE_PART" : 
"";
 
  100    std::vector<boost::shared_ptr<ParallelComm>> pcs(mesh_files.size());
 
  101    std::vector<EntityHandle> roots(mesh_files.size());
 
  103    for (
unsigned int i = 0; 
i < mesh_files.size(); 
i++) {
 
  104      pcs[
i] = boost::make_shared<ParallelComm>(&moab, MPI_COMM_WORLD);
 
  105      int index = pcs[
i]->get_id();
 
  106      std::string newread_opts;
 
  107      std::ostringstream extra_opt;
 
  108      extra_opt << 
";PARALLEL_COMM=" << index;
 
  109      newread_opts = read_opts + extra_opt.str();
 
  111      CHKERR moab.create_meshset(MESHSET_SET, roots[
i]);
 
  113      CHKERR moab.load_file(mesh_files[
i].c_str(), &roots[
i],
 
  114                            newread_opts.c_str());
 
  117    Range src_elems, targ_elems, targ_verts;
 
  118    CHKERR pcs[0]->get_part_entities(src_elems, 3);
 
  121    CHKERR moab.tag_get_handle(iterp_tag_name, interp_tag);
 
  124    CHKERR moab.tag_get_length(interp_tag, interp_tag_len);
 
  126    if (interp_tag_len != 1 && interp_tag_len != 3 && interp_tag_len != 9) {
 
  128              "Unsupported interpolation tag length: %d", interp_tag_len);
 
  131    Coupler mbc(&moab, pcs[0].get(), src_elems, 0, 
true);
 
  133    std::vector<double> vpos; 
 
  139    CHKERR pcs[1]->get_part_entities(targ_elems, 3);
 
  140    if (interp_order == 0) {
 
  141      targ_verts = targ_elems;
 
  143      CHKERR moab.get_adjacencies(targ_elems, 0, 
false, targ_verts,
 
  144                                  moab::Interface::UNION);
 
  148    CHKERR pcs[1]->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);
 
  149    targ_verts = subtract(targ_verts, tmp_verts);
 
  152    num_pts = (int)targ_verts.size();
 
  153    vpos.resize(3 * targ_verts.size());
 
  154    CHKERR moab.get_coords(targ_verts, &vpos[0]);
 
  157    boost::shared_ptr<TupleList> tl_ptr;
 
  158    tl_ptr = boost::make_shared<TupleList>();
 
  159    CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(), 
false);
 
  162    auto find_missing_points = [&](
Range &targ_verts, 
int &num_pts,
 
  163                                   std::vector<double> &vpos,
 
  164                                   Range &missing_verts) {
 
  166      int missing_pts_num = 0;
 
  168      auto vit = targ_verts.begin();
 
  169      for (; vit != targ_verts.end(); 
i++) {
 
  170        if (tl_ptr->vi_rd[3 * 
i + 1] == -1) {
 
  171          missing_verts.insert(*vit);
 
  172          vit = targ_verts.erase(vit);
 
  179      int missing_pts_num_global = 0;
 
  180      MPI_Allreduce(&missing_pts_num, &missing_pts_num_global, 1, MPI_INT,
 
  182      if (missing_pts_num_global) {
 
  184            << missing_pts_num_global
 
  185            << 
" points in target mesh were not located in source mesh. ";
 
  188      if (missing_pts_num) {
 
  189        num_pts = (int)targ_verts.size();
 
  190        vpos.resize(3 * targ_verts.size());
 
  191        CHKERR moab.get_coords(targ_verts, &vpos[0]);
 
  193        CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(),
 
  200    CHKERR find_missing_points(targ_verts, num_pts, vpos, missing_verts);
 
  202    std::vector<double> source_data(interp_tag_len * src_elems.size(), 0.0);
 
  203    std::vector<double> target_data(interp_tag_len * num_pts, 0.0);
 
  205    CHKERR moab.tag_get_data(interp_tag, src_elems, &source_data[0]);
 
  207    Tag scalar_tag, adj_count_tag;
 
  209    string scalar_tag_name = string(iterp_tag_name) + 
"_COMP";
 
  210    CHKERR moab.tag_get_handle(scalar_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
 
  211                               scalar_tag, MB_TAG_CREAT | MB_TAG_DENSE,
 
  214    string adj_count_tag_name = 
"ADJ_COUNT";
 
  216    CHKERR moab.tag_get_handle(adj_count_tag_name.c_str(), 1, MB_TYPE_DOUBLE,
 
  217                               adj_count_tag, MB_TAG_CREAT | MB_TAG_DENSE,
 
  222    auto create_scalar_tags = [&](
const Range &src_elems,
 
  223                                  const std::vector<double> &source_data,
 
  227      std::vector<double> source_data_scalar(src_elems.size());
 
  229      for (
int ielem = 0; ielem < src_elems.size(); ielem++) {
 
  230        source_data_scalar[ielem] = source_data[itag + ielem * interp_tag_len];
 
  234      CHKERR moab.tag_set_data(scalar_tag, src_elems, &source_data_scalar[0]);
 
  236      if (interp_order == 1) {
 
  239        CHKERR moab.get_connectivity(src_elems, src_verts, 
true);
 
  241        CHKERR moab.tag_clear_data(scalar_tag, src_verts, &def_scl);
 
  242        CHKERR moab.tag_clear_data(adj_count_tag, src_verts, &def_adj);
 
  244        for (
auto &tet : src_elems) {
 
  246          CHKERR moab.tag_get_data(scalar_tag, &tet, 1, &tet_data);
 
  249          CHKERR moab.get_connectivity(&tet, 1, adj_verts, 
true);
 
  251          std::vector<double> adj_vert_data(adj_verts.size(), 0.0);
 
  252          std::vector<double> adj_vert_count(adj_verts.size(), 0.0);
 
  254          CHKERR moab.tag_get_data(scalar_tag, adj_verts, &adj_vert_data[0]);
 
  255          CHKERR moab.tag_get_data(adj_count_tag, adj_verts,
 
  258          for (
int ivert = 0; ivert < adj_verts.size(); ivert++) {
 
  259            adj_vert_data[ivert] += tet_data;
 
  260            adj_vert_count[ivert] += 1;
 
  263          CHKERR moab.tag_set_data(scalar_tag, adj_verts, &adj_vert_data[0]);
 
  264          CHKERR moab.tag_set_data(adj_count_tag, adj_verts,
 
  269        std::vector<Tag> tags = {scalar_tag, adj_count_tag};
 
  270        pcs[0]->reduce_tags(tags, tags, MPI_SUM, src_verts);
 
  272        std::vector<double> src_vert_data(src_verts.size(), 0.0);
 
  273        std::vector<double> src_vert_adj_count(src_verts.size(), 0.0);
 
  275        CHKERR moab.tag_get_data(scalar_tag, src_verts, &src_vert_data[0]);
 
  276        CHKERR moab.tag_get_data(adj_count_tag, src_verts,
 
  277                                 &src_vert_adj_count[0]);
 
  279        for (
int ivert = 0; ivert < src_verts.size(); ivert++) {
 
  280          src_vert_data[ivert] /= src_vert_adj_count[ivert];
 
  282        CHKERR moab.tag_set_data(scalar_tag, src_verts, &src_vert_data[0]);
 
  287    for (
int itag = 0; itag < interp_tag_len; itag++) {
 
  289      CHKERR create_scalar_tags(src_elems, source_data, itag);
 
  291      std::vector<double> target_data_scalar(num_pts, 0.0);
 
  292      CHKERR mbc.interpolate(method, scalar_tag_name, &target_data_scalar[0],
 
  295      for (
int ielem = 0; ielem < num_pts; ielem++) {
 
  296        target_data[itag + ielem * interp_tag_len] = target_data_scalar[ielem];
 
  301    CHKERR moab.tag_set_data(interp_tag, targ_verts, &target_data[0]);
 
  303    if (missing_verts.size() && (interp_order == 1) && hybrid_interp) {
 
  304      MOFEM_LOG(
"WORLD", Sev::warning) << 
"Using hybrid interpolation for " 
  305                                          "missing points in the target mesh.";
 
  306      Range missing_adj_elems;
 
  307      CHKERR moab.get_adjacencies(missing_verts, 3, 
false, missing_adj_elems,
 
  308                                  moab::Interface::UNION);
 
  310      int num_adj_elems = (int)missing_adj_elems.size();
 
  311      std::vector<double> vpos_adj_elems;
 
  313      vpos_adj_elems.resize(3 * missing_adj_elems.size());
 
  314      CHKERR moab.get_coords(missing_adj_elems, &vpos_adj_elems[0]);
 
  318      CHKERR mbc.locate_points(&vpos_adj_elems[0], num_adj_elems, 0, toler,
 
  319                               tl_ptr.get(), 
false);
 
  322      CHKERR find_missing_points(missing_adj_elems, num_adj_elems,
 
  323                                 vpos_adj_elems, missing_tets);
 
  324      if (missing_tets.size()) {
 
  326            << missing_tets.size()
 
  327            << 
"points in target mesh were not located in source mesh. ";
 
  330      std::vector<double> target_data_adj_elems(interp_tag_len * num_adj_elems,
 
  333      for (
int itag = 0; itag < interp_tag_len; itag++) {
 
  334        CHKERR create_scalar_tags(src_elems, source_data, itag);
 
  336        std::vector<double> target_data_adj_elems_scalar(num_adj_elems, 0.0);
 
  337        CHKERR mbc.interpolate(method, scalar_tag_name,
 
  338                               &target_data_adj_elems_scalar[0], tl_ptr.get());
 
  340        for (
int ielem = 0; ielem < num_adj_elems; ielem++) {
 
  341          target_data_adj_elems[itag + ielem * interp_tag_len] =
 
  342              target_data_adj_elems_scalar[ielem];
 
  346      CHKERR moab.tag_set_data(interp_tag, missing_adj_elems,
 
  347                               &target_data_adj_elems[0]);
 
  350      for (
auto &vert : missing_verts) {
 
  352        CHKERR moab.get_adjacencies(&vert, 1, 3, 
false, adj_elems,
 
  353                                    moab::Interface::UNION);
 
  355        std::vector<double> adj_elems_data(adj_elems.size() * interp_tag_len,
 
  357        CHKERR moab.tag_get_data(interp_tag, adj_elems, &adj_elems_data[0]);
 
  359        std::vector<double> vert_data(interp_tag_len, 0.0);
 
  360        for (
int itag = 0; itag < interp_tag_len; itag++) {
 
  361          for (
int i = 0; 
i < adj_elems.size(); 
i++) {
 
  362            vert_data[itag] += adj_elems_data[
i * interp_tag_len + itag];
 
  364          vert_data[itag] /= adj_elems.size();
 
  366        CHKERR moab.tag_set_data(interp_tag, &vert, 1, &vert_data[0]);
 
  370    CHKERR moab.tag_delete(scalar_tag);
 
  371    CHKERR moab.tag_delete(adj_count_tag);
 
  376      if (interp_order == 1) {
 
  381        std::vector<Tag> tags;
 
  382        tags.push_back(interp_tag);
 
  383        pcs[1]->reduce_tags(tags, tags, MPI_SUM, targ_verts);
 
  385        for (
auto &tet : targ_elems) {
 
  387          CHKERR moab.get_connectivity(&tet, 1, adj_verts, 
true);
 
  389          std::vector<double> adj_vert_data(adj_verts.size() * interp_tag_len,
 
  391          CHKERR moab.tag_get_data(interp_tag, adj_verts, &adj_vert_data[0]);
 
  393          std::vector<double> tet_data(interp_tag_len, 0.0);
 
  394          for (
int itag = 0; itag < interp_tag_len; itag++) {
 
  395            for (
int i = 0; 
i < adj_verts.size(); 
i++) {
 
  396              tet_data[itag] += adj_vert_data[
i * interp_tag_len + itag];
 
  398            tet_data[itag] /= adj_verts.size();
 
  400          CHKERR moab.tag_set_data(interp_tag, &tet, 1, &tet_data[0]);
 
  404      std::vector<double> data_integ(interp_tag_len, 0.0);
 
  405      for (
auto &tet : targ_elems) {
 
  407        std::vector<double> tet_data(interp_tag_len, 0.0);
 
  408        CHKERR moab.tag_get_data(interp_tag, &tet, 1, &tet_data[0]);
 
  412        CHKERR moab.get_connectivity(tet, vert_conn, vert_num, 
true);
 
  413        std::vector<double> vpos(3 * vert_num);
 
  414        CHKERR moab.get_coords(vert_conn, vert_num, &vpos[0]);
 
  417        for (
int itag = 0; itag < interp_tag_len; itag++) {
 
  418          data_integ[itag] += tet_data[itag] * vol;
 
  422      std::vector<double> global_data_integ(interp_tag_len, 0.0);
 
  423      MPI_Allreduce(&data_integ[0], &global_data_integ[0], interp_tag_len,
 
  424                    MPI_DOUBLE, MPI_SUM, m_field.
get_comm());
 
  426      std::string 
temp = 
"Integrated stress for sslv116 test: ";
 
  427      for (
int i = 0; 
i < 9; ++
i) {
 
  428        std::stringstream ss;
 
  429        ss << std::scientific << std::setprecision(3) << global_data_integ[
i];
 
  430        temp += ss.str() + 
" ";
 
  434      double non_zero_val = 1.655e12;
 
  435      double non_zero_tol = 5.1e-2;
 
  436      if (interp_order == 1) {
 
  439      std::vector<int> non_zero_inds(3);
 
  440      std::vector<int> zero_inds(6);
 
  442        non_zero_inds = {0, 4, 8};
 
  443        zero_inds = {1, 2, 3, 5, 6, 7};
 
  445        non_zero_inds = {0, 1, 2};
 
  446        zero_inds = {3, 4, 5, 6, 7, 8};
 
  448      bool non_zero_check =
 
  449          all_of(begin(non_zero_inds), end(non_zero_inds), [&](
int i) {
 
  450            return abs(global_data_integ[
i] - non_zero_val) / non_zero_val <
 
  453      bool zero_check = all_of(begin(zero_inds), end(zero_inds), [&](
int i) {
 
  454        return abs(global_data_integ[
i]) < 1e-12;
 
  456      if (!non_zero_check || !zero_check) {
 
  458                "Wrong value of the integrated stress");
 
  467    part_sets.insert(roots[1]);
 
  468    std::string new_write_opts;
 
  469    std::ostringstream extra_opt;
 
  471      int index = pcs[1]->get_id();
 
  472      extra_opt << 
";PARALLEL_COMM=" << index;
 
  474    new_write_opts = write_opts + extra_opt.str();
 
  475    CHKERR moab.write_file(mesh_out_file, NULL, new_write_opts.c_str(),
 
  478      MOFEM_LOG(
"WORLD", Sev::inform) << 
"Wrote file " << mesh_out_file;