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_mesh_utils.cpp
8 ///
9 //-----------------------------------------------------------------------------
10
11 //-----------------------------------------------------------------------------
12 // std lib includes
13 #include <algorithm>
14 #include <cmath>
15 #include <deque>
16 #include <string>
17 #include <limits>
18 #include <map>
19 #include <memory>
20 #include <sstream>
21 #include <vector>
22
23 //-----------------------------------------------------------------------------
24 // conduit includes
25 //-----------------------------------------------------------------------------
26 #include "conduit_blueprint_mesh.hpp"
27 #include "conduit_blueprint_o2mrelation.hpp"
28 #include "conduit_blueprint_o2mrelation_iterator.hpp"
29 #include "conduit_blueprint_mesh_utils.hpp"
30
31 // access one-to-many index types
32 namespace O2MIndex = conduit::blueprint::o2mrelation;
33
34 //-----------------------------------------------------------------------------
35 // -- begin conduit --
36 //-----------------------------------------------------------------------------
37 namespace conduit
38 {
39
40 //-----------------------------------------------------------------------------
41 // -- begin conduit::blueprint --
42 //-----------------------------------------------------------------------------
43 namespace blueprint
44 {
45
46 //-----------------------------------------------------------------------------
47 // -- begin conduit::blueprint::mesh --
48 //-----------------------------------------------------------------------------
49 namespace mesh
50 {
51
52 //-----------------------------------------------------------------------------
53 // -- begin conduit::blueprint::mesh::utils --
54 //-----------------------------------------------------------------------------
55 namespace utils
56 {
57
58 //-----------------------------------------------------------------------------
59 /// blueprint mesh utility structures
60 //-----------------------------------------------------------------------------
61
62 //---------------------------------------------------------------------------//
ShapeType()63 ShapeType::ShapeType()
64 {
65 init(-1);
66 }
67
68 //---------------------------------------------------------------------------//
ShapeType(const index_t type_id)69 ShapeType::ShapeType(const index_t type_id)
70 {
71 init(type_id);
72 }
73
74 //---------------------------------------------------------------------------//
ShapeType(const std::string & type_name)75 ShapeType::ShapeType(const std::string &type_name)
76 {
77 init(type_name);
78 }
79
80 //---------------------------------------------------------------------------//
ShapeType(const conduit::Node & topology)81 ShapeType::ShapeType(const conduit::Node &topology)
82 {
83 init(-1);
84
85 if(topology["type"].as_string() == "unstructured" &&
86 topology["elements"].has_child("shape"))
87 {
88 init(topology["elements/shape"].as_string());
89 }
90 }
91
92
93 //---------------------------------------------------------------------------//
94 void
init(const std::string & type_name)95 ShapeType::init(const std::string &type_name)
96 {
97 init(-1);
98
99 for(index_t i = 0; i < (index_t)TOPO_SHAPES.size(); i++)
100 {
101 if(type_name == TOPO_SHAPES[i])
102 {
103 init(i);
104 }
105 }
106 }
107
108
109 //---------------------------------------------------------------------------//
110 void
init(const index_t type_id)111 ShapeType::init(const index_t type_id)
112 {
113 if(type_id < 0 || type_id >= (index_t)TOPO_SHAPES.size())
114 {
115 type = "";
116 id = dim = indices = embed_id = embed_count = -1;
117 embedding = NULL;
118 }
119 else
120 {
121 type = TOPO_SHAPES[type_id];
122 id = type_id;
123 dim = TOPO_SHAPE_DIMS[type_id];
124 indices = TOPO_SHAPE_INDEX_COUNTS[type_id];
125
126 embed_id = TOPO_SHAPE_EMBED_TYPES[type_id];
127 embed_count = TOPO_SHAPE_EMBED_COUNTS[type_id];
128 embedding = const_cast<index_t*>(TOPO_SHAPE_EMBEDDINGS[type_id]);
129 }
130 }
131
132
133 //---------------------------------------------------------------------------//
134 bool
is_poly() const135 ShapeType::is_poly() const
136 {
137 return embedding == NULL;
138 }
139
140
141 //---------------------------------------------------------------------------//
142 bool
is_polygonal() const143 ShapeType::is_polygonal() const
144 {
145 return embedding == NULL && dim == 2;
146 }
147
148
149 //---------------------------------------------------------------------------//
150 bool
is_polyhedral() const151 ShapeType::is_polyhedral() const
152 {
153 return embedding == NULL && dim == 3;
154 }
155
156
157 //---------------------------------------------------------------------------//
158 bool
is_valid() const159 ShapeType::is_valid() const
160 {
161 return id >= 0;
162 }
163
164
165 //---------------------------------------------------------------------------//
ShapeCascade(const conduit::Node & topology)166 ShapeCascade::ShapeCascade(const conduit::Node &topology)
167 {
168 ShapeType base_type(topology);
169 init(base_type);
170 }
171
172 //---------------------------------------------------------------------------//
ShapeCascade(const ShapeType & base_type)173 ShapeCascade::ShapeCascade(const ShapeType &base_type)
174 {
175 init(base_type);
176 }
177
178 //---------------------------------------------------------------------------//
179 index_t
get_num_embedded(const index_t level) const180 ShapeCascade::get_num_embedded(const index_t level) const
181 {
182 index_t num_embedded = -1;
183
184 if(!get_shape().is_poly())
185 {
186 num_embedded = 1;
187 for(index_t di = level + 1; di <= dim; di++)
188 {
189 num_embedded *= dim_types[di].embed_count;
190 }
191 }
192
193 return num_embedded;
194 }
195
196
197 //---------------------------------------------------------------------------//
198 const ShapeType&
get_shape(const index_t level) const199 ShapeCascade::get_shape(const index_t level) const
200 {
201 return dim_types[level < 0 ? dim : level];
202 }
203
204 //---------------------------------------------------------------------------//
205 void
init(const ShapeType & base_type)206 ShapeCascade::init(const ShapeType &base_type)
207 {
208 dim = base_type.dim;
209
210 dim_types[base_type.dim] = base_type;
211 for(index_t di = base_type.dim - 1; di >= 0; di--)
212 {
213 dim_types[di] = ShapeType(dim_types[di + 1].embed_id);
214 }
215 }
216
217 //---------------------------------------------------------------------------//
TopologyMetadata(const conduit::Node & topology,const conduit::Node & coordset)218 TopologyMetadata::TopologyMetadata(const conduit::Node &topology, const conduit::Node &coordset) :
219 topo(&topology), cset(&coordset),
220 int_dtype(find_widest_dtype(link_nodes(topology, coordset), DEFAULT_INT_DTYPES)),
221 float_dtype(find_widest_dtype(link_nodes(topology, coordset), DEFAULT_FLOAT_DTYPE)),
222 topo_cascade(topology), topo_shape(topology)
223 {
224 // NOTE(JRC): This type current only works at forming associations within
225 // an unstructured topology's hierarchy.
226 Node topo_offsets;
227 topology::unstructured::generate_offsets(*topo, topo_offsets);
228 const index_t topo_num_elems = topo_offsets.dtype().number_of_elements();
229 const index_t topo_num_coords = coordset::length(coordset);
230
231 // Allocate Data Templates for Outputs //
232
233 dim_topos.resize(topo_shape.dim + 1);
234 dim_geid_maps.resize(topo_shape.dim + 1);
235 dim_geassocs_maps.resize(topo_shape.dim + 1);
236 dim_leassocs_maps.resize(topo_shape.dim + 1);
237 dim_le2ge_maps.resize(topo_shape.dim + 1);
238
239 for(index_t di = 0; di < topo_shape.dim; di++)
240 {
241 Node &dim_topo = dim_topos[di];
242 dim_topo.reset();
243 dim_topo["type"].set("unstructured");
244 dim_topo["coordset"].set(topology["coordset"].as_string());
245 dim_topo["elements/shape"].set(topo_cascade.get_shape(di).type);
246 }
247 // NOTE: This is done so that the index values for the top-level entities
248 // can be extracted by the 'get_entity_data' function before DFS
249 // processing below.
250 dim_topos[topo_shape.dim].set_external(topology);
251 dim_topos[topo_shape.dim]["elements/offsets"].set(topo_offsets);
252 std::vector< std::vector<int64> > dim_buffers(topo_shape.dim + 1);
253
254 // Prepare Initial Values for Processing //
255
256 // Temporary nodes to manage pointers to important information (temp) and
257 // associated conversations (data).
258 Node temp, data;
259
260 // NOTE(JRC): A 'deque' is used so that queue behavior (FIFO)
261 // is responsible for ordering the identifiers in the cascade of entities,
262 // which more closely follows the average intuition.
263 const index_t bag_num_elems = topo_num_coords + topo_num_elems;
264 std::deque< std::vector<int64> > entity_index_bag(bag_num_elems);
265 std::deque< index_t > entity_dim_bag(bag_num_elems, -1);
266 std::deque< std::vector< std::pair<int64, int64> > > entity_parent_bag(bag_num_elems);
267
268 // NOTE(JRC): We start with processing the points of the topology followed
269 // by the top-level elements in order to ensure that order is preserved
270 // relative to the original topology for these entities.
271 for(index_t pi = 0; pi < topo_num_coords; pi++)
272 {
273 index_t bi = pi;
274 entity_index_bag[bi].push_back(bi);
275 entity_dim_bag[bi] = 0;
276 }
277 for(index_t ei = 0; ei < topo_num_elems; ei++)
278 {
279 index_t bi = topo_num_coords + ei;
280
281 temp.reset();
282 get_entity_data(TopologyMetadata::GLOBAL, ei, topo_shape.dim, temp);
283
284 std::vector<int64> &elem_indices = entity_index_bag[bi];
285 elem_indices.resize(temp.dtype().number_of_elements());
286 data.set_external(DataType::int64(elem_indices.size()), &elem_indices[0]);
287 temp.to_int64_array(data);
288
289 entity_dim_bag[bi] = topo_shape.dim;
290 }
291
292 while(!entity_index_bag.empty())
293 {
294 std::vector<int64> entity_indices = entity_index_bag.front();
295 entity_index_bag.pop_front();
296 index_t entity_dim = entity_dim_bag.front();
297 entity_dim_bag.pop_front();
298 std::vector< std::pair<int64, int64> > entity_parents = entity_parent_bag.front();
299 entity_parent_bag.pop_front();
300
301 std::vector<int64> &dim_buffer = dim_buffers[entity_dim];
302 std::map< std::set<index_t>, index_t > &dim_geid_map = dim_geid_maps[entity_dim];
303 auto &dim_geassocs = dim_geassocs_maps[entity_dim];
304 auto &dim_leassocs = dim_leassocs_maps[entity_dim];
305 std::vector<index_t> &dim_le2ge_map = dim_le2ge_maps[entity_dim];
306 ShapeType dim_shape = topo_cascade.get_shape(entity_dim);
307
308 // Add Element to Topology //
309
310 // NOTE: This code assumes that all entities can be uniquely
311 // identified by the list of coordinate indices of which they
312 // are comprised. This is certainly true of all implicit topologies
313 // and of 2D polygonal topologies, but it may not be always the
314 // case for 3D polygonal topologies.
315 std::set<int64> vert_ids;
316 if(!dim_shape.is_polyhedral())
317 {
318 vert_ids = std::set<int64>(entity_indices.begin(), entity_indices.end());
319 }
320 else // if(dim_shape.is_polyhedral())
321 {
322 const index_t elem_outer_count = entity_indices.size();
323 for(index_t oi = 0; oi < elem_outer_count; oi++)
324 {
325 temp.set_external((*topo)["subelements/offsets"]);
326 data.set_external(int_dtype, temp.element_ptr(entity_indices[oi]));
327 const index_t elem_inner_offset = data.to_index_t();
328
329 temp.set_external((*topo)["subelements/sizes"]);
330 data.set_external(int_dtype, temp.element_ptr(entity_indices[oi]));
331 const index_t elem_inner_count = data.to_index_t();
332
333 for(index_t ii = 0; ii < elem_inner_count; ii++)
334 {
335 temp.set_external((*topo)["subelements/connectivity"]);
336 data.set_external(int_dtype, temp.element_ptr(elem_inner_offset + ii));
337 const index_t vi = data.to_int64();
338 vert_ids.insert(vi);
339 }
340 }
341 }
342
343 const index_t local_id = dim_leassocs.size();
344 if(dim_geid_map.find(vert_ids) == dim_geid_map.end())
345 {
346 const index_t global_id = dim_geassocs.size();
347 dim_buffer.insert(dim_buffer.end(), entity_indices.begin(), entity_indices.end());
348 dim_geid_map[vert_ids] = global_id;
349 }
350 const index_t global_id = dim_geid_map.find(vert_ids)->second;
351
352 { // create_entity(global_id, local_id, entity_dim)
353 if((index_t)dim_geassocs.size() <= global_id)
354 {
355 dim_geassocs.resize(global_id + 1);
356 }
357 if((index_t)dim_leassocs.size() <= local_id)
358 {
359 dim_leassocs.resize(local_id + 1);
360 }
361 if((index_t)dim_le2ge_map.size() <= local_id)
362 {
363 dim_le2ge_map.resize(local_id + 1);
364 }
365 dim_le2ge_map[local_id] = global_id;
366 }
367
368 // Add Element to Associations //
369
370 add_entity_assoc(IndexType::GLOBAL, global_id, entity_dim, global_id, entity_dim);
371 add_entity_assoc(IndexType::LOCAL, local_id, entity_dim, local_id, entity_dim);
372 for(index_t pi = 0; pi < (index_t)entity_parents.size(); pi++)
373 {
374 index_t plevel = entity_parents.size() - pi - 1;
375 const index_t parent_global_id = entity_parents[plevel].first;
376 const index_t parent_local_id = entity_parents[plevel].second;
377 const index_t parent_dim = entity_dim + pi + 1;
378 add_entity_assoc(IndexType::GLOBAL, global_id, entity_dim, parent_global_id, parent_dim);
379 add_entity_assoc(IndexType::LOCAL, local_id, entity_dim, parent_local_id, parent_dim);
380 }
381
382 // Add Embedded Elements for Further Processing //
383
384 if(entity_dim > 0)
385 {
386 std::vector< std::pair<int64, int64> > embed_parents = entity_parents;
387 embed_parents.push_back(std::make_pair(global_id, local_id));
388 ShapeType embed_shape = topo_cascade.get_shape(entity_dim - 1);
389
390 index_t elem_outer_count = dim_shape.is_poly() ?
391 entity_indices.size() : dim_shape.embed_count;
392
393 // NOTE(JRC): This is horribly complicated for the poly case and needs
394 // to be refactored so that it's legible. There's a lot of overlap in
395 // used variables where it feels unnecessary (e.g. 'poly' being
396 // shoehorned into using 'implicit' variables), for example.
397 for(index_t oi = 0, ooff = 0; oi < elem_outer_count; oi++)
398 {
399 index_t elem_inner_count = embed_shape.indices;
400
401 if (dim_shape.is_polyhedral())
402 {
403 const Node &subelem_off_const = (*topo)["subelements/offsets"];
404 const Node &subelem_size_const = (*topo)["subelements/sizes"];
405
406 Node subelem_off; subelem_off.set_external(subelem_off_const);
407 Node subelem_size; subelem_size.set_external(subelem_size_const);
408
409 temp.set_external(int_dtype,
410 subelem_off.element_ptr(entity_indices[oi]));
411 ooff = temp.to_int64();
412 temp.set_external(int_dtype,
413 subelem_size.element_ptr(entity_indices[oi]));
414 elem_inner_count = temp.to_int64();
415 }
416
417 std::vector<int64> embed_indices;
418 for(index_t ii = 0; ii < elem_inner_count; ii++)
419 {
420 index_t ioff = ooff + (dim_shape.is_poly() ?
421 ii : dim_shape.embedding[oi * elem_inner_count + ii]);
422
423 if (dim_shape.is_polyhedral())
424 {
425 const Node &subelem_conn_const = (*topo)["subelements/connectivity"];
426 Node subelem_conn; subelem_conn.set_external(subelem_conn_const);
427
428 temp.set_external(int_dtype,
429 subelem_conn.element_ptr(ioff));
430 embed_indices.push_back(temp.to_int64());
431 }
432 else
433 {
434 embed_indices.push_back(
435 entity_indices[ioff % entity_indices.size()]);
436 }
437 }
438
439 ooff += dim_shape.is_polygonal() ? 1 : 0;
440
441 entity_index_bag.push_back(embed_indices);
442 entity_dim_bag.push_back(embed_shape.dim);
443 entity_parent_bag.push_back(embed_parents);
444 }
445 }
446 }
447
448 // Move Topological Data into Per-Dim Nodes //
449
450 for(index_t di = 0; di <= topo_shape.dim; di++)
451 {
452 Node &dim_conn = dim_topos[di]["elements/connectivity"];
453 dim_conn.set(DataType(int_dtype.id(), dim_buffers[di].size()));
454
455 data.reset();
456 data.set_external(DataType::int64(dim_buffers[di].size()),
457 &(dim_buffers[di][0]));
458 data.to_data_type(int_dtype.id(), dim_conn);
459
460 // Initialize element sizes for 2D polygonal mesh generating
461 // from 3D polyhedral mesh
462 if(di == 2 && topo_shape.is_polyhedral())
463 {
464 Node &poly_sizes = dim_topos[di]["elements/sizes"];
465 poly_sizes.set(DataType(int_dtype.id(), dim_geid_maps[di].size()));
466
467 temp.reset();
468 data.reset();
469
470 for(const auto &poly_pair : dim_geid_maps[di])
471 {
472 const std::set<index_t> &poly_verts = poly_pair.first;
473 const index_t &poly_geid = poly_pair.second;
474
475 temp.set_external(DataType(int_dtype.id(), 1),
476 poly_sizes.element_ptr(poly_geid));
477 data.set((index_t)poly_verts.size());
478 data.to_data_type(int_dtype.id(), temp);
479 }
480 }
481
482 Node dim_topo_offsets;
483 topology::unstructured::generate_offsets(dim_topos[di], dim_topo_offsets);
484 }
485 }
486
487 //---------------------------------------------------------------------------//
488 void
add_entity_assoc(IndexType type,index_t e0_id,index_t e0_dim,index_t e1_id,index_t e1_dim)489 TopologyMetadata::add_entity_assoc(IndexType type, index_t e0_id, index_t e0_dim, index_t e1_id, index_t e1_dim)
490 {
491 auto &assoc_maps = (type == IndexType::LOCAL) ? dim_leassocs_maps : dim_geassocs_maps;
492 std::vector< std::pair< std::vector<index_t>, std::set<index_t> > > *entity_assocs[2] = {
493 &assoc_maps[e0_dim][e0_id],
494 &assoc_maps[e1_dim][e1_id]
495 };
496
497 for(index_t ai = 0; ai < 2; ai++)
498 {
499 auto &curr_assocs = *entity_assocs[ai];
500 curr_assocs.resize(topo_shape.dim + 1);
501
502 const index_t cross_id = (ai == 0) ? e1_id : e0_id;
503 const index_t cross_dim = (ai == 0) ? e1_dim : e0_dim;
504 auto &cross_assocs = curr_assocs[cross_dim];
505 if(cross_assocs.second.find(cross_id) == cross_assocs.second.end())
506 {
507 cross_assocs.first.push_back(cross_id);
508 cross_assocs.second.insert(cross_id);
509 }
510 }
511 }
512
513
514 //---------------------------------------------------------------------------//
515 const std::vector<index_t>&
get_entity_assocs(IndexType type,index_t entity_id,index_t entity_dim,index_t assoc_dim) const516 TopologyMetadata::get_entity_assocs(IndexType type, index_t entity_id, index_t entity_dim, index_t assoc_dim) const
517 {
518 auto &dim_assocs = (type == IndexType::LOCAL) ? dim_leassocs_maps : dim_geassocs_maps;
519 return dim_assocs[entity_dim][entity_id][assoc_dim].first;
520 }
521
522
523 //---------------------------------------------------------------------------//
524 void
get_dim_map(IndexType type,index_t src_dim,index_t dst_dim,Node & map_node) const525 TopologyMetadata::get_dim_map(IndexType type, index_t src_dim, index_t dst_dim, Node &map_node) const
526 {
527 auto &dim_assocs = (type == IndexType::LOCAL) ? dim_leassocs_maps : dim_geassocs_maps;
528
529 std::vector<index_t> values, sizes, offsets;
530 for(index_t sdi = 0, so = 0; sdi < (index_t)dim_assocs[src_dim].size(); sdi++, so += sizes.back())
531 {
532 const std::vector<index_t> &src_assocs = get_entity_assocs(type, sdi, src_dim, dst_dim);
533 values.insert(values.end(), src_assocs.begin(), src_assocs.end());
534 sizes.push_back((index_t)src_assocs.size());
535 offsets.push_back(so);
536 }
537
538 std::vector<index_t>* path_data[] = { &values, &sizes, &offsets };
539 std::string path_names[] = { "values", "sizes", "offsets" };
540 const index_t path_count = sizeof(path_data) / sizeof(path_data[0]);
541 for(index_t pi = 0; pi < path_count; pi++)
542 {
543 Node data;
544 data.set(*path_data[pi]);
545 data.to_data_type(int_dtype.id(), map_node[path_names[pi]]);
546 }
547 }
548
549
550 //---------------------------------------------------------------------------//
551 void
get_entity_data(IndexType type,index_t entity_id,index_t entity_dim,Node & data) const552 TopologyMetadata::get_entity_data(IndexType type, index_t entity_id, index_t entity_dim, Node &data) const
553 {
554 Node temp;
555
556 // NOTE(JRC): This is done in order to get around 'const' casting for
557 // data pointers that won't be changed by the function anyway.
558 Node dim_conn; dim_conn.set_external(dim_topos[entity_dim]["elements/connectivity"]);
559 Node dim_off; dim_off.set_external(dim_topos[entity_dim]["elements/offsets"]);
560
561 const DataType conn_dtype(dim_conn.dtype().id(), 1);
562 const DataType off_dtype(dim_off.dtype().id(), 1);
563 const DataType data_dtype = data.dtype().is_number() ? data.dtype() : DataType::int64(1);
564
565 // FIXME(JRC): This code assumes that the per-element index data is packed
566 // in memory, which isn't guaranteed to be the case (could be stride between
567 // values, etc.).
568
569 const index_t entity_gid = (type == IndexType::LOCAL) ?
570 dim_le2ge_maps[entity_dim][entity_id] : entity_id;
571 temp.set_external(off_dtype, dim_off.element_ptr(entity_gid));
572 index_t entity_start_index = temp.to_int64();
573 temp.set_external(off_dtype, dim_off.element_ptr(entity_gid + 1));
574 index_t entity_end_index = (entity_gid < get_length(entity_dim) - 1) ?
575 temp.to_int64() : dim_conn.dtype().number_of_elements();
576
577 index_t entity_size = entity_end_index - entity_start_index;
578 temp.set_external(DataType(conn_dtype.id(), entity_size),
579 dim_conn.element_ptr(entity_start_index));
580 temp.to_data_type(data_dtype.id(), data);
581 }
582
583
584 //---------------------------------------------------------------------------//
585 void
get_point_data(IndexType type,index_t point_id,Node & data) const586 TopologyMetadata::get_point_data(IndexType type, index_t point_id, Node &data) const
587 {
588 const index_t point_gid = (type == IndexType::LOCAL) ?
589 dim_le2ge_maps[0][point_id] : point_id;
590
591 if(data.dtype().is_empty())
592 {
593 data.set(DataType::float64(3));
594 }
595 const DataType data_dtype(data.dtype().id(), 1);
596
597 Node temp1, temp2;
598 const std::vector<std::string> csys_axes = coordset::axes(*cset);
599 for(index_t di = 0; di < topo_shape.dim; di++)
600 {
601 temp1.set_external(float_dtype,
602 (void*)(*cset)["values"][csys_axes[di]].element_ptr(point_gid));
603 temp2.set_external(data_dtype, data.element_ptr(di));
604 temp1.to_data_type(data_dtype.id(), temp2);
605 }
606 }
607
608
609 //---------------------------------------------------------------------------//
610 index_t
get_length(index_t dim) const611 TopologyMetadata::get_length(index_t dim) const
612 {
613 // NOTE: The default version of 'get_length' gets the total length of all
614 // unique entities in the topology. The parameterized version fetches the
615 // length for just that parameter's dimension.
616
617 index_t start_dim = (dim >= 0) ? dim : 0;
618 index_t end_dim = (dim >= 0) ? dim : topo_shape.dim;
619
620 index_t topo_length = 0;
621 for(index_t di = start_dim; di <= end_dim; di++)
622 {
623 topo_length += topology::length(dim_topos[di]);
624 }
625
626 return topo_length;
627 }
628
629
630 //---------------------------------------------------------------------------//
631 index_t
get_embed_length(index_t entity_dim,index_t embed_dim) const632 TopologyMetadata::get_embed_length(index_t entity_dim, index_t embed_dim) const
633 {
634 // NOTE: The default version of 'get_embed_length' gets the total number of
635 // embeddings for each entity at the top level to the embedding level. The
636 // parameterized version just fetches the number of embeddings for one
637 // specific entity at the top level.
638
639 std::vector<index_t> entity_index_bag;
640 std::vector<index_t> entity_dim_bag;
641 for(index_t ei = 0; ei < this->get_length(entity_dim); ei++)
642 {
643 entity_index_bag.push_back(ei);
644 entity_dim_bag.push_back(entity_dim);
645 }
646
647 std::set<index_t> embed_set;
648 index_t embed_length = 0;
649 while(!entity_index_bag.empty())
650 {
651 index_t entity_index = entity_index_bag.back();
652 entity_index_bag.pop_back();
653 index_t entity_dim_back = entity_dim_bag.back();
654 entity_dim_bag.pop_back();
655
656 if(entity_dim_back == embed_dim)
657 {
658 if(embed_set.find(entity_index) == embed_set.end())
659 {
660 embed_length++;
661 }
662 embed_set.insert(entity_index);
663 }
664 else
665 {
666 const std::vector<index_t> &embed_ids = get_entity_assocs(
667 TopologyMetadata::LOCAL, entity_index, entity_dim_back, entity_dim_back - 1);
668 for(index_t ei = 0; ei < (index_t)embed_ids.size(); ei++)
669 {
670 entity_index_bag.push_back(embed_ids[ei]);
671 entity_dim_bag.push_back(entity_dim_back - 1);
672 }
673 }
674 }
675
676 return embed_length;
677 }
678
679
680 //---------------------------------------------------------------------------//
681 std::string
to_json() const682 TopologyMetadata::to_json() const
683 {
684 Node mesh;
685
686 Node &mesh_coords = mesh["coordsets"][(*topo)["coordset"].as_string()];
687 mesh_coords.set_external(*cset);
688
689 Node &mesh_topos = mesh["topologies"];
690 for(index_t di = 0; di <= topo_shape.dim; di++)
691 {
692 std::ostringstream oss;
693 oss << "d" << di;
694 mesh_topos[oss.str()].set_external(dim_topos[di]);
695 }
696
697 return mesh.to_json();
698 }
699
700 //-----------------------------------------------------------------------------
701 /// blueprint mesh utility query functions
702 //-----------------------------------------------------------------------------
703
704 //-----------------------------------------------------------------------------
link_nodes(const Node & lhs,const Node & rhs)705 Node link_nodes(const Node &lhs, const Node &rhs)
706 {
707 Node linker;
708 linker.append().set_external(lhs);
709 linker.append().set_external(rhs);
710 return linker;
711 }
712
713
714 //-----------------------------------------------------------------------------
715 DataType
find_widest_dtype(const Node & node,const std::vector<DataType> & default_dtypes)716 find_widest_dtype(const Node &node, const std::vector<DataType> &default_dtypes)
717 {
718 DataType widest_dtype(default_dtypes[0].id(), 0, 0, 0, 0, default_dtypes[0].endianness());
719
720 std::vector<const Node*> node_bag(1, &node);
721 while(!node_bag.empty())
722 {
723 const Node *curr_node = node_bag.back(); node_bag.pop_back();
724 const DataType curr_dtype = curr_node->dtype();
725 if( curr_dtype.is_list() || curr_dtype.is_object() )
726 {
727 NodeConstIterator curr_node_it = curr_node->children();
728 while(curr_node_it.has_next())
729 {
730 node_bag.push_back(&curr_node_it.next());
731 }
732 }
733 else
734 {
735 for(index_t ti = 0; ti < (index_t)default_dtypes.size(); ti++)
736 {
737 const DataType &valid_dtype = default_dtypes[ti];
738 bool is_valid_dtype =
739 (curr_dtype.is_floating_point() && valid_dtype.is_floating_point()) ||
740 (curr_dtype.is_signed_integer() && valid_dtype.is_signed_integer()) ||
741 (curr_dtype.is_unsigned_integer() && valid_dtype.is_unsigned_integer()) ||
742 (curr_dtype.is_string() && valid_dtype.is_string());
743 if(is_valid_dtype && (widest_dtype.element_bytes() < curr_dtype.element_bytes()))
744 {
745 widest_dtype.set(DataType(curr_dtype.id(), 1));
746 }
747 }
748 }
749 }
750
751 bool no_type_found = widest_dtype.element_bytes() == 0;
752 return no_type_found ? default_dtypes[0] : widest_dtype;
753 }
754
755
756 //-----------------------------------------------------------------------------
757 DataType
find_widest_dtype(const Node & node,const DataType & default_dtype)758 find_widest_dtype(const Node &node, const DataType &default_dtype)
759 {
760 return find_widest_dtype(node, std::vector<DataType>(1, default_dtype));
761 }
762
763
764 //-----------------------------------------------------------------------------
765 const Node *
find_reference_node(const Node & node,const std::string & ref_key)766 find_reference_node(const Node &node, const std::string &ref_key)
767 {
768 const Node *res = nullptr;
769
770 // NOTE: This segment of code is necessary to transform "topology" into
771 // "topologies" while keeping all other dependency names (e.g. "coordset")
772 // simply plural by just appending an "s" character.
773 const std::string ref_section = (ref_key[ref_key.length()-1] != 'y') ?
774 ref_key + "s" : ref_key.substr(0, ref_key.length()-1) + "ies";
775
776 if(node.has_child(ref_key))
777 {
778 const std::string &ref_value = node.fetch(ref_key).as_string();
779
780 const Node *traverse_node = node.parent();
781 while(traverse_node != NULL)
782 {
783 if(traverse_node->has_child(ref_section))
784 {
785 const Node &ref_parent = traverse_node->fetch(ref_section);
786 if(ref_parent.has_child(ref_value))
787 {
788 res = &ref_parent[ref_value];
789 }
790 break;
791 }
792 traverse_node = traverse_node->parent();
793 }
794 }
795
796 return res;
797 }
798
799
800 //-----------------------------------------------------------------------------
801 // NOTE: 'node' can be any subtree of a Blueprint-compliant mesh
802 index_t
find_domain_id(const Node & node)803 find_domain_id(const Node &node)
804 {
805 index_t domain_id = -1;
806
807 Node info;
808 const Node *curr_node = &node;
809 while(curr_node != NULL && domain_id == -1)
810 {
811 if(blueprint::mesh::verify(*curr_node, info))
812 {
813 const std::vector<const Node *> domains = blueprint::mesh::domains(*curr_node);
814 const Node &domain = *domains.front();
815 if(domain.has_path("state/domain_id"))
816 {
817 domain_id = domain["state/domain_id"].to_index_t();
818 }
819 }
820
821 curr_node = curr_node->parent();
822 }
823
824 return domain_id;
825 }
826
827 //-----------------------------------------------------------------------------
828 // -- begin conduit::blueprint::mesh::utils::connectivity --
829 //-----------------------------------------------------------------------------
830
831 //-----------------------------------------------------------------------------
832 void
make_element_2d(std::vector<index_t> & elem,index_t element,index_t iwidth)833 connectivity::make_element_2d(std::vector<index_t>& elem,
834 index_t element,
835 index_t iwidth)
836 {
837 index_t ilo = element % iwidth;
838 index_t jlo = element / iwidth;
839 index_t ihi = ilo + 1;
840 index_t jhi = jlo + 1;
841
842 index_t ilo_jlo = (iwidth+1)*jlo + ilo;
843 index_t ihi_jlo = (iwidth+1)*jlo + ihi;
844 index_t ihi_jhi = (iwidth+1)*jhi + ihi;
845 index_t ilo_jhi = (iwidth+1)*jhi + ilo;
846
847 elem.push_back(ilo_jlo);
848 elem.push_back(ihi_jlo);
849 elem.push_back(ihi_jhi);
850 elem.push_back(ilo_jhi);
851 }
852
853
854 void
make_element_3d(ElemType & connect,index_t element,index_t iwidth,index_t jwidth,index_t kwidth,SubelemMap & faces)855 connectivity::make_element_3d(ElemType& connect,
856 index_t element,
857 index_t iwidth,
858 index_t jwidth,
859 index_t kwidth,
860 SubelemMap& faces)
861 {
862 index_t ilo = element % iwidth;
863 index_t jlo = (element / iwidth) % jwidth;
864 index_t klo = element / (iwidth*jwidth);
865 index_t ihi = ilo + 1;
866 index_t jhi = jlo + 1;
867 index_t khi = klo + 1;
868
869 index_t jlo_offset = (iwidth+1)*jlo;
870 index_t jhi_offset = (iwidth+1)*jhi;
871 index_t klo_offset = (iwidth+1)*(jwidth+1)*klo;
872 index_t khi_offset = (iwidth+1)*(jwidth+1)*khi;
873
874
875 index_t iface_start = 0;
876 index_t jface_start = (iwidth+1)*jwidth*kwidth;
877 index_t kface_start = jface_start + iwidth*(jwidth+1)*kwidth;
878
879 //ifaces
880 {
881 index_t j_offset = jlo_offset;
882 index_t k_offset = (iwidth+1)*jwidth*klo;
883
884 index_t lo_face = iface_start + ilo + j_offset + k_offset;
885 index_t hi_face = iface_start + ihi + j_offset + k_offset;
886
887
888 //ilo face
889 if (faces.find(lo_face) == faces.end())
890 {
891 auto& ilo_face = faces[lo_face];
892 ilo_face.push_back(ilo + jlo_offset + klo_offset);
893 ilo_face.push_back(ilo + jhi_offset + klo_offset);
894 ilo_face.push_back(ilo + jhi_offset + khi_offset);
895 ilo_face.push_back(ilo + jlo_offset + khi_offset);
896 }
897 //ihi face
898 if (faces.find(hi_face) == faces.end())
899 {
900 auto& ihi_face = faces[hi_face];
901 ihi_face.push_back(ihi + jlo_offset + klo_offset);
902 ihi_face.push_back(ihi + jhi_offset + klo_offset);
903 ihi_face.push_back(ihi + jhi_offset + khi_offset);
904 ihi_face.push_back(ihi + jlo_offset + khi_offset);
905 }
906 connect.push_back(lo_face);
907 connect.push_back(hi_face);
908 }
909 //jfaces
910 {
911 index_t i_offset = ilo;
912 index_t jlo_face_offset = iwidth*jlo;
913 index_t jhi_face_offset = iwidth*jhi;
914 index_t k_offset = iwidth*(jwidth+1)*klo;
915
916 index_t lo_face = jface_start + i_offset + jlo_face_offset + k_offset;
917 index_t hi_face = jface_start + i_offset + jhi_face_offset + k_offset;
918
919 //jlo face
920 if (faces.find(lo_face) == faces.end())
921 {
922 auto& jlo_face = faces[lo_face];
923 jlo_face.push_back(ilo + jlo_offset + klo_offset);
924 jlo_face.push_back(ihi + jlo_offset + klo_offset);
925 jlo_face.push_back(ihi + jlo_offset + khi_offset);
926 jlo_face.push_back(ilo + jlo_offset + khi_offset);
927 }
928 //jhi face
929 if (faces.find(hi_face) == faces.end())
930 {
931 auto& jhi_face = faces[hi_face];
932 jhi_face.push_back(ilo + jhi_offset + klo_offset);
933 jhi_face.push_back(ihi + jhi_offset + klo_offset);
934 jhi_face.push_back(ihi + jhi_offset + khi_offset);
935 jhi_face.push_back(ilo + jhi_offset + khi_offset);
936 }
937 connect.push_back(lo_face);
938 connect.push_back(hi_face);
939 }
940 //kfaces
941 {
942 index_t i_offset = ilo;
943 index_t j_offset = iwidth*jlo;
944 index_t klo_face_offset = iwidth*jwidth*klo;
945 index_t khi_face_offset = iwidth*jwidth*khi;
946
947 index_t lo_face = kface_start + i_offset + j_offset + klo_face_offset;
948 index_t hi_face = kface_start + i_offset + j_offset + khi_face_offset;
949
950 //klo face
951 if (faces.find(lo_face) == faces.end())
952 {
953 auto& klo_face = faces[lo_face];
954 klo_face.push_back(ilo + jlo_offset + klo_offset);
955 klo_face.push_back(ihi + jlo_offset + klo_offset);
956 klo_face.push_back(ihi + jhi_offset + klo_offset);
957 klo_face.push_back(ilo + jhi_offset + klo_offset);
958 }
959 //khi face
960 if (faces.find(hi_face) == faces.end())
961 {
962 auto& khi_face = faces[hi_face];
963 khi_face.push_back(ilo + jlo_offset + khi_offset);
964 khi_face.push_back(ihi + jlo_offset + khi_offset);
965 khi_face.push_back(ihi + jhi_offset + khi_offset);
966 khi_face.push_back(ilo + jhi_offset + khi_offset);
967 }
968 connect.push_back(/*kface_start +*/lo_face);
969 connect.push_back(/*kface_start +*/hi_face);
970 }
971 }
972
973
974
975
976 //-----------------------------------------------------------------------------
977 void
create_elements_2d(const Node & ref_win,index_t i_lo,index_t j_lo,index_t iwidth,std::map<index_t,std::vector<index_t>> & elems)978 connectivity::create_elements_2d(const Node& ref_win,
979 index_t i_lo,
980 index_t j_lo,
981 index_t iwidth,
982 std::map<index_t, std::vector<index_t> >& elems)
983 {
984 index_t origin_iref = ref_win["origin/i"].to_index_t();
985 index_t origin_jref = ref_win["origin/j"].to_index_t();
986
987 index_t ref_size_i = ref_win["dims/i"].to_index_t();
988 index_t ref_size_j = ref_win["dims/j"].to_index_t();
989
990 if (ref_size_i == 1)
991 {
992 index_t jstart = origin_jref - j_lo;
993 index_t jend = origin_jref - j_lo + ref_size_j - 1;
994 if (origin_iref == i_lo)
995 {
996 for (index_t jidx = jstart; jidx < jend; ++jidx)
997 {
998 index_t offset = jidx * iwidth;
999 auto& elem_conn = elems[offset];
1000 if (elem_conn.empty())
1001 {
1002 connectivity::make_element_2d(elem_conn,
1003 offset,
1004 iwidth);
1005 }
1006 }
1007 }
1008 else
1009 {
1010 for (index_t jidx = jstart; jidx < jend; ++jidx)
1011 {
1012 index_t offset = jidx * iwidth + (origin_iref - i_lo - 1);
1013 auto& elem_conn = elems[offset];
1014 if (elem_conn.empty())
1015 {
1016 connectivity::make_element_2d(elem_conn,
1017 offset,
1018 iwidth);
1019 }
1020 }
1021 }
1022 }
1023 else if (ref_size_j == 1)
1024 {
1025 index_t istart = origin_iref - i_lo;
1026 index_t iend = origin_iref - i_lo + ref_size_i - 1;
1027 if (origin_jref == j_lo)
1028 {
1029 for (index_t iidx = istart; iidx < iend; ++iidx)
1030 {
1031 auto& elem_conn = elems[iidx];
1032 if (elem_conn.empty())
1033 {
1034 connectivity::make_element_2d(elem_conn,
1035 iidx,
1036 iwidth);
1037 }
1038 }
1039 }
1040 else
1041 {
1042 for (index_t iidx = istart; iidx < iend; ++iidx)
1043 {
1044 index_t offset = iidx + ((origin_jref - j_lo - 1) * iwidth);
1045 auto& elem_conn = elems[offset];
1046 if (elem_conn.empty())
1047 {
1048 connectivity::make_element_2d(elem_conn,
1049 offset,
1050 iwidth);
1051 }
1052 }
1053 }
1054 }
1055
1056 index_t istart = origin_iref - i_lo;
1057 index_t jstart = origin_jref - j_lo;
1058 index_t iend = istart + ref_size_i - 1;
1059 index_t jend = jstart + ref_size_j - 1;
1060
1061 if (ref_size_i == 1)
1062 {
1063 if (origin_iref != i_lo)
1064 {
1065 --istart;
1066 }
1067 iend = istart + 1;
1068 }
1069 if (ref_size_j == 1)
1070 {
1071 if (origin_jref != j_lo)
1072 {
1073 --jstart;
1074 }
1075 jend = jstart + 1;
1076 }
1077
1078 for (index_t jidx = jstart; jidx < jend; ++jidx)
1079 {
1080 index_t joffset = jidx * iwidth;
1081 for (index_t iidx = istart; iidx < iend; ++iidx)
1082 {
1083 index_t offset = joffset + iidx;
1084 auto& elem_conn = elems[offset];
1085 if (elem_conn.empty())
1086 {
1087 connectivity::make_element_2d(elem_conn,
1088 offset,
1089 iwidth);
1090 }
1091 }
1092 }
1093 }
1094
1095
1096 //-----------------------------------------------------------------------------
1097 void
create_elements_3d(const Node & ref_win,index_t i_lo,index_t j_lo,index_t k_lo,index_t iwidth,index_t jwidth,index_t kwidth,std::map<index_t,ElemType> & elems,SubelemMap & faces)1098 connectivity::create_elements_3d(const Node& ref_win,
1099 index_t i_lo,
1100 index_t j_lo,
1101 index_t k_lo,
1102 index_t iwidth,
1103 index_t jwidth,
1104 index_t kwidth,
1105 std::map<index_t, ElemType>& elems,
1106 SubelemMap& faces)
1107 {
1108 index_t origin_iref = ref_win["origin/i"].to_index_t();
1109 index_t origin_jref = ref_win["origin/j"].to_index_t();
1110 index_t origin_kref = ref_win["origin/k"].to_index_t();
1111
1112 index_t ref_size_i = ref_win["dims/i"].to_index_t();
1113 index_t ref_size_j = ref_win["dims/j"].to_index_t();
1114 index_t ref_size_k = ref_win["dims/k"].to_index_t();
1115
1116 index_t istart = origin_iref - i_lo;
1117 index_t jstart = origin_jref - j_lo;
1118 index_t kstart = origin_kref - k_lo;
1119 index_t iend = istart + ref_size_i - 1;
1120 index_t jend = jstart + ref_size_j - 1;
1121 index_t kend = kstart + ref_size_k - 1;
1122
1123 if (ref_size_i == 1)
1124 {
1125 iend = istart + 1;
1126 }
1127 if (ref_size_j == 1)
1128 {
1129 jend = jstart + 1;
1130 }
1131 if (ref_size_k == 1)
1132 {
1133 kend = kstart + 1;
1134 }
1135
1136 for (index_t kidx = kstart; kidx < kend; ++kidx)
1137 {
1138 index_t koffset = kidx * iwidth * jwidth;
1139 for (index_t jidx = jstart; jidx < jend; ++jidx)
1140 {
1141 index_t joffset = jidx * iwidth;
1142 for (index_t iidx = istart; iidx < iend; ++iidx)
1143 {
1144 index_t offset = koffset + joffset + iidx;
1145 auto& elem_conn = elems[offset];
1146 if (elem_conn.empty())
1147 {
1148 connectivity::make_element_3d(elem_conn,
1149 offset,
1150 iwidth,
1151 jwidth,
1152 kwidth,
1153 faces);
1154 }
1155 }
1156 }
1157 }
1158 }
1159
1160 void
connect_elements_3d(const Node & ref_win,index_t i_lo,index_t j_lo,index_t k_lo,index_t iwidth,index_t jwidth,index_t & new_vertex,std::map<index_t,ElemType> & elems)1161 connectivity::connect_elements_3d(const Node& ref_win,
1162 index_t i_lo,
1163 index_t j_lo,
1164 index_t k_lo,
1165 index_t iwidth,
1166 index_t jwidth,
1167 index_t& new_vertex,
1168 std::map<index_t, ElemType>& elems)
1169 {
1170 index_t origin_iref = ref_win["origin/i"].to_index_t();
1171 index_t origin_jref = ref_win["origin/j"].to_index_t();
1172 index_t origin_kref = ref_win["origin/k"].to_index_t();
1173
1174 index_t ref_size_i = ref_win["dims/i"].to_index_t();
1175 index_t ref_size_j = ref_win["dims/j"].to_index_t();
1176 index_t ref_size_k = ref_win["dims/k"].to_index_t();
1177
1178 index_t kstart = origin_kref - k_lo;
1179 index_t kend = origin_kref - k_lo + ref_size_k - 1;
1180 if (kstart == kend) kend = kstart + 1;
1181 index_t jstart = origin_jref - j_lo;
1182 index_t jend = origin_jref - j_lo + ref_size_j - 1;
1183 if (jstart == jend) jend = jstart + 1;
1184 index_t istart = origin_iref - i_lo;
1185 index_t iend = origin_iref - i_lo + ref_size_i - 1;
1186 if (istart == iend) iend = istart + 1;
1187
1188 for (index_t kidx = kstart; kidx < kend; ++kidx)
1189 {
1190 for (index_t jidx = jstart; jidx < jend; ++jidx)
1191 {
1192 for (index_t iidx = istart; iidx < iend; ++iidx)
1193 {
1194 index_t offset = kidx*iwidth*jwidth + jidx*iwidth + iidx;
1195 auto& elem_conn = elems[offset];
1196 elem_conn.push_back(new_vertex);
1197 ++new_vertex;
1198 }
1199 }
1200 }
1201 }
1202
1203 void
connect_elements_2d(const Node & ref_win,index_t i_lo,index_t j_lo,index_t iwidth,index_t ratio,index_t & new_vertex,std::map<index_t,std::vector<index_t>> & elems)1204 connectivity::connect_elements_2d(const Node& ref_win,
1205 index_t i_lo,
1206 index_t j_lo,
1207 index_t iwidth,
1208 index_t ratio,
1209 index_t& new_vertex,
1210 std::map<index_t, std::vector<index_t> >& elems)
1211 {
1212 index_t origin_iref = ref_win["origin/i"].to_index_t();
1213 index_t origin_jref = ref_win["origin/j"].to_index_t();
1214
1215 index_t ref_size_i = ref_win["dims/i"].to_index_t();
1216 index_t ref_size_j = ref_win["dims/j"].to_index_t();
1217
1218 if (ref_size_i == 1)
1219 {
1220 index_t jstart = origin_jref - j_lo;
1221 index_t jend = origin_jref - j_lo + ref_size_j - 1;
1222 if (origin_iref == i_lo)
1223 {
1224 for (index_t jidx = jstart; jidx < jend; ++jidx)
1225 {
1226 index_t offset = jidx * (iwidth);
1227 auto& elem_conn = elems[offset];
1228 if (ratio > 1)
1229 {
1230 for (index_t nr = ratio-1; nr > 0; --nr)
1231 {
1232 elem_conn.push_back(new_vertex+nr-1);
1233 }
1234 new_vertex += ratio - 1;
1235 }
1236 }
1237 }
1238 else
1239 {
1240 for (index_t jidx = jstart; jidx < jend; ++jidx)
1241 {
1242 index_t offset = jidx * iwidth + (origin_iref - i_lo - 1);
1243 auto& elem_conn = elems[offset];
1244 if (ratio > 1)
1245 {
1246 size_t new_size = elem_conn.size() + ratio - 1;
1247 elem_conn.resize(new_size);
1248 index_t corner = 1;
1249 if (elem_conn[1] - elem_conn[0] != 1)
1250 {
1251 index_t ioff = offset % iwidth;
1252 index_t joff = offset / iwidth;
1253 index_t target = (iwidth+1)*joff + ioff + 1;
1254 for (index_t nr = 1; nr < 1+ratio; ++nr)
1255 {
1256 if (elem_conn[nr] == target)
1257 {
1258 corner = nr;
1259 break;
1260 }
1261 }
1262 }
1263 for (index_t nr = new_size-1; nr > corner+ratio-1; --nr)
1264 {
1265 elem_conn[nr] = elem_conn[nr-ratio+1];
1266 }
1267 for (index_t nr = corner+1; nr < corner+ratio; ++nr)
1268 {
1269 elem_conn[nr] = new_vertex;
1270 ++new_vertex;
1271 }
1272 }
1273 }
1274 }
1275 }
1276 else if (ref_size_j == 1)
1277 {
1278 index_t istart = origin_iref - i_lo;
1279 index_t iend = origin_iref - i_lo + ref_size_i - 1;
1280 if (origin_jref == j_lo)
1281 {
1282 for (index_t iidx = istart; iidx < iend; ++iidx)
1283 {
1284 auto& elem_conn = elems[iidx];
1285 if (ratio > 1)
1286 {
1287 size_t new_size = elem_conn.size() + ratio - 1;
1288 elem_conn.resize(new_size);
1289 for (index_t nr = new_size-1; nr > 1; --nr)
1290 {
1291 elem_conn[nr] = elem_conn[nr-ratio+1];
1292 }
1293 for (index_t nr = 1; nr < ratio; ++nr)
1294 {
1295 elem_conn[nr] = (new_vertex+nr-1);
1296 }
1297 new_vertex += ratio - 1;
1298 }
1299 }
1300 }
1301 else
1302 {
1303 for (index_t iidx = istart; iidx < iend; ++iidx)
1304 {
1305 index_t offset = iidx + ((origin_jref - j_lo - 1) * iwidth);
1306 auto& elem_conn = elems[offset];
1307 if (ratio > 1)
1308 {
1309 size_t old_size = elem_conn.size();
1310 size_t new_size = old_size + ratio - 1;
1311 elem_conn.resize(new_size);
1312 index_t corner = 2;
1313 if (old_size != 4)
1314 {
1315 index_t ioff = offset % iwidth;
1316 index_t joff = offset / iwidth;
1317 index_t target = (iwidth+1)*(joff+1) + ioff + 1;
1318 for (index_t nr = 3; nr < 3+ratio; ++nr)
1319 {
1320 if (elem_conn[nr] == target) {
1321 corner = nr;
1322 break;
1323 }
1324 }
1325 }
1326 for (index_t nr = new_size-1; nr > corner+ratio-1; --nr)
1327 {
1328 elem_conn[nr] = elem_conn[nr-ratio+1];
1329 }
1330 for (index_t nr = corner+ratio-1; nr > corner; --nr) {
1331 elem_conn[nr] = new_vertex;
1332 ++new_vertex;
1333 }
1334 }
1335 }
1336 }
1337 }
1338 }
1339
1340 //-----------------------------------------------------------------------------
1341 // -- end conduit::blueprint::mesh::utils::connectivity --
1342 //-----------------------------------------------------------------------------
1343
1344 //-----------------------------------------------------------------------------
1345 // -- begin conduit::blueprint::mesh::utils::coordset --
1346 //-----------------------------------------------------------------------------
1347
1348 //-----------------------------------------------------------------------------
1349 std::pair<std::string, std::vector<std::string>>
get_coordset_info(const Node & n)1350 get_coordset_info(const Node &n)
1351 {
1352 std::pair<std::string, std::vector<std::string>> info;
1353 std::string &cset_coordsys = info.first;
1354 std::vector<std::string> &cset_axes = info.second;
1355
1356 std::string coords_path = "";
1357 if(n["type"].as_string() == "uniform")
1358 {
1359 if(n.has_child("origin"))
1360 {
1361 coords_path = "origin";
1362 }
1363 else if(n.has_child("spacing"))
1364 {
1365 coords_path = "spacing";
1366 }
1367 else
1368 {
1369 coords_path = "dims";
1370 }
1371 }
1372 else // if(n.has_child("values"))
1373 {
1374 coords_path = "values";
1375 }
1376
1377 Node axis_names;
1378 const Node &cset_coords = n[coords_path];
1379 NodeConstIterator itr = cset_coords.children();
1380 while(itr.has_next())
1381 {
1382 itr.next();
1383 const std::string axis_name = itr.name();
1384
1385 if(axis_name[0] == 'd' && axis_name.size() > 1)
1386 {
1387 axis_names[axis_name.substr(1, axis_name.length())];
1388 }
1389 else
1390 {
1391 axis_names[axis_name];
1392 }
1393 }
1394
1395 std::vector<std::string> cset_base_axes;
1396 cset_coordsys = "unknown";
1397 if(axis_names.has_child("theta") || axis_names.has_child("phi"))
1398 {
1399 cset_coordsys = "spherical";
1400 cset_base_axes = SPHERICAL_AXES;
1401 }
1402 else if(axis_names.has_child("r")) // rz, or r w/o theta, phi
1403 {
1404 cset_coordsys = "cylindrical";
1405 cset_base_axes = CYLINDRICAL_AXES;
1406 }
1407 else if(axis_names.has_child("x") || axis_names.has_child("y") || axis_names.has_child("z"))
1408 {
1409 cset_coordsys = "cartesian";
1410 cset_base_axes = CARTESIAN_AXES;
1411 }
1412 else if(axis_names.has_child("i") || axis_names.has_child("j") || axis_names.has_child("k"))
1413 {
1414 cset_coordsys = "logical";
1415 cset_base_axes = LOGICAL_AXES;
1416 }
1417
1418 // intersect 'cset_base_axes' and 'axis_names.child_names()'
1419 for(const std::string &base_axis : cset_base_axes)
1420 {
1421 for(const std::string &cset_axis : axis_names.child_names())
1422 {
1423 if(base_axis == cset_axis)
1424 {
1425 cset_axes.push_back(cset_axis);
1426 break;
1427 }
1428 }
1429 }
1430
1431 // TODO(JRC): The following is an alterate (though potentially more error-prone) solution:
1432 // std::vector<std::string>(cset_axes.begin(), cset_axes.begin() + cset_coords.number_of_children());
1433
1434 return info;
1435 }
1436
1437
1438 //-----------------------------------------------------------------------------
1439 index_t
dims(const Node & n)1440 coordset::dims(const Node &n)
1441 {
1442 const std::vector<std::string> csys_axes = coordset::axes(n);
1443 return (index_t)csys_axes.size();
1444 }
1445
1446
1447 //-----------------------------------------------------------------------------
1448 index_t
length(const Node & n)1449 coordset::length(const Node &n)
1450 {
1451 index_t coordset_length = 1;
1452
1453 const std::string csys_type = n["type"].as_string();
1454 const std::vector<std::string> csys_axes = coordset::axes(n);
1455 for(index_t i = 0; i < (index_t)csys_axes.size(); i++)
1456 {
1457 if(csys_type == "uniform")
1458 {
1459 coordset_length *= n["dims"][LOGICAL_AXES[i]].to_index_t();
1460 }
1461 else if(csys_type == "rectilinear")
1462 {
1463 coordset_length *= n["values"][csys_axes[i]].dtype().number_of_elements();
1464 }
1465 else // if(csys_type == "explicit")
1466 {
1467 coordset_length = n["values"][csys_axes[i]].dtype().number_of_elements();
1468 }
1469 }
1470
1471 return coordset_length;
1472 }
1473
1474
1475 //-----------------------------------------------------------------------------
1476 std::vector<std::string>
axes(const Node & n)1477 coordset::axes(const Node &n)
1478 {
1479 return get_coordset_info(n).second;
1480 }
1481
1482
1483 //-----------------------------------------------------------------------------
1484 std::string
coordsys(const Node & n)1485 coordset::coordsys(const Node &n)
1486 {
1487 return get_coordset_info(n).first;
1488 }
1489
1490
1491 //-----------------------------------------------------------------------------
1492 std::vector<float64>
coords(const Node & n,const index_t i)1493 coordset::_explicit::coords(const Node &n, const index_t i)
1494 {
1495 std::vector<float64> cvals;
1496
1497 Node temp;
1498 for(const std::string &axis : coordset::axes(n))
1499 {
1500 const Node &axis_node = n["values"][axis];
1501 temp.set_external(DataType(axis_node.dtype().id(), 1),
1502 (void*)axis_node.element_ptr(i));
1503 cvals.push_back(temp.to_float64());
1504 }
1505
1506 return std::vector<float64>(std::move(cvals));
1507 }
1508
1509 //-----------------------------------------------------------------------------
1510 void
logical_dims(const conduit::Node & n,index_t * d,index_t maxdims)1511 coordset::logical_dims(const conduit::Node &n, index_t *d, index_t maxdims)
1512 {
1513 for(index_t i = 0; i < maxdims; i++)
1514 {
1515 d[i] = 1;
1516 }
1517
1518 auto info = get_coordset_info(n);
1519 const std::string cset_type = n["type"].as_string();
1520 const std::vector<std::string> &cset_axes = info.second;
1521 if(cset_type == "uniform" || cset_type == "rectilinear")
1522 {
1523 const index_t dim = ((index_t)cset_axes.size() > maxdims) ? maxdims
1524 : (index_t)cset_axes.size();
1525 for(index_t i = 0; i < dim; i++)
1526 {
1527 if(cset_type == "uniform")
1528 {
1529 d[i] = n["dims"][LOGICAL_AXES[i]].to_index_t();
1530 }
1531 else // if(cset_type == "rectilinear")
1532 {
1533 d[i] = n["values"][cset_axes[i]].dtype().number_of_elements();
1534 }
1535 }
1536 }
1537 else // if(cset_type == "explicit")
1538 {
1539 d[0] = n["values"][cset_axes[0]].dtype().number_of_elements();
1540 }
1541 }
1542
1543 //-----------------------------------------------------------------------------
1544 template<typename data_array>
1545 static void
typed_minmax(const data_array & da,float64 & out_min,float64 & out_max)1546 typed_minmax(const data_array &da, float64 &out_min, float64 &out_max)
1547 {
1548 // Figure out what primitive type we are dealing with
1549 using T = decltype(std::declval<data_array>().min());
1550 const index_t nelem = da.number_of_elements();
1551 T min = std::numeric_limits<T>::max();
1552 T max = std::numeric_limits<T>::lowest();
1553 for(index_t i = 0; i < nelem; i++)
1554 {
1555 min = std::min(min, da[i]);
1556 max = std::max(max, da[i]);
1557 }
1558 out_min = (float64)min;
1559 out_max = (float64)max;
1560 }
1561
1562 //-----------------------------------------------------------------------------
1563 std::vector<float64>
extents(const Node & n)1564 coordset::extents(const Node &n)
1565 {
1566 std::vector<float64> cset_extents;
1567 const std::string csys_type = n["type"].as_string();
1568 const std::vector<std::string> csys_axes = coordset::axes(n);
1569 const index_t naxes = (index_t)csys_axes.size();
1570 cset_extents.reserve(naxes*2);
1571 for(index_t i = 0; i < naxes; i++)
1572 {
1573 float64 min, max;
1574 if(csys_type == "uniform")
1575 {
1576 index_t origin = 0;
1577 float64 spacing = 1.0;
1578 index_t dim = n["dims"][LOGICAL_AXES[i]].to_index_t();
1579 if(n.has_child("origin")
1580 && n["origin"].has_child(csys_axes[i]))
1581 {
1582 origin = n["origin"][csys_axes[i]].to_index_t();
1583 }
1584 if(n.has_child("spacing")
1585 && n["spacing"].has_child("d"+csys_axes[i]))
1586 {
1587 spacing = n["spacing"]["d" + csys_axes[i]].to_float64();
1588 }
1589 min = (float64)origin;
1590 max = (float64)origin + (spacing * ((float64)dim - 1.));
1591 if(spacing < 0.)
1592 {
1593 std::swap(min, max);
1594 }
1595 }
1596 else // csys_type == "rectilinear" || csys_type == "explicit"
1597 {
1598 const auto &axis = n["values"][csys_axes[i]];
1599 const auto id = axis.dtype().id();
1600 switch(id)
1601 {
1602 case conduit::DataType::INT8_ID:
1603 typed_minmax(axis.as_int8_array(), min, max);
1604 break;
1605 case conduit::DataType::INT16_ID:
1606 typed_minmax(axis.as_int16_array(), min, max);
1607 break;
1608 case conduit::DataType::INT32_ID:
1609 typed_minmax(axis.as_int32_array(), min, max);
1610 break;
1611 case conduit::DataType::INT64_ID:
1612 typed_minmax(axis.as_int64_array(), min, max);
1613 break;
1614 case conduit::DataType::UINT8_ID:
1615 typed_minmax(axis.as_uint8_array(), min, max);
1616 break;
1617 case conduit::DataType::UINT16_ID:
1618 typed_minmax(axis.as_uint16_array(), min, max);
1619 break;
1620 case conduit::DataType::UINT32_ID:
1621 typed_minmax(axis.as_uint32_array(), min, max);
1622 break;
1623 case conduit::DataType::UINT64_ID:
1624 typed_minmax(axis.as_uint64_array(), min, max);
1625 break;
1626 case conduit::DataType::FLOAT32_ID:
1627 typed_minmax(axis.as_float32_array(), min, max);
1628 break;
1629 case conduit::DataType::FLOAT64_ID:
1630 typed_minmax(axis.as_float64_array(), min, max);
1631 break;
1632 }
1633 }
1634 cset_extents.push_back(min);
1635 cset_extents.push_back(max);
1636 }
1637 return cset_extents;
1638 }
1639
1640 //-----------------------------------------------------------------------------
1641 // -- begin conduit::blueprint::mesh::utils::coordset::uniform --
1642 //-----------------------------------------------------------------------------
1643 //-----------------------------------------------------------------------------
1644 std::vector<double>
spacing(const Node & n)1645 coordset::uniform::spacing(const Node &n)
1646 {
1647 auto info = get_coordset_info(n);
1648 const auto &cset_axes = info.second;
1649 std::vector<double> retval(cset_axes.size(), 1);
1650 if(n.has_child("spacing"))
1651 {
1652 const Node &n_spacing = n["spacing"];
1653 for(index_t i = 0; i < (index_t)cset_axes.size(); i++)
1654 {
1655 const std::string child_name = "d"+cset_axes[i];
1656 if(n_spacing.has_child(child_name))
1657 {
1658 retval[i] = n_spacing[child_name].to_double();
1659 }
1660 }
1661 }
1662 return retval;
1663 }
1664
1665 //-----------------------------------------------------------------------------
1666 std::vector<index_t>
origin(const Node & n)1667 coordset::uniform::origin(const Node &n)
1668 {
1669 auto info = get_coordset_info(n);
1670 const auto &cset_axes = info.second;
1671 std::vector<index_t> retval(cset_axes.size(), 0);
1672 if(n.has_child("origin"))
1673 {
1674 const Node &n_spacing = n["origin"];
1675 for(index_t i = 0; i < (index_t)cset_axes.size(); i++)
1676 {
1677 const std::string child_name = cset_axes[i];
1678 if(n_spacing.has_child(child_name))
1679 {
1680 retval[i] = n_spacing[child_name].to_index_t();
1681 }
1682 }
1683 }
1684 return retval;
1685 }
1686
1687 //-----------------------------------------------------------------------------
1688 // -- end conduit::blueprint::mesh::utils::coordset::uniform --
1689 //-----------------------------------------------------------------------------
1690
1691 //-----------------------------------------------------------------------------
1692 // -- end conduit::blueprint::mesh::utils::coordset --
1693 //-----------------------------------------------------------------------------
1694
1695 //-----------------------------------------------------------------------------
1696 // -- begin conduit::blueprint::mesh::utils::topology --
1697 //-----------------------------------------------------------------------------
1698
1699 //-----------------------------------------------------------------------------
1700 index_t
dims(const Node & n)1701 topology::dims(const Node &n)
1702 {
1703 index_t topology_dims = 1;
1704
1705 const std::string type = n["type"].as_string();
1706 if(type != "unstructured")
1707 {
1708 const Node *coordset = find_reference_node(n, "coordset");
1709 topology_dims = coordset::dims(*coordset);
1710 }
1711 else // if(type == "unstructured")
1712 {
1713 ShapeType shape(n);
1714 topology_dims = (index_t)shape.dim;
1715 }
1716
1717 return topology_dims;
1718 }
1719
1720
1721 //-----------------------------------------------------------------------------
1722 void
logical_dims(const Node & n,index_t * d,index_t maxdims)1723 topology::logical_dims(const Node &n, index_t *d, index_t maxdims)
1724 {
1725 for(index_t i = 0; i < maxdims; i++)
1726 d[i] = 1;
1727
1728 const std::string type = n["type"].as_string();
1729 if(type == "uniform" || type == "rectilinear")
1730 {
1731 const Node *coordset = find_reference_node(n, "coordset");
1732 const std::vector<std::string> csys_axes = coordset::axes(*coordset);
1733 for(index_t i = 0; i < (index_t)csys_axes.size(); i++)
1734 {
1735 d[i] = ((type == "uniform") ?
1736 coordset->fetch_existing("dims")[LOGICAL_AXES[i]].to_index_t() :
1737 coordset->fetch_existing("values")[csys_axes[i]].dtype().number_of_elements()) - 1;
1738 }
1739 }
1740 else if(type == "structured")
1741 {
1742 const Node &dims = n["elements/dims"];
1743
1744 for(index_t i = 0; i < (index_t)dims.number_of_children(); i++)
1745 {
1746 d[i] = dims[LOGICAL_AXES[i]].to_index_t();
1747 }
1748 }
1749 else if(type == "points")
1750 {
1751 const Node *coordset = find_reference_node(n, "coordset");
1752 if(coordset)
1753 {
1754 coordset::logical_dims(*coordset, d, maxdims);
1755 }
1756 else
1757 {
1758 CONDUIT_ERROR("Unable to find reference node 'coordset'.");
1759 }
1760 }
1761 else // if(type == "unstructured")
1762 {
1763 // TODO(JRC): This is rather inefficient because the offsets array
1764 // is discarded after this calculation is complete.
1765 Node topo_offsets;
1766 topology::unstructured::generate_offsets(n, topo_offsets);
1767 d[0] = topo_offsets.dtype().number_of_elements();
1768 }
1769 }
1770
1771 //-----------------------------------------------------------------------------
1772 index_t
length(const Node & n)1773 topology::length(const Node &n)
1774 {
1775 index_t d[3]={1,1,1};
1776 logical_dims(n, d, 3);
1777 return d[0] * d[1] * d[2];
1778 }
1779
1780 //-----------------------------------------------------------------------------
1781 void
generate_offsets(Node & n,Node & dest)1782 topology::unstructured::generate_offsets(Node &n,
1783 Node &dest)
1784 {
1785 dest.reset();
1786
1787 // FIXME(JRC): There are weird cases wherein a polyhedral topology can have only
1788 // the 'elements/offsets' defined and not 'subelements/offsets', which isn't currently
1789 // properly handled by this function.
1790
1791 if(n["elements"].has_child("offsets") && !n["elements/offsets"].dtype().is_empty())
1792 {
1793 if(&dest != &n["elements/offsets"])
1794 {
1795 dest.set_external(n["elements/offsets"]);
1796 }
1797 }
1798 else
1799 {
1800 const Node &n_const = n;
1801 Node &offsets = n["elements/offsets"];
1802 blueprint::mesh::utils::topology::unstructured::generate_offsets(n_const, offsets);
1803 if(&dest != &offsets)
1804 {
1805 dest.set_external(offsets);
1806 }
1807 }
1808 }
1809
1810
1811 //-----------------------------------------------------------------------------
1812 void
generate_offsets(const Node & n,Node & dest)1813 topology::unstructured::generate_offsets(const Node &n,
1814 Node &dest)
1815 {
1816 const ShapeType topo_shape(n);
1817 const DataType int_dtype = find_widest_dtype(n, DEFAULT_INT_DTYPES);
1818 std::string key("elements/connectivity"), stream_key("elements/stream");
1819 if(!n.has_path(key))
1820 key = stream_key;
1821 const Node &topo_conn = n[key];
1822
1823 const DataType topo_dtype(topo_conn.dtype().id(), 1, 0, 0,
1824 topo_conn.dtype().element_bytes(), topo_conn.dtype().endianness());
1825
1826 if(n["elements"].has_child("offsets") && !n["elements/offsets"].dtype().is_empty())
1827 {
1828 if(&dest != &n["elements/offsets"])
1829 {
1830 dest.set_external(n["elements/offsets"]);
1831 }
1832 }
1833 else if(n.has_path(stream_key))
1834 {
1835 dest.reset();
1836 // Mixed element types
1837 std::map<int,int> stream_id_npts;
1838 const conduit::Node &n_element_types = n["elements/element_types"];
1839 for(index_t i = 0; i < n_element_types.number_of_children(); i++)
1840 {
1841 const Node &n_et = n_element_types[i];
1842 auto stream_id = n_et["stream_id"].to_int();
1843 std::string shape(n_et["shape"].as_string());
1844 for(size_t j = 0; j < TOPO_SHAPES.size(); j++)
1845 {
1846 if(shape == TOPO_SHAPES[j])
1847 {
1848 stream_id_npts[stream_id] = TOPO_SHAPE_EMBED_COUNTS[j];
1849 break;
1850 }
1851 }
1852 }
1853
1854 const Node &n_stream_ids = n["elements/element_index/stream_ids"];
1855 std::vector<index_t> offsets;
1856 if(n.has_path("elements/element_index/element_counts"))
1857 {
1858 const Node &n_element_counts = n["elements/element_index/element_counts"];
1859
1860 index_t offset = 0, elemid = 0;
1861 for(index_t j = 0; j < n_stream_ids.dtype().number_of_elements(); j++)
1862 {
1863 // Get the j'th elements from n_stream_ids, n_element_counts
1864 const Node n_elem_ct_j(int_dtype,
1865 const_cast<void*>(n_element_counts.element_ptr(j)), true);
1866 const Node n_stream_ids_j(int_dtype,
1867 const_cast<void*>(n_stream_ids.element_ptr(j)), true);
1868 auto ec = static_cast<index_t>(n_elem_ct_j.to_int64());
1869 auto sid = static_cast<index_t>(n_stream_ids_j.to_int64());
1870 auto npts = stream_id_npts[sid];
1871 for(index_t i = 0; i < ec; i++)
1872 {
1873 offsets.push_back(offset);
1874 offset += npts;
1875 elemid++;
1876 }
1877 }
1878 }
1879 else if(n.has_path("elements/element_index/offsets"))
1880 {
1881 const Node &n_stream = n["elements/stream"];
1882 const Node &n_element_offsets = n["elements/element_index/offsets"];
1883 index_t offset = 0, elemid = 0;
1884 for(index_t j = 0; j < n_stream_ids.dtype().number_of_elements(); j++)
1885 {
1886 // Get the j'th elements from n_stream_ids, n_element_offsets
1887 const Node n_stream_ids_j(int_dtype,
1888 const_cast<void*>(n_stream_ids.element_ptr(j)), true);
1889 const Node n_element_offsets_j(int_dtype,
1890 const_cast<void*>(n_element_offsets.element_ptr(j)), true);
1891 offset = n_element_offsets.to_index_t();
1892 index_t next_offset = offset;
1893 if(j == n_stream_ids.dtype().number_of_elements() - 1)
1894 {
1895 next_offset = n_stream.dtype().number_of_elements();
1896 }
1897 else
1898 {
1899 const Node n_element_offsets_j1(int_dtype,
1900 const_cast<void*>(n_element_offsets.element_ptr(j)), true);
1901 next_offset = n_element_offsets_j1.to_index_t();
1902 }
1903 const auto sid = static_cast<index_t>(n_stream_ids_j.to_int64());
1904 const auto npts = stream_id_npts[sid];
1905 while(offset < next_offset) {
1906 offsets.push_back(offset);
1907 offset += npts;
1908 elemid++;
1909 }
1910 }
1911 }
1912 else
1913 {
1914 CONDUIT_ERROR("Stream based mixed topology has no element_counts or offsets.")
1915 }
1916
1917 Node off_node;
1918 off_node.set_external(offsets);
1919 off_node.to_data_type(int_dtype.id(), dest);
1920 }
1921 else if(!topo_shape.is_poly())
1922 {
1923 dest.reset();
1924 // Single element type
1925 const index_t num_topo_shapes =
1926 topo_conn.dtype().number_of_elements() / topo_shape.indices;
1927
1928 Node shape_node(DataType::int64(num_topo_shapes));
1929 int64_array shape_array = shape_node.as_int64_array();
1930 for(index_t s = 0; s < num_topo_shapes; s++)
1931 {
1932 shape_array[s] = s * topo_shape.indices;
1933 }
1934 shape_node.to_data_type(int_dtype.id(), dest);
1935 }
1936 else if(topo_shape.type == "polygonal")
1937 {
1938 dest.reset();
1939
1940 const Node &topo_size = n["elements/sizes"];
1941 std::vector<int64> shape_array;
1942 index_t i = 0;
1943 index_t s = 0;
1944 while(i < topo_size.dtype().number_of_elements())
1945 {
1946 const Node index_node(int_dtype,
1947 const_cast<void*>(topo_size.element_ptr(i)), true);
1948 shape_array.push_back(s);
1949 s += index_node.to_int64();
1950 i++;
1951 }
1952
1953 Node shape_node;
1954 shape_node.set_external(shape_array);
1955 shape_node.to_data_type(int_dtype.id(), dest);
1956 }
1957 else if(topo_shape.type == "polyhedral")
1958 {
1959 Node &dest_elem_off = const_cast<Node &>(n)["elements/offsets"];
1960 Node &dest_subelem_off = const_cast<Node &>(n)["subelements/offsets"];
1961
1962 const Node& topo_elem_size = n["elements/sizes"];
1963 const Node& topo_subelem_size = n["subelements/sizes"];
1964
1965 Node elem_node;
1966 Node subelem_node;
1967
1968 std::vector<index_t> shape_array;
1969 index_t ei = 0;
1970 index_t es = 0;
1971 while(ei < topo_elem_size.dtype().number_of_elements())
1972 {
1973 const Node index_node(int_dtype,
1974 const_cast<void*>(topo_elem_size.element_ptr(ei)), true);
1975 shape_array.push_back(es);
1976 es += index_node.to_index_t();
1977 ei++;
1978 }
1979
1980 elem_node.set_external(shape_array);
1981 elem_node.to_data_type(int_dtype.id(), dest_elem_off);
1982 elem_node.to_data_type(int_dtype.id(), dest);
1983
1984 shape_array.clear();
1985 ei = 0;
1986 es = 0;
1987 while(ei < topo_subelem_size.dtype().number_of_elements())
1988 {
1989 const Node index_node(int_dtype,
1990 const_cast<void*>(topo_subelem_size.element_ptr(ei)), true);
1991 shape_array.push_back(es);
1992 es += index_node.to_index_t();
1993 ei++;
1994 }
1995
1996 subelem_node.set_external(shape_array);
1997 subelem_node.to_data_type(int_dtype.id(), dest_subelem_off);
1998 }
1999 }
2000
2001
2002 //-----------------------------------------------------------------------------
2003 std::vector<index_t>
points(const Node & n,const index_t ei)2004 topology::unstructured::points(const Node &n,
2005 const index_t ei)
2006 {
2007 // NOTE(JRC): This is a workaround to ensure offsets are generated up-front
2008 // if they don't exist and aren't regenerated for each subcall that needs them.
2009 Node ntemp;
2010 ntemp.set_external(n);
2011 generate_offsets(ntemp, ntemp["elements/offsets"]);
2012
2013 Node temp;
2014 const ShapeType topo_shape(ntemp);
2015
2016 std::set<index_t> pidxs;
2017 if(!topo_shape.is_poly())
2018 {
2019 const Node &poffs_node = ntemp["elements/offsets"];
2020 temp.set_external(DataType(poffs_node.dtype().id(), 1),
2021 (void*)poffs_node.element_ptr(ei));
2022 const index_t eoff = temp.to_index_t();
2023
2024 const Node &pidxs_node = ntemp["elements/connectivity"];
2025 for(index_t pi = 0; pi < topo_shape.indices; pi++)
2026 {
2027 temp.set_external(DataType(pidxs_node.dtype().id(), 1),
2028 (void*)pidxs_node.element_ptr(eoff + pi));
2029 pidxs.insert(temp.to_index_t());
2030 }
2031 }
2032 else // if(topo_shape.is_poly())
2033 {
2034 Node enode;
2035 std::set<index_t> eidxs;
2036 if(topo_shape.is_polygonal())
2037 {
2038 enode.set_external(ntemp["elements"]);
2039
2040 eidxs.insert(ei);
2041 }
2042 else // if(topo_shape.is_polyhedral())
2043 {
2044 enode.set_external(ntemp["subelements"]);
2045
2046 const Node &eidxs_node = ntemp["elements/connectivity"];
2047 o2mrelation::O2MIterator eiter(ntemp["elements"]);
2048 eiter.to(ei, O2MIndex::ONE);
2049 eiter.to_front(O2MIndex::MANY);
2050 while(eiter.has_next(O2MIndex::MANY))
2051 {
2052 eiter.next(O2MIndex::MANY);
2053 const index_t ii = eiter.index(O2MIndex::DATA);
2054 temp.set_external(DataType(eidxs_node.dtype().id(), 1),
2055 (void*)eidxs_node.element_ptr(ii));
2056 eidxs.insert(temp.to_index_t());
2057 }
2058 }
2059
2060 for(const index_t eidx : eidxs)
2061 {
2062 const Node &pidxs_node = enode["connectivity"];
2063 o2mrelation::O2MIterator piter(enode);
2064 piter.to(eidx, O2MIndex::ONE);
2065 piter.to_front(O2MIndex::MANY);
2066 while(piter.has_next(O2MIndex::MANY))
2067 {
2068 piter.next(O2MIndex::MANY);
2069 const index_t pi = piter.index(O2MIndex::DATA);
2070 temp.set_external(DataType(pidxs_node.dtype().id(), 1),
2071 (void*)pidxs_node.element_ptr(pi));
2072 pidxs.insert(temp.to_index_t());
2073 }
2074 }
2075 }
2076
2077 return std::vector<index_t>(pidxs.begin(), pidxs.end());
2078 }
2079
2080 //-----------------------------------------------------------------------------
2081 // -- end conduit::blueprint::mesh::utils::topology --
2082 //-----------------------------------------------------------------------------
2083
2084 //-----------------------------------------------------------------------------
2085 // -- begin conduit::blueprint::mesh::utils::adjset --
2086 //-----------------------------------------------------------------------------
2087
2088 //-----------------------------------------------------------------------------
2089 void
canonicalize(Node & adjset)2090 adjset::canonicalize(Node &adjset)
2091 {
2092 const index_t domain_id = find_domain_id(adjset);
2093
2094 const std::vector<std::string> &adjset_group_names = adjset["groups"].child_names();
2095 for(const std::string &old_group_name : adjset_group_names)
2096 {
2097 const Node &group_node = adjset["groups"][old_group_name];
2098 const Node &neighbors_node = group_node["neighbors"];
2099
2100 std::string new_group_name;
2101 {
2102 std::ostringstream oss;
2103 oss << "group";
2104
2105 Node temp;
2106 DataType temp_dtype(neighbors_node.dtype().id(), 1);
2107
2108 // NOTE(JRC): Need to use a vector instead of direct 'Node::to_index_t'
2109 // because the local node ID isn't included in the neighbor list and
2110 // 'DataArray' uses a static array size.
2111 std::vector<index_t> group_neighbors(1, domain_id);
2112 for(index_t ni = 0; ni < neighbors_node.dtype().number_of_elements(); ni++)
2113 {
2114 temp.set_external(temp_dtype, (void*)neighbors_node.element_ptr(ni));
2115 group_neighbors.push_back(temp.to_index_t());
2116 }
2117 std::sort(group_neighbors.begin(), group_neighbors.end());
2118
2119 for(const index_t &neighbor_id : group_neighbors)
2120 {
2121 oss << "_" << neighbor_id;
2122 }
2123
2124 new_group_name = oss.str();
2125 }
2126
2127 adjset["groups"].rename_child(old_group_name, new_group_name);
2128 }
2129 }
2130
2131 //-----------------------------------------------------------------------------
2132 // -- end conduit::blueprint::mesh::utils::adjset --
2133 //-----------------------------------------------------------------------------
2134
2135 }
2136 //-----------------------------------------------------------------------------
2137 // -- end conduit::blueprint::mesh::utils --
2138 //-----------------------------------------------------------------------------
2139
2140
2141 }
2142 //-----------------------------------------------------------------------------
2143 // -- end conduit::blueprint::mesh --
2144 //-----------------------------------------------------------------------------
2145
2146
2147 }
2148 //-----------------------------------------------------------------------------
2149 // -- end conduit::blueprint --
2150 //-----------------------------------------------------------------------------
2151
2152 }
2153 //-----------------------------------------------------------------------------
2154 // -- end conduit:: --
2155 //-----------------------------------------------------------------------------
2156
2157