MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
conduitdatacollection.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_CONDUIT
15 
16 #include "fem.hpp"
17 #include "../general/text.hpp"
18 #include <conduit_relay.hpp>
19 #include <conduit_blueprint.hpp>
20 
21 #include <string>
22 #include <sstream>
23 
24 using namespace conduit;
25 
26 namespace mfem
27 {
28 
29 //---------------------------------------------------------------------------//
30 // class ConduitDataCollection implementation
31 //---------------------------------------------------------------------------//
32 
33 //------------------------------
34 // begin public methods
35 //------------------------------
36 
37 //---------------------------------------------------------------------------//
38 ConduitDataCollection::ConduitDataCollection(const std::string& coll_name,
39  Mesh *mesh)
40  : DataCollection(coll_name, mesh),
41  relay_protocol("hdf5")
42 {
43  appendRankToFileName = true; // always include rank in file names
44  cycle = 0; // always include cycle in directory names
45 }
46 
47 #ifdef MFEM_USE_MPI
48 //---------------------------------------------------------------------------//
50  const std::string& coll_name,
51  Mesh *mesh)
52  : DataCollection(coll_name, mesh),
53  relay_protocol("hdf5")
54 {
55  m_comm = comm;
56  MPI_Comm_rank(comm, &myid);
57  MPI_Comm_size(comm, &num_procs);
58  appendRankToFileName = true; // always include rank in file names
59  cycle = 0; // always include cycle in directory names
60 }
61 #endif
62 
63 //---------------------------------------------------------------------------//
65 {
66  // empty
67 }
68 
69 //---------------------------------------------------------------------------//
71 {
72  std::string dir_name = MeshDirectoryName();
73  int err = create_directory(dir_name, mesh, myid);
74  if (err)
75  {
76  MFEM_ABORT("Error creating directory: " << dir_name);
77  }
78 
79  Node n_mesh;
80  // future? If moved into Mesh class
81  // mesh->toConduitBlueprint(n_mesh);
82  MeshToBlueprintMesh(mesh,n_mesh);
83 
84  Node verify_info;
85  if (!blueprint::mesh::verify(n_mesh,verify_info))
86  {
87  MFEM_ABORT("Conduit Mesh Blueprint Verify Failed:\n"
88  << verify_info.to_json());
89  }
90 
92  for ( itr = field_map.begin(); itr != field_map.end(); itr++)
93  {
94  std::string name = itr->first;
95  GridFunction *gf = itr->second;
96  // don't save mesh nodes twice ...
97  if ( gf != mesh->GetNodes())
98  {
99  // future? If moved into GridFunction class
100  //gf->toConduitBlueprint(n_mesh["fields"][it->first]);
102  n_mesh["fields"][name]);
103  }
104  }
105 
106  // save mesh data
108  n_mesh,
110 
111  if (myid == 0)
112  {
113  // save root file
115  n_mesh,
117  }
118 }
119 
120 //---------------------------------------------------------------------------//
122 {
123  DeleteAll();
124  this->cycle = cycle;
125 
126  // Note: We aren't currently using much info from the root file ...
127  // with cycle, we can use implicit mfem conduit file layout
128 
129  Node n_root;
130  LoadRootFile(n_root);
131  relay_protocol = n_root["protocol/name"].as_string();
132 
133  // for MPI case, we assume that we have # of mpi tasks
134  // == number of domains
135 
136  int num_domains = n_root["number_of_trees"].to_int();
137 
138  if (num_procs != num_domains)
139  {
140  error = READ_ERROR;
141  MFEM_WARNING("num_procs must equal num_domains");
142  return;
143  }
144 
145  // load the mesh and fields
147 
148  // TODO: am I properly wielding this?
149  own_data = true;
150 
151 }
152 
153 //---------------------------------------------------------------------------//
154 void
155 ConduitDataCollection::SetProtocol(const std::string &protocol)
156 {
157  relay_protocol = protocol;
158 }
159 
160 //------------------------------
161 // begin static public methods
162 //------------------------------
163 
164 //---------------------------------------------------------------------------//
165 mfem::Mesh *
167  const std::string &main_toplogy_name,
168  bool zero_copy)
169 {
170  // n_conv holds converted data (when necessary for mfem api)
171  // if n_conv is used ( !n_conv.dtype().empty() ) we
172  // now that some data allocation was necessary, so we
173  // can't return a mesh that zero copies the conduit data
174  Node n_conv;
175 
176  //
177  // we need to find the topology and its coordset.
178  //
179 
180  std::string topo_name = main_toplogy_name;
181  // if topo name is not set, look for first topology
182  if (topo_name == "")
183  {
184  topo_name = n_mesh["topologies"].schema().child_name(0);
185  }
186 
187  MFEM_ASSERT(n_mesh.has_path("topologies/" + topo_name),
188  "Expected topology named \"" + topo_name + "\" "
189  "(node is missing path \"topologies/" + topo_name + "\")");
190 
191  // find the coord set
192  std::string coords_name =
193  n_mesh["topologies"][topo_name]["coordset"].as_string();
194 
195 
196  MFEM_ASSERT(n_mesh.has_path("coordsets/" + coords_name),
197  "Expected topology named \"" + coords_name + "\" "
198  "(node is missing path \"coordsets/" + coords_name + "\")");
199 
200  const Node &n_coordset = n_mesh["coordsets"][coords_name];
201  const Node &n_coordset_vals = n_coordset["values"];
202 
203  // get the number of dims of the coordset
204  int ndims = n_coordset_vals.number_of_children();
205 
206  // get the number of points
207  int num_verts = n_coordset_vals[0].dtype().number_of_elements();
208  // get vals for points
209  const double *verts_ptr = NULL;
210 
211  // the mfem mesh constructor needs coords with interleaved (aos) type
212  // ordering, even for 1d + 2d we always need 3 doubles b/c it uses
213  // Array<Vertex> and Vertex is a pod of 3 doubles. we check for this
214  // case, if we don't have it we convert the data
215 
216  if (ndims == 3 &&
217  n_coordset_vals[0].dtype().is_double() &&
218  blueprint::mcarray::is_interleaved(n_coordset_vals) )
219  {
220  // already interleaved mcarray of 3 doubles,
221  // return ptr to beginning
222  verts_ptr = n_coordset_vals[0].value();
223  }
224  else
225  {
226  Node n_tmp;
227  // check all vals, if we don't have doubles convert
228  // to doubles
229  NodeConstIterator itr = n_coordset_vals.children();
230  while (itr.has_next())
231  {
232  const Node &c_vals = itr.next();
233  std::string c_name = itr.name();
234 
235  if ( c_vals.dtype().is_double() )
236  {
237  // zero copy current coords
238  n_tmp[c_name].set_external(c_vals);
239 
240  }
241  else
242  {
243  // convert
244  c_vals.to_double_array(n_tmp[c_name]);
245  }
246  }
247 
248  // check if we need to add extra dims to get
249  // proper interleaved array
250  if (ndims < 3)
251  {
252  // add dummy z
253  n_tmp["z"].set(DataType::c_double(num_verts));
254  }
255 
256  if (ndims < 2)
257  {
258  // add dummy y
259  n_tmp["y"].set(DataType::c_double(num_verts));
260  }
261 
262  Node &n_conv_coords_vals = n_conv["coordsets"][coords_name]["values"];
263  blueprint::mcarray::to_interleaved(n_tmp,
264  n_conv_coords_vals);
265  verts_ptr = n_conv_coords_vals[0].value();
266  }
267 
268 
269 
270  const Node &n_mesh_topo = n_mesh["topologies"][topo_name];
271  std::string mesh_ele_shape = n_mesh_topo["elements/shape"].as_string();
272 
273  mfem::Geometry::Type mesh_geo = ShapeNameToGeomType(mesh_ele_shape);
274  int num_idxs_per_ele = Geometry::NumVerts[mesh_geo];
275 
276  const Node &n_mesh_conn = n_mesh_topo["elements/connectivity"];
277 
278  const int *elem_indices = NULL;
279  // mfem requires ints, we could have int64s, etc convert if necessary
280  if (n_mesh_conn.dtype().is_int() &&
281  n_mesh_conn.is_compact() )
282  {
283  elem_indices = n_mesh_topo["elements/connectivity"].value();
284  }
285  else
286  {
287  Node &n_mesh_conn_conv=
288  n_conv["topologies"][topo_name]["elements/connectivity"];
289  n_mesh_conn.to_int_array(n_mesh_conn_conv);
290  elem_indices = n_mesh_conn_conv.value();
291  }
292 
293  int num_mesh_ele =
294  n_mesh_topo["elements/connectivity"].dtype().number_of_elements();
295  num_mesh_ele = num_mesh_ele / num_idxs_per_ele;
296 
297 
298  const int *bndry_indices = NULL;
299  int num_bndry_ele = 0;
300  // init to something b/c the mesh constructor will use this for a
301  // table lookup, even if we don't have boundary info.
303 
304  if ( n_mesh_topo.has_child("boundary_topology") )
305  {
306  std::string bndry_topo_name = n_mesh_topo["boundary_topology"].as_string();
307 
308  // In VisIt, we encountered a case were a mesh specified a boundary
309  // topology, but the boundary topology was omitted from the blueprint
310  // index, so it's data could not be obtained.
311  //
312  // This guard prevents an error in that case, allowing the mesh to be
313  // created without boundary info
314 
315  if (n_mesh["topologies"].has_child(bndry_topo_name))
316  {
317  const Node &n_bndry_topo = n_mesh["topologies"][bndry_topo_name];
318  std::string bndry_ele_shape = n_bndry_topo["elements/shape"].as_string();
319 
320  bndry_geo = ShapeNameToGeomType(bndry_ele_shape);
321  int num_idxs_per_bndry_ele = Geometry::NumVerts[mesh_geo];
322 
323  const Node &n_bndry_conn = n_bndry_topo["elements/connectivity"];
324 
325  // mfem requires ints, we could have int64s, etc convert if necessary
326  if ( n_bndry_conn.dtype().is_int() &&
327  n_bndry_conn.is_compact())
328  {
329  bndry_indices = n_bndry_conn.value();
330  }
331  else
332  {
333  Node &(n_bndry_conn_conv) =
334  n_conv["topologies"][bndry_topo_name]["elements/connectivity"];
335  n_bndry_conn.to_int_array(n_bndry_conn_conv);
336  bndry_indices = (n_bndry_conn_conv).value();
337 
338  }
339 
340  num_bndry_ele =
341  n_bndry_topo["elements/connectivity"].dtype().number_of_elements();
342  num_bndry_ele = num_bndry_ele / num_idxs_per_bndry_ele;
343  }
344  }
345  else
346  {
347  // Skipping Boundary Element Data
348  }
349 
350  const int *mesh_atts = NULL;
351  const int *bndry_atts = NULL;
352 
353  // These variables are used in debug code below.
354  // int num_mesh_atts_entires = 0;
355  // int num_bndry_atts_entires = 0;
356 
357  // the attribute fields could have several names
358  // for the element attributes check for first occurrence of field with
359  // name containing "_attribute", that doesn't contain "boundary"
360  std::string main_att_name = "";
361 
362  const Node &n_fields = n_mesh["fields"];
363  NodeConstIterator itr = n_fields.children();
364 
365  while ( itr.has_next() && main_att_name == "" )
366  {
367  itr.next();
368  std::string fld_name = itr.name();
369  if ( fld_name.find("boundary") == std::string::npos &&
370  fld_name.find("_attribute") != std::string::npos )
371  {
372  main_att_name = fld_name;
373  }
374  }
375 
376  if ( main_att_name != "" )
377  {
378  const Node &n_mesh_atts_vals = n_fields[main_att_name]["values"];
379 
380  // mfem requires ints, we could have int64s, etc convert if necessary
381  if (n_mesh_atts_vals.dtype().is_int() &&
382  n_mesh_atts_vals.is_compact() )
383  {
384  mesh_atts = n_mesh_atts_vals.value();
385  }
386  else
387  {
388  Node &n_mesh_atts_vals_conv = n_conv["fields"][main_att_name]["values"];
389  n_mesh_atts_vals.to_int_array(n_mesh_atts_vals_conv);
390  mesh_atts = n_mesh_atts_vals_conv.value();
391  }
392 
393  // num_mesh_atts_entires = n_mesh_atts_vals.dtype().number_of_elements();
394  }
395  else
396  {
397  // Skipping Mesh Attribute Data
398  }
399 
400  // for the boundary attributes check for first occurrence of field with
401  // name containing "_attribute", that also contains "boundary"
402  std::string bnd_att_name = "";
403  itr = n_fields.children();
404 
405  while ( itr.has_next() && bnd_att_name == "" )
406  {
407  itr.next();
408  std::string fld_name = itr.name();
409  if ( fld_name.find("boundary") != std::string::npos &&
410  fld_name.find("_attribute") != std::string::npos )
411  {
412  bnd_att_name = fld_name;
413  }
414  }
415 
416  if ( bnd_att_name != "" )
417  {
418  // Info: "Getting Boundary Attribute Data"
419  const Node &n_bndry_atts_vals =n_fields[bnd_att_name]["values"];
420 
421  // mfem requires ints, we could have int64s, etc convert if necessary
422  if ( n_bndry_atts_vals.dtype().is_int() &&
423  n_bndry_atts_vals.is_compact())
424  {
425  bndry_atts = n_bndry_atts_vals.value();
426  }
427  else
428  {
429  Node &n_bndry_atts_vals_conv = n_conv["fields"][bnd_att_name]["values"];
430  n_bndry_atts_vals.to_int_array(n_bndry_atts_vals_conv);
431  bndry_atts = n_bndry_atts_vals_conv.value();
432  }
433 
434  // num_bndry_atts_entires = n_bndry_atts_vals.dtype().number_of_elements();
435 
436  }
437  else
438  {
439  // Skipping Boundary Attribute Data
440  }
441 
442  // Info: "Number of Vertices: " << num_verts << endl
443  // << "Number of Mesh Elements: " << num_mesh_ele << endl
444  // << "Number of Boundary Elements: " << num_bndry_ele << endl
445  // << "Number of Mesh Attribute Entries: "
446  // << num_mesh_atts_entires << endl
447  // << "Number of Boundary Attribute Entries: "
448  // << num_bndry_atts_entires << endl);
449 
450  // Construct MFEM Mesh Object with externally owned data
451  // Note: if we don't have a gf, we need to provide the proper space dim
452  // if nodes gf is attached later, it resets the space dim based
453  // on the gf's fes.
454  Mesh *mesh = new Mesh(// from coordset
455  const_cast<double*>(verts_ptr),
456  num_verts,
457  // from topology
458  const_cast<int*>(elem_indices),
459  mesh_geo,
460  // from mesh_attribute field
461  const_cast<int*>(mesh_atts),
462  num_mesh_ele,
463  // from boundary topology
464  const_cast<int*>(bndry_indices),
465  bndry_geo,
466  // from boundary_attribute field
467  const_cast<int*>(bndry_atts),
468  num_bndry_ele,
469  ndims, // dim
470  ndims); // space dim
471 
472  // Attach Nodes Grid Function, if it exists
473  if (n_mesh_topo.has_child("grid_function"))
474  {
475  std::string nodes_gf_name = n_mesh_topo["grid_function"].as_string();
476 
477  // fetch blueprint field for the nodes gf
478  const Node &n_mesh_gf = n_mesh["fields"][nodes_gf_name];
479  // create gf
481  n_mesh_gf);
482  // attach to mesh
483  mesh->NewNodes(*nodes,true);
484  }
485 
486 
487  if (zero_copy && !n_conv.dtype().is_empty())
488  {
489  //Info: "Cannot zero-copy since data conversions were necessary"
490  zero_copy = false;
491  }
492 
493  Mesh *res = NULL;
494 
495  if (zero_copy)
496  {
497  res = mesh;
498  }
499  else
500  {
501  // the mesh above contains references to external data, to get a
502  // copy independent of the conduit data, we use:
503  res = new Mesh(*mesh,true);
504  delete mesh;
505  }
506 
507  return res;
508 }
509 
510 //---------------------------------------------------------------------------//
513  const Node &n_field,
514  bool zero_copy)
515 {
516  // n_conv holds converted data (when necessary for mfem api)
517  // if n_conv is used ( !n_conv.dtype().empty() ) we
518  // know that some data allocation was necessary, so we
519  // can't return a gf that zero copies the conduit data
520  Node n_conv;
521 
522  const double *vals_ptr = NULL;
523 
524  int vdim = 1;
525 
527 
528  if (n_field["values"].dtype().is_object())
529  {
530  vdim = n_field["values"].number_of_children();
531 
532  // need to check that we have doubles and
533  // cover supported layouts
534 
535  if ( n_field["values"][0].dtype().is_double() )
536  {
537  // check for contig
538  if (n_field["values"].is_contiguous())
539  {
540  // conduit mcarray contig == mfem byNODES
541  vals_ptr = n_field["values"].child(0).value();
542  }
543  // check for interleaved
544  else if (blueprint::mcarray::is_interleaved(n_field["values"]))
545  {
546  // conduit mcarray interleaved == mfem byVDIM
547  ordering = Ordering::byVDIM;
548  vals_ptr = n_field["values"].child(0).value();
549  }
550  else
551  {
552  // for mcarray generic case -- default to byNODES
553  // and provide values w/ contiguous (soa) ordering
554  blueprint::mcarray::to_contiguous(n_field["values"],
555  n_conv["values"]);
556  vals_ptr = n_conv["values"].child(0).value();
557  }
558  }
559  else // convert to doubles and use contig
560  {
561  Node n_tmp;
562  // check all vals, if we don't have doubles convert
563  // to doubles
564  NodeConstIterator itr = n_field["values"].children();
565  while (itr.has_next())
566  {
567  const Node &c_vals = itr.next();
568  std::string c_name = itr.name();
569 
570  if ( c_vals.dtype().is_double() )
571  {
572  // zero copy current coords
573  n_tmp[c_name].set_external(c_vals);
574 
575  }
576  else
577  {
578  // convert
579  c_vals.to_double_array(n_tmp[c_name]);
580  }
581  }
582 
583  // for mcarray generic case -- default to byNODES
584  // and provide values w/ contiguous (soa) ordering
585  blueprint::mcarray::to_contiguous(n_tmp,
586  n_conv["values"]);
587  vals_ptr = n_conv["values"].child(0).value();
588  }
589  }
590  else
591  {
592  if (n_field["values"].dtype().is_double() &&
593  n_field["values"].is_compact())
594  {
595  vals_ptr = n_field["values"].value();
596  }
597  else
598  {
599  n_field["values"].to_double_array(n_conv["values"]);
600  vals_ptr = n_conv["values"].value();
601  }
602  }
603 
604  if (zero_copy && !n_conv.dtype().is_empty())
605  {
606  //Info: "Cannot zero-copy since data conversions were necessary"
607  zero_copy = false;
608  }
609 
610  // we need basis name to create the proper mfem fec
611  std::string fec_name = n_field["basis"].as_string();
612 
613  GridFunction *res = NULL;
615  fec_name.c_str());
617  fec,
618  vdim,
619  ordering);
620 
621  if (zero_copy)
622  {
623  res = new GridFunction(fes,const_cast<double*>(vals_ptr));
624  }
625  else
626  {
627  // copy case, this constructor will alloc the space for the GF data
628  res = new GridFunction(fes);
629  // create an mfem vector that wraps the conduit data
630  Vector vals_vec(const_cast<double*>(vals_ptr),fes->GetVSize());
631  // copy values into the result
632  (*res) = vals_vec;
633  }
634 
635  // TODO: I believe the GF already has ownership of fes, so this should be all
636  // we need to do to avoid leaking objs created here?
637  res->MakeOwner(fec);
638 
639  return res;
640 }
641 
642 //---------------------------------------------------------------------------//
643 void
645  Node &n_mesh,
646  const std::string &coordset_name,
647  const std::string &main_topology_name,
648  const std::string &boundary_topology_name)
649 {
650  int dim = mesh->SpaceDimension();
651 
652  MFEM_ASSERT(dim >= 1 && dim <= 3, "invalid mesh dimension");
653 
654  ////////////////////////////////////////////
655  // Setup main coordset
656  ////////////////////////////////////////////
657 
658  // Assumes mfem::Vertex has the layout of a double array.
659 
660  // this logic assumes an mfem vertex is always 3 doubles wide
661  int stride = sizeof(mfem::Vertex);
662  int num_vertices = mesh->GetNV();
663 
664  MFEM_ASSERT( ( stride == 3 * sizeof(double) ),
665  "Unexpected stride for Vertex");
666 
667  Node &n_mesh_coords = n_mesh["coordsets"][coordset_name];
668  n_mesh_coords["type"] = "explicit";
669 
670 
671  double *coords_ptr = mesh->GetVertex(0);
672 
673  n_mesh_coords["values/x"].set_external(coords_ptr,
674  num_vertices,
675  0,
676  stride);
677 
678  if (dim >= 2)
679  {
680  n_mesh_coords["values/y"].set_external(coords_ptr,
681  num_vertices,
682  sizeof(double),
683  stride);
684  }
685  if (dim >= 3)
686  {
687  n_mesh_coords["values/z"].set_external(coords_ptr,
688  num_vertices,
689  sizeof(double) * 2,
690  stride);
691  }
692 
693  ////////////////////////////////////////////
694  // Setup main topo
695  ////////////////////////////////////////////
696 
697  Node &n_topo = n_mesh["topologies"][main_topology_name];
698 
699  n_topo["type"] = "unstructured";
700  n_topo["coordset"] = coordset_name;
701 
702  Element::Type ele_type = mesh->GetElementType(0);
703 
704  std::string ele_shape = ElementTypeToShapeName(ele_type);
705 
706  n_topo["elements/shape"] = ele_shape;
707 
708  GridFunction *gf_mesh_nodes = mesh->GetNodes();
709 
710  if (gf_mesh_nodes != NULL)
711  {
712  n_topo["grid_function"] = "mesh_nodes";
713  }
714 
715  // connectivity
716  // TODO: generic case, i don't think we can zero-copy (mfem allocs
717  // an array per element) so we alloc our own contig array and
718  // copy out. Some other cases (sidre) may actually have contig
719  // allocation but I am not sure how to detect this case from mfem
720  int num_ele = mesh->GetNE();
721  int geom = mesh->GetElementBaseGeometry(0);
722  int idxs_per_ele = Geometry::NumVerts[geom];
723  int num_conn_idxs = num_ele * idxs_per_ele;
724 
725  n_topo["elements/connectivity"].set(DataType::c_int(num_conn_idxs));
726 
727  int *conn_ptr = n_topo["elements/connectivity"].value();
728 
729  for (int i=0; i < num_ele; i++)
730  {
731  const Element *ele = mesh->GetElement(i);
732  const int *ele_verts = ele->GetVertices();
733 
734  memcpy(conn_ptr, ele_verts, idxs_per_ele * sizeof(int));
735 
736  conn_ptr += idxs_per_ele;
737  }
738 
739  if (gf_mesh_nodes != NULL)
740  {
741  GridFunctionToBlueprintField(gf_mesh_nodes,
742  n_mesh["fields/mesh_nodes"],
743  main_topology_name);
744  }
745 
746  ////////////////////////////////////////////
747  // Setup mesh attribute
748  ////////////////////////////////////////////
749 
750  Node &n_mesh_att = n_mesh["fields/element_attribute"];
751 
752  n_mesh_att["association"] = "element";
753  n_mesh_att["topology"] = main_topology_name;
754  n_mesh_att["values"].set(DataType::c_int(num_ele));
755 
756  int_array att_vals = n_mesh_att["values"].value();
757  for (int i = 0; i < num_ele; i++)
758  {
759  att_vals[i] = mesh->GetAttribute(i);
760  }
761 
762  ////////////////////////////////////////////
763  // Setup bndry topo "boundary"
764  ////////////////////////////////////////////
765 
766  // guard vs if we have boundary elements
767  if (mesh->GetNBE() > 0)
768  {
769  n_topo["boundary_topology"] = boundary_topology_name;
770 
771  Node &n_bndry_topo = n_mesh["topologies"][boundary_topology_name];
772 
773  n_bndry_topo["type"] = "unstructured";
774  n_bndry_topo["coordset"] = coordset_name;
775 
776  Element::Type bndry_ele_type = mesh->GetBdrElementType(0);
777 
778  std::string bndry_ele_shape = ElementTypeToShapeName(bndry_ele_type);
779 
780  n_bndry_topo["elements/shape"] = bndry_ele_shape;
781 
782 
783  int num_bndry_ele = mesh->GetNBE();
784  int bndry_geom = mesh->GetBdrElementBaseGeometry(0);
785  int bndry_idxs_per_ele = Geometry::NumVerts[bndry_geom];
786  int num_bndry_conn_idxs = num_bndry_ele * bndry_idxs_per_ele;
787 
788  n_bndry_topo["elements/connectivity"].set(DataType::c_int(num_bndry_conn_idxs));
789 
790  int *bndry_conn_ptr = n_bndry_topo["elements/connectivity"].value();
791 
792  for (int i=0; i < num_bndry_ele; i++)
793  {
794  const Element *bndry_ele = mesh->GetBdrElement(i);
795  const int *bndry_ele_verts = bndry_ele->GetVertices();
796 
797  memcpy(bndry_conn_ptr, bndry_ele_verts, bndry_idxs_per_ele * sizeof(int));
798 
799  bndry_conn_ptr += bndry_idxs_per_ele;
800  }
801 
802  ////////////////////////////////////////////
803  // Setup bndry mesh attribute
804  ////////////////////////////////////////////
805 
806  Node &n_bndry_mesh_att = n_mesh["fields/boundary_attribute"];
807 
808  n_bndry_mesh_att["association"] = "element";
809  n_bndry_mesh_att["topology"] = boundary_topology_name;
810  n_bndry_mesh_att["values"].set(DataType::c_int(num_bndry_ele));
811 
812  int_array bndry_att_vals = n_bndry_mesh_att["values"].value();
813  for (int i = 0; i < num_bndry_ele; i++)
814  {
815  bndry_att_vals[i] = mesh->GetBdrAttribute(i);
816  }
817  }
818 }
819 
820 //---------------------------------------------------------------------------//
821 void
823  Node &n_field,
824  const std::string &main_topology_name)
825 {
826  n_field["basis"] = gf->FESpace()->FEColl()->Name();
827  n_field["topology"] = main_topology_name;
828 
829  int vdim = gf->FESpace()->GetVDim();
830  int ndofs = gf->FESpace()->GetNDofs();
831 
832  if (vdim == 1) // scalar case
833  {
834  n_field["values"].set_external(gf->GetData(),
835  ndofs);
836  }
837  else // vector case
838  {
839  // deal with striding of all components
840 
841  Ordering::Type ordering = gf->FESpace()->GetOrdering();
842 
843  int entry_stride = (ordering == Ordering::byNODES ? 1 : vdim);
844  int vdim_stride = (ordering == Ordering::byNODES ? ndofs : 1);
845 
846  index_t offset = 0;
847  index_t stride = sizeof(double) * entry_stride;
848 
849  for (int d = 0; d < vdim; d++)
850  {
851  std::ostringstream oss;
852  oss << "v" << d;
853  std::string comp_name = oss.str();
854  n_field["values"][comp_name].set_external(gf->GetData(),
855  ndofs,
856  offset,
857  stride);
858  offset += sizeof(double) * vdim_stride;
859  }
860  }
861 
862 }
863 
864 //------------------------------
865 // end static public methods
866 //------------------------------
867 
868 //------------------------------
869 // end public methods
870 //------------------------------
871 
872 //------------------------------
873 // begin protected methods
874 //------------------------------
875 
876 //---------------------------------------------------------------------------//
877 std::string
879 {
880  std::string res = prefix_path + name + "_" +
882  ".root";
883  return res;
884 }
885 
886 //---------------------------------------------------------------------------//
887 std::string
889  const std::string &relay_protocol)
890 {
891  std::string res = prefix_path +
892  name +
893  "_" +
895  "/domain_" +
896  to_padded_string(domain_id, pad_digits_rank) +
897  "." +
899 
900  return res;
901 }
902 
903 //---------------------------------------------------------------------------//
904 std::string
906 {
907  std::string res = prefix_path +
908  name +
909  "_" +
911  return res;
912 }
913 
914 //---------------------------------------------------------------------------//
915 std::string
916 ConduitDataCollection::MeshFilePattern(const std::string &relay_protocol)
917 {
918  std::ostringstream oss;
919  oss << prefix_path
920  << name
921  << "_"
923  << "/domain_%0"
924  << pad_digits_rank
925  << "d."
926  << relay_protocol;
927 
928  return oss.str();
929 }
930 
931 
932 //---------------------------------------------------------------------------//
933 void
935  const Node &n_mesh,
936  const std::string &relay_protocol)
937 {
938  // default to json root file, except for hdf5 case
939  std::string root_proto = "json";
940 
941  if (relay_protocol == "hdf5")
942  {
943  root_proto = relay_protocol;
944  }
945 
946  Node n_root;
947  // create blueprint index
948  Node &n_bp_idx = n_root["blueprint_index"];
949 
950  blueprint::mesh::generate_index(n_mesh,
951  "",
952  num_domains,
953  n_bp_idx["mesh"]);
954 
955  // there are cases where the data backing the gf fields doesn't
956  // accurately represent the number of components in physical space,
957  // so we loop over all gfs and fix those that are incorrect
958 
960  for ( itr = field_map.begin(); itr != field_map.end(); itr++)
961  {
962  std::string gf_name = itr->first;
963  GridFunction *gf = itr->second;
964 
965  Node &idx_gf_ncomps = n_bp_idx["mesh/fields"][gf_name]["number_of_components"];
966  // check that the number_of_components in the index matches what we expect
967  // correct if necessary
968  if ( idx_gf_ncomps.to_int() != gf->VectorDim() )
969  {
970  idx_gf_ncomps = gf->VectorDim();
971  }
972  }
973  // add extra header info
974  n_root["protocol/name"] = relay_protocol;
975  n_root["protocol/version"] = "0.3.1";
976 
977 
978  // we will save one file per domain, so trees == files
979  n_root["number_of_files"] = num_domains;
980  n_root["number_of_trees"] = num_domains;
981  n_root["file_pattern"] = MeshFilePattern(relay_protocol);
982  n_root["tree_pattern"] = "";
983 
984  relay::io::save(n_root, RootFileName(), root_proto);
985 }
986 
987 //---------------------------------------------------------------------------//
988 void
990  const Node &n_mesh,
991  const std::string &relay_protocol)
992 {
993  relay::io::save(n_mesh, MeshFileName(domain_id, relay_protocol));
994 }
995 
996 //---------------------------------------------------------------------------//
997 void
999 {
1000  if (myid == 0)
1001  {
1002  // assume root file is json, unless hdf5 is specified
1003  std::string root_protocol = "json";
1004 
1005  if ( relay_protocol.find("hdf5") != std::string::npos )
1006  {
1007  root_protocol = "hdf5";
1008  }
1009 
1010 
1011  relay::io::load(RootFileName(), root_protocol, root_out);
1012 #ifdef MFEM_USE_MPI
1013  // broadcast contents of root file other ranks
1014  // (conduit relay mpi would simplify, but we would need to link another
1015  // lib for mpi case)
1016 
1017  // create json string
1018  std::string root_json = root_out.to_json();
1019  // string size +1 for null term
1020  int json_str_size = root_json.size() + 1;
1021 
1022  // broadcast json string buffer size
1023  int mpi_status = MPI_Bcast((void*)&json_str_size, // ptr
1024  1, // size
1025  MPI_INT, // type
1026  0, // root
1027  m_comm); // comm
1028 
1029  if (mpi_status != MPI_SUCCESS)
1030  {
1031  MFEM_ABORT("Broadcast of root file json string size failed");
1032  }
1033 
1034  // broadcast json string
1035  mpi_status = MPI_Bcast((void*)root_json.c_str(), // ptr
1036  json_str_size, // size
1037  MPI_CHAR, // type
1038  0, // root
1039  m_comm); // comm
1040 
1041  if (mpi_status != MPI_SUCCESS)
1042  {
1043  MFEM_ABORT("Broadcast of root file json string failed");
1044  }
1045 
1046 #endif
1047  }
1048 
1049 #ifdef MFEM_USE_MPI
1050  else
1051  {
1052  // recv json string buffer size via broadcast
1053  int json_str_size = -1;
1054  int mpi_status = MPI_Bcast(&json_str_size, // ptr
1055  1, // size
1056  MPI_INT, // type
1057  0, // root
1058  m_comm); // comm
1059 
1060  if (mpi_status != MPI_SUCCESS)
1061  {
1062  MFEM_ABORT("Broadcast of root file json string size failed");
1063  }
1064 
1065  // recv json string buffer via broadcast
1066  char *json_buff = new char[json_str_size];
1067  mpi_status = MPI_Bcast(json_buff, // ptr
1068  json_str_size, // size
1069  MPI_CHAR, // type
1070  0, // root
1071  m_comm); // comm
1072 
1073  if (mpi_status != MPI_SUCCESS)
1074  {
1075  MFEM_ABORT("Broadcast of root file json string failed");
1076  }
1077 
1078  // reconstruct root file contents
1079  Generator g(std::string(json_buff),"json");
1080  g.walk(root_out);
1081  // cleanup temp buffer
1082  delete [] json_buff;
1083  }
1084 #endif
1085 }
1086 
1087 //---------------------------------------------------------------------------//
1088 void
1090  const std::string &relay_protocol)
1091 {
1092  // Note: This path doesn't use any info from the root file
1093  // it uses the implicit mfem ConduitDataCollection layout
1094 
1095  Node n_mesh;
1096  relay::io::load( MeshFileName(domain_id, relay_protocol), n_mesh);
1097 
1098 
1099  Node verify_info;
1100  if (!blueprint::mesh::verify(n_mesh,verify_info))
1101  {
1102  MFEM_ABORT("Conduit Mesh Blueprint Verify Failed:\n"
1103  << verify_info.to_json());
1104  }
1105 
1106  mesh = BlueprintMeshToMesh(n_mesh);
1107 
1108  field_map.clear();
1109 
1110  NodeConstIterator itr = n_mesh["fields"].children();
1111 
1112  std::string nodes_gf_name = "";
1113 
1114  const Node &n_topo = n_mesh["topologies/main"];
1115  if (n_topo.has_child("grid_function"))
1116  {
1117  nodes_gf_name = n_topo["grid_function"].as_string();
1118  }
1119 
1120  while (itr.has_next())
1121  {
1122  const Node &n_field = itr.next();
1123  std::string field_name = itr.name();
1124 
1125  // skip mesh nodes gf since they are already processed
1126  // skip attribute fields, they aren't grid functions
1127  if ( field_name != nodes_gf_name &&
1128  field_name.find("_attribute") == std::string::npos
1129  )
1130  {
1132  field_map.Register(field_name, gf, true);
1133  }
1134  }
1135 }
1136 
1137 //------------------------------
1138 // end protected methods
1139 //------------------------------
1140 
1141 //------------------------------
1142 // begin static private methods
1143 //------------------------------
1144 
1145 //---------------------------------------------------------------------------//
1146 std::string
1147 ConduitDataCollection::ElementTypeToShapeName(Element::Type element_type)
1148 {
1149  // Adapted from SidreDataCollection
1150 
1151  // Note -- the mapping from Element::Type to string is based on
1152  // enum Element::Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
1153  // TETRAHEDRON, HEXAHEDRON };
1154  // Note: -- the string names are from conduit's blueprint
1155 
1156  switch (element_type)
1157  {
1158  case Element::POINT: return "point";
1159  case Element::SEGMENT: return "line";
1160  case Element::TRIANGLE: return "tri";
1161  case Element::QUADRILATERAL: return "quad";
1162  case Element::TETRAHEDRON: return "tet";
1163  case Element::HEXAHEDRON: return "hex";
1164  case Element::WEDGE:
1165  default: ;
1166  }
1167 
1168  return "unknown";
1169 }
1170 
1171 //---------------------------------------------------------------------------//
1173 ConduitDataCollection::ShapeNameToGeomType(const std::string &shape_name)
1174 {
1175  // Note: must init to something to avoid invalid memory access
1176  // in the mfem mesh constructor
1178 
1179  if (shape_name == "point")
1180  {
1181  res = mfem::Geometry::POINT;
1182  }
1183  else if (shape_name == "line")
1184  {
1186  }
1187  else if (shape_name == "tri")
1188  {
1190  }
1191  else if (shape_name == "quad")
1192  {
1193  res = mfem::Geometry::SQUARE;
1194  }
1195  else if (shape_name == "tet")
1196  {
1198  }
1199  else if (shape_name == "hex")
1200  {
1201  res = mfem::Geometry::CUBE;
1202  }
1203  else
1204  {
1205  MFEM_ABORT("Unsupported Element Shape: " << shape_name);
1206  }
1207 
1208  return res;
1209 }
1210 
1211 //------------------------------
1212 // end static private methods
1213 //------------------------------
1214 
1215 } // end namespace mfem
1216 
1217 #endif
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:412
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:397
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1053
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:794
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void SaveMeshAndFields(int domain_id, const conduit::Node &n_mesh, const std::string &file_protocol)
Saves all meshes and fields for the current cycle.
int error
Error state.
bool appendRankToFileName
Append rank to any output file names.
iterator begin()
Returns a begin iterator to the registered fields.
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
std::string RootFileName()
Returns blueprint root file name for the current cycle.
void DeleteAll()
Delete data owned by the DataCollection including field information.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:740
void MakeOwner(FiniteElementCollection *_fec)
Make the GridFunction the owner of fec and fes.
Definition: gridfunc.hpp:115
iterator end()
Returns an end iterator to the registered fields.
const Geometry::Type geom
Definition: ex1.cpp:40
int pad_digits_cycle
Number of digits used for the cycle and MPI rank in filenames.
void NewNodes(GridFunction &nodes, bool make_owner=false)
Replace the internal node GridFunction with the given GridFunction.
Definition: mesh.cpp:6627
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:5030
Element::Type GetBdrElementType(int i) const
Returns the type of boundary element i.
Definition: mesh.cpp:5035
int VectorDim() const
Definition: gridfunc.cpp:300
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
Data type for vertex.
Definition: vertex.hpp:22
bool own_data
Should the collection delete its mesh and fields.
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:169
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:834
virtual void Load(int cycle=0)
Load the collection based blueprint data.
ConduitDataCollection(const std::string &collection_name, Mesh *mesh=NULL)
Constructor. The collection name is used when saving the data.
MPI_Comm m_comm
Associated MPI communicator.
static Mesh * BlueprintMeshToMesh(const conduit::Node &n_mesh, const std::string &main_toplogy_name="", bool zero_copy=false)
Constructs and MFEM mesh from a Conduit Blueprint Description.
Mesh * mesh
The (common) mesh for the collected fields.
virtual ~ConduitDataCollection()
We will delete the mesh and fields if we own them.
int cycle
Time cycle; for time-dependent simulations cycle &gt;= 0, otherwise = -1.
std::string MeshFileName(int domain_id, const std::string &file_protocol="hdf5")
Returns the mesh file name for a given domain at the current cycle.
std::string to_padded_string(int i, int digits)
Convert an integer to a 0-padded string with the given number of digits.
Definition: text.hpp:63
Type
Ordering methods:
Definition: fespace.hpp:32
static const int NumVerts[NumGeom]
Definition: geom.hpp:48
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
static int create_directory(const std::string &dir_name, const Mesh *mesh, int myid)
virtual void Save()
Save the collection and a Conduit blueprint root file.
void Register(const std::string &fname, T *field, bool own_data)
Register field field with name fname.
const Element * GetElement(int i) const
Definition: mesh.hpp:819
void LoadRootFile(conduit::Node &n_root_out)
Loads contents of the root field for the current cycle into n_root_out.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:389
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:595
int SpaceDimension() const
Definition: mesh.hpp:789
void SaveRootFile(int num_domains, const conduit::Node &n_mesh, const std::string &file_protocol)
Saves root file for the current cycle.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
void LoadMeshAndFields(int domain_id, const std::string &file_protocol)
Loads all meshes and fields of a given domain id for the current cycle.
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
int myid
MPI rank (in parallel)
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:47
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:734
virtual const char * Name() const
Definition: fe_coll.hpp:61
GFieldMap::const_iterator FieldMapConstIterator
void clear()
Clears the map of registered fields without reclaiming memory.
std::string MeshDirectoryName()
Returns the mesh output directory for the current cycle.
int dim
Definition: ex24.cpp:53
std::string prefix_path
A path where the directory with results is saved. If not empty, it has &#39;/&#39; at the end...
Geometry::Type GetBdrElementBaseGeometry(int i) const
Definition: mesh.hpp:839
Vector data type.
Definition: vector.hpp:51
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:414
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6603
static void MeshToBlueprintMesh(Mesh *m, conduit::Node &out, const std::string &coordset_name="coords", const std::string &main_topology_name="main", const std::string &boundary_topology_name="boundary")
Describes a MFEM mesh using the mesh blueprint.
std::string name
Name of the collection, used as a directory name when saving.
std::string MeshFilePattern(const std::string &file_protocol="hdf5")
Returns the mesh file pattern for the current cycle.
static void GridFunctionToBlueprintField(GridFunction *gf, conduit::Node &out, const std::string &main_topology_name="main")
Describes a MFEM grid function using the mesh blueprint.
int num_procs
Number of MPI ranks (in parallel)
static GridFunction * BlueprintFieldToGridFunction(Mesh *mesh, const conduit::Node &n_field, bool zero_copy=false)
Constructs and MFEM Grid Function from a Conduit Blueprint Description.
Abstract data type element.
Definition: element.hpp:28
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1047
void SetProtocol(const std::string &protocol)
Set the Conduit relay i/o protocol to use.
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:823