21                                 {
   23 
   24  try {
   25 
   26    moab::Core mb_instance;
   27    moab::Interface &moab = mb_instance;
   28 
   29    
   32 
   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";
   37 
   38    int interp_order = 1;
   39    PetscBool hybrid_interp = PETSC_TRUE;
   40 
   42 
   43    double toler = 5.e-10;
   44 
   45    PetscOptionsBegin(m_field.
get_comm(), 
"", 
"mesh data transfer tool",
 
   46                      "none");
   47    CHKERR PetscOptionsString(
"-source_file", 
"source mesh file name", 
"",
 
   48                              "source.h5m", mesh_source_file, 255,
   49                              PETSC_NULLPTR);
   50    CHKERR PetscOptionsString(
"-target_file", 
"target mesh file name", 
"",
 
   51                              "target.h5m", mesh_target_file, 255,
   52                              PETSC_NULLPTR);
   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,
   57                              PETSC_NULLPTR);
   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,
 
   63                           PETSC_NULLPTR);
   64 
   65    PetscOptionsEnd();
   66 
   69 
   70    Coupler::Method method;
   71    switch (interp_order) {
   72    case 0:
   73      method = Coupler::CONSTANT;
   74      break;
   75    case 1:
   76      method = Coupler::LINEAR_FE;
   77      break;
   78    default:
   80              "Unsupported interpolation order");
   81    }
   82 
   83    std::vector<std::string> mesh_files(2);
   84    mesh_files[0] = string(mesh_source_file);
   85    mesh_files[1] = string(mesh_target_file);
   86 
   87    int nprocs, rank;
   88    ierr = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
 
   90    ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
   92 
   93    std::string read_opts, write_opts;
   94    read_opts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_"
   95                "DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS";
   96    if (nprocs > 1)
   97      read_opts += ";PARALLEL_GHOSTS=3.0.1";
   98    write_opts = (nprocs > 1) ? "PARALLEL=WRITE_PART" : "";
   99 
  100    std::vector<boost::shared_ptr<ParallelComm>> pcs(mesh_files.size());
  101    std::vector<EntityHandle> roots(mesh_files.size());
  102 
  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();
  110 
  111      CHKERR moab.create_meshset(MESHSET_SET, roots[
i]);
 
  112 
  113      CHKERR moab.load_file(mesh_files[
i].c_str(), &roots[
i],
 
  114                            newread_opts.c_str());
  115    }
  116 
  117    Range src_elems, targ_elems, targ_verts;
 
  118    CHKERR pcs[0]->get_part_entities(src_elems, 3);
 
  119 
  121    CHKERR moab.tag_get_handle(iterp_tag_name, interp_tag);
 
  122 
  123    int interp_tag_len;
  124    CHKERR moab.tag_get_length(interp_tag, interp_tag_len);
 
  125 
  126    if (interp_tag_len != 1 && interp_tag_len != 3 && interp_tag_len != 9) {
  128              "Unsupported interpolation tag length: %d", interp_tag_len);
  129    }
  130 
  131    Coupler mbc(&moab, pcs[0].get(), src_elems, 0, true);
  132 
  133    std::vector<double> vpos; 
  134    int num_pts = 0;
  135 
  137 
  138    
  139    CHKERR pcs[1]->get_part_entities(targ_elems, 3);
 
  140    if (interp_order == 0) {
  141      targ_verts = targ_elems;
  142    } else {
  143      CHKERR moab.get_adjacencies(targ_elems, 0, 
false, targ_verts,
 
  144                                  moab::Interface::UNION);
  145    }
  146 
  147    
  148    CHKERR pcs[1]->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);
 
  149    targ_verts = subtract(targ_verts, tmp_verts);
  150 
  151    
  152    num_pts = (
int)targ_verts.size();
 
  153    vpos.resize(3 * targ_verts.size());
  154    CHKERR moab.get_coords(targ_verts, &vpos[0]);
 
  155 
  156    
  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);
 
  160 
  161    
  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);
  173          missing_pts_num++;
  174        } else {
  175          vit++;
  176        }
  177      }
  178 
  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. ";
  186      }
  187 
  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]);
 
  192        tl_ptr->reset();
  193        CHKERR mbc.locate_points(&vpos[0], num_pts, 0, toler, tl_ptr.get(),
 
  194                                 false);
  195      }
  197    };
  198 
  200    CHKERR find_missing_points(targ_verts, num_pts, vpos, missing_verts);
 
  201 
  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);
  204 
  205    CHKERR moab.tag_get_data(interp_tag, src_elems, &source_data[0]);
 
  206 
  207    Tag scalar_tag, adj_count_tag;
 
  208    double def_scl = 0;
  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,
  212                               &def_scl);
  213 
  214    string adj_count_tag_name = "ADJ_COUNT";
  215    double def_adj = 0;
  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,
  218                               &def_adj);
  219 
  220    
  221    
  222    auto create_scalar_tags = [&](
const Range &src_elems,
 
  223                                  const std::vector<double> &source_data,
  224                                  int itag) {
  226 
  227      std::vector<double> source_data_scalar(src_elems.size());
  228      
  229      for (int ielem = 0; ielem < src_elems.size(); ielem++) {
  230        source_data_scalar[ielem] = source_data[itag + ielem * interp_tag_len];
  231      }
  232 
  233      
  234      CHKERR moab.tag_set_data(scalar_tag, src_elems, &source_data_scalar[0]);
 
  235 
  236      if (interp_order == 1) {
  237        
  239        CHKERR moab.get_connectivity(src_elems, src_verts, 
true);
 
  240 
  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);
 
  243 
  244        for (auto &tet : src_elems) {
  245          double tet_data = 0;
  246          CHKERR moab.tag_get_data(scalar_tag, &tet, 1, &tet_data);
 
  247 
  249          CHKERR moab.get_connectivity(&tet, 1, adj_verts, 
true);
 
  250 
  251          std::vector<double> adj_vert_data(adj_verts.size(), 0.0);
  252          std::vector<double> adj_vert_count(adj_verts.size(), 0.0);
  253 
  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,
 
  256                                   &adj_vert_count[0]);
  257 
  258          for (int ivert = 0; ivert < adj_verts.size(); ivert++) {
  259            adj_vert_data[ivert] += tet_data;
  260            adj_vert_count[ivert] += 1;
  261          }
  262 
  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,
 
  265                                   &adj_vert_count[0]);
  266        }
  267 
  268        
  269        std::vector<Tag> tags = {scalar_tag, adj_count_tag};
  270        pcs[0]->reduce_tags(tags, tags, MPI_SUM, src_verts);
  271 
  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);
  274 
  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]);
  278 
  279        for (int ivert = 0; ivert < src_verts.size(); ivert++) {
  280          src_vert_data[ivert] /= src_vert_adj_count[ivert];
  281        }
  282        CHKERR moab.tag_set_data(scalar_tag, src_verts, &src_vert_data[0]);
 
  283      }
  285    };
  286 
  287    for (int itag = 0; itag < interp_tag_len; itag++) {
  288 
  289      CHKERR create_scalar_tags(src_elems, source_data, itag);
 
  290 
  291      std::vector<double> target_data_scalar(num_pts, 0.0);
  292      CHKERR mbc.interpolate(method, scalar_tag_name, &target_data_scalar[0],
 
  293                             tl_ptr.get());
  294 
  295      for (int ielem = 0; ielem < num_pts; ielem++) {
  296        target_data[itag + ielem * interp_tag_len] = target_data_scalar[ielem];
  297      }
  298    }
  299 
  300    
  301    CHKERR moab.tag_set_data(interp_tag, targ_verts, &target_data[0]);
 
  302 
  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);
  309 
  310      int num_adj_elems = (
int)missing_adj_elems.size();
 
  311      std::vector<double> vpos_adj_elems;
  312 
  313      vpos_adj_elems.resize(3 * missing_adj_elems.size());
  314      CHKERR moab.get_coords(missing_adj_elems, &vpos_adj_elems[0]);
 
  315 
  316      
  317      tl_ptr->reset();
  318      CHKERR mbc.locate_points(&vpos_adj_elems[0], num_adj_elems, 0, toler,
 
  319                               tl_ptr.get(), false);
  320 
  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. ";
  328      }
  329 
  330      std::vector<double> target_data_adj_elems(interp_tag_len * num_adj_elems,
  331                                                0.0);
  332 
  333      for (int itag = 0; itag < interp_tag_len; itag++) {
  334        CHKERR create_scalar_tags(src_elems, source_data, itag);
 
  335 
  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());
  339 
  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];
  343        }
  344      }
  345 
  346      CHKERR moab.tag_set_data(interp_tag, missing_adj_elems,
 
  347                               &target_data_adj_elems[0]);
  348 
  349      
  350      for (auto &vert : missing_verts) {
  352        CHKERR moab.get_adjacencies(&vert, 1, 3, 
false, adj_elems,
 
  353                                    moab::Interface::UNION);
  354 
  355        std::vector<double> adj_elems_data(adj_elems.size() * interp_tag_len,
  356                                           0.0);
  357        CHKERR moab.tag_get_data(interp_tag, adj_elems, &adj_elems_data[0]);
 
  358 
  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];
 
  363          }
  364          vert_data[itag] /= adj_elems.size();
  365        }
  366        CHKERR moab.tag_set_data(interp_tag, &vert, 1, &vert_data[0]);
 
  367      }
  368    }
  369 
  370    CHKERR moab.tag_delete(scalar_tag);
 
  371    CHKERR moab.tag_delete(adj_count_tag);
 
  372 
  374      
  375 
  376      if (interp_order == 1) {
  377        
  378        
  379 
  380        
  381        std::vector<Tag> tags;
  382        tags.push_back(interp_tag);
  383        pcs[1]->reduce_tags(tags, tags, MPI_SUM, targ_verts);
  384 
  385        for (auto &tet : targ_elems) {
  387          CHKERR moab.get_connectivity(&tet, 1, adj_verts, 
true);
 
  388 
  389          std::vector<double> adj_vert_data(adj_verts.size() * interp_tag_len,
  390                                            0.0);
  391          CHKERR moab.tag_get_data(interp_tag, adj_verts, &adj_vert_data[0]);
 
  392 
  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];
 
  397            }
  398            tet_data[itag] /= adj_verts.size();
  399          }
  400          CHKERR moab.tag_set_data(interp_tag, &tet, 1, &tet_data[0]);
 
  401        }
  402      }
  403 
  404      std::vector<double> data_integ(interp_tag_len, 0.0);
  405      for (auto &tet : targ_elems) {
  406 
  407        std::vector<double> tet_data(interp_tag_len, 0.0);
  408        CHKERR moab.tag_get_data(interp_tag, &tet, 1, &tet_data[0]);
 
  409 
  411        int vert_num;
  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]);
 
  416 
  417        for (int itag = 0; itag < interp_tag_len; itag++) {
  418          data_integ[itag] += tet_data[itag] * vol;
  419        }
  420      }
  421 
  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());
 
  425 
  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() + 
