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.