MFEM  v4.6.0
Finite element discretization library
conduitdatacollection.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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  const std::string &main_adjset_name)
650 {
651  int dim = mesh->SpaceDimension();
652 
653  MFEM_ASSERT(dim >= 1 && dim <= 3, "invalid mesh dimension");
654 
655  ////////////////////////////////////////////
656  // Setup main coordset
657  ////////////////////////////////////////////
658 
659  // Assumes mfem::Vertex has the layout of a double array.
660 
661  // this logic assumes an mfem vertex is always 3 doubles wide
662  int stride = sizeof(mfem::Vertex);
663  int num_vertices = mesh->GetNV();
664 
665  MFEM_ASSERT( ( stride == 3 * sizeof(double) ),
666  "Unexpected stride for Vertex");
667 
668  Node &n_mesh_coords = n_mesh["coordsets"][coordset_name];
669  n_mesh_coords["type"] = "explicit";
670 
671 
672  double *coords_ptr = mesh->GetVertex(0);
673 
674  n_mesh_coords["values/x"].set_external(coords_ptr,
675  num_vertices,
676  0,
677  stride);
678 
679  if (dim >= 2)
680  {
681  n_mesh_coords["values/y"].set_external(coords_ptr,
682  num_vertices,
683  sizeof(double),
684  stride);
685  }
686  if (dim >= 3)
687  {
688  n_mesh_coords["values/z"].set_external(coords_ptr,
689  num_vertices,
690  sizeof(double) * 2,
691  stride);
692  }
693 
694  ////////////////////////////////////////////
695  // Setup main topo
696  ////////////////////////////////////////////
697 
698  Node &n_topo = n_mesh["topologies"][main_topology_name];
699 
700  n_topo["type"] = "unstructured";
701  n_topo["coordset"] = coordset_name;
702 
703  Element::Type ele_type = mesh->GetElementType(0);
704 
705  std::string ele_shape = ElementTypeToShapeName(ele_type);
706 
707  n_topo["elements/shape"] = ele_shape;
708 
709  GridFunction *gf_mesh_nodes = mesh->GetNodes();
710 
711  if (gf_mesh_nodes != NULL)
712  {
713  n_topo["grid_function"] = "mesh_nodes";
714  }
715 
716  // connectivity
717  // TODO: generic case, i don't think we can zero-copy (mfem allocs
718  // an array per element) so we alloc our own contig array and
719  // copy out. Some other cases (sidre) may actually have contig
720  // allocation but I am not sure how to detect this case from mfem
721  int num_ele = mesh->GetNE();
722  int geom = mesh->GetElementBaseGeometry(0);
723  int idxs_per_ele = Geometry::NumVerts[geom];
724  int num_conn_idxs = num_ele * idxs_per_ele;
725 
726  n_topo["elements/connectivity"].set(DataType::c_int(num_conn_idxs));
727 
728  int *conn_ptr = n_topo["elements/connectivity"].value();
729 
730  for (int i=0; i < num_ele; i++)
731  {
732  const Element *ele = mesh->GetElement(i);
733  const int *ele_verts = ele->GetVertices();
734 
735  memcpy(conn_ptr, ele_verts, idxs_per_ele * sizeof(int));
736 
737  conn_ptr += idxs_per_ele;
738  }
739 
740  if (gf_mesh_nodes != NULL)
741  {
742  GridFunctionToBlueprintField(gf_mesh_nodes,
743  n_mesh["fields/mesh_nodes"],
744  main_topology_name);
745  }
746 
747  ////////////////////////////////////////////
748  // Setup mesh attribute
749  ////////////////////////////////////////////
750 
751  Node &n_mesh_att = n_mesh["fields/element_attribute"];
752 
753  n_mesh_att["association"] = "element";
754  n_mesh_att["topology"] = main_topology_name;
755  n_mesh_att["values"].set(DataType::c_int(num_ele));
756 
757  int_array att_vals = n_mesh_att["values"].value();
758  for (int i = 0; i < num_ele; i++)
759  {
760  att_vals[i] = mesh->GetAttribute(i);
761  }
762 
763  ////////////////////////////////////////////
764  // Setup bndry topo "boundary"
765  ////////////////////////////////////////////
766 
767  // guard vs if we have boundary elements
768  if (mesh->HasBoundaryElements())
769  {
770  n_topo["boundary_topology"] = boundary_topology_name;
771 
772  Node &n_bndry_topo = n_mesh["topologies"][boundary_topology_name];
773 
774  n_bndry_topo["type"] = "unstructured";
775  n_bndry_topo["coordset"] = coordset_name;
776 
777  int num_bndry_ele = mesh->GetNBE();
778 
779  Element *BE0 = NULL; // representative boundary element
780  if (num_bndry_ele > 0) { BE0 = mesh->GetBdrElement(0); }
781 
782  // must initialize this to something, pick POINT if no boundary elements
783  Element::Type bndry_ele_type = (BE0) ? BE0->GetType() : Element::POINT;
784  std::string bndry_ele_shape = ElementTypeToShapeName(bndry_ele_type);
785  n_bndry_topo["elements/shape"] = bndry_ele_shape;
786 
787  // must initialize this to something, pick POINT if no boundary elements
788  int bndry_geom = (BE0) ? BE0->GetGeometryType() : Geometry::POINT;
789  int bndry_idxs_per_ele = Geometry::NumVerts[bndry_geom];
790  int num_bndry_conn_idxs = num_bndry_ele * bndry_idxs_per_ele;
791 
792  n_bndry_topo["elements/connectivity"].set(DataType::c_int(num_bndry_conn_idxs));
793 
794  int *bndry_conn_ptr = n_bndry_topo["elements/connectivity"].value();
795 
796  for (int i=0; i < num_bndry_ele; i++)
797  {
798  const Element *bndry_ele = mesh->GetBdrElement(i);
799  const int *bndry_ele_verts = bndry_ele->GetVertices();
800 
801  memcpy(bndry_conn_ptr, bndry_ele_verts, bndry_idxs_per_ele * sizeof(int));
802 
803  bndry_conn_ptr += bndry_idxs_per_ele;
804  }
805 
806  ////////////////////////////////////////////
807  // Setup bndry mesh attribute
808  ////////////////////////////////////////////
809 
810  Node &n_bndry_mesh_att = n_mesh["fields/boundary_attribute"];
811 
812  n_bndry_mesh_att["association"] = "element";
813  n_bndry_mesh_att["topology"] = boundary_topology_name;
814  n_bndry_mesh_att["values"].set(DataType::c_int(num_bndry_ele));
815 
816  int_array bndry_att_vals = n_bndry_mesh_att["values"].value();
817  for (int i = 0; i < num_bndry_ele; i++)
818  {
819  bndry_att_vals[i] = mesh->GetBdrAttribute(i);
820  }
821  }
822 
823  ////////////////////////////////////////////
824  // Setup adjsets
825  ////////////////////////////////////////////
826 
827 #ifdef MFEM_USE_MPI
828  ParMesh *pmesh = dynamic_cast<ParMesh*>(mesh);
829  if (pmesh)
830  {
831  ////////////////////////////////////////////
832  // Setup main adjset
833  ////////////////////////////////////////////
834 
835  Node &n_adjset = n_mesh["adjsets"][main_adjset_name];
836 
837  n_adjset["association"] = "vertex";
838  n_adjset["topology"] = main_topology_name;
839  n_adjset["groups"].set(DataType::object());
840 
841  const GroupTopology &pmesh_gtopo = pmesh->gtopo;
842  const int local_rank = pmesh->GetMyRank();
843  const int num_groups = pmesh_gtopo.NGroups();
844  // NOTE: skip the first group since its the local-only group
845  for (int i = 1; i < num_groups; i++)
846  {
847  const int num_group_nbrs = pmesh_gtopo.GetGroupSize(i);
848  const int *group_nbrs = pmesh_gtopo.GetGroup(i);
849  const int num_group_verts = pmesh->GroupNVertices(i);
850 
851  // NOTE: 'neighbor' values are local to this processor, but Blueprint
852  // expects global domain identifiers, so we collapse this layer of
853  // indirection
854  Array<int> group_ranks(num_group_nbrs);
855  std::string group_name = "group";
856  {
857  for (int j = 0; j < num_group_nbrs; j++)
858  {
859  group_ranks[j] = pmesh_gtopo.GetNeighborRank(group_nbrs[j]);
860  }
861  group_ranks.Sort();
862  for (int j = 0; j < num_group_nbrs; j++)
863  {
864  group_name += "_" + std::to_string(group_ranks[j]);
865  }
866 
867  // NOTE: Blueprint only wants remote ranks in its neighbor list,
868  // so we remove the local rank after the canonicalized Blueprint
869  // group name is formed
870  group_ranks.DeleteFirst(local_rank);
871  }
872  Node &n_group = n_adjset["groups"][group_name];
873 
874  n_group["neighbors"].set(group_ranks.GetData(), group_ranks.Size());
875  n_group["values"].set(DataType::c_int(num_group_verts));
876 
877  int_array group_vals = n_group["values"].value();
878  for (int j = 0; j < num_group_verts; j++)
879  {
880  group_vals[j] = pmesh->GroupVertex(i, j);
881  }
882  }
883 
884  // NOTE: We don't create an adjset for face neighbor data because
885  // these faces aren't listed in the 'boundary_topology_name' topology
886  // (this topology only covers the faces between 'main_topology_name'
887  // elements and void). To include a face neighbor data adjset, this
888  // function would need to export a topology with either (1) all faces
889  // in the mesh topology or (2) all boundary faces, including neighbors.
890 
891  ////////////////////////////////////////////
892  // Setup distributed state
893  ////////////////////////////////////////////
894 
895  Node &n_domid = n_mesh["state/domain_id"];
896  n_domid.set(local_rank);
897  }
898 #endif
899 }
900 
901 //---------------------------------------------------------------------------//
902 void
904  Node &n_field,
905  const std::string &main_topology_name)
906 {
907  n_field["basis"] = gf->FESpace()->FEColl()->Name();
908  n_field["topology"] = main_topology_name;
909 
910  int vdim = gf->FESpace()->GetVDim();
911  int ndofs = gf->FESpace()->GetNDofs();
912 
913  if (vdim == 1) // scalar case
914  {
915  n_field["values"].set_external(gf->GetData(),
916  ndofs);
917  }
918  else // vector case
919  {
920  // deal with striding of all components
921 
922  Ordering::Type ordering = gf->FESpace()->GetOrdering();
923 
924  int entry_stride = (ordering == Ordering::byNODES ? 1 : vdim);
925  int vdim_stride = (ordering == Ordering::byNODES ? ndofs : 1);
926 
927  index_t offset = 0;
928  index_t stride = sizeof(double) * entry_stride;
929 
930  for (int d = 0; d < vdim; d++)
931  {
932  std::ostringstream oss;
933  oss << "v" << d;
934  std::string comp_name = oss.str();
935  n_field["values"][comp_name].set_external(gf->GetData(),
936  ndofs,
937  offset,
938  stride);
939  offset += sizeof(double) * vdim_stride;
940  }
941  }
942 
943 }
944 
945 //------------------------------
946 // end static public methods
947 //------------------------------
948 
949 //------------------------------
950 // end public methods
951 //------------------------------
952 
953 //------------------------------
954 // begin protected methods
955 //------------------------------
956 
957 //---------------------------------------------------------------------------//
958 std::string
960 {
961  std::string res = prefix_path + name + "_" +
963  ".root";
964  return res;
965 }
966 
967 //---------------------------------------------------------------------------//
968 std::string
970  const std::string &relay_protocol)
971 {
972  std::string res = prefix_path +
973  name +
974  "_" +
976  "/domain_" +
977  to_padded_string(domain_id, pad_digits_rank) +
978  "." +
980 
981  return res;
982 }
983 
984 //---------------------------------------------------------------------------//
985 std::string
987 {
988  std::string res = prefix_path +
989  name +
990  "_" +
992  return res;
993 }
994 
995 //---------------------------------------------------------------------------//
996 std::string
997 ConduitDataCollection::MeshFilePattern(const std::string &relay_protocol)
998 {
999  std::ostringstream oss;
1000  oss << prefix_path
1001  << name
1002  << "_"
1004  << "/domain_%0"
1005  << pad_digits_rank
1006  << "d."
1007  << relay_protocol;
1008 
1009  return oss.str();
1010 }
1011 
1012 
1013 //---------------------------------------------------------------------------//
1014 void
1016  const Node &n_mesh,
1017  const std::string &relay_protocol)
1018 {
1019  // default to json root file, except for hdf5 case
1020  std::string root_proto = "json";
1021 
1022  if (relay_protocol == "hdf5")
1023  {
1024  root_proto = relay_protocol;
1025  }
1026 
1027  Node n_root;
1028  // create blueprint index
1029  Node &n_bp_idx = n_root["blueprint_index"];
1030 
1031  blueprint::mesh::generate_index(n_mesh,
1032  "",
1033  num_domains,
1034  n_bp_idx["mesh"]);
1035 
1036  // there are cases where the data backing the gf fields doesn't
1037  // accurately represent the number of components in physical space,
1038  // so we loop over all gfs and fix those that are incorrect
1039 
1041  for ( itr = field_map.begin(); itr != field_map.end(); itr++)
1042  {
1043  std::string gf_name = itr->first;
1044  GridFunction *gf = itr->second;
1045 
1046  Node &idx_gf_ncomps = n_bp_idx["mesh/fields"][gf_name]["number_of_components"];
1047  // check that the number_of_components in the index matches what we expect
1048  // correct if necessary
1049  if ( idx_gf_ncomps.to_int() != gf->VectorDim() )
1050  {
1051  idx_gf_ncomps = gf->VectorDim();
1052  }
1053  }
1054  // add extra header info
1055  n_root["protocol/name"] = relay_protocol;
1056  n_root["protocol/version"] = "0.3.1";
1057 
1058 
1059  // we will save one file per domain, so trees == files
1060  n_root["number_of_files"] = num_domains;
1061  n_root["number_of_trees"] = num_domains;
1062  n_root["file_pattern"] = MeshFilePattern(relay_protocol);
1063  n_root["tree_pattern"] = "";
1064 
1065  // Add the time, time step, and cycle
1066  n_root["blueprint_index/mesh/state/time"] = time;
1067  n_root["blueprint_index/mesh/state/time_step"] = time_step;
1068  n_root["blueprint_index/mesh/state/cycle"] = cycle;
1069 
1070  relay::io::save(n_root, RootFileName(), root_proto);
1071 }
1072 
1073 //---------------------------------------------------------------------------//
1074 void
1076  const Node &n_mesh,
1077  const std::string &relay_protocol)
1078 {
1079  relay::io::save(n_mesh, MeshFileName(domain_id, relay_protocol));
1080 }
1081 
1082 //---------------------------------------------------------------------------//
1083 void
1085 {
1086  if (myid == 0)
1087  {
1088  // assume root file is json, unless hdf5 is specified
1089  std::string root_protocol = "json";
1090 
1091  if ( relay_protocol.find("hdf5") != std::string::npos )
1092  {
1093  root_protocol = "hdf5";
1094  }
1095 
1096 
1097  relay::io::load(RootFileName(), root_protocol, root_out);
1098 #ifdef MFEM_USE_MPI
1099  // broadcast contents of root file other ranks
1100  // (conduit relay mpi would simplify, but we would need to link another
1101  // lib for mpi case)
1102 
1103  // create json string
1104  std::string root_json = root_out.to_json();
1105  // string size +1 for null term
1106  int json_str_size = root_json.size() + 1;
1107 
1108  // broadcast json string buffer size
1109  int mpi_status = MPI_Bcast((void*)&json_str_size, // ptr
1110  1, // size
1111  MPI_INT, // type
1112  0, // root
1113  m_comm); // comm
1114 
1115  if (mpi_status != MPI_SUCCESS)
1116  {
1117  MFEM_ABORT("Broadcast of root file json string size failed");
1118  }
1119 
1120  // broadcast json string
1121  mpi_status = MPI_Bcast((void*)root_json.c_str(), // ptr
1122  json_str_size, // size
1123  MPI_CHAR, // type
1124  0, // root
1125  m_comm); // comm
1126 
1127  if (mpi_status != MPI_SUCCESS)
1128  {
1129  MFEM_ABORT("Broadcast of root file json string failed");
1130  }
1131 
1132 #endif
1133  }
1134 
1135 #ifdef MFEM_USE_MPI
1136  else
1137  {
1138  // recv json string buffer size via broadcast
1139  int json_str_size = -1;
1140  int mpi_status = MPI_Bcast(&json_str_size, // ptr
1141  1, // size
1142  MPI_INT, // type
1143  0, // root
1144  m_comm); // comm
1145 
1146  if (mpi_status != MPI_SUCCESS)
1147  {
1148  MFEM_ABORT("Broadcast of root file json string size failed");
1149  }
1150 
1151  // recv json string buffer via broadcast
1152  char *json_buff = new char[json_str_size];
1153  mpi_status = MPI_Bcast(json_buff, // ptr
1154  json_str_size, // size
1155  MPI_CHAR, // type
1156  0, // root
1157  m_comm); // comm
1158 
1159  if (mpi_status != MPI_SUCCESS)
1160  {
1161  MFEM_ABORT("Broadcast of root file json string failed");
1162  }
1163 
1164  // reconstruct root file contents
1165  Generator g(std::string(json_buff),"json");
1166  g.walk(root_out);
1167  // cleanup temp buffer
1168  delete [] json_buff;
1169  }
1170 #endif
1171 }
1172 
1173 //---------------------------------------------------------------------------//
1174 void
1176  const std::string &relay_protocol)
1177 {
1178  // Note: This path doesn't use any info from the root file
1179  // it uses the implicit mfem ConduitDataCollection layout
1180 
1181  Node n_mesh;
1182  relay::io::load( MeshFileName(domain_id, relay_protocol), n_mesh);
1183 
1184 
1185  Node verify_info;
1186  if (!blueprint::mesh::verify(n_mesh,verify_info))
1187  {
1188  MFEM_ABORT("Conduit Mesh Blueprint Verify Failed:\n"
1189  << verify_info.to_json());
1190  }
1191 
1192  mesh = BlueprintMeshToMesh(n_mesh);
1193 
1194  field_map.clear();
1195 
1196  NodeConstIterator itr = n_mesh["fields"].children();
1197 
1198  std::string nodes_gf_name = "";
1199 
1200  const Node &n_topo = n_mesh["topologies/main"];
1201  if (n_topo.has_child("grid_function"))
1202  {
1203  nodes_gf_name = n_topo["grid_function"].as_string();
1204  }
1205 
1206  while (itr.has_next())
1207  {
1208  const Node &n_field = itr.next();
1209  std::string field_name = itr.name();
1210 
1211  // skip mesh nodes gf since they are already processed
1212  // skip attribute fields, they aren't grid functions
1213  if ( field_name != nodes_gf_name &&
1214  field_name.find("_attribute") == std::string::npos
1215  )
1216  {
1218  field_map.Register(field_name, gf, true);
1219  }
1220  }
1221 }
1222 
1223 //------------------------------
1224 // end protected methods
1225 //------------------------------
1226 
1227 //------------------------------
1228 // begin static private methods
1229 //------------------------------
1230 
1231 //---------------------------------------------------------------------------//
1232 std::string
1233 ConduitDataCollection::ElementTypeToShapeName(Element::Type element_type)
1234 {
1235  // Adapted from SidreDataCollection
1236 
1237  // Note -- the mapping from Element::Type to string is based on
1238  // enum Element::Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
1239  // TETRAHEDRON, HEXAHEDRON };
1240  // Note: -- the string names are from conduit's blueprint
1241 
1242  switch (element_type)
1243  {
1244  case Element::POINT: return "point";
1245  case Element::SEGMENT: return "line";
1246  case Element::TRIANGLE: return "tri";
1247  case Element::QUADRILATERAL: return "quad";
1248  case Element::TETRAHEDRON: return "tet";
1249  case Element::HEXAHEDRON: return "hex";
1250  case Element::WEDGE:
1251  default: ;
1252  }
1253 
1254  return "unknown";
1255 }
1256 
1257 //---------------------------------------------------------------------------//
1259 ConduitDataCollection::ShapeNameToGeomType(const std::string &shape_name)
1260 {
1261  // Note: must init to something to avoid invalid memory access
1262  // in the mfem mesh constructor
1264 
1265  if (shape_name == "point")
1266  {
1267  res = mfem::Geometry::POINT;
1268  }
1269  else if (shape_name == "line")
1270  {
1272  }
1273  else if (shape_name == "tri")
1274  {
1276  }
1277  else if (shape_name == "quad")
1278  {
1279  res = mfem::Geometry::SQUARE;
1280  }
1281  else if (shape_name == "tet")
1282  {
1284  }
1285  else if (shape_name == "hex")
1286  {
1287  res = mfem::Geometry::CUBE;
1288  }
1289  else
1290  {
1291  MFEM_ABORT("Unsupported Element Shape: " << shape_name);
1292  }
1293 
1294  return res;
1295 }
1296 
1297 //------------------------------
1298 // end static private methods
1299 //------------------------------
1300 
1301 } // end namespace mfem
1302 
1303 #endif
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.
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
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.
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1242
std::string RootFileName()
Returns blueprint root file name for the current cycle.
void DeleteAll()
Delete data owned by the DataCollection including field information.
iterator end()
Returns an end iterator to the registered fields.
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:8329
T * GetData()
Returns the data.
Definition: array.hpp:115
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:706
int NGroups() const
Return the number of groups.
int GetNeighborRank(int i) const
Return the MPI rank of neighbor &#39;i&#39;.
const Element * GetElement(int i) const
Return pointer to the i&#39;th element object.
Definition: mesh.hpp:1143
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", const std::string &main_adjset_name="main_adjset")
Describes a MFEM mesh using the mesh blueprint.
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
Definition: gridfunc.hpp:122
Data type for vertex.
Definition: vertex.hpp:22
bool own_data
Should the collection delete its mesh and fields.
void DeleteFirst(const T &el)
Delete the first entry with value == &#39;el&#39;.
Definition: array.hpp:837
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
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.
double time
Physical time (for time-dependent simulations)
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.
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:1083
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1190
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 >= 0, otherwise = -1.
Geometry::Type GetGeometryType() const
Definition: element.hpp:52
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.
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
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:54
int GroupNVertices(int group) const
Definition: pmesh.hpp:395
Type
Ordering methods:
Definition: fespace.hpp:33
static const int NumVerts[NumGeom]
Definition: geom.hpp:49
int GroupVertex(int group, int i) const
Definition: pmesh.hpp:400
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1125
virtual const char * Name() const
Definition: fe_coll.hpp:80
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.
const int * GetGroup(int g) const
Return a pointer to a list of neighbors for a given group. Neighbor 0 is the local processor...
void Register(const std::string &fname, T *field, bool own_data)
Register field field with name fname.
void LoadRootFile(conduit::Node &n_root_out)
Loads contents of the root field for the current cycle into n_root_out.
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
Definition: array.hpp:251
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
int GetMyRank() const
Definition: pmesh.hpp:353
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1196
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:219
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)
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:116
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:6644
int VectorDim() const
Definition: gridfunc.cpp:323
int GetGroupSize(int g) const
Get the number of processors in a group.
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.
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:721
int dim
Definition: ex24.cpp:53
virtual bool HasBoundaryElements() const
Checks if the mesh has boundary elements.
Definition: mesh.hpp:1050
std::string prefix_path
A path where the directory with results is saved. If not empty, it has &#39;/&#39; at the end...
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Vector data type.
Definition: vector.hpp:58
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
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.
const Element * GetBdrElement(int i) const
Return pointer to the i&#39;th boundary element object.
Definition: mesh.hpp:1158
static void GridFunctionToBlueprintField(GridFunction *gf, conduit::Node &out, const std::string &main_topology_name="main")
Describes a MFEM grid function using the mesh blueprint.
double time_step
Time step i.e. delta_t (for time-dependent simulations)
GroupTopology gtopo
Definition: pmesh.hpp:375
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.
Class for parallel meshes.
Definition: pmesh.hpp:32
Abstract data type element.
Definition: element.hpp:28
virtual Type GetType() const =0
Returns element&#39;s type.
void SetProtocol(const std::string &protocol)
Set the Conduit relay i/o protocol to use.