" ";
 
  431      }
  433 
  434      double non_zero_val = 1.655e12;
  435      double non_zero_tol = 5.1e-2;
  436      if (interp_order == 1) {
  437        non_zero_tol = 1e-3;
  438      }
  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};
  444      } else {
  445        non_zero_inds = {0, 1, 2};
  446        zero_inds = {3, 4, 5, 6, 7, 8};
  447      }
  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 <
 
  451                   non_zero_tol;
  452          });
  453      bool zero_check = all_of(begin(zero_inds), end(zero_inds), [&](
int i) {
 
  454        return abs(global_data_integ[
i]) < 1e-12;
 
  455      });
  456      if (!non_zero_check || !zero_check) {
  458                "Wrong value of the integrated stress");
  459      }
  463    }
  464 
  466    
  467    part_sets.insert(roots[1]);
  468    std::string new_write_opts;
  469    std::ostringstream extra_opt;
  470    if (nprocs > 1) {
  471      int index = pcs[1]->get_id();
  472      extra_opt << ";PARALLEL_COMM=" << index;
  473    }
  474    new_write_opts = write_opts + extra_opt.str();
  475    CHKERR moab.write_file(mesh_out_file, NULL, new_write_opts.c_str(),
 
  476                           part_sets);
  477    if (0 == rank) {
  478      MOFEM_LOG(
"WORLD", Sev::inform) << 
"Wrote file " << mesh_out_file;
 
  479    }
  480  }
  482 
  484 
  485  return 0;
  486}
#define CATCH_ERRORS
Catch errors.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
void temp(int x, int y=10)
virtual MPI_Comm & get_comm() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.