v0.15.0
Loading...
Searching...
No Matches
mesh_data_transfer.cpp File Reference

Transfer data from source mesh to target mesh. More...

#include <MoFEM.hpp>
#include <mbcoupler/Coupler.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Detailed Description

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

Definition in file mesh_data_transfer.cpp.

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 21 of file mesh_data_transfer.cpp.

21 {
22 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
23
24 try {
25
26 moab::Core mb_instance;
27 moab::Interface &moab = mb_instance;
28
29 // Create MoFEM database
30 MoFEM::Core core(moab);
31 MoFEM::Interface &m_field = core;
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
41 int atom_test = 0;
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
67 MOFEM_LOG_CHANNEL("WORLD");
68 MOFEM_LOG_TAG("WORLD", "mesh_data_transfer");
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:
79 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
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);
89 CHKERRQ(ierr);
90 ierr = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
91 CHKERRQ(ierr);
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
120 Tag interp_tag;
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) {
127 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
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; // the positions we are interested in
134 int num_pts = 0;
135
136 Range tmp_verts;
137
138 // First get all vertices adj to partition entities in target mesh
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 // Then get non-owned verts and subtract
148 CHKERR pcs[1]->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);
149 targ_verts = subtract(targ_verts, tmp_verts);
150
151 // get position of these entities; these are the target points
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 // Locate those points in the source mesh
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 // If some points were not located, we need to process them
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;
167 int i = 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,
181 MPI_SUM, m_field.get_comm());
182 if (missing_pts_num_global) {
183 MOFEM_LOG("WORLD", Sev::warning)
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
199 Range missing_verts;
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 // MBCoupler functionality supports only scalar tags. For the case of
221 // vector or tensor tags we need to save each component as a scalar tag
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 // Populate source_data_scalar
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 // Set data on the scalar tag
234 CHKERR moab.tag_set_data(scalar_tag, src_elems, &source_data_scalar[0]);
235
236 if (interp_order == 1) {
237 // Linear interpolation: compute average value of data on vertices
238 Range src_verts;
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
248 Range adj_verts;
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 // Reduce tags for the parallel case
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 // Use original tag
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 // Locate those points in the source mesh
317 tl_ptr->reset();
318 CHKERR mbc.locate_points(&vpos_adj_elems[0], num_adj_elems, 0, toler,
319 tl_ptr.get(), false);
320
321 Range missing_tets;
322 CHKERR find_missing_points(missing_adj_elems, num_adj_elems,
323 vpos_adj_elems, missing_tets);
324 if (missing_tets.size()) {
325 MOFEM_LOG("WORLD", Sev::warning)
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 // FIXME: add implementation for parallel case
350 for (auto &vert : missing_verts) {
351 Range adj_elems;
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
373 if (atom_test == 1 || atom_test == 2) {
374 // Check the integrated stress for the SSLV116 test
375
376 if (interp_order == 1) {
377 // If the data was interpolated to target mesh vertices, we need to
378 // average the vertex values for each 3D element and save on its tag
379
380 // We need to reduce tags for the parallel case
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) {
386 Range adj_verts;
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
410 const EntityHandle *vert_conn;
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]);
415 double vol = Tools::tetVolume(vpos.data());
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 }
432 MOFEM_LOG("WORLD", Sev::inform) << temp;
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);
441 if (atom_test == 1) {
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) {
457 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
458 "Wrong value of the integrated stress");
459 }
460 } else if (atom_test) {
461 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED,
462 "Unknown test number: %d", atom_test);
463 }
464
465 Range part_sets;
466 // Only save the target mesh
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
Definition definitions.h:40
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#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 char help[]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
int atom_test
Atom test.
Definition plastic.cpp:121
void temp(int x, int y=10)
Definition simple.cpp:4
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition Tools.cpp:30

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 19 of file mesh_data_transfer.cpp.