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_examples_julia.cpp
8 ///
9 //-----------------------------------------------------------------------------
10
11 #if defined(CONDUIT_PLATFORM_WINDOWS)
12 #define NOMINMAX
13 #undef min
14 #undef max
15 #include "windows.h"
16 #endif
17
18 //-----------------------------------------------------------------------------
19 // std lib includes
20 //-----------------------------------------------------------------------------
21 #include <string.h>
22 #include <stdio.h>
23 #include <math.h>
24 #include <algorithm>
25 #include <cassert>
26 #include <map>
27 #include <set>
28 #include <vector>
29 #include <queue>
30
31 //-----------------------------------------------------------------------------
32 // conduit includes
33 //-----------------------------------------------------------------------------
34 #include "conduit_blueprint_mesh_examples.hpp"
35 #include "conduit_blueprint_mesh.hpp"
36
37
38 //-----------------------------------------------------------------------------
39 // -- begin conduit:: --
40 //-----------------------------------------------------------------------------
41 namespace conduit
42 {
43
44
45 //-----------------------------------------------------------------------------
46 // -- begin conduit::blueprint:: --
47 //-----------------------------------------------------------------------------
48 namespace blueprint
49 {
50
51
52 //-----------------------------------------------------------------------------
53 // -- begin conduit::blueprint::mesh --
54 //-----------------------------------------------------------------------------
55 namespace mesh
56 {
57
58
59 //-----------------------------------------------------------------------------
60 // -- begin conduit::blueprint::mesh::examples --
61 //-----------------------------------------------------------------------------
62 namespace examples
63 {
64
65 //---------------------------------------------------------------------------//
venn_full_matset(Node & res)66 void venn_full_matset(Node &res)
67 {
68 // create the material sets
69 index_t nx = res["coordsets/coords/params/nx"].value();
70 index_t ny = res["coordsets/coords/params/ny"].value();
71
72 index_t elements = nx * ny;
73
74 // Nodes are evenly spaced from 0 through 1.
75 float64 dx = 1.0 / float64(nx);
76 float64 dy = 1.0 / float64(ny);
77 float64 element_area = dx * dy;
78
79 // create importance reference
80 Node &mat_importance = res["meta/importance"];
81
82 res["matsets/matset/topology"] = "topo";
83 res["matsets/matset/volume_fractions/background"] = DataType::float64(elements);
84 res["matsets/matset/volume_fractions/circle_a"] = DataType::float64(elements);
85 res["matsets/matset/volume_fractions/circle_b"] = DataType::float64(elements);
86 res["matsets/matset/volume_fractions/circle_c"] = DataType::float64(elements);
87
88 res["fields/area/matset_values/background"] = DataType::float64(elements);
89 res["fields/area/matset_values/circle_a"] = DataType::float64(elements);
90 res["fields/area/matset_values/circle_b"] = DataType::float64(elements);
91 res["fields/area/matset_values/circle_c"] = DataType::float64(elements);
92
93 res["fields/importance/matset_values/background"] = DataType::float64(elements);
94 res["fields/importance/matset_values/circle_a"] = DataType::float64(elements);
95 res["fields/importance/matset_values/circle_b"] = DataType::float64(elements);
96 res["fields/importance/matset_values/circle_c"] = DataType::float64(elements);
97
98 res["fields/mat_check/matset_values/background"] = DataType::int64(elements);
99 res["fields/mat_check/matset_values/circle_a"] = DataType::int64(elements);
100 res["fields/mat_check/matset_values/circle_b"] = DataType::int64(elements);
101 res["fields/mat_check/matset_values/circle_c"] = DataType::int64(elements);
102
103 float64_array cir_a = res["fields/circle_a/values"].value();
104 float64_array cir_b = res["fields/circle_b/values"].value();
105 float64_array cir_c = res["fields/circle_c/values"].value();
106 float64_array bg = res["fields/background/values"].value();
107
108 float64_array area = res["fields/area/values"].value();
109 float64_array matset_area_bg = res["fields/area/matset_values/background"].value();
110 float64_array matset_area_cir_a = res["fields/area/matset_values/circle_a"].value();
111 float64_array matset_area_cir_b = res["fields/area/matset_values/circle_b"].value();
112 float64_array matset_area_cir_c = res["fields/area/matset_values/circle_c"].value();
113
114 float64_array importance = res["fields/importance/values"].value();
115 float64_array matset_importance_bg = res["fields/importance/matset_values/background"].value();
116 float64_array matset_importance_cir_a = res["fields/importance/matset_values/circle_a"].value();
117 float64_array matset_importance_cir_b = res["fields/importance/matset_values/circle_b"].value();
118 float64_array matset_importance_cir_c = res["fields/importance/matset_values/circle_c"].value();
119
120 int64_array mat_check = res["fields/mat_check/values"].value();
121 int64_array mat_check_bg = res["fields/mat_check/matset_values/background"].value();
122 int64_array mat_check_bg_cir_a = res["fields/mat_check/matset_values/circle_a"].value();
123 int64_array mat_check_bg_cir_b = res["fields/mat_check/matset_values/circle_b"].value();
124 int64_array mat_check_bg_cir_c = res["fields/mat_check/matset_values/circle_c"].value();
125
126 float64_array mat_bg = res["matsets/matset/volume_fractions/background"].value();
127 float64_array mat_ca = res["matsets/matset/volume_fractions/circle_a"].value();
128 float64_array mat_cb = res["matsets/matset/volume_fractions/circle_b"].value();
129 float64_array mat_cc = res["matsets/matset/volume_fractions/circle_c"].value();
130
131 for(index_t idx = 0; idx < elements; idx++)
132 {
133 mat_ca[idx] = cir_a[idx];
134 mat_cb[idx] = cir_b[idx];
135 mat_cc[idx] = cir_c[idx];
136 mat_bg[idx] = bg[idx];
137
138 if (mat_ca[idx] > 0.) { matset_area_cir_a[idx] = element_area; }
139 if (mat_cb[idx] > 0.) { matset_area_cir_b[idx] = element_area; }
140 if (mat_cc[idx] > 0.) { matset_area_cir_c[idx] = element_area; }
141 if (mat_bg[idx] > 0.) { matset_area_bg[idx] = element_area; }
142
143 // the overall sum provides a value that can id all mats
144 // in a zone, these are the per-material components
145 if (mat_bg[idx] > 0.) { mat_check_bg[idx] = 1; }
146 if (mat_ca[idx] > 0.) { mat_check_bg_cir_a[idx] = 20; }
147 if (mat_cb[idx] > 0.) { mat_check_bg_cir_b[idx] = 300; }
148 if (mat_cc[idx] > 0.) { mat_check_bg_cir_c[idx] = 4000; }
149
150 area[idx] = mat_ca[idx] * matset_area_cir_a[idx] +
151 mat_cb[idx] * matset_area_cir_b[idx] +
152 mat_cc[idx] * matset_area_cir_c[idx] +
153 mat_bg[idx] * matset_area_bg[idx];
154
155 if (mat_ca[idx] > 0.) { matset_importance_cir_a[idx] = mat_importance["a"].value(); }
156 if (mat_cb[idx] > 0.) { matset_importance_cir_b[idx] = mat_importance["b"].value(); }
157 if (mat_cc[idx] > 0.) { matset_importance_cir_c[idx] = mat_importance["c"].value(); }
158 float64 x_pos = ((float64)(idx % nx)) / nx;
159 float64 y_pos = ((float64)(idx / nx)) / ny;
160 if (mat_bg[idx] > 0.) { matset_importance_bg[idx] = x_pos + y_pos; }
161
162 importance[idx] = mat_ca[idx] * matset_importance_cir_a[idx] +
163 mat_cb[idx] * matset_importance_cir_b[idx] +
164 mat_cc[idx] * matset_importance_cir_c[idx] +
165 mat_bg[idx] * matset_importance_bg[idx];
166 }
167 }
168
169 //---------------------------------------------------------------------------//
build_material_sparse(Node & src,index_t len,const std::string & mat_name,float64 element_area,float64 material_importance,Node & matset_area,Node & matset_importance,Node & matset)170 void build_material_sparse(Node & src, index_t len,
171 const std::string & mat_name,
172 float64 element_area, float64 material_importance,
173 Node & matset_area, Node & matset_importance, Node & matset)
174 {
175 float64_array src_val = src.value();
176
177 index_t nsparse = 0;
178 for (index_t idx = 0; idx < len; ++idx)
179 {
180 if (src_val[idx] > 0)
181 {
182 nsparse += 1;
183 }
184 }
185
186 matset["volume_fractions/" + mat_name].set(DataType::float64(nsparse));
187 matset["element_ids/" + mat_name].set(DataType::int32(nsparse));
188 float64_array sparse_val = matset["volume_fractions/" + mat_name].value();
189 int32_array sparse_element_ids = matset["element_ids/" + mat_name].value();
190
191 matset_area.set(DataType::float64(nsparse));
192 float64_array matset_area_val = matset_area.value();
193 matset_importance.set(DataType::float64(nsparse));
194 float64_array matset_importance_val = matset_importance.value();
195
196 index_t sparse_idx = 0;
197 for (index_t idx = 0; idx < len; ++idx)
198 {
199 if (src_val[idx] > 0)
200 {
201 sparse_element_ids[sparse_idx] = (int32)idx;
202 sparse_val[sparse_idx] = src_val[idx];
203
204 matset_area_val[sparse_idx] = element_area;
205 matset_importance_val[sparse_idx] = material_importance;
206
207 sparse_idx += 1;
208 }
209 }
210 }
211
212 //---------------------------------------------------------------------------//
compute_material_sparse_matset_field(Node & res,const std::string & field_name)213 void compute_material_sparse_matset_field(Node &res,
214 const std::string & field_name)
215 {
216 index_t nx = res["coordsets/coords/params/nx"].value();
217 index_t ny = res["coordsets/coords/params/ny"].value();
218 index_t elements = nx * ny;
219
220 Node & n = res["fields/" + field_name + "/values"];
221 n.set(DataType::float64(elements));
222 float64_array n_val = n.value();
223
224 Node & matset_values = res["fields/" + field_name + "/matset_values"];
225
226 NodeIterator itr = matset_values.children();
227 while (itr.has_next())
228 {
229 Node &cld = itr.next();
230 const std::string & cld_name = itr.name();
231 float64_array matset_vals = cld.value();
232
233 float64_array vf_vals = res["matsets/matset/volume_fractions/" + cld_name].value();
234 int32_array vf_elt_ids = res["matsets/matset/element_ids/" + cld_name].value();
235 index_t sparse_elements = vf_elt_ids.number_of_elements();
236
237 index_t sparse_index = 0;
238 for (index_t elt = 0; elt < elements && sparse_index < sparse_elements; ++elt)
239 {
240 if (vf_elt_ids[sparse_index] == elt)
241 {
242 n_val[elt] += matset_vals[sparse_index] * vf_vals[sparse_index];
243 sparse_index += 1;
244 }
245 }
246 }
247 }
248
249 //---------------------------------------------------------------------------//
venn_sparse_by_material_matset(Node & res)250 void venn_sparse_by_material_matset(Node &res)
251 {
252 // create the materials
253 index_t nx = res["coordsets/coords/params/nx"].value();
254 index_t ny = res["coordsets/coords/params/ny"].value();
255
256 // make sure our materials appear in the correct order,
257 // by pre-creating the nodes.
258 res["matsets/matset/volume_fractions/background"];
259 res["matsets/matset/volume_fractions/circle_a"];
260 res["matsets/matset/volume_fractions/circle_b"];
261 res["matsets/matset/volume_fractions/circle_c"];
262
263 res["matsets/matset/element_ids/background"];
264 res["matsets/matset/element_ids/circle_a"];
265 res["matsets/matset/element_ids/circle_b"];
266 res["matsets/matset/element_ids/circle_c"];
267
268 // we could also use the following material map:
269 //
270 // res["matsets/matset/material_map/background"] = 0;
271 // res["matsets/matset/material_map/circle_a"] = 1;
272 // res["matsets/matset/material_map/circle_b"] = 2;
273 // res["matsets/matset/material_map/circle_c"] = 3;
274 //
275 // however, this current setup is a good test of the
276 // non explicit map case
277
278
279 float64_array cir_a = res["fields/circle_a/values"].value();
280 float64_array cir_b = res["fields/circle_b/values"].value();
281 float64_array cir_c = res["fields/circle_c/values"].value();
282
283 res["matsets/matset/topology"] = "topo";
284
285 // Nodes are evenly spaced from 0 through 1.
286 float64 dx = 1.0 / float64(nx);
287 float64 dy = 1.0 / float64(ny);
288 float64 element_area = dx * dy;
289 index_t elements = nx * ny;
290
291 // create importance reference
292 Node &mat_importance = res["meta/importance"];
293
294 build_material_sparse(res["fields/circle_a/values"],
295 elements,
296 "circle_a",
297 element_area,
298 mat_importance["a"].value(),
299 res["fields/area/matset_values/circle_a"],
300 res["fields/importance/matset_values/circle_a"],
301 res["matsets/matset"]);
302
303 build_material_sparse(res["fields/circle_b/values"],
304 elements,
305 "circle_b",
306 element_area,
307 mat_importance["b"].value(),
308 res["fields/area/matset_values/circle_b"],
309 res["fields/importance/matset_values/circle_b"],
310 res["matsets/matset"]);
311
312 build_material_sparse(res["fields/circle_c/values"],
313 elements,
314 "circle_c",
315 element_area,
316 mat_importance["c"].value(),
317 res["fields/area/matset_values/circle_c"],
318 res["fields/importance/matset_values/circle_c"],
319 res["matsets/matset"]);
320
321 // The background material volume fraction depends on the other three
322 // materials, so we deal with it in a custom loop.
323 index_t bgcount = 0;
324 for (index_t idx = 0; idx < elements; ++idx)
325 {
326 if (cir_a[idx] + cir_b[idx] + cir_c[idx] < 1.) bgcount += 1;
327 }
328
329 res["matsets/matset/volume_fractions/background"].set(DataType::float64(bgcount));
330 res["matsets/matset/element_ids/background"].set(DataType::int32(bgcount));
331 float64_array bg_val = res["matsets/matset/volume_fractions/background"].value();
332 int32_array bg_idx = res["matsets/matset/element_ids/background"].value();
333
334 Node &matset_area_bg = res["fields/area/matset_values/background"];
335 matset_area_bg.set(DataType::float64(bgcount));
336 float64_array matset_area_bg_value = matset_area_bg.value();
337
338 Node &matset_importance_bg = res["fields/importance/matset_values/background"];
339 matset_importance_bg.set(DataType::float64(bgcount));
340 float64_array matset_importance_bg_value = matset_importance_bg.value();
341
342 index_t nidx = 0;
343 for (index_t idx = 0; idx < elements; ++idx)
344 {
345 float64 x_pos = ((float64)(idx % nx)) / nx;
346 float64 y_pos = ((float64)(idx / nx)) / ny;
347
348 float64 fgvf = cir_a[idx] + cir_b[idx] + cir_c[idx];
349 if (fgvf < 1.)
350 {
351 bg_idx[nidx] = (int32)idx;
352
353 bg_val[nidx] = 1. - fgvf;
354
355 matset_area_bg_value[nidx] = element_area;
356 matset_importance_bg_value[nidx] = x_pos + y_pos;
357
358 nidx += 1;
359 }
360 }
361
362 // --------- //
363 // mat_check //
364 // --------- //
365
366 // the matset now has element id's for each material computed.
367 // to create mat_check, we simply need to fill arrays of the correct size
368 // with our unique value for each mat
369
370 // Note: we already have bgcount
371 index_t cir_a_count = res["matsets/matset/element_ids/circle_a"].dtype().number_of_elements();
372 index_t cir_b_count = res["matsets/matset/element_ids/circle_b"].dtype().number_of_elements();
373 index_t cir_c_count = res["matsets/matset/element_ids/circle_c"].dtype().number_of_elements();
374
375 // init our matset_values arrays
376 res["fields/mat_check/matset_values/background"].set(DataType::int64(bgcount));
377 res["fields/mat_check/matset_values/circle_a"].set(DataType::int64(cir_a_count));
378 res["fields/mat_check/matset_values/circle_b"].set(DataType::int64(cir_b_count));
379 res["fields/mat_check/matset_values/circle_c"].set(DataType::int64(cir_c_count));
380
381 // set with unique values
382
383 // NOTE: this is good use case for adding DataArray.fill ...
384 int64_array m_chk_bg_vals = res["fields/mat_check/matset_values/background"].value();
385
386 for(index_t idx=0;idx < bgcount; idx++)
387 {
388 m_chk_bg_vals[idx] = 1;
389 }
390
391 int64_array m_chk_cir_a_vals = res["fields/mat_check/matset_values/circle_a"].value();
392 for(index_t idx=0;idx < cir_a_count; idx++)
393 {
394 m_chk_cir_a_vals[idx] = 20;
395 }
396
397 int64_array m_chk_cir_b_vals = res["fields/mat_check/matset_values/circle_b"].value();
398 for(index_t idx=0;idx < cir_b_count; idx++)
399 {
400 m_chk_cir_b_vals[idx] = 300;
401 }
402
403 int64_array m_chk_cir_c_vals = res["fields/mat_check/matset_values/circle_c"].value();
404 for(index_t idx=0;idx < cir_c_count; idx++)
405 {
406 m_chk_cir_c_vals[idx] = 4000;
407 }
408
409 // Now we've computed the matset values for the fields area and
410 // importance, sum the product of the volume fraction and the matset
411 // value for each element to compute the field itself.
412
413 compute_material_sparse_matset_field(res, "area");
414 compute_material_sparse_matset_field(res, "importance");
415 }
416
venn_sparse_by_element_matset(Node & res)417 void venn_sparse_by_element_matset(Node &res)
418 {
419 // create the materials
420 index_t nx = res["coordsets/coords/params/nx"].value();
421 index_t ny = res["coordsets/coords/params/ny"].value();
422
423 float64_array cir_a = res["fields/circle_a/values"].value();
424 float64_array cir_b = res["fields/circle_b/values"].value();
425 float64_array cir_c = res["fields/circle_c/values"].value();
426
427 // Nodes are evenly spaced from 0 through 1.
428 float64 dx = 1.0 / float64(nx);
429 float64 dy = 1.0 / float64(ny);
430 float64 element_area = dx * dy;
431 index_t elements = nx * ny;
432 index_t vfcount = 0;
433
434 // create importance reference
435 Node &mat_importance = res["meta/importance"];
436
437 // Count up all the non-zero volume fragments
438 // (so we can allocate a correctly-sized buffer)
439 for (index_t idx = 0; idx < elements; ++idx)
440 {
441 if (cir_a[idx] > 0.) vfcount += 1;
442 if (cir_b[idx] > 0.) vfcount += 1;
443 if (cir_c[idx] > 0.) vfcount += 1;
444 if (cir_a[idx] + cir_b[idx] + cir_c[idx] < 1.) vfcount += 1;
445 }
446
447 // Build the rest of the single-buffer matset
448 res["matsets/matset/topology"] = "topo";
449 // This is the "key" that tells what material each volume fraction refers to
450 res["matsets/matset/material_map/circle_a"] = 1;
451 res["matsets/matset/material_map/circle_b"] = 2;
452 res["matsets/matset/material_map/circle_c"] = 3;
453 res["matsets/matset/material_map/background"] = 0;
454
455 // All the volume fractions go here ("one big buffer")
456 res["matsets/matset/volume_fractions"].set(DataType::float64(vfcount));
457 // The material measured by each volume fraction
458 res["matsets/matset/material_ids"].set(DataType::int32(vfcount));
459 // The number of volume fractions in an element
460 res["matsets/matset/sizes"].set(DataType::int32(elements));
461 // The offset of the first vf in an element
462 res["matsets/matset/offsets"].set(DataType::int32(elements));
463
464 // The matset values (for fields that have them) get built up
465 // in the same way as the sparse fields, into "one big array."
466 res["fields/area/matset_values"].set(DataType::float64(vfcount));
467 res["fields/importance/matset_values"].set(DataType::float64(vfcount));
468 res["fields/mat_check/matset_values"].set(DataType::float64(vfcount));
469
470 // The actual fields
471 res["fields/area/values"].set(DataType::float64(elements));
472 res["fields/importance/values"].set(DataType::float64(elements));
473
474 float64_array vf = res["matsets/matset/volume_fractions"].value();
475 int32_array id = res["matsets/matset/material_ids"].value();
476 int32_array sizes = res["matsets/matset/sizes"].value();
477 int32_array offsets = res["matsets/matset/offsets"].value();
478
479 float64_array matset_area = res["fields/area/matset_values"].value();
480 float64_array matset_impt = res["fields/importance/matset_values"].value();
481
482 float64_array field_area = res["fields/area/values"].value();
483 float64_array field_impt = res["fields/importance/values"].value();
484
485 float64_array matset_check_vals = res["fields/mat_check/matset_values"].value();
486
487 // Build up the arrays!
488 index_t vfidx = 0;
489 for (index_t idx = 0; idx < elements; ++idx)
490 {
491 int size = 0;
492 float64 ca = cir_a[idx];
493 float64 cb = cir_b[idx];
494 float64 cc = cir_c[idx];
495
496 field_area[idx] = 0.;
497 field_impt[idx] = 0.;
498
499 auto fill_in = [&](float64 frac, float64 imp, int64 mcheck, int32 mat_id)
500 {
501 vf[vfidx + size] = frac;
502 matset_area[vfidx + size] = element_area;
503 matset_impt[vfidx + size] = imp;
504 matset_check_vals[vfidx + size] = mcheck;
505 id[vfidx + size] = mat_id;
506
507 field_area[idx] += matset_area[vfidx + size] * vf[vfidx + size];
508 field_impt[idx] += matset_impt[vfidx + size] * vf[vfidx + size];
509
510 size += 1;
511 };
512
513 if (ca > 0.)
514 {
515 fill_in(ca, mat_importance["a"].value(), 20, 1);
516 }
517 if (cb > 0.)
518 {
519 fill_in(cb, mat_importance["b"].value(), 300, 2);
520 }
521 if (cc > 0.)
522 {
523 fill_in(cc, mat_importance["c"].value(), 4000, 3);
524 }
525 if (ca + cb + cc < 1.)
526 {
527 float64 x_pos = ((float64)(idx % nx)) / nx;
528 float64 y_pos = ((float64)(idx / nx)) / ny;
529 fill_in(1 - (ca + cb + cc), x_pos + y_pos, 1, 0);
530 }
531
532 sizes[idx] = size;
533 offsets[idx] = (int32) vfidx;
534 vfidx += size;
535 }
536 }
537
venn(const std::string & matset_type,index_t nx,index_t ny,float64 radius,Node & res)538 void venn(const std::string &matset_type,
539 index_t nx,
540 index_t ny,
541 float64 radius,
542 Node &res)
543 {
544 res.reset();
545 // create a rectilinear coordset
546 res["coordsets/coords/type"] = "rectilinear";
547 res["coordsets/coords/values/x"] = DataType::float64(nx+1);
548 res["coordsets/coords/values/y"] = DataType::float64(ny+1);
549
550 // Not part of the blueprint, but I want these values handy
551 res["coordsets/coords/params/nx"] = nx;
552 res["coordsets/coords/params/ny"] = ny;
553
554 float64_array x_coords = res["coordsets/coords/values/x"].value();
555 float64_array y_coords = res["coordsets/coords/values/y"].value();
556
557 // 0 <-> 1
558 float64 dx = 1.0/float64(nx);
559 float64 dy = 1.0/float64(ny);
560
561 float64 vx = 0;
562 for(index_t i =0; i< nx+1; i++)
563 {
564 x_coords[i] = vx;
565 vx+=dx;
566 }
567
568 float64 vy = 0;
569 for(index_t i =0; i< ny+1; i++)
570 {
571 y_coords[i] = vy;
572 vy+=dy;
573 }
574
575 // create the topology
576
577 res["topologies/topo/type"] = "rectilinear";
578 res["topologies/topo/coordset"] = "coords";
579
580 // create the fields
581
582 // circle a distance field
583 res["fields/radius_a/association"] = "element";
584 res["fields/radius_a/topology"] = "topo";
585 res["fields/radius_a/values"] = DataType::float64(nx * ny);
586 // circle a vf
587 res["fields/circle_a/association"] = "element";
588 res["fields/circle_a/topology"] = "topo";
589 res["fields/circle_a/values"] = DataType::float64(nx * ny);
590
591 // circle b distance field
592 res["fields/radius_b/association"] = "element";
593 res["fields/radius_b/topology"] = "topo";
594 res["fields/radius_b/values"] = DataType::float64(nx * ny);
595 // circle b vf
596 res["fields/circle_b/association"] = "element";
597 res["fields/circle_b/topology"] = "topo";
598 res["fields/circle_b/values"] = DataType::float64(nx * ny);
599
600 // circle c distance field
601 res["fields/radius_c/association"] = "element";
602 res["fields/radius_c/topology"] = "topo";
603 res["fields/radius_c/values"] = DataType::float64(nx * ny);
604 // circle b vf
605 res["fields/circle_c/association"] = "element";
606 res["fields/circle_c/topology"] = "topo";
607 res["fields/circle_c/values"] = DataType::float64(nx * ny);
608
609 // per element how many circles overlap
610 res["fields/overlap/association"] = "element";
611 res["fields/overlap/topology"] = "topo";
612 res["fields/overlap/values"] = DataType::float64(nx * ny);
613
614 // per element background
615 res["fields/background/association"] = "element";
616 res["fields/background/topology"] = "topo";
617 res["fields/background/values"] = DataType::float64(nx * ny);
618
619 // per element field with matset values.
620 //
621 // For a field with matset values, each element will have a value equal
622 // to the sum of each material's volume fraction times the material's
623 // matset value.
624 //
625 // As with ordinary fields, fields/field/values holds the overall value
626 // for each cell. Additionally, fields/field/matset_values holds the
627 // matset value for each contributing material. The matset values are
628 // stored in the same way as matsets/matset/volume_fractions (full,
629 // sparse-by-material, one-buffer-sparse-by-element).
630 //
631 // Area is trivial to compute and easy to verify.
632 res["fields/area/association"] = "element";
633 res["fields/area/topology"] = "topo";
634 res["fields/area/matset"] = "matset";
635 res["fields/area/values"] = DataType::float64(nx * ny);
636 // "Importance" is a made-up field.
637 // Circles a, b, and c have differing importance.
638 // Background material has importance that varies with position.
639 res["fields/importance/association"] = "element";
640 res["fields/importance/topology"] = "topo";
641 res["fields/importance/matset"] = "matset";
642 res["fields/importance/values"] = DataType::float64(nx * ny);
643
644 // "mat_check" is a made-up field.
645 // It computes unique sum that encodes the material ids
646 // for verification, and provides per-material components
647 // via matset_values
648 res["fields/mat_check/association"] = "element";
649 res["fields/mat_check/topology"] = "topo";
650 res["fields/mat_check/matset"] = "matset";
651 res["fields/mat_check/values"] = DataType::int64(nx * ny);
652
653 float64_array rad_a = res["fields/radius_a/values"].value();
654 float64_array cir_a = res["fields/circle_a/values"].value();
655
656 float64_array rad_b = res["fields/radius_b/values"].value();
657 float64_array cir_b = res["fields/circle_b/values"].value();
658
659 float64_array rad_c = res["fields/radius_c/values"].value();
660 float64_array cir_c = res["fields/circle_c/values"].value();
661
662 float64_array bg = res["fields/background/values"].value();
663
664 float64_array area = res["fields/area/values"].value();
665 float64_array importance = res["fields/importance/values"].value();
666
667 float64_array olap = res["fields/overlap/values"].value();
668
669 int64_array mat_check = res["fields/mat_check/values"].value();
670
671
672 // circle a
673 // centered at:
674 // x = 2/3 of width
675 // y = 2/3 of width
676
677 float64 a_center_x = 0.3333333333;
678 float64 a_center_y = 0.6666666666;
679 float64 a_importance = 0.1f;
680
681 // circle b
682 // centered at:
683 // x = 2/3 of width
684 // y = 2/3 of width
685
686 float64 b_center_x = 0.6666666666;
687 float64 b_center_y = 0.6666666666;
688 float64 b_importance = 0.2f;
689
690 // circle c
691 // centered at:
692 // x = 1/2 of width
693 // y = 2/3 of width
694
695 float64 c_center_x = 0.5;
696 float64 c_center_y = 0.3333333333;
697 float64 c_importance = 0.6f;
698
699 // Fill in fields
700
701 float64 y = 0;
702
703 index_t idx = 0;
704 for(index_t j = 0; j < ny; j++)
705 {
706 float64 x = 0;
707 for(index_t i = 0; i < nx; i++)
708 {
709 // dist to circle a
710 rad_a[idx] = sqrt( (x - a_center_x)*(x - a_center_x) +
711 (y - a_center_y)*(y - a_center_y));
712 if(rad_a[idx] > radius)
713 rad_a[idx] = 0.0; // clamp outside of radius
714
715 // dist to circle b
716 rad_b[idx] = sqrt( (x - b_center_x)*(x - b_center_x) +
717 (y - b_center_y)*(y - b_center_y));
718 if(rad_b[idx] > radius)
719 rad_b[idx] = 0.0; // clamp outside of radius
720
721 // dist to circle c
722 rad_c[idx] = sqrt( (x - c_center_x)*(x - c_center_x) +
723 (y - c_center_y)*(y - c_center_y));
724 if(rad_c[idx] > radius)
725 rad_c[idx] = 0.0; // clamp outside of radius
726
727 // populate overlap with count
728 if(rad_a[idx] > 0)
729 olap[idx]++;
730
731 if(rad_b[idx] > 0)
732 olap[idx]++;
733
734 if(rad_c[idx] > 0)
735 olap[idx]++;
736
737 // circle vfs
738 if(rad_a[idx] > 0)
739 {
740 mat_check[idx] += 20;
741 cir_a[idx] = 1.0/olap[idx];
742 }
743 else
744 {
745 cir_a[idx] = 0.0;
746 }
747
748 if(rad_b[idx] > 0)
749 {
750 mat_check[idx] += 300;
751 cir_b[idx] = 1.0/olap[idx];
752 }
753 else
754 {
755 cir_b[idx] = 0.0;
756 }
757
758 if(rad_c[idx] > 0)
759 {
760 mat_check[idx] += 4000;
761 cir_c[idx] = 1.0/olap[idx];
762 }
763 else
764 {
765 cir_c[idx] = 0.0;
766 }
767
768 // bg vf
769 bg[idx] = 1. - (cir_a[idx] + cir_b[idx] + cir_c[idx]);
770
771 if(bg[idx] > 0.0 )
772 {
773 mat_check[idx] += 1;
774 }
775
776 // initialize area and importance to 0.
777 area[idx] = 0.;
778 importance[idx] = 0.;
779
780 x+=dx;
781 idx++;
782 }
783 y+=dy;
784 }
785
786 // Fill in metadata (used by helper functions)
787
788 res["meta/importance/a"] = a_importance;
789 res["meta/importance/b"] = b_importance;
790 res["meta/importance/c"] = c_importance;
791
792 // Shape in materials; compute fields with matset values.
793
794 if (matset_type == "full")
795 {
796 venn_full_matset(res);
797 }
798 else if (matset_type == "sparse_by_material")
799 {
800 venn_sparse_by_material_matset(res);
801 }
802 else if (matset_type == "sparse_by_element")
803 {
804 venn_sparse_by_element_matset(res);
805 }
806 else
807 {
808 CONDUIT_ERROR("unknown matset_type = " << matset_type);
809 }
810
811 // remove temp tree used during construction
812 res.remove("meta");
813 }
814
815 }
816 //-----------------------------------------------------------------------------
817 // -- end conduit::blueprint::mesh::examples --
818 //-----------------------------------------------------------------------------
819
820
821 }
822 //-----------------------------------------------------------------------------
823 // -- end conduit::blueprint::mesh --
824 //-----------------------------------------------------------------------------
825
826
827 }
828 //-----------------------------------------------------------------------------
829 // -- end conduit::blueprint:: --
830 //-----------------------------------------------------------------------------
831
832
833 }
834 //-----------------------------------------------------------------------------
835 // -- end conduit:: --
836 //-----------------------------------------------------------------------------
837