1 // Copyright (c) Lawrence Livermore National Security, LLC and other Conduit
2 // Project developers. See top-level LICENSE AND COPYRIGHT files for dates and
3 // other details. No copyright assignment is required to contribute to Conduit.
4
5 //-----------------------------------------------------------------------------
6 ///
7 /// file: conduit_blueprint_mpi_mesh.cpp
8 ///
9 //-----------------------------------------------------------------------------
10
11 //-----------------------------------------------------------------------------
12 // std lib includes
13 //-----------------------------------------------------------------------------
14 #include <algorithm>
15 #include <tuple>
16 #include <vector>
17 #include <cmath>
18
19 //-----------------------------------------------------------------------------
20 // conduit includes
21 //-----------------------------------------------------------------------------
22 #include "conduit_blueprint_mesh.hpp"
23 #include "conduit_blueprint_mesh_utils.hpp"
24 #include "conduit_blueprint_mpi_mesh.hpp"
25 #include "conduit_blueprint_mpi_mesh_flatten.hpp"
26 #include "conduit_blueprint_mpi_mesh_partition.hpp"
27 #include "conduit_blueprint_mesh_utils.hpp"
28 #include "conduit_blueprint_o2mrelation.hpp"
29 #include "conduit_blueprint_o2mrelation_iterator.hpp"
30 #include "conduit_relay_mpi.hpp"
31 #include <assert.h>
32 #include <cmath>
33 #include <limits>
34 #include <list>
35 // access conduit blueprint mesh utilities
36 namespace bputils = conduit::blueprint::mesh::utils;
37
38 #include "conduit_blueprint_mpi_mesh.hpp"
39
40 // TODO(JRC): Consider moving an improved version of this type to a more accessible location.
41 struct ffloat64
42 {
43 conduit::float64 data;
44
ffloat64ffloat6445 ffloat64(conduit::float64 input = 0.0)
46 {
47 data = input;
48 }
49
operator conduit::float64ffloat6450 operator conduit::float64() const
51 {
52 return this->data;
53 }
54
operator <ffloat6455 bool operator<(const ffloat64 &other) const
56 {
57 return this->data < other.data && std::abs(this->data - other.data) > CONDUIT_EPSILON;
58 }
59 };
60
61 // access conduit blueprint mesh utilities
62 namespace bputils = conduit::blueprint::mesh::utils;
63 // access one-to-many index types
64 namespace O2MIndex = conduit::blueprint::o2mrelation;
65
66 // typedefs for verbose but commonly used types
67 typedef std::tuple<conduit::Node*, conduit::Node*, conduit::Node*> DomMapsTuple;
68 typedef std::tuple<ffloat64, ffloat64, ffloat64> PointTuple;
69 // typedefs to enable passing around function pointers
70 typedef void (*GenDerivedFun)(const conduit::Node&, conduit::Node&, conduit::Node&, conduit::Node&);
71 typedef void (*GenDecomposedFun)(const conduit::Node&, conduit::Node&, conduit::Node&, conduit::Node&, conduit::Node&);
72 typedef conduit::index_t (*IdDecomposedFun)(const bputils::TopologyMetadata&, const conduit::index_t /*ei*/, const conduit::index_t /*di*/);
73 typedef std::vector<conduit::index_t> (*CalcDimDecomposedFun)(const bputils::ShapeType&);
74
75 //-----------------------------------------------------------------------------
76 // -- begin conduit --
77 //-----------------------------------------------------------------------------
78 namespace conduit
79 {
80
81 //-----------------------------------------------------------------------------
82 // -- begin conduit::blueprint --
83 //-----------------------------------------------------------------------------
84 namespace blueprint
85 {
86
87 //-----------------------------------------------------------------------------
88 // -- begin conduit::blueprint::mpi --
89 //-----------------------------------------------------------------------------
90 namespace mpi
91 {
92
93 //-----------------------------------------------------------------------------
94 // -- begin conduit::blueprint::mesh --
95 //-----------------------------------------------------------------------------
96
97 namespace mesh
98 {
99 //-----------------------------------------------------------------------------
100
101 // TODO(JRC): Consider moving these structures somewhere else.
102
103 //-------------------------------------------------------------------------
104 struct SharedFace
105 {
106 index_t m_face_id;
107 std::vector<index_t> m_fine_subelem;
108 };
109
110
111 struct PolyBndry
112 {
113 index_t side; //which 3D side 0-5
114 index_t m_nbr_rank;
115 index_t m_nbr_id;
116 size_t m_nbrs_per_face;
117 std::vector<index_t> m_elems; //elems of nbr domain that touch side
118 std::map<index_t, index_t> m_bface; //map from nbr elem to face of nbr elem
119 std::map<index_t, std::vector<index_t> > m_nbr_elems; //map from local
120 //elem to all
121 //nbr elems that
122 //touch it
123 //outer map: local elem, inner map: nbr elem to face
124 std::map<index_t, std::map<index_t, SharedFace> > m_nbr_faces;
125 };
126
127 //-----------------------------------------------------------------------------
128 // blueprint::mesh::index protocol interface
129 //-----------------------------------------------------------------------------
130
131 //-----------------------------------------------------------------------------
132 bool
verify(const conduit::Node & n,conduit::Node & info,MPI_Comm comm)133 verify(const conduit::Node &n,
134 conduit::Node &info,
135 MPI_Comm comm)
136 {
137 int par_size = relay::mpi::size(comm);
138
139 // NOTE(JRC): MPI tasks without any domains should use a multi-domain
140 // format with empty contents (i.e. an empty object or list node).
141 int local_verify_ok = conduit::blueprint::mesh::verify(n, info) ? 1 : 0;
142 int global_verify_ok = 0;
143
144 Node n_snd, n_reduce;
145 // make sure some MPI task actually had bp data
146 n_snd.set_external(&local_verify_ok,1);
147 n_reduce.set_external(&global_verify_ok,1);
148
149 relay::mpi::sum_all_reduce(n_snd, n_reduce, comm);
150 return global_verify_ok == par_size;
151 }
152
153 //-------------------------------------------------------------------------
154 void
generate_index(const conduit::Node & mesh,const std::string & ref_path,Node & index_out,MPI_Comm comm)155 generate_index(const conduit::Node &mesh,
156 const std::string &ref_path,
157 Node &index_out,
158 MPI_Comm comm)
159 {
160 // we need a list of all possible topos, coordsets, etc
161 // for the blueprint index in the root file.
162 //
163 // across ranks, domains may be sparse
164 // for example: a topo may only exist in one domain
165 // so we union all local mesh indices, and then
166 // se an all gather and union the results together
167 // to create an accurate global index.
168
169 index_t local_num_domains = blueprint::mesh::number_of_domains(mesh);
170 // note:
171 // find global # of domains w/o conduit_blueprint_mpi for now
172 // since we aren't yet linking conduit_blueprint_mpi
173 Node n_src, n_reduce;
174 n_src = local_num_domains;
175
176 relay::mpi::sum_all_reduce(n_src,
177 n_reduce,
178 comm);
179
180 index_t global_num_domains = n_reduce.to_int();
181
182 index_out.reset();
183
184 Node local_idx, gather_idx;
185
186 if(local_num_domains > 0)
187 {
188 ::conduit::blueprint::mesh::generate_index(mesh,
189 ref_path,
190 global_num_domains,
191 local_idx);
192 }
193
194 relay::mpi::all_gather_using_schema(local_idx,
195 gather_idx,
196 comm);
197
198 // union all entries into final index that reps
199 // all domains
200 NodeConstIterator itr = gather_idx.children();
201 while(itr.has_next())
202 {
203 const Node &curr = itr.next();
204 index_out.update(curr);
205 }
206 }
207
208
209 //-----------------------------------------------------------------------------
generate_domain_to_rank_map(const conduit::Node & mesh,Node & domain_to_rank_map,MPI_Comm comm)210 void generate_domain_to_rank_map(const conduit::Node &mesh,
211 Node &domain_to_rank_map,
212 MPI_Comm comm)
213 {
214 int64 par_rank = relay::mpi::rank(comm);
215 int64 max_local_id = -1;
216
217 std::vector<const Node *> domains = ::conduit::blueprint::mesh::domains(mesh);
218 std::vector<int64> local_domains;
219 for(index_t di = 0; di < (index_t)domains.size(); di++)
220 {
221 const conduit::Node &domain = *domains[di];
222
223 int64 domain_id = par_rank;
224 if(domain.has_child("state") && domain["state"].has_child("domain_id"))
225 {
226 domain_id = domain["state/domain_id"].to_int64();
227 }
228 local_domains.push_back(domain_id);
229
230 max_local_id = (domain_id > max_local_id) ? domain_id : max_local_id;
231 }
232
233 Node max_local, max_global;
234 max_local.set_int64(max_local_id);
235 max_global.set_int64(-1);
236 relay::mpi::max_all_reduce(max_local, max_global, comm);
237
238 std::vector<int64> local_map(max_global.as_int64() + 1, -1);
239 for(auto m_itr = local_domains.begin(); m_itr != local_domains.end(); ++m_itr)
240 {
241 local_map[*m_itr] = par_rank;
242 }
243
244 Node local_par;
245 local_par.set_external(&local_map[0], local_map.size());
246
247 relay::mpi::max_all_reduce(local_par, domain_to_rank_map, comm);
248 }
249
250
251 //-----------------------------------------------------------------------------
252 index_t
number_of_domains(const conduit::Node & n,MPI_Comm comm)253 number_of_domains(const conduit::Node &n,
254 MPI_Comm comm)
255 {
256 // called only when mesh bp very is true, simplifies logic here
257 index_t local_num_domains = 0;
258 if(!n.dtype().is_empty())
259 {
260 local_num_domains = ::conduit::blueprint::mesh::number_of_domains(n);
261 }
262
263 index_t global_num_domains = 0;
264
265 Node n_snd, n_reduce;
266 // count all domains with mpi
267 n_snd.set_external(&local_num_domains,1);
268 n_reduce.set_external(&global_num_domains,1);
269
270 relay::mpi::all_reduce(n_snd, n_reduce, MPI_SUM, comm);
271 return global_num_domains;
272 }
273
274 //-------------------------------------------------------------------------
to_polytopal(const Node & n,Node & dest,const std::string & name,MPI_Comm comm)275 void to_polytopal(const Node &n,
276 Node &dest,
277 const std::string& name,
278 MPI_Comm comm)
279 {
280
281 const std::vector<const conduit::Node *> doms = ::conduit::blueprint::mesh::domains(n);
282
283 // make sure all topos match
284 index_t ok = 1;
285 index_t num_doms = (index_t) doms.size();
286 index_t dims = -1;
287
288 for (const auto& dom_ptr : doms)
289 {
290 if(dom_ptr->fetch("topologies").has_child(name))
291 {
292 const Node &topo = dom_ptr->fetch("topologies")[name];
293 if(topo["type"].as_string() == "structured")
294 {
295 dims = topo["elements/dims"].number_of_children();
296 }
297 else
298 {
299 ok = 0;
300 }
301 }else
302 {
303 ok = 0;
304 }
305 }
306
307 // reduce and check for consistency (all ok, and all doms are either 2d or 3d)
308 Node local, gather;
309 local.set(DataType::index_t(3));
310 index_t_array local_vals = local.value();
311 local_vals[0] = ok;
312 local_vals[1] = num_doms;
313 local_vals[2] = dims;
314
315 // Note: this might be more efficient as
316 // a set of flat gathers into separate arrays
317 relay::mpi::all_gather_using_schema(local,
318 gather,
319 comm);
320
321 NodeConstIterator gitr = gather.children();
322 index_t gather_dims = -1;
323 while(gitr.has_next() && (ok == 1))
324 {
325 const Node &curr = gitr.next();
326 index_t_array gather_vals = curr.value();
327 if(gather_vals[0] != 1)
328 {
329 ok = 0;
330 }
331 else
332 {
333 // this proc has domains and we haven't inited dims
334 if( gather_vals[1] > 0 && gather_dims == -1)
335 {
336 gather_dims = gather_vals[2];
337 }
338 else if(gather_vals[1] > 0) // this proc has domains
339 {
340 if(gather_dims != gather_vals[2])
341 {
342 ok = 0;
343 }
344 }
345 }
346 }
347
348 if(ok == 1)
349 {
350 if(dims == 2)
351 {
352 to_polygonal(n,dest,name,comm);
353 }
354 else if(dims == 3)
355 {
356 to_polyhedral(n,dest,name,comm);
357 }
358 else
359 {
360 CONDUIT_ERROR("to_polytopal only supports 2d or 3d structured toplogies"
361 " (passed mesh has dims = " << dims << ")");
362 }
363 }
364 else
365 {
366 CONDUIT_ERROR("to_polytopal only supports structured toplogies");
367 }
368 }
369
370 //-------------------------------------------------------------------------
to_polygonal(const Node & n,Node & dest,const std::string & name,MPI_Comm comm)371 void to_polygonal(const Node &n,
372 Node &dest,
373 const std::string& name,
374 MPI_Comm comm)
375 {
376 // Helper Functions //
377
378 const static auto gen_default_name = [] (const std::string &prefix, const index_t id)
379 {
380 std::ostringstream oss;
381 oss << prefix << "_" << std::setw(6) << std::setfill('0') << id;
382 return oss.str();
383 };
384
385 const static auto is_domain_local = [] (const conduit::Node &root, const std::string &name)
386 {
387 return ::conduit::blueprint::mesh::is_multi_domain(root) && root.has_child(name);
388 };
389
390 // Implementation //
391
392 // TODO(JRC): Use widest data type available in given Node for both ints and doubles
393 // TODO(JRC): Communicate this across the network?
394
395 dest.reset();
396
397 Node temp;
398
399 std::map<index_t, std::map<index_t, std::vector<index_t> > > poly_elems_map;
400 std::map<index_t, std::map<index_t, std::vector<double> > > dom_to_nbr_to_xbuffer;
401 std::map<index_t, std::map<index_t, std::vector<double> > > dom_to_nbr_to_ybuffer;
402
403 const std::vector<const conduit::Node *> doms = ::conduit::blueprint::mesh::domains(n);
404 for (index_t si = 0; si < 3; si++)
405 {
406 // TODO(JRC): Figure out better names for all of these stages.
407 // TODO(JRC): Consider using a local enumeration to track this
408 // instead, to have compile-time consistency.
409 const std::string stage_name = (
410 (si == 0) ? "send_missing_windows" : (
411 (si == 1) ? "fill_all_windows" : (
412 (si == 2) ? "generate_polygons" : (""))));
413
414 for (const auto& dom_ptr : doms)
415 {
416 const Node &dom = *dom_ptr;
417 Node &dest_dom = ::conduit::blueprint::mesh::is_multi_domain(n) ?
418 dest[dom.name()] : dest;
419
420 const Node &in_topo = dom["topologies"][name];
421 const Node *in_cset =
422 bputils::find_reference_node(in_topo, "coordset");
423 Node &out_cset = dest_dom["coordsets"][in_topo["coordset"].as_string()];
424
425 const index_t domain_id = dom["state/domain_id"].to_index_t();
426 const index_t level_id = dom["state/level_id"].to_index_t();
427 const std::string ref_name = gen_default_name("window", domain_id);
428
429 const index_t i_lo = in_topo["elements/origin/i0"].to_index_t();
430 const index_t j_lo = in_topo["elements/origin/j0"].to_index_t();
431 const index_t iwidth = in_topo["elements/dims/i"].to_index_t();
432 const index_t jwidth = in_topo["elements/dims/j"].to_index_t();
433 const index_t niwidth = in_topo["elements/dims/i"].to_index_t() + 1;
434
435 // TODO(JRC): There's an implicit assumption here about each
436 // domain having a 'state' path, which is optional in the general
437 // case. This needs to be checked and warned about before any
438 // processing occurs.
439 if (si == 0)
440 {
441 if (level_id == 0)
442 {
443 continue;
444 }
445 }
446 else if (si == 2)
447 {
448 if (blueprint::mesh::coordset::uniform::verify(*in_cset, temp))
449 {
450 blueprint::mesh::coordset::uniform::to_explicit(*in_cset, out_cset);
451 }
452 else if (blueprint::mesh::coordset::rectilinear::verify(*in_cset, temp))
453 {
454 blueprint::mesh::coordset::rectilinear::to_explicit(*in_cset, out_cset);
455 }
456 else
457 {
458 out_cset.set(*in_cset);
459 }
460 }
461
462 auto &poly_elems = poly_elems_map[domain_id];
463 auto &nbr_to_xbuffer = dom_to_nbr_to_xbuffer[domain_id];
464 auto &nbr_to_ybuffer = dom_to_nbr_to_ybuffer[domain_id];
465
466 if (dom.has_path("adjsets/adjset/groups"))
467 {
468 NodeConstIterator grp_itr = dom["adjsets/adjset/groups"].children();
469 while(grp_itr.has_next())
470 {
471 const Node& group = grp_itr.next();
472 if (group.has_child("neighbors") && group.has_child("windows"))
473 {
474 temp.reset();
475 temp.set_external(DataType(group["neighbors"].dtype().id(), 1),
476 (void*)group["neighbors"].element_ptr(1));
477 const index_t nbr_id = temp.to_index_t();
478 const std::string nbr_name = gen_default_name("domain", nbr_id);
479
480 const Node &in_windows = group["windows"];
481 const std::string nbr_win_name = gen_default_name("window", nbr_id);
482
483 const Node &ref_win = in_windows[ref_name];
484 const Node &nbr_win = in_windows[nbr_win_name];
485
486 const index_t nbr_level = nbr_win["level_id"].to_index_t();
487 const index_t ref_level = ref_win["level_id"].to_index_t();
488
489 if ((si == 0 && nbr_level < ref_level) ||
490 (si != 0 && nbr_level > ref_level))
491 {
492 const index_t ref_size_i = ref_win["dims/i"].to_index_t();
493 const index_t ref_size_j = ref_win["dims/j"].to_index_t();
494 const index_t ref_size = ref_size_i * ref_size_j;
495
496 const index_t nbr_size_i = nbr_win["dims/i"].to_index_t();
497 const index_t nbr_size_j = nbr_win["dims/j"].to_index_t();
498 const index_t nbr_size = nbr_size_i * nbr_size_j;
499
500 const index_t ratio_i = nbr_win["ratio/i"].to_index_t();
501 const index_t ratio_j = nbr_win["ratio/j"].to_index_t();
502
503 if (si == 0 && nbr_size < ref_size && !is_domain_local(n, nbr_name))
504 {
505 std::vector<double> xbuffer, ybuffer;
506 const Node &fcoords = (*in_cset)["values"];
507
508 const index_t origin_i = ref_win["origin/i"].to_index_t();
509 const index_t origin_j = ref_win["origin/j"].to_index_t();
510
511 temp.reset();
512
513 if (ref_size_i == 1)
514 {
515 const index_t icnst = origin_i - i_lo;
516 const index_t jstart = origin_j - j_lo;
517 const index_t jend = jstart + ref_size_j;
518 for (index_t jidx = jstart; jidx < jend; ++jidx)
519 {
520 const index_t offset = jidx * niwidth + icnst;
521 temp.set_external(DataType(fcoords["x"].dtype().id(), 1),
522 (void*)fcoords["x"].element_ptr(offset));
523 xbuffer.push_back(temp.to_double());
524 temp.set_external(DataType(fcoords["y"].dtype().id(), 1),
525 (void*)fcoords["y"].element_ptr(offset));
526 ybuffer.push_back(temp.to_double());
527 }
528 }
529 else if (ref_size_j == 1)
530 {
531 const index_t jcnst = origin_j - j_lo;
532 const index_t istart = origin_i - i_lo;
533 const index_t iend = istart + ref_size_i;
534 for (index_t iidx = istart; iidx < iend; ++iidx)
535 {
536 const index_t offset = jcnst * niwidth + iidx;
537 temp.set_external(DataType(fcoords["x"].dtype().id(), 1),
538 (void*)fcoords["x"].element_ptr(offset));
539 xbuffer.push_back(temp.to_double());
540 temp.set_external(DataType(fcoords["y"].dtype().id(), 1),
541 (void*)fcoords["y"].element_ptr(offset));
542 ybuffer.push_back(temp.to_double());
543 }
544 }
545 const index_t nbr_rank = group["rank"].to_index_t();
546 MPI_Send(&xbuffer[0],
547 xbuffer.size(),
548 MPI_DOUBLE,
549 nbr_rank,
550 domain_id,
551 comm);
552 MPI_Send(&ybuffer[0],
553 ybuffer.size(),
554 MPI_DOUBLE,
555 nbr_rank,
556 domain_id,
557 comm);
558 }
559 else if (si == 1 && nbr_size > ref_size)
560 {
561 auto& xbuffer = nbr_to_xbuffer[nbr_id];
562 auto& ybuffer = nbr_to_ybuffer[nbr_id];
563
564 if (!is_domain_local(n, nbr_name))
565 {
566 if (nbr_size_i == 1)
567 {
568 xbuffer.resize(nbr_size_j);
569 ybuffer.resize(nbr_size_j);
570 }
571 else if (nbr_size_j == 1)
572 {
573 xbuffer.resize(nbr_size_i);
574 ybuffer.resize(nbr_size_i);
575 }
576
577 index_t nbr_rank = group["rank"].to_index_t();
578 MPI_Recv(&xbuffer[0],
579 xbuffer.size(),
580 MPI_DOUBLE,
581 nbr_rank,
582 nbr_id,
583 comm,
584 MPI_STATUS_IGNORE);
585 MPI_Recv(&ybuffer[0],
586 ybuffer.size(),
587 MPI_DOUBLE,
588 nbr_rank,
589 nbr_id,
590 comm,
591 MPI_STATUS_IGNORE);
592
593 }
594 else
595 {
596 const Node& nbr_dom = n[nbr_name];
597 const Node& nbr_coords = nbr_dom["coordsets/coords"];
598
599 const Node& ntopo = nbr_dom["topologies"][name];
600 index_t ni_lo = ntopo["elements/origin/i0"].to_index_t();
601 index_t nj_lo = ntopo["elements/origin/j0"].to_index_t();
602 index_t nbr_iwidth = ntopo["elements/dims/i"].to_index_t() + 1;
603
604 // FIXME(JRC): These arrays aren't guaranteed to be double arrays;
605 // this code must be changed to accomodate arbitrary types.
606 const Node& fcoords = nbr_coords["values"];
607 const double_array& xarray = fcoords["x"].as_double_array();
608 const double_array& yarray = fcoords["y"].as_double_array();
609
610 index_t origin_i = nbr_win["origin/i"].to_index_t();
611 index_t origin_j = nbr_win["origin/j"].to_index_t();
612
613 index_t istart = origin_i - ni_lo;
614 index_t jstart = origin_j - nj_lo;
615 index_t iend = istart + nbr_size_i;
616 index_t jend = jstart + nbr_size_j;
617 for (index_t jidx = jstart; jidx < jend; ++jidx)
618 {
619 index_t joffset = jidx*nbr_iwidth;
620 for (index_t iidx = istart; iidx < iend; ++iidx)
621 {
622 index_t offset = joffset+iidx;
623 xbuffer.push_back(xarray[offset]);
624 ybuffer.push_back(yarray[offset]);
625 }
626 }
627 }
628 }
629 else if (si == 2 && ref_size < nbr_size)
630 {
631 bputils::connectivity::create_elements_2d(ref_win,
632 i_lo,
633 j_lo,
634 iwidth,
635 poly_elems);
636
637 const index_t use_ratio = (
638 (nbr_size_j == 1) ? ratio_i : (
639 (nbr_size_i == 1) ? ratio_j : (0)));
640
641 auto& xbuffer = nbr_to_xbuffer[nbr_id];
642 auto& ybuffer = nbr_to_ybuffer[nbr_id];
643
644 const index_t added = (
645 (nbr_size_j == 1) ? (xbuffer.size() - ref_size_i) : (
646 (nbr_size_i == 1) ? (ybuffer.size() - ref_size_j) : (0)));
647
648 // FIXME(JRC): These arrays aren't guaranteed to be double arrays;
649 // this code must be changed to accomodate arbitrary types.
650 const auto& out_x = out_cset["values"]["x"].as_double_array();
651 const auto& out_y = out_cset["values"]["y"].as_double_array();
652 index_t new_vertex = out_x.number_of_elements();
653
654 const index_t out_x_size = out_x.number_of_elements();
655 const index_t out_y_size = out_y.number_of_elements();
656
657 std::vector<double> new_x;
658 std::vector<double> new_y;
659 new_x.reserve(out_x_size + added);
660 new_y.reserve(out_y_size + added);
661 const double* out_x_ptr = static_cast<const double*>(out_x.element_ptr(0));
662 const double* out_y_ptr = static_cast<const double*>(out_y.element_ptr(0));
663
664 new_x.insert(new_x.end(), out_x_ptr, out_x_ptr + out_x_size);
665 new_y.insert(new_y.end(), out_y_ptr, out_y_ptr + out_y_size);
666
667 if ((xbuffer.size()-1)%use_ratio)
668 {
669 new_x.reserve(out_x_size + added*2);
670 }
671 for (index_t ni = 0; ni < (index_t)xbuffer.size(); ++ni)
672 {
673 if (ni % use_ratio)
674 {
675 new_x.push_back(xbuffer[ni]);
676 new_y.push_back(ybuffer[ni]);
677 }
678 }
679
680 out_cset["values"]["x"].set(new_x);
681 out_cset["values"]["y"].set(new_y);
682
683 bputils::connectivity::connect_elements_2d(ref_win,
684 i_lo,
685 j_lo,
686 iwidth,
687 use_ratio,
688 new_vertex,
689 poly_elems);
690 }
691 }
692 }
693 }
694 }
695
696 if (si == 2)
697 {
698 dest_dom["state"].set(dom["state"]);
699 dest_dom["topologies"][name]["coordset"].set(dom["topologies"][name]["coordset"]);
700
701 Node& out_topo = dest_dom["topologies"][name];
702 out_topo["type"].set("unstructured");
703 out_topo["elements/shape"].set("polygonal");
704
705 const index_t elemsize = iwidth*jwidth;
706
707 std::vector<index_t> connect;
708 std::vector<index_t> num_vertices;
709 std::vector<index_t> elem_offsets;
710 index_t offset_sum = 0;
711 for (index_t elem = 0; elem < elemsize; ++elem)
712 {
713 auto elem_itr = poly_elems.find(elem);
714 if (elem_itr == poly_elems.end())
715 {
716 bputils::connectivity::make_element_2d(connect, elem, iwidth);
717 num_vertices.push_back(4);
718 }
719 else
720 {
721 std::vector<index_t>& poly_elem = elem_itr->second;
722 connect.insert(connect.end(), poly_elem.begin(), poly_elem.end());
723 num_vertices.push_back(poly_elem.size());
724 }
725 elem_offsets.push_back(offset_sum);
726 offset_sum += num_vertices.back();
727 }
728
729 // TODO(JRC): Remove extra copy here.
730 out_topo["elements/connectivity"].set(connect);
731 out_topo["elements/sizes"].set(num_vertices);
732 out_topo["elements/offsets"].set(elem_offsets);
733
734 //TODO: zero copy
735 if (dom.has_child("fields"))
736 {
737 dest_dom["fields"].set(dom["fields"]);
738 }
739 }
740 }
741 }
742 }
743
744 void
partition(const conduit::Node & n_mesh,const conduit::Node & options,conduit::Node & output,MPI_Comm comm)745 partition(const conduit::Node &n_mesh,
746 const conduit::Node &options,
747 conduit::Node &output,
748 MPI_Comm comm)
749 {
750 ParallelPartitioner p(comm);
751 output.reset();
752
753 // Partitioners on different ranks ought to return the same value but
754 // perhaps some did not when they examined their own domains against
755 // selection.
756 int iinit, ginit;
757 iinit = p.initialize(n_mesh, options) ? 1 : 0;
758 MPI_Allreduce(&iinit, &ginit, 1, MPI_INT, MPI_MAX, comm);
759 if(ginit > 0)
760 {
761 p.split_selections();
762 p.execute(output);
763 }
764 }
765
766 //-----------------------------------------------------------------------------
767 void
find_delegate_domain(const conduit::Node & n,conduit::Node & domain,MPI_Comm comm)768 find_delegate_domain(const conduit::Node &n,
769 conduit::Node &domain,
770 MPI_Comm comm)
771 {
772 const index_t par_rank = relay::mpi::rank(comm);
773 const index_t par_size = relay::mpi::size(comm);
774
775 const std::vector<const Node *> domains = ::conduit::blueprint::mesh::domains(n);
776 const index_t local_ind = domains.empty() ? 0 : 1;
777
778 std::vector<int64> local_ind_list(par_size, 0);
779 local_ind_list[par_rank] = local_ind;
780
781 Node local_ind_node, global_ind_node;
782 local_ind_node.set_external(&local_ind_list[0], local_ind_list.size());
783 relay::mpi::sum_all_reduce(local_ind_node, global_ind_node, comm);
784
785 Node temp;
786 index_t min_delegate_rank = -1;
787 for(index_t ri = 0; ri < par_size && min_delegate_rank == -1; ri++)
788 {
789 temp.set_external(DataType(global_ind_node.dtype().id(), 1),
790 global_ind_node.element_ptr(ri));
791
792 const index_t rank_ind = temp.to_index_t();
793 if(rank_ind == 1)
794 {
795 min_delegate_rank = ri;
796 }
797 }
798
799 if(min_delegate_rank != -1)
800 {
801 if(par_rank == min_delegate_rank)
802 {
803 // TODO(JRC): Consider using a more consistent evaluation for the
804 // "first" domain (e.g. the domain with the lowest "state/domain_id"
805 // value).
806 domain.set(*domains[0]);
807 }
808 relay::mpi::broadcast_using_schema(domain, min_delegate_rank, comm);
809 }
810 }
811
812
813 //-----------------------------------------------------------------------------
match_nbr_elems(PolyBndry & pbnd,std::map<index_t,bputils::connectivity::ElemType> & nbr_elems,std::map<index_t,bputils::connectivity::SubelemMap> & allfaces_map,const Node & ref_topo,const Node & ref_win,const Node & nbr_win,index_t nbr_iwidth,index_t nbr_jwidth,index_t nbr_kwidth,index_t ni_lo,index_t nj_lo,index_t nk_lo,index_t ratio_i,index_t ratio_j,index_t ratio_k,bool local)814 void match_nbr_elems(PolyBndry& pbnd,
815 std::map<index_t, bputils::connectivity::ElemType>& nbr_elems,
816 std::map<index_t, bputils::connectivity::SubelemMap>& allfaces_map,
817 const Node& ref_topo,
818 const Node& ref_win,
819 const Node& nbr_win,
820 index_t nbr_iwidth, index_t nbr_jwidth,
821 index_t nbr_kwidth,
822 index_t ni_lo, index_t nj_lo, index_t nk_lo,
823 index_t ratio_i, index_t ratio_j, index_t ratio_k,
824 bool local)
825 {
826 index_t nbr_size_i = nbr_win["dims/i"].to_index_t();
827 index_t nbr_size_j = nbr_win["dims/j"].to_index_t();
828 index_t nbr_size_k = nbr_win["dims/k"].to_index_t();
829
830 index_t ri_lo = ref_topo["elements/origin/i0"].to_index_t();
831 index_t rj_lo = ref_topo["elements/origin/j0"].to_index_t();
832 index_t rk_lo = ref_topo["elements/origin/k0"].to_index_t();
833
834 index_t ref_iwidth = ref_topo["elements/dims/i"].to_index_t();
835 index_t ref_jwidth = ref_topo["elements/dims/j"].to_index_t();
836
837 index_t ref_origin_i = ref_win["origin/i"].to_index_t();
838 index_t ref_origin_j = ref_win["origin/j"].to_index_t();
839 index_t ref_origin_k = ref_win["origin/k"].to_index_t();
840
841 index_t origin_i = nbr_win["origin/i"].to_index_t();
842 index_t origin_j = nbr_win["origin/j"].to_index_t();
843 index_t origin_k = nbr_win["origin/k"].to_index_t();
844
845 index_t side = pbnd.side;
846
847 if (side == 0 || side == 1)
848 {
849 pbnd.m_nbrs_per_face = static_cast<size_t>(ratio_j*ratio_k);
850
851 index_t nside = (side+1)%2; //nbr side counterpart to ref side
852 index_t shift = -(nside%2); //0 on low side, -1 on high side
853 index_t icnst = origin_i - ni_lo + shift;
854 index_t jstart = origin_j;
855 index_t jend = jstart + nbr_size_j - 1;
856 index_t kstart = origin_k;
857 index_t kend = kstart + nbr_size_k - 1;
858
859 index_t rshift = -(side%2);
860 index_t ricnst = ref_origin_i - ri_lo + rshift;
861
862 for (index_t kidx = kstart; kidx < kend; ++kidx)
863 {
864 index_t rkidx = kidx/ratio_k;
865 index_t nk = kidx - nk_lo;
866 index_t rk = rkidx - rk_lo;
867 index_t nkoffset = icnst + nk*nbr_iwidth*nbr_jwidth;
868 index_t rkoffset = ricnst + rk*ref_iwidth*ref_jwidth;
869 for (index_t jidx = jstart; jidx < jend; ++jidx)
870 {
871 index_t rjidx = jidx/ratio_j;
872 index_t nj = jidx - nj_lo;
873 index_t rj = rjidx - rj_lo;
874 index_t noffset = nkoffset + nj*nbr_iwidth;
875 index_t roffset = rkoffset + rj*ref_iwidth;
876
877 pbnd.m_nbr_elems[roffset].push_back(noffset);
878
879 if (local)
880 {
881 auto& nbr_elem = nbr_elems[noffset];
882 index_t face_id =
883 nside == 0 ?
884 nbr_elem[0] :
885 nbr_elem[1];
886 pbnd.m_nbr_faces[roffset][noffset].m_face_id = face_id;
887 pbnd.m_nbr_faces[roffset][noffset].m_fine_subelem =
888 allfaces_map[pbnd.m_nbr_id][face_id];
889 }
890 else
891 {
892 std::vector<index_t> nbr_elem;
893 std::map< index_t, std::vector<index_t> > elem_faces;
894 bputils::connectivity::make_element_3d(nbr_elem,
895 noffset,
896 nbr_iwidth,
897 nbr_jwidth,
898 nbr_kwidth,
899 elem_faces);
900
901 index_t jfoffset = nj*(nbr_iwidth+1);
902 index_t kfoffset = nk*nbr_jwidth*(nbr_iwidth+1);
903 index_t face_id;
904 if (nside == 0)
905 {
906 face_id =
907 icnst + jfoffset + kfoffset;
908 }
909 else
910 {
911 face_id =
912 icnst + jfoffset + kfoffset + 1;
913 }
914 pbnd.m_nbr_faces[roffset][noffset].m_face_id =
915 face_id;
916 pbnd.m_nbr_faces[roffset][noffset].m_fine_subelem =
917 elem_faces[face_id];
918 }
919 }
920 }
921 }
922 else if (side == 2 || side == 3)
923 {
924 pbnd.m_nbrs_per_face = static_cast<size_t>(ratio_i*ratio_k);
925
926 index_t jface_start = (nbr_iwidth+1)*nbr_jwidth*nbr_kwidth;
927 index_t nside = 2 + (side+1)%2; //nbr side counterpart to ref side
928 index_t shift = -(nside%2); //0 on low side, -1 on high side
929 index_t jcnst = origin_j - nj_lo + shift;
930 index_t istart = origin_i;
931 index_t iend = istart + nbr_size_i - 1;
932 index_t kstart = origin_k;
933 index_t kend = kstart + nbr_size_k - 1;
934
935 index_t rshift = -(side%2);
936 index_t rjcnst = ref_origin_j - rj_lo + rshift;
937
938 for (index_t kidx = kstart; kidx < kend; ++kidx)
939 {
940 index_t rkidx = kidx/ratio_k;
941 index_t nk = kidx - nk_lo;
942 index_t rk = rkidx - rk_lo;
943 index_t nkoffset = nk*nbr_iwidth*nbr_jwidth;
944 index_t rkoffset = rk*ref_iwidth*ref_jwidth;
945 for (index_t iidx = istart; iidx < iend; ++iidx)
946 {
947 index_t riidx = iidx/ratio_i;
948 index_t ni = iidx - ni_lo;
949 index_t ri = riidx - ri_lo;
950 index_t noffset = nkoffset + jcnst*nbr_iwidth + ni;
951 index_t roffset = rkoffset + rjcnst*ref_iwidth + ri;
952
953 pbnd.m_nbr_elems[roffset].push_back(noffset);
954
955 if (local)
956 {
957 auto& nbr_elem = nbr_elems[noffset];
958 index_t face_id =
959 nside == 2 ?
960 nbr_elem[2] :
961 nbr_elem[3];
962 pbnd.m_nbr_faces[roffset][noffset].m_face_id = face_id;
963 pbnd.m_nbr_faces[roffset][noffset].m_fine_subelem =
964 allfaces_map[pbnd.m_nbr_id][face_id];
965 }
966 else
967 {
968 std::vector<index_t> nbr_elem;
969 std::map< index_t, std::vector<index_t> > elem_faces;
970 bputils::connectivity::make_element_3d(nbr_elem,
971 noffset,
972 nbr_iwidth,
973 nbr_jwidth,
974 nbr_kwidth,
975 elem_faces);
976
977 index_t ifoffset = ni;
978 index_t kfoffset = nk*(nbr_jwidth+1)*nbr_iwidth;
979 index_t face_id;
980 if (nside == 2)
981 {
982 face_id =
983 jface_start +
984 nbr_iwidth*jcnst + ifoffset + kfoffset;
985 }
986 else
987 {
988 face_id =
989 jface_start +
990 nbr_iwidth*(jcnst+1) + ifoffset + kfoffset;
991 }
992 pbnd.m_nbr_faces[roffset][noffset].m_face_id =
993 face_id;
994 pbnd.m_nbr_faces[roffset][noffset].m_fine_subelem =
995 elem_faces[face_id];
996 }
997 }
998 }
999 }
1000 else if (side == 4 || side == 5)
1001 {
1002 pbnd.m_nbrs_per_face = static_cast<size_t>(ratio_i*ratio_j);
1003
1004 index_t jface_start = (nbr_iwidth+1)*nbr_jwidth*nbr_kwidth;
1005 index_t kface_start = jface_start +
1006 nbr_iwidth*(nbr_jwidth+1)*nbr_kwidth;
1007 index_t nside = 4 + (side+1)%2; //nbr side counterpart to ref side
1008 index_t shift = -(nside%2); //0 on low side, -1 on high side
1009 index_t kcnst = origin_k - nk_lo + shift;
1010 index_t jstart = origin_j;
1011 index_t jend = jstart + nbr_size_j - 1;
1012 index_t istart = origin_i;
1013 index_t iend = istart + nbr_size_i - 1;
1014
1015 index_t rshift = -(side%2);
1016 index_t rkcnst = ref_origin_k - rk_lo + rshift;
1017
1018 for (index_t jidx = jstart; jidx < jend; ++jidx)
1019 {
1020 index_t rjidx = jidx/ratio_j;
1021 index_t nj = jidx - nj_lo;
1022 index_t rj = rjidx - rj_lo;
1023 index_t njoffset = kcnst*nbr_iwidth*nbr_jwidth + nj*nbr_iwidth;
1024 index_t rjoffset = rkcnst*ref_iwidth*ref_jwidth + rj*ref_iwidth;
1025 for (index_t iidx = istart; iidx < iend; ++iidx)
1026 {
1027 index_t riidx = iidx/ratio_i;
1028 index_t ni = iidx - ni_lo;
1029 index_t ri = riidx - ri_lo;
1030 index_t noffset = njoffset + ni;
1031 index_t roffset = rjoffset + ri;
1032
1033 pbnd.m_nbr_elems[roffset].push_back(noffset);
1034
1035 if (local)
1036 {
1037 auto& nbr_elem = nbr_elems[noffset];
1038 index_t face_id =
1039 nside == 4 ?
1040 nbr_elem[4] :
1041 nbr_elem[5];
1042 pbnd.m_nbr_faces[roffset][noffset].m_face_id = face_id;
1043 pbnd.m_nbr_faces[roffset][noffset].m_fine_subelem =
1044 allfaces_map[pbnd.m_nbr_id][face_id];
1045 }
1046 else
1047 {
1048 std::vector<index_t> nbr_elem;
1049 std::map< index_t, std::vector<index_t> > elem_faces;
1050 bputils::connectivity::make_element_3d(nbr_elem,
1051 noffset,
1052 nbr_iwidth,
1053 nbr_jwidth,
1054 nbr_kwidth,
1055 elem_faces);
1056 index_t ifoffset = ni;
1057 index_t jfoffset = nj*nbr_iwidth;
1058 index_t face_id;
1059 if (nside == 4)
1060 {
1061 face_id =
1062 kface_start +
1063 nbr_iwidth*nbr_jwidth*kcnst + ifoffset + jfoffset;
1064 }
1065 else
1066 {
1067 face_id =
1068 kface_start +
1069 nbr_iwidth*nbr_jwidth*(kcnst+1) + ifoffset + jfoffset;
1070 }
1071 pbnd.m_nbr_faces[roffset][noffset].m_face_id =
1072 face_id;
1073 pbnd.m_nbr_faces[roffset][noffset].m_fine_subelem =
1074 elem_faces[face_id];
1075 }
1076 }
1077 }
1078 }
1079 else
1080 {
1081 //not a side
1082 }
1083
1084 }
1085
1086
1087 //-------------------------------------------------------------------------
fix_duplicated_vertices(std::set<index_t> & check_verts,std::vector<index_t> & elem,bputils::connectivity::SubelemMap & allfaces,const double_array & xarray,const double_array & yarray,const double_array & zarray)1088 void fix_duplicated_vertices(std::set<index_t>& check_verts,
1089 std::vector<index_t>& elem,
1090 bputils::connectivity::SubelemMap& allfaces,
1091 const double_array& xarray,
1092 const double_array& yarray,
1093 const double_array& zarray)
1094 {
1095
1096 double eps = sqrt(std::numeric_limits<double>::epsilon());
1097
1098 std::set<index_t> erase_verts;
1099 for (auto it0 = check_verts.begin(); it0 != check_verts.end(); ++it0)
1100 {
1101 const index_t& v0 = *it0;
1102 double x0 = xarray[v0];
1103 double y0 = yarray[v0];
1104 double z0 = zarray[v0];
1105
1106 auto it1 = it0;
1107 ++it1;
1108
1109 for ( ; it1 != check_verts.end(); ++it1)
1110 {
1111 const index_t& v1 = *it1;
1112 double x1 = xarray[v1];
1113 double y1 = yarray[v1];
1114 double z1 = zarray[v1];
1115 double xdiff = x0-x1;
1116 double ydiff = y0-y1;
1117 double zdiff = z0-z1;
1118 double distsqr = xdiff*xdiff+ydiff*ydiff+zdiff*zdiff;
1119
1120 if (distsqr < eps)
1121 {
1122 for (auto eitr = elem.begin(); eitr != elem.end(); ++eitr)
1123 {
1124 auto& face = allfaces[*eitr];
1125 for (auto fitr = face.begin(); fitr != face.end(); ++fitr)
1126 {
1127 index_t& vert = *fitr;
1128 if (vert == v1)
1129 {
1130 vert = v0;
1131 }
1132 }
1133 }
1134 erase_verts.insert(v1);
1135 }
1136 }
1137 }
1138
1139 for (auto ev = erase_verts.begin(); ev != erase_verts.end(); ++ev)
1140 {
1141 check_verts.erase(*ev);
1142 }
1143 }
1144
1145 //-------------------------------------------------------------------------
to_polyhedral(const Node & n,Node & dest,const std::string & name,MPI_Comm comm)1146 void to_polyhedral(const Node &n,
1147 Node &dest,
1148 const std::string& name,
1149 MPI_Comm comm)
1150 {
1151 dest.reset();
1152
1153 index_t par_rank = relay::mpi::rank(comm);
1154
1155 NodeConstIterator itr = n.children();
1156
1157 std::map<index_t, std::map<index_t, bputils::connectivity::ElemType> > poly_elems_map;
1158 std::map<index_t, std::map<index_t, std::set<index_t> > > new_vert_map;
1159 std::map<index_t, std::map<index_t, std::set<index_t> > > new_face_map;
1160 std::map<index_t, bputils::connectivity::SubelemMap> allfaces_map;
1161
1162 std::map<index_t, std::vector<index_t> > elem_connect;
1163 std::map<index_t, std::vector<index_t> > elem_sizes;
1164 std::map<index_t, std::vector<index_t> > elem_offsets;
1165 std::map<index_t, std::vector<index_t> > subelem_connect;
1166 std::map<index_t, std::vector<index_t> > subelem_sizes;
1167 std::map<index_t, std::vector<index_t> > subelem_offsets;
1168
1169 while(itr.has_next())
1170 {
1171 const Node& chld = itr.next();
1172 std::string domain_name = itr.name();
1173 dest[domain_name]["state"] = chld["state"];
1174 const Node& in_coords = chld["coordsets/coords"];
1175 const Node& in_topo = chld["topologies"][name];
1176
1177 index_t iwidth = in_topo["elements/dims/i"].to_index_t();
1178 index_t jwidth = in_topo["elements/dims/j"].to_index_t();
1179 index_t kwidth = in_topo["elements/dims/k"].to_index_t();
1180
1181 Node& out_coords = dest[domain_name]["coordsets/coords"];
1182
1183 index_t domain_id = chld["state/domain_id"].to_index_t();
1184 Node& out_values = out_coords["values"];
1185 if (in_coords["type"].as_string() == "uniform")
1186 {
1187 blueprint::mesh::coordset::uniform::to_explicit(in_coords, out_coords);
1188 }
1189 else
1190 {
1191 out_coords["type"] = in_coords["type"];
1192 const Node& in_values = in_coords["values"];
1193 out_values = in_values;
1194 }
1195
1196 auto& poly_elems = poly_elems_map[domain_id];
1197 auto& allfaces = allfaces_map[domain_id];
1198
1199 index_t elemsize = iwidth*jwidth*kwidth;
1200
1201 for (index_t elem = 0; elem < elemsize; ++elem)
1202 {
1203 bputils::connectivity::make_element_3d(poly_elems[elem],
1204 elem,
1205 iwidth,
1206 jwidth,
1207 kwidth,
1208 allfaces);
1209 }
1210 }
1211
1212 itr = n.children();
1213 while(itr.has_next())
1214 {
1215 const Node& chld = itr.next();
1216 std::string domain_name = itr.name();
1217
1218 index_t domain_id = chld["state/domain_id"].to_index_t();
1219
1220 const Node& in_topo = chld["topologies"][name];
1221
1222 const Node *in_cset = bputils::find_reference_node(in_topo, "coordset");
1223
1224 const index_t i_lo = in_topo["elements/origin/i0"].to_index_t();
1225 const index_t j_lo = in_topo["elements/origin/j0"].to_index_t();
1226 const index_t k_lo = in_topo["elements/origin/k0"].to_index_t();
1227 const index_t niwidth = in_topo["elements/dims/i"].to_index_t() + 1;
1228 const index_t njwidth = in_topo["elements/dims/j"].to_index_t() + 1;
1229
1230 const Node* in_parent = chld.parent();
1231
1232 std::ostringstream win_oss;
1233 win_oss << "window_" << std::setw(6) << std::setfill('0') << domain_id;
1234 std::string win_name = win_oss.str();
1235
1236 if (chld.has_path("adjsets/adjset/groups"))
1237 {
1238 const Node& in_groups = chld["adjsets/adjset/groups"];
1239 NodeConstIterator grp_itr = in_groups.children();
1240 while(grp_itr.has_next())
1241 {
1242 const Node& group = grp_itr.next();
1243 if (group.has_child("neighbors") && group.has_child("windows"))
1244 {
1245 int64_array neighbors = group["neighbors"].as_int64_array();
1246
1247 index_t nbr_id = neighbors[1];
1248 const Node& in_windows = group["windows"];
1249 std::ostringstream nw_oss;
1250 nw_oss << "window_" << std::setw(6)
1251 << std::setfill('0') << nbr_id;
1252 std::string nbr_win_name = nw_oss.str();
1253
1254 const Node& ref_win = in_windows[win_name];
1255 const Node& nbr_win = in_windows[nbr_win_name];
1256
1257 if (nbr_win["level_id"].to_index_t() < ref_win["level_id"].to_index_t())
1258 {
1259 index_t ref_size_i = ref_win["dims/i"].to_index_t();
1260 index_t ref_size_j = ref_win["dims/j"].to_index_t();
1261 index_t ref_size_k = ref_win["dims/k"].to_index_t();
1262 index_t ref_size = ref_size_i*ref_size_j*ref_size_k;
1263
1264 index_t nbr_size_i = nbr_win["dims/i"].to_index_t();
1265 index_t nbr_size_j = nbr_win["dims/j"].to_index_t();
1266 index_t nbr_size_k = nbr_win["dims/k"].to_index_t();
1267 index_t nbr_size = nbr_size_i*nbr_size_j*nbr_size_k;
1268
1269 if (nbr_size < ref_size)
1270 {
1271 std::ostringstream nbr_oss;
1272 nbr_oss << "domain_" << std::setw(6)
1273 << std::setfill('0') << nbr_id;
1274 std::string nbr_name = nbr_oss.str();
1275
1276 bool is_side = true;
1277 if (nbr_size_i * nbr_size_j == 1 ||
1278 nbr_size_i * nbr_size_k == 1 ||
1279 nbr_size_j * nbr_size_k == 1)
1280 {
1281 is_side = false;
1282 }
1283
1284 if (is_side && !in_parent->has_child(nbr_name))
1285 {
1286 index_t buffer[6];
1287 buffer[0] = in_topo["elements/dims/i"].to_index_t();
1288 buffer[1] = in_topo["elements/dims/j"].to_index_t();
1289 buffer[2] = in_topo["elements/dims/k"].to_index_t();
1290 buffer[3] = in_topo["elements/origin/i0"].to_index_t();
1291 buffer[4] = in_topo["elements/origin/j0"].to_index_t();
1292 buffer[5] = in_topo["elements/origin/k0"].to_index_t();
1293
1294 index_t nbr_rank = group["rank"].to_index_t();
1295 MPI_Send(buffer,
1296 6,
1297 MPI_INT64_T,
1298 nbr_rank,
1299 domain_id,
1300 comm);
1301
1302 std::vector<index_t> vertices;
1303 std::vector<double> xbuffer, ybuffer, zbuffer;
1304 const Node &fcoords = (*in_cset)["values"];
1305
1306 const index_t origin_i = ref_win["origin/i"].to_index_t();
1307 const index_t origin_j = ref_win["origin/j"].to_index_t();
1308 const index_t origin_k = ref_win["origin/k"].to_index_t();
1309
1310 Node temp;
1311
1312 if (ref_size_i == 1)
1313 {
1314 const index_t icnst = origin_i - i_lo;
1315 const index_t jstart = origin_j - j_lo;
1316 const index_t jend = jstart + ref_size_j;
1317 const index_t kstart = origin_k - k_lo;
1318 const index_t kend = kstart + ref_size_k;
1319 for (index_t kidx = kstart; kidx < kend; ++kidx)
1320 {
1321 const index_t koffset = kidx * niwidth * njwidth;
1322 for (index_t jidx = jstart; jidx < jend; ++jidx)
1323 {
1324 const index_t offset = koffset + jidx * niwidth + icnst;
1325 vertices.push_back(offset);
1326 temp.set_external(DataType(fcoords["x"].dtype().id(), 1),
1327 (void*)fcoords["x"].element_ptr(offset));
1328 xbuffer.push_back(temp.to_double());
1329 temp.set_external(DataType(fcoords["y"].dtype().id(), 1),
1330 (void*)fcoords["y"].element_ptr(offset));
1331 ybuffer.push_back(temp.to_double());
1332 temp.set_external(DataType(fcoords["z"].dtype().id(), 1),
1333 (void*)fcoords["z"].element_ptr(offset));
1334 zbuffer.push_back(temp.to_double());
1335 }
1336 }
1337 }
1338 else if (ref_size_j == 1)
1339 {
1340 const index_t istart = origin_i - i_lo;
1341 const index_t iend = istart + ref_size_i;
1342 const index_t jcnst = origin_j - j_lo;
1343 const index_t kstart = origin_k - k_lo;
1344 const index_t kend = kstart + ref_size_k;
1345 for (index_t kidx = kstart; kidx < kend; ++kidx)
1346 {
1347 const index_t koffset = kidx * niwidth * njwidth;
1348 for (index_t iidx = istart; iidx < iend; ++iidx)
1349 {
1350 const index_t offset = koffset + jcnst * niwidth + iidx;
1351 vertices.push_back(offset);
1352 temp.set_external(DataType(fcoords["x"].dtype().id(), 1),
1353 (void*)fcoords["x"].element_ptr(offset));
1354 xbuffer.push_back(temp.to_double());
1355 temp.set_external(DataType(fcoords["y"].dtype().id(), 1),
1356 (void*)fcoords["y"].element_ptr(offset));
1357 ybuffer.push_back(temp.to_double());
1358 temp.set_external(DataType(fcoords["z"].dtype().id(), 1),
1359 (void*)fcoords["z"].element_ptr(offset));
1360 zbuffer.push_back(temp.to_double());
1361 }
1362 }
1363 }
1364 else
1365 {
1366 const index_t istart = origin_i - i_lo;
1367 const index_t iend = istart + ref_size_i;
1368 const index_t jstart = origin_j - j_lo;
1369 const index_t jend = jstart + ref_size_j;
1370 const index_t kcnst = origin_k - k_lo;
1371 const index_t koffset = kcnst * niwidth * njwidth;
1372 for (index_t jidx = jstart; jidx < jend; ++jidx)
1373 {
1374 const index_t joffset = jidx * niwidth;
1375 for (index_t iidx = istart; iidx < iend; ++iidx)
1376 {
1377 const index_t offset = koffset + joffset + iidx;
1378 vertices.push_back(offset);
1379 temp.set_external(DataType(fcoords["x"].dtype().id(), 1),
1380 (void*)fcoords["x"].element_ptr(offset));
1381 xbuffer.push_back(temp.to_double());
1382 temp.set_external(DataType(fcoords["y"].dtype().id(), 1),
1383 (void*)fcoords["y"].element_ptr(offset));
1384 ybuffer.push_back(temp.to_double());
1385 temp.set_external(DataType(fcoords["z"].dtype().id(), 1),
1386 (void*)fcoords["z"].element_ptr(offset));
1387 zbuffer.push_back(temp.to_double());
1388 }
1389 }
1390 }
1391
1392
1393 MPI_Send(&vertices[0],
1394 vertices.size(),
1395 MPI_INT64_T,
1396 nbr_rank,
1397 domain_id,
1398 comm);
1399 MPI_Send(&xbuffer[0],
1400 xbuffer.size(),
1401 MPI_DOUBLE,
1402 nbr_rank,
1403 domain_id,
1404 comm);
1405 MPI_Send(&ybuffer[0],
1406 ybuffer.size(),
1407 MPI_DOUBLE,
1408 nbr_rank,
1409 domain_id,
1410 comm);
1411 MPI_Send(&zbuffer[0],
1412 zbuffer.size(),
1413 MPI_DOUBLE,
1414 nbr_rank,
1415 domain_id,
1416 comm);
1417
1418
1419 }
1420 }
1421 }
1422 }
1423 }
1424 }
1425 }
1426
1427 //Outer map: domain_id, inner map: nbr_id to PolyBndry
1428 std::map<index_t, std::map<index_t, PolyBndry> > poly_bndry_map;
1429
1430 std::map<index_t, std::vector<index_t> > nbr_ratio;
1431
1432 itr = n.children();
1433 while(itr.has_next())
1434 {
1435 const Node& chld = itr.next();
1436 std::string domain_name = itr.name();
1437
1438 index_t domain_id = chld["state/domain_id"].to_index_t();
1439
1440 auto& poly_elems = poly_elems_map[domain_id];
1441
1442 const Node& in_topo = chld["topologies"][name];
1443
1444 index_t iwidth = in_topo["elements/dims/i"].to_index_t();
1445 index_t jwidth = in_topo["elements/dims/j"].to_index_t();
1446
1447 index_t i_lo = in_topo["elements/origin/i0"].to_index_t();
1448 index_t j_lo = in_topo["elements/origin/j0"].to_index_t();
1449 index_t k_lo = in_topo["elements/origin/k0"].to_index_t();
1450
1451 const Node* in_parent = chld.parent();
1452
1453 std::ostringstream win_oss;
1454 win_oss << "window_" << std::setw(6) << std::setfill('0') << domain_id;
1455 std::string win_name = win_oss.str();
1456
1457 if (chld.has_path("adjsets/adjset/groups"))
1458 {
1459 const Node& in_groups = chld["adjsets/adjset/groups"];
1460 NodeConstIterator grp_itr = in_groups.children();
1461 while(grp_itr.has_next())
1462 {
1463 const Node& group = grp_itr.next();
1464
1465 if (group.has_child("neighbors") && group.has_child("windows"))
1466 {
1467 int64_array neighbors = group["neighbors"].as_int64_array();
1468
1469 index_t nbr_id = neighbors[1];
1470 const Node& in_windows = group["windows"];
1471 std::ostringstream nw_oss;
1472 nw_oss << "window_" << std::setw(6)
1473 << std::setfill('0') << nbr_id;
1474 std::string nbr_win_name = nw_oss.str();
1475
1476 const Node& ref_win = in_windows[win_name];
1477 const Node& nbr_win = in_windows[nbr_win_name];
1478 if (nbr_win["level_id"].to_index_t() > ref_win["level_id"].to_index_t())
1479 {
1480 index_t ratio_i = nbr_win["ratio/i"].to_index_t();
1481 index_t ratio_j = nbr_win["ratio/j"].to_index_t();
1482 index_t ratio_k = nbr_win["ratio/k"].to_index_t();
1483
1484 nbr_ratio[domain_id] = {ratio_i, ratio_j, ratio_k};
1485
1486 index_t ref_size_i = ref_win["dims/i"].to_index_t();
1487 index_t ref_size_j = ref_win["dims/j"].to_index_t();
1488 index_t ref_size_k = ref_win["dims/k"].to_index_t();
1489 index_t ref_size = ref_size_i*ref_size_j*ref_size_k;
1490
1491 index_t nbr_size_i = nbr_win["dims/i"].to_index_t();
1492 index_t nbr_size_j = nbr_win["dims/j"].to_index_t();
1493 index_t nbr_size_k = nbr_win["dims/k"].to_index_t();
1494 index_t nbr_size = nbr_size_i*nbr_size_j*nbr_size_k;
1495
1496 if (nbr_size > ref_size)
1497 {
1498 index_t origin_i = ref_win["origin/i"].to_index_t();
1499 index_t origin_j = ref_win["origin/j"].to_index_t();
1500 index_t origin_k = ref_win["origin/k"].to_index_t();
1501
1502 PolyBndry& pbnd = poly_bndry_map[domain_id][nbr_id];
1503
1504 if (ref_size_i == 1 && ref_size_j != 1 &&
1505 ref_size_k != 1)
1506 {
1507 index_t shift = 0;
1508 if (origin_i == i_lo)
1509 {
1510 pbnd.side = 0;
1511 }
1512 else
1513 {
1514 pbnd.side = 1;
1515 shift = -1;
1516 }
1517 index_t icnst = origin_i - i_lo + shift;
1518 index_t jstart = origin_j - j_lo;
1519 index_t jend = jstart + ref_size_j - 1;
1520 index_t kstart = origin_k - k_lo;
1521 index_t kend = kstart + ref_size_k - 1;
1522 for (index_t kidx = kstart; kidx < kend; ++kidx)
1523 {
1524 index_t koffset = icnst + kidx*iwidth*jwidth;
1525 for (index_t jidx = jstart; jidx < jend; ++jidx)
1526 {
1527 index_t offset = koffset + jidx*iwidth;
1528 pbnd.m_elems.push_back(offset);
1529 pbnd.m_bface[offset] =
1530 pbnd.side == 0 ?
1531 poly_elems[offset][0] :
1532 poly_elems[offset][1];
1533 }
1534 }
1535 }
1536 else if (ref_size_j == 1 && ref_size_i != 1 &&
1537 ref_size_k != 1)
1538 {
1539 index_t shift = 0;
1540 if (origin_j == j_lo)
1541 {
1542 pbnd.side = 2;
1543 }
1544 else
1545 {
1546 pbnd.side = 3;
1547 shift = -1;
1548 }
1549 index_t jcnst = origin_j - j_lo + shift;
1550 index_t istart = origin_i - i_lo;
1551 index_t iend = istart + ref_size_i - 1;
1552 index_t kstart = origin_k - k_lo;
1553 index_t kend = kstart + ref_size_k - 1;
1554 for (index_t kidx = kstart; kidx < kend; ++kidx)
1555 {
1556 index_t koffset = jcnst*iwidth + kidx*iwidth*jwidth;
1557 for (index_t iidx = istart; iidx < iend; ++iidx)
1558 {
1559 index_t offset = koffset + iidx;
1560 pbnd.m_elems.push_back(offset);
1561 pbnd.m_bface[offset] =
1562 pbnd.side == 2 ?
1563 poly_elems[offset][2] :
1564 poly_elems[offset][3];
1565 }
1566 }
1567 }
1568 else if (ref_size_k == 1 && ref_size_i != 1 &&
1569 ref_size_j != 1)
1570 {
1571 index_t shift = 0;
1572 if (origin_k == k_lo)
1573 {
1574 pbnd.side = 4;
1575 }
1576 else
1577 {
1578 pbnd.side = 5;
1579 shift = -1;
1580 }
1581 index_t kcnst = origin_k - k_lo + shift;
1582 index_t istart = origin_i - i_lo;
1583 index_t iend = istart + ref_size_i - 1;
1584 index_t jstart = origin_j - j_lo;
1585 index_t jend = jstart + ref_size_j - 1;
1586 for (index_t jidx = jstart; jidx < jend; ++jidx)
1587 {
1588 index_t joffset = jidx*iwidth + kcnst*iwidth*jwidth;
1589 for (index_t iidx = istart; iidx < iend; ++iidx)
1590 {
1591 index_t offset = joffset + iidx;
1592 pbnd.m_elems.push_back(offset);
1593 pbnd.m_bface[offset] =
1594 pbnd.side == 4 ?
1595 poly_elems[offset][4] :
1596 poly_elems[offset][5];
1597 }
1598 }
1599 }
1600 else
1601 {
1602 pbnd.side = -1;
1603 }
1604
1605
1606 std::ostringstream nbr_oss;
1607 nbr_oss << "domain_" << std::setw(6)
1608 << std::setfill('0') << nbr_id;
1609 std::string nbr_name = nbr_oss.str();
1610
1611 if (in_parent->has_child(nbr_name))
1612 {
1613 const Node& nbr_dom =
1614 (*in_parent)[nbr_name];
1615
1616 const Node& ntopo =
1617 nbr_dom["topologies"][name];
1618
1619 index_t ni_lo = ntopo["elements/origin/i0"].to_index_t();
1620 index_t nj_lo = ntopo["elements/origin/j0"].to_index_t();
1621 index_t nk_lo = ntopo["elements/origin/k0"].to_index_t();
1622
1623 index_t nbr_iwidth =
1624 ntopo["elements/dims/i"].to_index_t();
1625 index_t nbr_jwidth =
1626 ntopo["elements/dims/j"].to_index_t();
1627 index_t nbr_kwidth =
1628 ntopo["elements/dims/k"].to_index_t();
1629
1630 auto& nbr_elems = poly_elems_map[nbr_id];
1631
1632 if (pbnd.side >= 0)
1633 {
1634 pbnd.m_nbr_rank = par_rank;
1635 pbnd.m_nbr_id = nbr_id;
1636 match_nbr_elems(pbnd, nbr_elems, allfaces_map, in_topo,
1637 ref_win, nbr_win,
1638 nbr_iwidth, nbr_jwidth,
1639 nbr_kwidth,
1640 ni_lo, nj_lo, nk_lo,
1641 ratio_i,ratio_j,ratio_k,
1642 true);
1643 }
1644 }
1645 }
1646 }
1647 }
1648 }
1649 }
1650 }
1651
1652 std::map<index_t, std::map<index_t, std::map<index_t,index_t> > > dom_to_nbr_to_buffidx;
1653 std::map<index_t, std::map<index_t, std::vector<double> > > dom_to_nbr_to_xbuffer;
1654 std::map<index_t, std::map<index_t, std::vector<double> > > dom_to_nbr_to_ybuffer;
1655 std::map<index_t, std::map<index_t, std::vector<double> > > dom_to_nbr_to_zbuffer;
1656
1657 itr = n.children();
1658 while(itr.has_next())
1659 {
1660 const Node& chld = itr.next();
1661 std::string domain_name = itr.name();
1662
1663 index_t domain_id = chld["state/domain_id"].to_index_t();
1664 const Node& in_topo = chld["topologies"][name];
1665
1666 std::ostringstream win_oss;
1667 win_oss << "window_" << std::setw(6) << std::setfill('0') << domain_id;
1668 std::string win_name = win_oss.str();
1669
1670 const Node* in_parent = chld.parent();
1671
1672 auto& nbr_to_buffidx = dom_to_nbr_to_buffidx[domain_id];
1673 auto& nbr_to_xbuffer = dom_to_nbr_to_xbuffer[domain_id];
1674 auto& nbr_to_ybuffer = dom_to_nbr_to_ybuffer[domain_id];
1675 auto& nbr_to_zbuffer = dom_to_nbr_to_zbuffer[domain_id];
1676
1677 if (chld.has_path("adjsets/adjset/groups"))
1678 {
1679 const Node& in_groups = chld["adjsets/adjset/groups"];
1680 NodeConstIterator grp_itr = in_groups.children();
1681 while(grp_itr.has_next())
1682 {
1683 const Node& group = grp_itr.next();
1684
1685 if (group.has_child("neighbors") && group.has_child("windows"))
1686 {
1687 int64_array neighbors = group["neighbors"].as_int64_array();
1688
1689 index_t nbr_id = neighbors[1];
1690 const Node& in_windows = group["windows"];
1691 std::ostringstream nw_oss;
1692 nw_oss << "window_" << std::setw(6)
1693 << std::setfill('0') << nbr_id;
1694 std::string nbr_win_name = nw_oss.str();
1695
1696 const Node& ref_win = in_windows[win_name];
1697 const Node& nbr_win = in_windows[nbr_win_name];
1698 if (nbr_win["level_id"].to_index_t() > ref_win["level_id"].to_index_t())
1699 {
1700 index_t ratio_i = nbr_win["ratio/i"].to_index_t();
1701 index_t ratio_j = nbr_win["ratio/j"].to_index_t();
1702 index_t ratio_k = nbr_win["ratio/k"].to_index_t();
1703
1704 index_t ref_size_i = ref_win["dims/i"].to_index_t();
1705 index_t ref_size_j = ref_win["dims/j"].to_index_t();
1706 index_t ref_size_k = ref_win["dims/k"].to_index_t();
1707 index_t ref_size = ref_size_i*ref_size_j*ref_size_k;
1708
1709 index_t nbr_size_i = nbr_win["dims/i"].to_index_t();
1710 index_t nbr_size_j = nbr_win["dims/j"].to_index_t();
1711 index_t nbr_size_k = nbr_win["dims/k"].to_index_t();
1712 index_t nbr_size = nbr_size_i*nbr_size_j*nbr_size_k;
1713
1714 if (nbr_size > ref_size)
1715 {
1716 std::vector<index_t> vertices;
1717 auto& buffidx = nbr_to_buffidx[nbr_id];
1718 auto& xbuffer = nbr_to_xbuffer[nbr_id];
1719 auto& ybuffer = nbr_to_ybuffer[nbr_id];
1720 auto& zbuffer = nbr_to_zbuffer[nbr_id];
1721
1722 std::ostringstream nbr_oss;
1723 nbr_oss << "domain_" << std::setw(6)
1724 << std::setfill('0') << nbr_id;
1725 std::string nbr_name = nbr_oss.str();
1726
1727 if (!in_parent->has_child(nbr_name))
1728 {
1729 PolyBndry& pbnd = poly_bndry_map[domain_id][nbr_id];
1730
1731 if (pbnd.side >= 0)
1732 {
1733 index_t nbr_rank = group["rank"].to_index_t();
1734 index_t buffer[6];
1735
1736 MPI_Recv(buffer,
1737 6,
1738 MPI_INT64_T,
1739 nbr_rank,
1740 nbr_id,
1741 comm,
1742 MPI_STATUS_IGNORE);
1743
1744 index_t nbr_iwidth = buffer[0];
1745 index_t nbr_jwidth = buffer[1];
1746 index_t nbr_kwidth = buffer[2];
1747 index_t ni_lo = buffer[3];
1748 index_t nj_lo = buffer[4];
1749 index_t nk_lo = buffer[5];
1750
1751 auto& nbr_elems = poly_elems_map[nbr_id];
1752
1753 pbnd.m_nbr_rank = nbr_rank;
1754 pbnd.m_nbr_id = nbr_id;
1755 match_nbr_elems(pbnd, nbr_elems, allfaces_map, in_topo,
1756 ref_win, nbr_win,
1757 nbr_iwidth, nbr_jwidth,
1758 nbr_kwidth,
1759 ni_lo, nj_lo, nk_lo,
1760 ratio_i,ratio_j,ratio_k,
1761 false);
1762 if (nbr_size_i == 1)
1763 {
1764 vertices.resize(nbr_size_j*nbr_size_k);
1765 xbuffer.resize(nbr_size_j*nbr_size_k);
1766 ybuffer.resize(nbr_size_j*nbr_size_k);
1767 zbuffer.resize(nbr_size_j*nbr_size_k);
1768 }
1769 else if (nbr_size_j == 1)
1770 {
1771 vertices.resize(nbr_size_i*nbr_size_k);
1772 xbuffer.resize(nbr_size_i*nbr_size_k);
1773 ybuffer.resize(nbr_size_i*nbr_size_k);
1774 zbuffer.resize(nbr_size_i*nbr_size_k);
1775 }
1776 else
1777 {
1778 vertices.resize(nbr_size_i*nbr_size_j);
1779 xbuffer.resize(nbr_size_i*nbr_size_j);
1780 ybuffer.resize(nbr_size_i*nbr_size_j);
1781 zbuffer.resize(nbr_size_i*nbr_size_j);
1782 }
1783
1784 MPI_Recv(&vertices[0],
1785 vertices.size(),
1786 MPI_INT64_T,
1787 nbr_rank,
1788 nbr_id,
1789 comm,
1790 MPI_STATUS_IGNORE);
1791 MPI_Recv(&xbuffer[0],
1792 xbuffer.size(),
1793 MPI_DOUBLE,
1794 nbr_rank,
1795 nbr_id,
1796 comm,
1797 MPI_STATUS_IGNORE);
1798 MPI_Recv(&ybuffer[0],
1799 ybuffer.size(),
1800 MPI_DOUBLE,
1801 nbr_rank,
1802 nbr_id,
1803 comm,
1804 MPI_STATUS_IGNORE);
1805 MPI_Recv(&zbuffer[0],
1806 zbuffer.size(),
1807 MPI_DOUBLE,
1808 nbr_rank,
1809 nbr_id,
1810 comm,
1811 MPI_STATUS_IGNORE);
1812 index_t v = 0;
1813 for (auto vitr = vertices.begin();
1814 vitr != vertices.end(); ++vitr)
1815 {
1816 buffidx[*vitr] = v;
1817 ++v;
1818 }
1819 }
1820 }
1821 else
1822 {
1823 const Node& nbr_dom =
1824 (*in_parent)[nbr_name];
1825 const Node& nbr_coords =
1826 nbr_dom["coordsets/coords"];
1827
1828 const Node& ntopo =
1829 nbr_dom["topologies"][name];
1830 index_t ni_lo =
1831 ntopo["elements/origin/i0"].to_index_t();
1832 index_t nj_lo =
1833 ntopo["elements/origin/j0"].to_index_t();
1834 index_t nk_lo =
1835 ntopo["elements/origin/k0"].to_index_t();
1836 index_t nbr_iwidth =
1837 ntopo["elements/dims/i"].to_index_t() + 1;
1838 index_t nbr_jwidth =
1839 ntopo["elements/dims/j"].to_index_t() + 1;
1840
1841 const Node& fcoords =
1842 nbr_coords["values"];
1843 const double_array& xarray =
1844 fcoords["x"].as_double_array();
1845 const double_array& yarray =
1846 fcoords["y"].as_double_array();
1847 const double_array& zarray =
1848 fcoords["z"].as_double_array();
1849
1850 index_t origin_i = nbr_win["origin/i"].to_index_t();
1851 index_t origin_j = nbr_win["origin/j"].to_index_t();
1852 index_t origin_k = nbr_win["origin/k"].to_index_t();
1853
1854 index_t istart = origin_i - ni_lo;
1855 index_t jstart = origin_j - nj_lo;
1856 index_t kstart = origin_k - nk_lo;
1857 index_t iend = istart + nbr_size_i;
1858 index_t jend = jstart + nbr_size_j;
1859 index_t kend = kstart + nbr_size_k;
1860
1861 for (index_t kidx = kstart; kidx < kend; ++kidx)
1862 {
1863 index_t koffset = kidx*nbr_iwidth*nbr_jwidth;
1864 for (index_t jidx = jstart; jidx < jend; ++jidx)
1865 {
1866 index_t joffset = jidx*nbr_iwidth;
1867 for (index_t iidx = istart; iidx < iend; ++iidx)
1868 {
1869 index_t offset = koffset+joffset+iidx;
1870 index_t new_idx = xbuffer.size();
1871 buffidx[offset] = new_idx;
1872 xbuffer.push_back(xarray[offset]);
1873 ybuffer.push_back(yarray[offset]);
1874 zbuffer.push_back(zarray[offset]);
1875 }
1876 }
1877 }
1878 }
1879 }
1880 }
1881 }
1882 }
1883 }
1884 }
1885
1886 std::map<index_t, std::map<index_t, std::map<index_t,index_t> > > dom_to_nbr_to_sharedmap;
1887
1888 itr = n.children();
1889 while(itr.has_next())
1890 {
1891 const Node& chld = itr.next();
1892 std::string domain_name = itr.name();
1893
1894 index_t domain_id = chld["state/domain_id"].to_index_t();
1895
1896 const Node& dom_coords = chld["coordsets/coords"];
1897
1898 const double_array& ref_xarray =
1899 dom_coords["values/x"].as_double_array();
1900 const double_array& ref_yarray =
1901 dom_coords["values/y"].as_double_array();
1902 const double_array& ref_zarray =
1903 dom_coords["values/z"].as_double_array();
1904
1905 Node& out_coords = dest[domain_name]["coordsets/coords"];
1906
1907 std::vector<double> out_xvec;
1908 std::vector<double> out_yvec;
1909 std::vector<double> out_zvec;
1910 if (dom_coords["type"].as_string() == "uniform")
1911 {
1912 Node tmp_coords;
1913 blueprint::mesh::coordset::uniform::to_explicit(dom_coords, tmp_coords);
1914 const Node& tmp_x = tmp_coords["values/x"];
1915 const Node& tmp_y = tmp_coords["values/y"];
1916 const Node& tmp_z = tmp_coords["values/z"];
1917 const double* tmp_x_ptr = static_cast< const double*>(tmp_x.data_ptr());
1918 const double* tmp_y_ptr = static_cast< const double*>(tmp_y.data_ptr());
1919 const double* tmp_z_ptr = static_cast< const double*>(tmp_z.data_ptr());
1920 index_t num_coords = tmp_x.dtype().number_of_elements();
1921 out_xvec.insert(out_xvec.end(), tmp_x_ptr, tmp_x_ptr + num_coords);
1922 out_yvec.insert(out_yvec.end(), tmp_y_ptr, tmp_y_ptr + num_coords);
1923 out_zvec.insert(out_zvec.end(), tmp_z_ptr, tmp_z_ptr + num_coords);
1924 out_coords["type"] = "explicit";
1925 }
1926 else
1927 {
1928 out_coords["type"] = dom_coords["type"];
1929 const double* tmp_x_ptr = static_cast< const double*>(ref_xarray.data_ptr());
1930 const double* tmp_y_ptr = static_cast< const double*>(ref_yarray.data_ptr());
1931 const double* tmp_z_ptr = static_cast< const double*>(ref_zarray.data_ptr());
1932
1933 index_t num_coords = ref_xarray.dtype().number_of_elements();
1934 out_xvec.insert(out_xvec.end(), tmp_x_ptr, tmp_x_ptr + num_coords);
1935 out_yvec.insert(out_yvec.end(), tmp_y_ptr, tmp_y_ptr + num_coords);
1936 out_zvec.insert(out_zvec.end(), tmp_z_ptr, tmp_z_ptr + num_coords);
1937 }
1938
1939
1940 auto& poly_elems = poly_elems_map[domain_id];
1941 auto& new_vertices = new_vert_map[domain_id];
1942 auto& new_subelems = new_face_map[domain_id];
1943 auto& allfaces = allfaces_map[domain_id];
1944
1945 auto& nbr_to_sharedmap = dom_to_nbr_to_sharedmap[domain_id];
1946 auto& nbr_to_buffidx = dom_to_nbr_to_buffidx[domain_id];
1947 auto& nbr_to_xbuffer = dom_to_nbr_to_xbuffer[domain_id];
1948 auto& nbr_to_ybuffer = dom_to_nbr_to_ybuffer[domain_id];
1949 auto& nbr_to_zbuffer = dom_to_nbr_to_zbuffer[domain_id];
1950
1951 std::set<index_t> elems_on_bdry;
1952 std::set<index_t> multi_nbrs;
1953
1954 if (poly_bndry_map.find(domain_id) != poly_bndry_map.end())
1955 {
1956 std::map<index_t, std::map<index_t, size_t> > elem_to_new_faces;
1957
1958 //std::map<index_t, PolyBndry> bndries
1959 auto& bndries = poly_bndry_map[domain_id];
1960 for (auto bitr = bndries.begin(); bitr != bndries.end(); ++bitr)
1961 {
1962 //One PolyBndry for each fine neighbor domain
1963 PolyBndry& pbnd = bitr->second;
1964 index_t nbr_id = pbnd.m_nbr_id; //domain id of nbr
1965
1966 auto& sharedmap = nbr_to_sharedmap[nbr_id];
1967 auto& buffidx = nbr_to_buffidx[nbr_id];
1968 auto& xbuffer = nbr_to_xbuffer[nbr_id];
1969 auto& ybuffer = nbr_to_ybuffer[nbr_id];
1970 auto& zbuffer = nbr_to_zbuffer[nbr_id];
1971
1972 index_t orig_num_vertices = out_xvec.size();
1973
1974 std::map<index_t, std::map<index_t, std::vector<index_t> > >fv_map;
1975
1976 auto& nbr_elems = pbnd.m_nbr_elems;
1977 for (auto eitr = nbr_elems.begin(); eitr != nbr_elems.end(); ++eitr)
1978 {
1979 //offset for a single domain elem
1980 index_t ref_offset = eitr->first;
1981
1982 auto& num_new_faces = elem_to_new_faces[ref_offset];
1983 if (num_new_faces.find(pbnd.side) == num_new_faces.end())
1984 {
1985 num_new_faces[pbnd.side] = 0;
1986 }
1987
1988 //Holds neighbor elem offsets
1989 std::vector<index_t>& nbrs = eitr->second;
1990 size_t num_nbrs = nbrs.size();
1991 for (size_t n = 0; n < num_nbrs; ++n)
1992 {
1993 index_t nbr_elem = nbrs[n];
1994
1995 //nbr_face is subelem offset for the face of
1996 //nbr_elem touching the boundary
1997 index_t nbr_face = pbnd.m_nbr_faces[ref_offset][nbr_elem].m_face_id;
1998 //face_subelem holds the vertices for the nbr_face
1999 auto& face_subelem = pbnd.m_nbr_faces[ref_offset][nbr_elem].m_fine_subelem;
2000
2001 //subelem at ref/nbr interface
2002 fv_map[ref_offset][nbr_face] = face_subelem;
2003 }
2004 }
2005
2006 std::map<index_t, std::set<index_t> > new_nbr_vertices;
2007
2008 //fv_map maps a coarse subelem to a vector of its
2009 //fine subelems;
2010 for (auto fitr = fv_map.begin(); fitr != fv_map.end(); ++fitr)
2011 {
2012 index_t ref_offset = fitr->first;
2013 bputils::connectivity::ElemType& ref_elem = poly_elems[ref_offset];
2014 index_t ref_face = ref_elem[pbnd.side];
2015 std::vector<index_t>& ref_subelem = allfaces[ref_face];
2016
2017 auto& new_nbr_verts = new_nbr_vertices[ref_offset];
2018 std::map<index_t, index_t> face_verts;
2019
2020 auto nbrs = fitr->second;
2021 std::set<index_t> sh_verts;
2022 std::set<index_t> others;
2023
2024 // This loop causes sh_verts to be filled with the
2025 // vertices that only exist once in the neighbor
2026 // subelems. Those vertices are the ones shared by
2027 // the reference subelem.
2028 for (auto nitr = nbrs.begin(); nitr != nbrs.end(); ++nitr)
2029 {
2030 std::vector<index_t>& nbr_subelem = nitr->second;
2031 for (auto nv = nbr_subelem.begin(); nv != nbr_subelem.end(); ++nv)
2032 {
2033 index_t nvert = *nv;
2034 if (sh_verts.find(nvert) == sh_verts.end())
2035 {
2036 if (others.find(nvert) == others.end())
2037 {
2038 sh_verts.insert(nvert);
2039 }
2040 else
2041 {
2042 ++face_verts[nvert];
2043 }
2044 }
2045 else
2046 {
2047 sh_verts.erase(nvert);
2048 others.insert(nvert);
2049 face_verts[nvert] = 1;
2050 }
2051 }
2052 }
2053 for (auto fv = face_verts.begin(); fv != face_verts.end();
2054 ++fv)
2055 {
2056 if (fv->second == 1)
2057 {
2058 new_nbr_verts.insert(fv->first);
2059 }
2060 }
2061
2062 double eps = sqrt(std::numeric_limits<double>::epsilon());
2063
2064 // Determine which reference vertex is shared with
2065 // each neighbor vertex.
2066 for (auto rv = ref_subelem.begin(); rv != ref_subelem.end(); ++rv)
2067 {
2068 double min_dist_sqr = std::numeric_limits<double>::max();
2069 index_t shared_vert = -1;
2070 for (auto sv = sh_verts.begin(); sv != sh_verts.end();
2071 ++sv)
2072 {
2073 index_t nvert = *sv;
2074
2075 index_t buf_offset = buffidx[nvert];
2076 double nbr_x = xbuffer[buf_offset];
2077 double nbr_y = ybuffer[buf_offset];
2078 double nbr_z = zbuffer[buf_offset];
2079
2080 double xdiff = nbr_x-ref_xarray[*rv];
2081 double ydiff = nbr_y-ref_yarray[*rv];
2082 double zdiff = nbr_z-ref_zarray[*rv];
2083
2084 double dist_sqr = xdiff*xdiff + ydiff*ydiff +
2085 zdiff*zdiff;
2086
2087 if (dist_sqr < min_dist_sqr)
2088 {
2089 min_dist_sqr = dist_sqr;
2090 shared_vert = nvert;
2091 }
2092 }
2093
2094 if (min_dist_sqr < eps)
2095 {
2096 sharedmap[shared_vert] = *rv;
2097 sh_verts.erase(shared_vert);
2098 }
2099 }
2100 }
2101
2102 index_t num_vertices = out_xvec.size();
2103
2104 //Replace the coarse subelem with new fine subelems
2105 for (auto eitr = nbr_elems.begin(); eitr != nbr_elems.end(); ++eitr)
2106 {
2107 index_t ref_offset = eitr->first;
2108 auto& ref_elem = poly_elems[ref_offset];
2109 auto& new_nbr_verts = new_nbr_vertices[ref_offset];
2110 auto& new_verts = new_vertices[ref_offset];
2111
2112 index_t ref_face = ref_elem[pbnd.side];
2113
2114 auto& new_subs = new_subelems[ref_offset];
2115
2116 auto& nbr_subelems = fv_map[ref_offset];
2117
2118 auto& num_added_faces =
2119 elem_to_new_faces[ref_offset][pbnd.side];
2120
2121 if (elems_on_bdry.find(ref_offset) == elems_on_bdry.end())
2122 {
2123 elems_on_bdry.insert(ref_offset);
2124 }
2125 else
2126 {
2127 multi_nbrs.insert(ref_offset);
2128 }
2129
2130 // Set to true if the current neighbor domain
2131 // has faces that cover only part of the current coarse
2132 // face. Another neighbor domain will provide faces
2133 // to cover the remainder, and we need to keep track
2134 // of all new vertices added.
2135 bool partial_nbr = false;
2136 if (nbr_subelems.size() != pbnd.m_nbrs_per_face)
2137 {
2138 partial_nbr = true;
2139 }
2140
2141 index_t last_added_face = allfaces.rbegin()->first;
2142 for (auto nb = nbr_subelems.begin();
2143 nb != nbr_subelems.end(); ++nb)
2144 {
2145 std::vector<index_t> new_face;
2146 auto& n_subelem = nb->second;
2147
2148 for (unsigned int v = 0; v < n_subelem.size(); ++v)
2149 {
2150 index_t n_vtx = n_subelem[v];
2151 if (sharedmap.find(n_vtx) == sharedmap.end())
2152 {
2153 new_face.push_back(num_vertices);
2154 sharedmap[n_vtx] = num_vertices;
2155
2156 ++num_vertices;
2157
2158 index_t buf_offset = buffidx[n_vtx];
2159 out_xvec.push_back(xbuffer[buf_offset]);
2160 out_yvec.push_back(ybuffer[buf_offset]);
2161 out_zvec.push_back(zbuffer[buf_offset]);
2162 if (partial_nbr)
2163 {
2164 new_verts.insert(new_face.back());
2165 }
2166 }
2167 else
2168 {
2169 new_face.push_back(sharedmap[n_vtx]);
2170 if (partial_nbr &&
2171 new_face.back() > orig_num_vertices)
2172 {
2173 new_verts.insert(new_face.back());
2174 }
2175 }
2176 if (!partial_nbr &&
2177 new_nbr_verts.find(n_vtx) !=
2178 new_nbr_verts.end())
2179 {
2180 new_verts.insert(new_face.back());
2181 }
2182 }
2183 ++num_added_faces;
2184 ++last_added_face;
2185 allfaces[last_added_face] = new_face;
2186 ref_elem.push_back(last_added_face);
2187 new_subs.insert(last_added_face);
2188 }
2189
2190 if (num_added_faces == pbnd.m_nbrs_per_face)
2191 {
2192 allfaces[ref_face] = allfaces[last_added_face];
2193 allfaces.erase(last_added_face);
2194 ref_elem.pop_back();
2195 new_subs.erase(last_added_face);
2196 new_subs.insert(ref_face);
2197 }
2198 }
2199 }
2200 }
2201
2202 out_coords["values/x"].set(out_xvec);
2203 out_coords["values/y"].set(out_yvec);
2204 out_coords["values/z"].set(out_zvec);
2205 const double_array& out_xarray =
2206 out_coords["values/x"].as_double_array();
2207 const double_array& out_yarray =
2208 out_coords["values/y"].as_double_array();
2209 const double_array& out_zarray =
2210 out_coords["values/z"].as_double_array();
2211
2212 // Elements that had vertices contributed by multiple neighbor
2213 // domains need a fix-up step.
2214 for (auto mn = multi_nbrs.begin(); mn != multi_nbrs.end(); ++mn)
2215 {
2216 const index_t& offset = *mn;
2217 fix_duplicated_vertices(new_vertices[offset], poly_elems[offset],
2218 allfaces, out_xarray, out_yarray,
2219 out_zarray);
2220 }
2221 }
2222
2223 itr = n.children();
2224 while(itr.has_next())
2225 {
2226 const Node& chld = itr.next();
2227 std::string domain_name = itr.name();
2228
2229 index_t domain_id = chld["state/domain_id"].to_index_t();
2230
2231 const Node& in_topo = chld["topologies"][name];
2232
2233 index_t iwidth = in_topo["elements/dims/i"].to_index_t();
2234 index_t jwidth = in_topo["elements/dims/j"].to_index_t();
2235 index_t kwidth = in_topo["elements/dims/k"].to_index_t();
2236
2237 index_t istride = 1;
2238 index_t jstride = iwidth+1;
2239 index_t kstride = (iwidth+1)*(jwidth+1);
2240
2241 const auto& ratio = nbr_ratio[domain_id];
2242
2243 Node& out_coords = dest[domain_name]["coordsets/coords"];
2244
2245 const double_array& xarray =
2246 out_coords["values/x"].as_double_array();
2247 const double_array& yarray =
2248 out_coords["values/y"].as_double_array();
2249 const double_array& zarray =
2250 out_coords["values/z"].as_double_array();
2251
2252 auto& poly_elems = poly_elems_map[domain_id];
2253 auto& new_vertices = new_vert_map[domain_id];
2254 auto& new_subelems = new_face_map[domain_id];
2255 auto& allfaces = allfaces_map[domain_id];
2256
2257 index_t elemsize = iwidth*jwidth*kwidth;
2258
2259 for (index_t elem = 0; elem < elemsize; ++elem)
2260 {
2261 auto& poly_elem = poly_elems[elem];
2262
2263 //Unchanged elems have 6 sides, no need to do more
2264 if (poly_elem.size() > 6)
2265 {
2266 auto& new_faces = new_subelems[elem];
2267
2268 auto& new_verts = new_vertices[elem];
2269
2270 for (auto pi = poly_elem.begin(); pi != poly_elem.end();
2271 ++pi)
2272 {
2273 index_t pface = *pi;
2274 //only if "pface" is an old (coarse) subelem
2275 if (new_faces.find(pface) == new_faces.end())
2276 {
2277 auto& subelem = allfaces[pface];
2278 if (subelem.size() > 4)
2279 {
2280 //This subelem has already added new vertices
2281 continue;
2282 }
2283
2284 //new subelems that are adjacent to pface.
2285 std::set<index_t> edge_verts;
2286 for (auto vi = subelem.begin(); vi != subelem.end();
2287 ++vi)
2288 {
2289 for (auto si = new_faces.begin();
2290 si != new_faces.end(); ++si)
2291 {
2292 auto& nface = allfaces[*si];
2293 for (auto ni = nface.begin();
2294 ni != nface.end(); ++ni)
2295 {
2296 if (*ni == *vi)
2297 {
2298 edge_verts.insert(*vi);
2299 }
2300 }
2301 }
2302 }
2303
2304 if (edge_verts.empty())
2305 {
2306 continue;
2307 }
2308
2309
2310 size_t num_edges;
2311 if (edge_verts.size() == 4)
2312 {
2313 //Special case when all 4 vertices of subelem
2314 //are selected. There may be only 3 edges that
2315 //get new vertices but we have to check them all.
2316 num_edges = edge_verts.size();
2317 }
2318 else
2319 {
2320 num_edges = edge_verts.size() - 1;
2321 }
2322
2323 size_t edge_count = 0;
2324 std::vector<std::vector<index_t> > edges(num_edges);
2325
2326 for (auto vi = subelem.begin(); vi != subelem.end();
2327 ++vi)
2328 {
2329 if (edge_count == num_edges)
2330 {
2331 break;
2332 }
2333 auto next = vi + 1;
2334 if (next == subelem.end())
2335 {
2336 next = subelem.begin();
2337 }
2338 if (edge_verts.find(*vi) != edge_verts.end() &&
2339 edge_verts.find(*next) != edge_verts.end())
2340 {
2341 edges[edge_count].push_back(*vi);
2342 edges[edge_count].push_back(*next);
2343 ++edge_count;
2344 }
2345 }
2346
2347
2348 //The edges are going to get more vertices.
2349 //Figure out which of the "new verts" fit on these
2350 //edges.
2351
2352 for (auto ei = edges.begin(); ei != edges.end();
2353 ++ei)
2354 {
2355 auto this_edge = *ei;
2356
2357 index_t stride =
2358 std::abs(this_edge[1]-this_edge[0]);
2359
2360 index_t edge_ratio;
2361 if (stride == istride)
2362 {
2363 edge_ratio = ratio[0];
2364 }
2365 else if (stride == jstride)
2366 {
2367 edge_ratio = ratio[1];
2368 }
2369 else
2370 {
2371 assert(stride == kstride);
2372 edge_ratio = ratio[2];
2373 }
2374
2375 //Test each new_vert point (xn,yn,zn) to see
2376 //if it's colinear with (x0,y0,z0) and (x1,y1,z1).
2377 //
2378 //We need something better than these
2379 //epsilon-based floating point comparisons
2380 double x0 = xarray[this_edge[0]];
2381 double y0 = yarray[this_edge[0]];
2382 double z0 = zarray[this_edge[0]];
2383 double x1 = xarray[this_edge[1]];
2384 double y1 = yarray[this_edge[1]];
2385 double z1 = zarray[this_edge[1]];
2386
2387 double eps = sqrt(std::numeric_limits<double>::epsilon());
2388 double xedge = x1-x0+eps;
2389 double yedge = y1-y0+eps;
2390 double zedge = z1-z0+eps;
2391
2392 std::list<index_t> add_verts;
2393 std::multimap<double,index_t> test_verts;
2394 size_t verts_needed =
2395 static_cast<size_t>(edge_ratio-1);
2396 for (auto nv = new_verts.begin(); nv != new_verts.end();
2397 ++nv)
2398 {
2399 index_t nvert = *nv;
2400 double xn = xarray[nvert];
2401 double yn = yarray[nvert];
2402 double zn = zarray[nvert];
2403 double xnedge = (xn-x0)*(xn-x0) > eps ? xn-x0 : 0.0;
2404 double ynedge = (yn-y0)*(yn-y0) > eps ? yn-y0 : 0.0;
2405 double znedge = (zn-z0)*(zn-z0) > eps ? zn-z0 : 0.0;
2406 double edgesum = xnedge+ynedge+znedge;
2407
2408 double xlam = xnedge/xedge;
2409 double ylam = ynedge/yedge;
2410 double zlam = znedge/zedge;
2411 double lam = (xnedge+ynedge+znedge) /
2412 (xedge+yedge+zedge);
2413 double testsum = 0.0;
2414 if (xlam*xlam > eps)
2415 {
2416 testsum += (xlam-lam)*(xlam-lam);
2417 }
2418 if (ylam*ylam > eps)
2419 {
2420 testsum += (ylam-lam)*(ylam-lam);
2421 }
2422 if (zlam*zlam > eps)
2423 {
2424 testsum += (zlam-lam)*(zlam-lam);
2425 }
2426 if (testsum > eps && edgesum/testsum < eps)
2427 {
2428 continue;
2429 }
2430 else
2431 {
2432 test_verts.emplace(testsum, nvert);
2433 }
2434 }
2435
2436 assert(test_verts.size() >= verts_needed ||
2437 num_edges == 4);
2438
2439 for (auto tv = test_verts.begin();
2440 tv != test_verts.end(); ++tv)
2441 {
2442 add_verts.push_back(tv->second);
2443 if (add_verts.size() == verts_needed)
2444 {
2445 break;
2446 }
2447 }
2448
2449 if (add_verts.size() != 1)
2450 {
2451 std::multimap<double,index_t> map_verts;
2452
2453 std::list<index_t> tmp_list;
2454 tmp_list.swap(add_verts);
2455 for (auto av = tmp_list.begin();
2456 av != tmp_list.end(); ++av)
2457 {
2458 double xa = xarray[*av];
2459 double ya = yarray[*av];
2460 double za = zarray[*av];
2461
2462 double dstsq = (xa-x0)*(xa-x0) +
2463 (ya-y0)*(ya-y0) +
2464 (za-z0)*(za-z0);
2465 map_verts.emplace(dstsq, *av);
2466 }
2467 for (auto mv = map_verts.begin();
2468 mv != map_verts.end(); ++mv)
2469 {
2470 add_verts.push_back(mv->second);
2471 new_verts.erase(mv->second);
2472 }
2473 }
2474 else
2475 {
2476 new_verts.erase(add_verts.front());
2477 }
2478
2479 if (this_edge[0] == subelem.back())
2480 {
2481 subelem.insert(subelem.end(),
2482 add_verts.begin(), add_verts.end());
2483 }
2484 else
2485 {
2486 index_t ctr = 0;
2487 for (auto vi = subelem.begin(); vi != subelem.end();
2488 ++vi)
2489 {
2490 if (this_edge[0] == *vi)
2491 {
2492 subelem.insert(vi+1,
2493 add_verts.begin(), add_verts.end());
2494 break;
2495 }
2496 ++ctr;
2497 }
2498 }
2499 }
2500 }
2501 }
2502 }
2503 }
2504 }
2505
2506 itr = n.children();
2507 while(itr.has_next())
2508 {
2509 const Node& chld = itr.next();
2510 std::string domain_name = itr.name();
2511
2512 index_t domain_id = chld["state/domain_id"].to_index_t();
2513
2514 const Node& in_topo = chld["topologies"][name];
2515
2516 index_t iwidth = in_topo["elements/dims/i"].to_index_t();
2517 index_t jwidth = in_topo["elements/dims/j"].to_index_t();
2518 index_t kwidth = in_topo["elements/dims/k"].to_index_t();
2519
2520 auto& poly_elems = poly_elems_map[domain_id];
2521 auto& allfaces = allfaces_map[domain_id];
2522
2523 index_t elemsize = iwidth*jwidth*kwidth;
2524
2525 std::vector<index_t>& e_connect = elem_connect[domain_id];
2526 std::vector<index_t>& e_sizes = elem_sizes[domain_id];
2527 std::vector<index_t>& e_offsets = elem_offsets[domain_id];
2528 std::vector<index_t>& sub_connect = subelem_connect[domain_id];
2529 std::vector<index_t>& sub_sizes = subelem_sizes[domain_id];
2530 std::vector<index_t>& sub_offsets = subelem_offsets[domain_id];
2531 index_t elem_offset_sum = 0;
2532 index_t subelem_offset_sum = 0;
2533 for (index_t elem = 0; elem < elemsize; ++elem)
2534 {
2535 auto& poly_elem = poly_elems[elem];
2536 e_connect.insert(e_connect.end(), poly_elem.begin(), poly_elem.end());
2537 e_sizes.push_back(poly_elem.size());
2538 e_offsets.push_back(elem_offset_sum);
2539 elem_offset_sum += e_sizes.back();
2540 }
2541 for (auto if_itr = allfaces.begin(); if_itr != allfaces.end(); ++if_itr)
2542 {
2543 auto& if_elem = if_itr->second;
2544 sub_connect.insert(sub_connect.end(), if_elem.begin(), if_elem.end());
2545 sub_sizes.push_back(if_elem.size());
2546 sub_offsets.push_back(subelem_offset_sum);
2547 subelem_offset_sum += sub_sizes.back();
2548 }
2549 }
2550
2551 itr = n.children();
2552 while(itr.has_next())
2553 {
2554 const Node& chld = itr.next();
2555 std::string domain_name = itr.name();
2556
2557 index_t domain_id = chld["state/domain_id"].to_index_t();
2558
2559 std::vector<index_t>& e_connect = elem_connect[domain_id];
2560 std::vector<index_t>& e_sizes = elem_sizes[domain_id];
2561 std::vector<index_t>& e_offsets = elem_offsets[domain_id];
2562 std::vector<index_t>& sub_connect = subelem_connect[domain_id];
2563 std::vector<index_t>& sub_sizes = subelem_sizes[domain_id];
2564 std::vector<index_t>& sub_offsets = subelem_offsets[domain_id];
2565
2566 Node& topo = dest[domain_name]["topologies"][name];
2567 const Node& in_topo = chld["topologies"][name];
2568
2569 topo["coordset"] = in_topo["coordset"];
2570
2571 topo["type"] = "unstructured";
2572 topo["elements/shape"] = "polyhedral";
2573 topo["elements/shape"].set_string("polyhedral");
2574 topo["elements/connectivity"].set(e_connect);
2575 topo["elements/sizes"].set(e_sizes);
2576 topo["elements/offsets"].set(e_offsets);
2577 topo["subelements/shape"].set_string("polygonal");
2578 topo["subelements/connectivity"].set(sub_connect);
2579 topo["subelements/sizes"].set(sub_sizes);
2580 topo["subelements/offsets"].set(sub_offsets);
2581
2582 }
2583 }
2584
2585 //
2586 // This function is responsible for collecting the domains within the given mesh
2587 // with subnodes of the given maps based on the domain's path, e.g.:
2588 //
2589 // input:
2590 // mesh: {"domain0": {/*1*/}, "domain1": {/*2*/}}
2591 // s2dmap: {}
2592 // d2smap: {}
2593 //
2594 // output:
2595 // mesh: {"domain0": {/*1*/}, "domain1": {/*2*/}}
2596 // s2dmap: {"domain0": {/*A*/}, "domain1": {/*B*/}}
2597 // d2smap: {"domain0": {/*a*/}, "domain1": {/*b*/}}
2598 // return: [<&1, &A, &a>, <&2, &B, &b>]
2599 //
2600 std::vector<DomMapsTuple>
group_domains_and_maps(conduit::Node & mesh,conduit::Node & s2dmap,conduit::Node & d2smap)2601 group_domains_and_maps(conduit::Node &mesh, conduit::Node &s2dmap, conduit::Node &d2smap)
2602 {
2603 std::vector<DomMapsTuple> doms_and_maps;
2604
2605 s2dmap.reset();
2606 d2smap.reset();
2607
2608 if(!conduit::blueprint::mesh::is_multi_domain(mesh))
2609 {
2610 doms_and_maps.emplace_back(&mesh, &s2dmap, &d2smap);
2611 }
2612 else
2613 {
2614 NodeIterator domains_it = mesh.children();
2615 while(domains_it.has_next())
2616 {
2617 conduit::Node& domain = domains_it.next();
2618 if(mesh.dtype().is_object())
2619 {
2620 doms_and_maps.emplace_back(&domain,
2621 &s2dmap[domain.name()],
2622 &d2smap[domain.name()]);
2623 }
2624 else
2625 {
2626 doms_and_maps.emplace_back(&domain,
2627 &s2dmap.append(),
2628 &d2smap.append());
2629 }
2630 }
2631 }
2632
2633 return std::vector<DomMapsTuple>(std::move(doms_and_maps));
2634 }
2635
2636
2637 //-----------------------------------------------------------------------------
2638 std::vector<index_t>
calculate_decomposed_dims(const conduit::Node & mesh,const std::string & adjset_name,CalcDimDecomposedFun calc_dims)2639 calculate_decomposed_dims(const conduit::Node &mesh, const std::string &adjset_name, CalcDimDecomposedFun calc_dims)
2640 {
2641 // NOTE(JRC): This strategy works even if some ranks have empty meshes because
2642 // the empty ranks won't use the resulting 'dims' array to compute indices and
2643 // thus it doesn't matter if the variable has the technically incorrect value.
2644 const std::vector<const Node *> domains = ::conduit::blueprint::mesh::domains(mesh);
2645 if(domains.empty())
2646 {
2647 return std::vector<index_t>();
2648 }
2649 else // if(!domains.empty())
2650 {
2651 const Node &domain = *domains.front();
2652
2653 const Node &adjset = domain["adjsets"][adjset_name];
2654 const Node *topo_ptr = bputils::find_reference_node(adjset, "topology");
2655 const Node &topo = *topo_ptr;
2656
2657 const bputils::ShapeType shape(topo);
2658 return calc_dims(shape);
2659 }
2660 }
2661
2662
2663 //-----------------------------------------------------------------------------
2664 void
verify_generate_mesh(const conduit::Node & mesh,const std::string & adjset_name)2665 verify_generate_mesh(const conduit::Node &mesh,
2666 const std::string &adjset_name)
2667 {
2668 const std::vector<const Node *> domains = blueprint::mesh::domains(mesh);
2669 for(index_t di = 0; di < (index_t)domains.size(); di++)
2670 {
2671 const Node &domain = *domains[di];
2672 Node info;
2673
2674 if(!domain["adjsets"].has_child(adjset_name))
2675 {
2676 CONDUIT_ERROR("<blueprint::mpi::mesh::generate_*> " <<
2677 "Requested source adjacency set '" << adjset_name << "' " <<
2678 "doesn't exist on domain '" << domain.name() << ".'");
2679 }
2680
2681 if(domain["adjsets"][adjset_name]["association"].as_string() != "vertex")
2682 {
2683 CONDUIT_ERROR("<blueprint::mpi::mesh::generate_*> " <<
2684 "Given adjacency set has an unsupported association type 'element.'\n" <<
2685 "Supported associations:\n" <<
2686 " 'vertex'");
2687 }
2688
2689 const Node &adjset = domain["adjsets"][adjset_name];
2690 const Node *topo_ptr = bputils::find_reference_node(adjset, "topology");
2691 const Node &topo = *topo_ptr;
2692 if(!conduit::blueprint::mesh::topology::unstructured::verify(topo, info))
2693 {
2694 CONDUIT_ERROR("<blueprint::mpi::mesh::generate_*> " <<
2695 "Requested source topology '" << topo.name() << "' " <<
2696 "is of unsupported type '" << topo["type"].as_string() << ".'\n" <<
2697 "Supported types:\n" <<
2698 " 'unstructured'");
2699 }
2700 }
2701 }
2702
2703
2704 //-----------------------------------------------------------------------------
2705 void
generate_derived_entities(conduit::Node & mesh,const std::string & src_adjset_name,const std::string & dst_adjset_name,const std::string & dst_topo_name,conduit::Node & s2dmap,conduit::Node & d2smap,MPI_Comm,GenDerivedFun generate_derived)2706 generate_derived_entities(conduit::Node &mesh,
2707 const std::string &src_adjset_name,
2708 const std::string &dst_adjset_name,
2709 const std::string &dst_topo_name,
2710 conduit::Node &s2dmap,
2711 conduit::Node &d2smap,
2712 MPI_Comm /*comm*/,
2713 GenDerivedFun generate_derived)
2714 {
2715 const std::vector<DomMapsTuple> doms_and_maps = group_domains_and_maps(mesh, s2dmap, d2smap);
2716 for(index_t di = 0; di < (index_t)doms_and_maps.size(); di++)
2717 {
2718 conduit::Node &domain = *std::get<0>(doms_and_maps[di]);
2719 conduit::Node &domain_s2dmap = *std::get<1>(doms_and_maps[di]);
2720 conduit::Node &domain_d2smap = *std::get<2>(doms_and_maps[di]);
2721
2722 const conduit::Node &src_adjset = domain["adjsets"][src_adjset_name];
2723 const Node *src_topo_ptr = bputils::find_reference_node(src_adjset, "topology");
2724 const Node &src_topo = *src_topo_ptr;
2725
2726 conduit::Node &dst_topo = domain["topologies"][dst_topo_name];
2727 generate_derived(src_topo, dst_topo, domain_s2dmap, domain_d2smap);
2728
2729 conduit::Node &dst_adjset = domain["adjsets"][dst_adjset_name];
2730 dst_adjset.reset();
2731 dst_adjset["association"].set("element");
2732 dst_adjset["topology"].set(dst_topo_name);
2733 }
2734
2735 Node src_data, dst_data;
2736 for(index_t di = 0; di < (index_t)doms_and_maps.size(); di++)
2737 {
2738 conduit::Node &domain = *std::get<0>(doms_and_maps[di]);
2739 const index_t domain_id = domain["state/domain_id"].to_index_t();
2740
2741 const Node *src_topo_ptr = bputils::find_reference_node(domain["adjsets"][src_adjset_name], "topology");
2742 const Node &src_topo = *src_topo_ptr;
2743 const Node *src_cset_ptr = bputils::find_reference_node(src_topo, "coordset");
2744 const Node &src_cset = *src_cset_ptr;
2745
2746 const Node &dst_topo = domain["topologies"][dst_topo_name];
2747 const index_t dst_topo_len = bputils::topology::length(dst_topo);
2748
2749 const conduit::Node &src_adjset_groups = domain["adjsets"][src_adjset_name]["groups"];
2750 conduit::Node &dst_adjset_groups = domain["adjsets"][dst_adjset_name]["groups"];
2751
2752 // Organize Adjset Points into Interfaces (Pair-Wise Groups) //
2753
2754 // {(neighbor domain id): <(participating points for domain interface)>}
2755 std::map<index_t, std::set<index_t>> neighbor_pidxs_map;
2756 for(const std::string &group_name : src_adjset_groups.child_names())
2757 {
2758 const conduit::Node &src_group = src_adjset_groups[group_name];
2759 const conduit::Node &src_neighbors = src_group["neighbors"];
2760 const conduit::Node &src_values = src_group["values"];
2761
2762 for(index_t ni = 0; ni < src_neighbors.dtype().number_of_elements(); ni++)
2763 {
2764 src_data.set_external(DataType(src_neighbors.dtype().id(), 1),
2765 (void*)src_neighbors.element_ptr(ni));
2766 std::set<index_t> &neighbor_pidxs = neighbor_pidxs_map[src_data.to_index_t()];
2767 for(index_t pi = 0; pi < src_values.dtype().number_of_elements(); pi++)
2768 {
2769 src_data.set_external(DataType(src_values.dtype().id(), 1),
2770 (void*)src_values.element_ptr(pi));
2771 neighbor_pidxs.insert(src_data.to_index_t());
2772 }
2773 }
2774 }
2775
2776 // Collect Viable Entities for All Interfaces //
2777
2778 // {(entity id in topology): <(neighbor domain ids that contain this entity)>}
2779 std::map<index_t, std::set<index_t>> entity_neighbor_map;
2780 for(index_t ei = 0; ei < dst_topo_len; ei++)
2781 {
2782 std::vector<index_t> entity_pidxs = bputils::topology::unstructured::points(dst_topo, ei);
2783 for(const auto &neighbor_pair : neighbor_pidxs_map)
2784 {
2785 const index_t &ni = neighbor_pair.first;
2786 const std::set<index_t> &neighbor_pidxs = neighbor_pair.second;
2787
2788 bool entity_in_neighbor = true;
2789 for(index_t pi = 0; pi < (index_t)entity_pidxs.size() && entity_in_neighbor; pi++)
2790 {
2791 entity_in_neighbor &= neighbor_pidxs.find(entity_pidxs[pi]) != neighbor_pidxs.end();
2792 }
2793
2794 if(entity_in_neighbor)
2795 {
2796 entity_neighbor_map[ei].insert(ni);
2797 }
2798 }
2799 }
2800
2801 // Use Entity Interfaces to Construct Group Entity Lists //
2802
2803 std::map<std::set<index_t>, std::vector<std::tuple<std::set<PointTuple>, index_t>>> group_entity_map;
2804 for(const auto &entity_neighbor_pair : entity_neighbor_map)
2805 {
2806 const index_t &ei = entity_neighbor_pair.first;
2807 const std::set<index_t> &entity_neighbors = entity_neighbor_pair.second;
2808 std::tuple<std::set<PointTuple>, index_t> entity;
2809
2810 std::vector<index_t> entity_pidxs = bputils::topology::unstructured::points(dst_topo, ei);
2811 std::set<PointTuple> &entity_points = std::get<0>(entity);
2812 for(const index_t &entity_pidx : entity_pidxs)
2813 {
2814 const std::vector<float64> point_coords = bputils::coordset::_explicit::coords(
2815 src_cset, entity_pidx);
2816 entity_points.emplace(
2817 point_coords[0],
2818 (point_coords.size() > 1) ? point_coords[1] : 0.0,
2819 (point_coords.size() > 2) ? point_coords[2] : 0.0);
2820 }
2821
2822 index_t &entity_id = std::get<1>(entity);
2823 entity_id = ei;
2824
2825 // NOTE(JRC): Inserting with this method allows this algorithm to sort new
2826 // elements as they're generated, rather than as a separate process at the
2827 // end (slight optimization overall).
2828 std::vector<std::tuple<std::set<PointTuple>, index_t>> &group_entities =
2829 group_entity_map[entity_neighbors];
2830 auto entity_itr = std::upper_bound(group_entities.begin(), group_entities.end(), entity);
2831 group_entities.insert(entity_itr, entity);
2832 }
2833
2834 for(const auto &group_pair : group_entity_map)
2835 {
2836 // NOTE(JRC): It's possible for the 'src_adjset_groups' node to be empty,
2837 // so we only want to query child data types if we know there is at least
2838 // 1 non-empty group.
2839 const conduit::DataType src_neighbors_dtype = src_adjset_groups.child(0)["neighbors"].dtype();
2840 const conduit::DataType src_values_dtype = src_adjset_groups.child(0)["values"].dtype();
2841
2842 const std::set<index_t> &group_nidxs = group_pair.first;
2843 const std::vector<std::tuple<std::set<PointTuple>, index_t>> &group_entities = group_pair.second;
2844 std::string group_name;
2845 {
2846 // NOTE(JRC): The current domain is included in the domain name so that
2847 // it matches across all domains and processors (also, using std::set
2848 // ensures that items are sorted and the order is the same across ranks).
2849 std::set<index_t> group_all_nidxs = group_nidxs;
2850 group_all_nidxs.insert(domain_id);
2851
2852 std::ostringstream oss;
2853 oss << "group";
2854 for(const index_t &group_nidx : group_all_nidxs)
2855 {
2856 oss << "_" << group_nidx;
2857 }
2858 group_name = oss.str();
2859 }
2860
2861 conduit::Node &dst_group = dst_adjset_groups[group_name];
2862 conduit::Node &dst_neighbors = dst_group["neighbors"];
2863 conduit::Node &dst_values = dst_group["values"];
2864
2865 dst_neighbors.set(DataType(src_neighbors_dtype.id(), group_nidxs.size()));
2866 index_t ni = 0;
2867 for(auto nitr = group_nidxs.begin(); nitr != group_nidxs.end(); ++nitr)
2868 {
2869 src_data.set_external(DataType::index_t(1),
2870 (void*)&(*nitr));
2871 dst_data.set_external(DataType(src_neighbors_dtype.id(), 1),
2872 (void*)dst_neighbors.element_ptr(ni++));
2873 src_data.to_data_type(dst_data.dtype().id(), dst_data);
2874 }
2875
2876 dst_values.set(DataType(src_values_dtype.id(), group_entities.size()));
2877 for(index_t ei = 0; ei < (index_t)group_entities.size(); ei++)
2878 {
2879 src_data.set_external(DataType::index_t(1),
2880 (void*)&std::get<1>(group_entities[ei]));
2881 dst_data.set_external(DataType(src_values_dtype.id(), 1),
2882 (void*)dst_values.element_ptr(ei));
2883 src_data.to_data_type(dst_data.dtype().id(), dst_data);
2884 }
2885 }
2886 }
2887
2888 // TODO(JRC): Waitall?
2889 }
2890
2891
2892 //-----------------------------------------------------------------------------
2893 void
generate_decomposed_entities(conduit::Node & mesh,const std::string & src_adjset_name,const std::string & dst_adjset_name,const std::string & dst_topo_name,const std::string & dst_cset_name,conduit::Node & s2dmap,conduit::Node & d2smap,MPI_Comm,GenDecomposedFun generate_decomposed,IdDecomposedFun identify_decomposed,const std::vector<index_t> & decomposed_centroid_dims)2894 generate_decomposed_entities(conduit::Node &mesh,
2895 const std::string &src_adjset_name,
2896 const std::string &dst_adjset_name,
2897 const std::string &dst_topo_name,
2898 const std::string &dst_cset_name,
2899 conduit::Node &s2dmap,
2900 conduit::Node &d2smap,
2901 MPI_Comm /*comm*/,
2902 GenDecomposedFun generate_decomposed,
2903 IdDecomposedFun identify_decomposed,
2904 const std::vector<index_t> &decomposed_centroid_dims)
2905 {
2906 const std::vector<DomMapsTuple> doms_and_maps = group_domains_and_maps(mesh, s2dmap, d2smap);
2907 for(index_t di = 0; di < (index_t)doms_and_maps.size(); di++)
2908 {
2909 Node &domain = *std::get<0>(doms_and_maps[di]);
2910 Node &domain_s2dmap = *std::get<1>(doms_and_maps[di]);
2911 Node &domain_d2smap = *std::get<2>(doms_and_maps[di]);
2912
2913 const Node &src_adjset = domain["adjsets"][src_adjset_name];
2914 const Node *src_topo_ptr = bputils::find_reference_node(src_adjset, "topology");
2915 const Node &src_topo = *src_topo_ptr;
2916
2917 // NOTE(JRC): Diff- new code below
2918 Node &dst_topo = domain["topologies"][dst_topo_name];
2919 Node &dst_cset = domain["coordsets"][dst_cset_name];
2920 generate_decomposed(src_topo, dst_topo, dst_cset, domain_s2dmap, domain_d2smap);
2921
2922 Node &dst_adjset = domain["adjsets"][dst_adjset_name];
2923 dst_adjset.reset();
2924 // NOTE(JRC): Diff- different association (decomposed entity -> assoc: vertex)
2925 dst_adjset["association"].set("vertex");
2926 dst_adjset["topology"].set(dst_topo_name);
2927 }
2928
2929 Node src_data, dst_data;
2930 src_data.reset();
2931 dst_data.reset();
2932 for(index_t di = 0; di < (index_t)doms_and_maps.size(); di++)
2933 {
2934 Node &domain = *std::get<0>(doms_and_maps[di]);
2935 const index_t domain_id = domain["state/domain_id"].to_index_t();
2936
2937 const Node *src_topo_ptr = bputils::find_reference_node(domain["adjsets"][src_adjset_name], "topology");
2938 const Node &src_topo = *src_topo_ptr;
2939 const Node *src_cset_ptr = bputils::find_reference_node(src_topo, "coordset");
2940 const Node &src_cset = *src_cset_ptr;
2941 // NOTE(JRC): Diff- generate topology metadata for source topology to find
2942 // centroids that may exist within an adjset group.
2943 const bputils::TopologyMetadata src_topo_data(src_topo, src_cset);
2944
2945 // const Node &dst_topo = domain["topologies"][dst_topo_name];
2946 const Node &dst_cset = domain["coordsets"][dst_cset_name];
2947
2948 const Node &src_adjset_groups = domain["adjsets"][src_adjset_name]["groups"];
2949 Node &dst_adjset_groups = domain["adjsets"][dst_adjset_name]["groups"];
2950
2951 // Organize Adjset Points into Interfaces (Pair-Wise Groups) //
2952
2953 // {(neighbor domain id): <(participating points for domain interface)>}
2954 std::map<index_t, std::set<index_t>> neighbor_pidxs_map;
2955 for(const std::string &group_name : src_adjset_groups.child_names())
2956 {
2957 const Node &src_group = src_adjset_groups[group_name];
2958 const Node &src_neighbors = src_group["neighbors"];
2959 const Node &src_values = src_group["values"];
2960
2961 for(index_t ni = 0; ni < src_neighbors.dtype().number_of_elements(); ni++)
2962 {
2963 src_data.set_external(DataType(src_neighbors.dtype().id(), 1),
2964 (void*)src_neighbors.element_ptr(ni));
2965 std::set<index_t> &neighbor_pidxs = neighbor_pidxs_map[src_data.to_index_t()];
2966 for(index_t pi = 0; pi < src_values.dtype().number_of_elements(); pi++)
2967 {
2968 src_data.set_external(DataType(src_values.dtype().id(), 1),
2969 (void*)src_values.element_ptr(pi));
2970 neighbor_pidxs.insert(src_data.to_index_t());
2971 }
2972 }
2973 }
2974
2975 // Collect Viable Entities for All Interfaces //
2976
2977 // {(entity centroid id): <(neighbor domain ids that contain this entity)>}
2978 std::map<index_t, std::set<index_t>> entity_neighbor_map;
2979 // NOTE(JRC): Diff, entirely different iteration strategy for finding entities
2980 // to consider on individual adjset interfaces.
2981 for(const index_t &di : decomposed_centroid_dims)
2982 {
2983 const Node &dim_topo = src_topo_data.dim_topos[di];
2984 for(index_t ei = 0; ei < src_topo_data.get_length(di); ei++)
2985 {
2986 std::vector<index_t> entity_pidxs = bputils::topology::unstructured::points(dim_topo, ei);
2987 for(const auto &neighbor_pair : neighbor_pidxs_map)
2988 {
2989 const index_t &ni = neighbor_pair.first;
2990 const std::set<index_t> &neighbor_pidxs = neighbor_pair.second;
2991
2992 bool entity_in_neighbor = true;
2993 for(index_t pi = 0; pi < (index_t)entity_pidxs.size() && entity_in_neighbor; pi++)
2994 {
2995 entity_in_neighbor &= neighbor_pidxs.find(entity_pidxs[pi]) != neighbor_pidxs.end();
2996 }
2997
2998 if(entity_in_neighbor)
2999 {
3000 const index_t entity_cidx = identify_decomposed(src_topo_data, ei, di);
3001 entity_neighbor_map[entity_cidx].insert(ni);
3002 }
3003 }
3004 }
3005 }
3006
3007 // Use Entity Interfaces to Construct Group Entity Lists //
3008
3009 std::map<std::set<index_t>, std::vector<std::tuple<std::set<PointTuple>, index_t>>> group_entity_map;
3010 for(const auto &entity_neighbor_pair : entity_neighbor_map)
3011 {
3012 const index_t &entity_cidx = entity_neighbor_pair.first;
3013 const std::set<index_t> &entity_neighbors = entity_neighbor_pair.second;
3014 std::tuple<std::set<PointTuple>, index_t> entity;
3015
3016 std::set<PointTuple> &entity_points = std::get<0>(entity);
3017 // NOTE(JRC): Diff: Substitute entity for centroid point at the end here.
3018 const std::vector<float64> point_coords = bputils::coordset::_explicit::coords(
3019 dst_cset, entity_cidx);
3020 entity_points.emplace(
3021 point_coords[0],
3022 (point_coords.size() > 1) ? point_coords[1] : 0.0,
3023 (point_coords.size() > 2) ? point_coords[2] : 0.0);
3024
3025 index_t &entity_id = std::get<1>(entity);
3026 entity_id = entity_cidx;
3027
3028 // NOTE(JRC): Inserting with this method allows this algorithm to sort new
3029 // elements as they're generated, rather than as a separate process at the
3030 // end (slight optimization overall).
3031 std::vector<std::tuple<std::set<PointTuple>, index_t>> &group_entities =
3032 group_entity_map[entity_neighbors];
3033 auto entity_itr = std::upper_bound(group_entities.begin(), group_entities.end(), entity);
3034 group_entities.insert(entity_itr, entity);
3035 }
3036
3037 for(const auto &group_pair : group_entity_map)
3038 {
3039 // NOTE(JRC): It's possible for the 'src_adjset_groups' node to be empty,
3040 // so we only want to query child data types if we know there is at least
3041 // 1 non-empty group.
3042 const DataType src_neighbors_dtype = src_adjset_groups.child(0)["neighbors"].dtype();
3043 const DataType src_values_dtype = src_adjset_groups.child(0)["values"].dtype();
3044
3045 const std::set<index_t> &group_nidxs = group_pair.first;
3046 const std::vector<std::tuple<std::set<PointTuple>, index_t>> &group_entities = group_pair.second;
3047 std::string group_name;
3048 {
3049 // NOTE(JRC): The current domain is included in the domain name so that
3050 // it matches across all domains and processors (also, using std::set
3051 // ensures that items are sorted and the order is the same across ranks).
3052 std::set<index_t> group_all_nidxs = group_nidxs;
3053 group_all_nidxs.insert(domain_id);
3054
3055 std::ostringstream oss;
3056 oss << "group";
3057 for(const index_t &group_nidx : group_all_nidxs)
3058 {
3059 oss << "_" << group_nidx;
3060 }
3061 group_name = oss.str();
3062 }
3063
3064 Node &dst_group = dst_adjset_groups[group_name];
3065 Node &dst_neighbors = dst_group["neighbors"];
3066 Node &dst_values = dst_group["values"];
3067
3068 dst_neighbors.set(DataType(src_neighbors_dtype.id(), group_nidxs.size()));
3069 index_t ni = 0;
3070 for(auto nitr = group_nidxs.begin(); nitr != group_nidxs.end(); ++nitr)
3071 {
3072 src_data.set_external(DataType::index_t(1),
3073 (void*)&(*nitr));
3074 dst_data.set_external(DataType(src_neighbors_dtype.id(), 1),
3075 (void*)dst_neighbors.element_ptr(ni++));
3076 src_data.to_data_type(dst_data.dtype().id(), dst_data);
3077 }
3078
3079 dst_values.set(DataType(src_values_dtype.id(), group_entities.size()));
3080 for(index_t ei = 0; ei < (index_t)group_entities.size(); ei++)
3081 {
3082 src_data.set_external(DataType::index_t(1),
3083 (void*)&std::get<1>(group_entities[ei]));
3084 dst_data.set_external(DataType(src_values_dtype.id(), 1),
3085 (void*)dst_values.element_ptr(ei));
3086 src_data.to_data_type(dst_data.dtype().id(), dst_data);
3087 }
3088 }
3089 }
3090
3091 // TODO(JRC): Waitall?
3092 }
3093
3094
3095 //-----------------------------------------------------------------------------
3096 void
generate_points(conduit::Node & mesh,const std::string & src_adjset_name,const std::string & dst_adjset_name,const std::string & dst_topo_name,conduit::Node & s2dmap,conduit::Node & d2smap,MPI_Comm comm)3097 generate_points(conduit::Node &mesh,
3098 const std::string &src_adjset_name,
3099 const std::string &dst_adjset_name,
3100 const std::string &dst_topo_name,
3101 conduit::Node &s2dmap,
3102 conduit::Node &d2smap,
3103 MPI_Comm comm)
3104 {
3105 verify_generate_mesh(mesh, src_adjset_name);
3106 generate_derived_entities(
3107 mesh, src_adjset_name, dst_adjset_name, dst_topo_name, s2dmap, d2smap, comm,
3108 conduit::blueprint::mesh::topology::unstructured::generate_points);
3109 }
3110
3111
3112 //-----------------------------------------------------------------------------
3113 void
generate_lines(conduit::Node & mesh,const std::string & src_adjset_name,const std::string & dst_adjset_name,const std::string & dst_topo_name,conduit::Node & s2dmap,conduit::Node & d2smap,MPI_Comm comm)3114 generate_lines(conduit::Node &mesh,
3115 const std::string &src_adjset_name,
3116 const std::string &dst_adjset_name,
3117 const std::string &dst_topo_name,
3118 conduit::Node &s2dmap,
3119 conduit::Node &d2smap,
3120 MPI_Comm comm)
3121 {
3122 verify_generate_mesh(mesh, src_adjset_name);
3123 generate_derived_entities(
3124 mesh, src_adjset_name, dst_adjset_name, dst_topo_name, s2dmap, d2smap, comm,
3125 conduit::blueprint::mesh::topology::unstructured::generate_lines);
3126 }
3127
3128
3129 //-----------------------------------------------------------------------------
3130 void
generate_faces(conduit::Node & mesh,const std::string & src_adjset_name,const std::string & dst_adjset_name,const std::string & dst_topo_name,conduit::Node & s2dmap,conduit::Node & d2smap,MPI_Comm comm)3131 generate_faces(conduit::Node &mesh,
3132 const std::string& src_adjset_name,
3133 const std::string& dst_adjset_name,
3134 const std::string& dst_topo_name,
3135 conduit::Node& s2dmap,
3136 conduit::Node& d2smap,
3137 MPI_Comm comm)
3138 {
3139 verify_generate_mesh(mesh, src_adjset_name);
3140 generate_derived_entities(
3141 mesh, src_adjset_name, dst_adjset_name, dst_topo_name, s2dmap, d2smap, comm,
3142 conduit::blueprint::mesh::topology::unstructured::generate_faces);
3143 }
3144
3145
3146 //-----------------------------------------------------------------------------
3147 void
generate_centroids(conduit::Node & mesh,const std::string & src_adjset_name,const std::string & dst_adjset_name,const std::string & dst_topo_name,const std::string & dst_cset_name,conduit::Node & s2dmap,conduit::Node & d2smap,MPI_Comm comm)3148 generate_centroids(conduit::Node& mesh,
3149 const std::string& src_adjset_name,
3150 const std::string& dst_adjset_name,
3151 const std::string& dst_topo_name,
3152 const std::string& dst_cset_name,
3153 conduit::Node& s2dmap,
3154 conduit::Node& d2smap,
3155 MPI_Comm comm)
3156 {
3157 const static auto identify_centroid = []
3158 (const bputils::TopologyMetadata &/*topo_data*/, const index_t ei, const index_t /*di*/)
3159 {
3160 return ei;
3161 };
3162
3163 const static auto calculate_centroid_dims = [] (const bputils::ShapeType &topo_shape)
3164 {
3165 return std::vector<index_t>(1, topo_shape.dim);
3166 };
3167
3168 verify_generate_mesh(mesh, src_adjset_name);
3169
3170 const std::vector<index_t> centroid_dims = calculate_decomposed_dims(
3171 mesh, src_adjset_name, calculate_centroid_dims);
3172
3173 generate_decomposed_entities(
3174 mesh, src_adjset_name, dst_adjset_name, dst_topo_name, dst_cset_name, s2dmap, d2smap, comm,
3175 conduit::blueprint::mesh::topology::unstructured::generate_centroids, identify_centroid, centroid_dims);
3176 }
3177
3178
3179 //-----------------------------------------------------------------------------
3180 void
generate_sides(conduit::Node & mesh,const std::string & src_adjset_name,const std::string & dst_adjset_name,const std::string & dst_topo_name,const std::string & dst_cset_name,conduit::Node & s2dmap,conduit::Node & d2smap,MPI_Comm comm)3181 generate_sides(conduit::Node& mesh,
3182 const std::string& src_adjset_name,
3183 const std::string& dst_adjset_name,
3184 const std::string& dst_topo_name,
3185 const std::string& dst_cset_name,
3186 conduit::Node& s2dmap,
3187 conduit::Node& d2smap,
3188 MPI_Comm comm)
3189 {
3190 const static auto identify_side = []
3191 (const bputils::TopologyMetadata &topo_data, const index_t ei, const index_t di)
3192 {
3193 index_t doffset = 0;
3194 for(index_t dii = 0; dii < di; dii++)
3195 {
3196 if(dii != 1)
3197 {
3198 doffset += topo_data.get_length(dii);
3199 }
3200 }
3201
3202 return doffset + ei;
3203 };
3204
3205 const static auto calculate_side_dims = [] (const bputils::ShapeType &topo_shape)
3206 {
3207 std::vector<index_t> side_dims;
3208
3209 side_dims.push_back(0);
3210 if(topo_shape.dim == 3)
3211 {
3212 side_dims.push_back(2);
3213 }
3214
3215 return std::vector<index_t>(std::move(side_dims));
3216 };
3217
3218 verify_generate_mesh(mesh, src_adjset_name);
3219
3220 const std::vector<index_t> side_dims = calculate_decomposed_dims(
3221 mesh, src_adjset_name, calculate_side_dims);
3222
3223 generate_decomposed_entities(
3224 mesh, src_adjset_name, dst_adjset_name, dst_topo_name, dst_cset_name, s2dmap, d2smap, comm,
3225 conduit::blueprint::mesh::topology::unstructured::generate_sides, identify_side, side_dims);
3226 }
3227
3228
3229 //-----------------------------------------------------------------------------
3230 void
generate_corners(conduit::Node & mesh,const std::string & src_adjset_name,const std::string & dst_adjset_name,const std::string & dst_topo_name,const std::string & dst_cset_name,conduit::Node & s2dmap,conduit::Node & d2smap,MPI_Comm comm)3231 generate_corners(conduit::Node& mesh,
3232 const std::string& src_adjset_name,
3233 const std::string& dst_adjset_name,
3234 const std::string& dst_topo_name,
3235 const std::string& dst_cset_name,
3236 conduit::Node& s2dmap,
3237 conduit::Node& d2smap,
3238 MPI_Comm comm)
3239 {
3240 const static auto identify_corner = []
3241 (const bputils::TopologyMetadata &topo_data, const index_t ei, const index_t di)
3242 {
3243 index_t doffset = 0;
3244 for(index_t dii = 0; dii < di; dii++)
3245 {
3246 doffset += topo_data.get_length(dii);
3247 }
3248
3249 return doffset + ei;
3250 };
3251
3252 const static auto calculate_corner_dims = [] (const bputils::ShapeType &topo_shape)
3253 {
3254 std::vector<index_t> corner_dims;
3255
3256 for(index_t di = 0; di < topo_shape.dim; di++)
3257 {
3258 corner_dims.push_back(di);
3259 }
3260
3261 return std::vector<index_t>(std::move(corner_dims));
3262 };
3263
3264 verify_generate_mesh(mesh, src_adjset_name);
3265
3266 const std::vector<index_t> corner_dims = calculate_decomposed_dims(
3267 mesh, src_adjset_name, calculate_corner_dims);
3268
3269 generate_decomposed_entities(
3270 mesh, src_adjset_name, dst_adjset_name, dst_topo_name, dst_cset_name, s2dmap, d2smap, comm,
3271 conduit::blueprint::mesh::topology::unstructured::generate_corners, identify_corner, corner_dims);
3272 }
3273
3274 //-------------------------------------------------------------------------
3275 void
flatten(const conduit::Node & mesh,const conduit::Node & options,conduit::Node & output,MPI_Comm comm)3276 flatten(const conduit::Node &mesh, const conduit::Node &options,
3277 conduit::Node &output, MPI_Comm comm)
3278 {
3279 output.reset();
3280
3281 ParallelMeshFlattener do_flatten(comm);
3282 do_flatten.set_options(options);
3283 do_flatten.execute(mesh, output);
3284 }
3285
3286 //-----------------------------------------------------------------------------
3287 }
3288 //-----------------------------------------------------------------------------
3289 // -- end conduit::blueprint::mpi::mesh --
3290 //-----------------------------------------------------------------------------
3291
3292
3293 }
3294 //-----------------------------------------------------------------------------
3295 // -- end conduit::blueprint::mpi --
3296 //-----------------------------------------------------------------------------
3297
3298
3299 }
3300 //-----------------------------------------------------------------------------
3301 // -- end conduit::blueprint --
3302 //-----------------------------------------------------------------------------
3303
3304 }
3305 //-----------------------------------------------------------------------------
3306 // -- end conduit:: --
3307 //-----------------------------------------------------------------------------
3308