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