MFEM  v3.3
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
sidredatacollection.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 
13 #include "../config/config.hpp"
14 
15 #ifdef MFEM_USE_SIDRE
16 
17 #include "fem.hpp"
18 
19 #ifdef MFEM_USE_MPI
20 #include <spio/IOManager.hpp>
21 #endif
22 
23 #include <string>
24 #include <iomanip> // for setw, setfill
25 #include <cstdio> // for snprintf()
26 
27 namespace sidre = asctoolkit::sidre;
28 
29 namespace mfem
30 {
31 
32 // Constructor that will automatically create the sidre data store and necessary
33 // data groups for domain and global data.
34 SidreDataCollection::SidreDataCollection(const std::string& collection_name,
35  Mesh * the_mesh, bool own_mesh_data)
36  : mfem::DataCollection(collection_name, the_mesh),
37  m_owns_datastore(true),
38  m_owns_mesh_data(own_mesh_data),
39  m_meshNodesGFName("mesh_nodes")
40 {
41  m_datastore_ptr = new sidre::DataStore();
42 
43  sidre::DataGroup * global_grp =
44  m_datastore_ptr->getRoot()->createGroup(collection_name + "_global");
45  sidre::DataGroup * domain_grp =
46  m_datastore_ptr->getRoot()->createGroup(collection_name);
47 
48  bp_grp = domain_grp->createGroup("blueprint");
49  // Currently only rank 0 adds anything to bp_index.
50  bp_index_grp = global_grp->createGroup("blueprint_index/" + name);
51 
52  named_bufs_grp = domain_grp->createGroup("named_buffers");
53 
54  if (the_mesh)
55  {
56  SetMesh(the_mesh);
57  }
58 #ifdef MFEM_USE_MPI
59  else
60  {
61  m_comm = MPI_COMM_NULL;
62  }
63 #endif
64 }
65 
66 // Second constructor that allows external code to specify data groups to place
67 // domain and global data in.
68 
69 // TODO - Conduit will have the capability to generate a blueprint index group
70 // in the future. When this is available, all the blueprint index code can be
71 // removed from the data collection class.
72 SidreDataCollection::SidreDataCollection(const std::string& collection_name,
73  asctoolkit::sidre::DataGroup* global_grp,
74  asctoolkit::sidre::DataGroup* domain_grp,
75  bool own_mesh_data)
76  : mfem::DataCollection(collection_name),
77  m_owns_datastore(false),
78  m_owns_mesh_data(own_mesh_data),
79  m_meshNodesGFName("mesh_nodes"),
80  m_datastore_ptr(NULL)
81 {
82  bp_grp = domain_grp->createGroup("blueprint");
83 
84  // Currently only rank 0 adds anything to bp_index.
85  bp_index_grp = global_grp->createGroup("blueprint_index/" + name);
86 
87  named_bufs_grp = domain_grp->createGroup("named_buffers");
88 
89 #ifdef MFEM_USE_MPI
90  m_comm = MPI_COMM_NULL;
91 #endif
92 }
93 
95 {
96  if (m_owns_datastore)
97  {
98  delete m_datastore_ptr;
99  }
100 }
101 
102 #ifdef MFEM_USE_MPI
103 void SidreDataCollection::SetComm(MPI_Comm comm)
104 {
105  m_comm = comm;
106  serial = false;
107  appendRankToFileName = true;
108  MPI_Comm_rank(m_comm, &myid);
109  MPI_Comm_size(m_comm, &num_procs);
110 }
111 #endif
112 
113 // protected method
114 sidre::DataGroup *SidreDataCollection::named_buffers_grp() const
115 {
116  MFEM_ASSERT(named_bufs_grp != NULL,
117  "No group 'named_buffers' in data collection. Verify that"
118  " SetMesh was called to set the mesh in the data collection.");
119  return named_bufs_grp;
120 }
121 
122 // protected method
123 asctoolkit::sidre::DataView *
124 SidreDataCollection::alloc_view(asctoolkit::sidre::DataGroup *grp,
125  const std::string &view_name)
126 {
127  MFEM_ASSERT(grp, "DataGroup pointer is NULL");
128  sidre::DataView *v = grp->getView(view_name);
129  if (!v)
130  {
131  v = grp->createView(view_name);
132  MFEM_ASSERT(v, "error allocating DataView " << view_name
133  << " in group " << grp->getPathName());
134  }
135  return v;
136 }
137 
138 // protected method
139 asctoolkit::sidre::DataView *
140 SidreDataCollection::alloc_view(asctoolkit::sidre::DataGroup *grp,
141  const std::string &view_name,
142  const asctoolkit::sidre::DataType &dtype)
143 {
144  MFEM_ASSERT(grp, "DataGroup pointer is NULL");
145  sidre::DataView *v = grp->getView(view_name);
146  if (!v)
147  {
148  v = grp->createView(view_name, dtype);
149  MFEM_ASSERT(v, "error allocating DataView " << view_name
150  << " in group " << grp->getPathName());
151  }
152  else
153  {
154  MFEM_ASSERT(v->getSchema().dtype().equals(dtype), "");
155  }
156  return v;
157 }
158 
159 // protected method
160 asctoolkit::sidre::DataGroup *
161 SidreDataCollection::alloc_group(asctoolkit::sidre::DataGroup *grp,
162  const std::string &group_name)
163 {
164  MFEM_ASSERT(grp, "DataGroup pointer is NULL");
165  sidre::DataGroup *g = grp->getGroup(group_name);
166  if (!g)
167  {
168  g = grp->createGroup(group_name);
169  MFEM_ASSERT(g, "error allocating DataGroup " << group_name
170  << " in group " << grp->getPathName());
171  }
172  return g;
173 }
174 
175 // protected method
176 std::string
177 SidreDataCollection::get_file_path(const std::string &filename) const
178 {
179  std::stringstream fNameSstr;
180 
181  // Note: If non-empty, prefix_path has a separator ('/') at the end
182  fNameSstr << prefix_path << filename;
183 
184  if (GetCycle() >= 0)
185  {
186  fNameSstr << "_" << std::setfill('0') << std::setw(pad_digits)
187  << GetCycle();
188  }
189 
190  return fNameSstr.str();
191 }
192 
193 asctoolkit::sidre::DataView *
194 SidreDataCollection::AllocNamedBuffer(const std::string& buffer_name,
195  asctoolkit::sidre::SidreLength sz,
196  asctoolkit::sidre::TypeID type)
197 {
198  sz = std::max(sz, sidre::SidreLength(0));
199  sidre::DataGroup *f = named_buffers_grp();
200  sidre::DataView *v = f->getView(buffer_name);
201  if ( v == NULL )
202  {
203  // create a buffer view
204  v = f->createViewAndAllocate(buffer_name, type, sz);
205  }
206  else
207  {
208  MFEM_ASSERT(v->getTypeID() == type, "type does not match existing type");
209 
210  // Here v is the view holding the buffer in the named_buffers group, so
211  // its size is the full size of the buffer.
212 
213  // check if we need to resize.
214  if (!v->isApplied() || v->getNumElements() < sz)
215  {
216  // resize, even if the buffer has more than 1 DataView.
217  // v->reallocate(sz); // this will not work for more than 1 view.
218  sidre::DataType dtype(v->getSchema().dtype());
219  dtype.set_number_of_elements(sz);
220  f->destroyViewAndData(buffer_name);
221  v = f->createViewAndAllocate(buffer_name, dtype);
222  }
223  }
224  MFEM_ASSERT(v && v->isApplied(), "allocation failed");
225  return v;
226 }
227 
228 // private method
229 void SidreDataCollection::createMeshBlueprintStubs(bool hasBP)
230 {
231  if (!hasBP)
232  {
233  bp_grp->createGroup("state");
234  bp_grp->createGroup("coordsets");
235  bp_grp->createGroup("topologies");
236  bp_grp->createGroup("fields");
237  }
238 
239  // If rank is 0, set up blueprint index state group.
240  if (myid == 0)
241  {
242  bp_index_grp->createGroup("state");
243  bp_index_grp->createGroup("coordsets");
244  bp_index_grp->createGroup("topologies");
245  bp_index_grp->createGroup("fields");
246  }
247 }
248 
249 // private method
250 void SidreDataCollection::createMeshBlueprintState(bool hasBP)
251 {
252  if (!hasBP)
253  {
254  // Set up blueprint state group.
255  bp_grp->createViewScalar("state/cycle", 0);
256  bp_grp->createViewScalar("state/time", 0.);
257  bp_grp->createViewScalar("state/domain", myid);
258  bp_grp->createViewScalar("state/time_step", 0.);
259  }
260 
261  // If rank is 0, set up blueprint index state group.
262  if (myid == 0)
263  {
264  bp_index_grp->createViewScalar("state/cycle", 0);
265  bp_index_grp->createViewScalar("state/time", 0.);
266  bp_index_grp->createViewScalar("state/number_of_domains", num_procs);
267  }
268 }
269 
270 // private method
271 void SidreDataCollection::createMeshBlueprintCoordset(bool hasBP)
272 {
273  int dim = mesh->SpaceDimension();
274  MFEM_ASSERT(dim >= 1 && dim <= 3, "invalid mesh dimension");
275 
276  // Assuming mfem::Vertex has the layout of a double array.
277  const int NUM_COORDS = sizeof(mfem::Vertex) / sizeof(double);
278 
279  const int num_vertices = mesh->GetNV();
280  const int coordset_len = NUM_COORDS * num_vertices;
281 
282  // Add blueprint if not present
283  if ( !hasBP )
284  {
285  bp_grp->createViewString("coordsets/coords/type", "explicit");
286 
287  sidre::DataType dtype =
288  sidre::DataType::c_double(num_vertices);
289  const size_t stride = dtype.stride();
290  dtype.set_stride(stride*NUM_COORDS);
291 
292  // Set up views for x, y, z values
293  sidre::DataView *vx, *vy = NULL, *vz = NULL;
294  vx = bp_grp->createView("coordsets/coords/values/x", dtype);
295 
296  if (dim >= 2)
297  {
298  dtype.set_offset(dtype.offset() + stride);
299  vy = bp_grp->createView("coordsets/coords/values/y", dtype);
300  }
301  if (dim >= 3)
302  {
303  dtype.set_offset(dtype.offset() + stride);
304  vz = bp_grp->createView("coordsets/coords/values/z", dtype);
305  }
306 
307  if (m_owns_mesh_data)
308  {
309  // Allocate buffer for coord values.
310  sidre::DataBuffer* coordbuf =
311  AllocNamedBuffer("vertex_coords", coordset_len)->getBuffer();
312 
313  vx->attachBuffer(coordbuf);
314  if (dim >= 2) { vy->attachBuffer(coordbuf); }
315  if (dim >= 3) { vz->attachBuffer(coordbuf); }
316  }
317  else
318  {
319  double *coordbuf = mesh->GetVertex(0);
320 
321  vx->setExternalDataPtr(coordbuf);
322  if (dim >= 2) { vy->setExternalDataPtr(coordbuf); }
323  if (dim >= 3) { vz->setExternalDataPtr(coordbuf); }
324  }
325 
326  }
327 
328  // If rank 0, set up blueprint index for coordinate set.
329  if (myid == 0)
330  {
331  bp_index_grp->createViewString(
332  "coordsets/coords/path", bp_grp->getPathName() + "/coordsets/coords");
333 
334  bp_index_grp->getGroup("coordsets/coords")->copyView(
335  bp_grp->getView("coordsets/coords/type") );
336 
337  bp_index_grp->createViewString(
338  "coordsets/coords/coord_system/type", "cartesian");
339 
340  // These are empty views, their existence in the group tree is used to
341  // define the number of dims
342  bp_index_grp->createView("coordsets/coords/coord_system/axes/x");
343 
344  if (dim >= 2)
345  {
346  bp_index_grp->createView("coordsets/coords/coord_system/axes/y");
347  }
348 
349  if (dim == 3)
350  {
351  bp_index_grp->createView("coordsets/coords/coord_system/axes/z");
352  }
353  }
354 
355  if (m_owns_mesh_data)
356  {
357  double *coord_values = GetNamedBuffer("vertex_coords")->getData();
358  // Change ownership of the mesh vertex data to sidre
359  mesh->ChangeVertexDataOwnership(coord_values, coordset_len, hasBP);
360  }
361 }
362 
363 // private method
364 void SidreDataCollection::
365 createMeshBlueprintTopologies(bool hasBP, const std::string& mesh_name)
366 {
367  const bool isBdry = (mesh_name == "boundary");
368 
369  const int num_elements = !isBdry
370  ? mesh->GetNE()
371  : mesh->GetNBE();
372 
373  MFEM_VERIFY(num_elements > 0,
374  "TODO: processors with 0 " << mesh_name << " elements");
375 
376  const int element_size = !isBdry
377  ? mesh->GetElement(0)->GetNVertices()
379 
380  const int num_indices = num_elements * element_size;
381 
382  // Find the element shape
383  // Note: Assumes homogeneous elements, so only check the first element
384  const int geom =
385  isBdry ?
388  const std::string eltTypeStr =
389  !isBdry
390  ? getElementName( static_cast<Element::Type>(
391  mesh->GetElement(0)->GetType() ) )
392  : getElementName( static_cast<Element::Type>(
393  mesh->GetBdrElement(0)->GetType() ) );
394 
395  const std::string mesh_topo_str = "topologies/" + mesh_name;
396  const std::string mesh_attr_str = "fields/"+mesh_name+"_material_attribute";
397 
398  if ( !hasBP )
399  {
400  sidre::DataGroup* topology_grp = bp_grp->createGroup(mesh_topo_str);
401 
402  // Add mesh topology
403  topology_grp->createViewString("type", "unstructured");
404  // Note: eltTypeStr comes form the mesh
405  topology_grp->createViewString("elements/shape", eltTypeStr);
406  topology_grp->createViewAndAllocate(
407  "elements/connectivity", sidre::INT_ID, num_indices);
408  topology_grp->createViewString("coordset", "coords");
409 
410  // If the mesh has nodes, set the name of the GridFunction holding the
411  // mesh nodes in the blueprint group.
412  if (!isBdry && mesh->GetNodes() != NULL)
413  {
414  topology_grp->createViewString("grid_function",m_meshNodesGFName);
415  }
416 
417  // Add material attribute field to blueprint
418  sidre::DataGroup* attr_grp = bp_grp->createGroup(mesh_attr_str);
419  attr_grp->createViewString("association", "element");
420  attr_grp->createViewAndAllocate("values", sidre::INT_ID, num_elements);
421  attr_grp->createViewString("topology", mesh_name);
422  }
423 
424  // If rank 0, set up blueprint index for topologies group and material
425  // attribute field.
426  if (myid == 0)
427  {
428  const std::string bp_grp_path = bp_grp->getPathName();
429 
430  // Create blueprint index for topologies.
431  if (isBdry)
432  {
433  // "Shallow" copy the bp_grp view into the bp_index_grp sub-group.
434  // Note that the "topologies/mesh" sub-group has to exist, i.e. this
435  // method should be called first with mesh_name = "mesh".
436  bp_index_grp->getGroup("topologies/mesh")
437  ->copyView( bp_grp->getView("topologies/mesh/boundary_topology") );
438  }
439 
440  sidre::DataGroup *bp_index_topo_grp =
441  bp_index_grp->createGroup(mesh_topo_str);
442  sidre::DataGroup *topology_grp = bp_grp->getGroup(mesh_topo_str);
443 
444  bp_index_topo_grp->createViewString(
445  "path", bp_grp_path + "/" + mesh_topo_str);
446  bp_index_topo_grp->copyView( topology_grp->getView("type") );
447  bp_index_topo_grp->copyView( topology_grp->getView("coordset") );
448 
449  // If the mesh has nodes, set the name of the GridFunction holding the
450  // mesh nodes in the blueprint_index group.
451  if (!isBdry && mesh->GetNodes() != NULL)
452  {
453  bp_index_topo_grp->copyView(topology_grp->getView("grid_function"));
454  }
455 
456  // Create blueprint index for material attributes.
457  sidre::DataGroup *bp_index_attr_grp =
458  bp_index_grp->createGroup(mesh_attr_str);
459  sidre::DataGroup *attr_grp = bp_grp->getGroup(mesh_attr_str);
460 
461  bp_index_attr_grp->createViewString(
462  "path", bp_grp_path + "/" + mesh_attr_str );
463  bp_index_attr_grp->copyView( attr_grp->getView("association") );
464  bp_index_attr_grp->copyView( attr_grp->getView("topology") );
465 
466  int number_of_components = 1;
467  bp_index_attr_grp->createViewScalar("number_of_components",
468  number_of_components);
469  }
470 
471  // Finally, change ownership or copy the element arrays into Sidre
472  sidre::DataView* conn_view =
473  bp_grp->getGroup(mesh_topo_str)->getView("elements/connectivity");
474  sidre::DataView* attr_view =
475  bp_grp->getGroup(mesh_attr_str)->getView("values");
476  // The SidreDataCollection always owns these arrays:
477  Array<int> conn_array(conn_view->getData<int*>(), num_indices);
478  Array<int> attr_array(attr_view->getData<int*>(), num_elements);
479  if (!isBdry)
480  {
481  mesh->GetElementData(geom, conn_array, attr_array);
482  }
483  else
484  {
485  mesh->GetBdrElementData(geom, conn_array, attr_array);
486  }
487  MFEM_ASSERT(!conn_array.OwnsData(), "");
488  MFEM_ASSERT(!attr_array.OwnsData(), "");
489 }
490 
491 // private method
492 void SidreDataCollection::verifyMeshBlueprint()
493 {
494  // Conduit will have a verify mesh blueprint capability in the future.
495  // Add call to that when it's available to check actual contents in sidre.
496 }
497 
499 {
500  DataCollection::SetMesh(new_mesh);
501 
502 #ifdef MFEM_USE_MPI
503  ParMesh *pmesh = dynamic_cast<ParMesh*>(new_mesh);
504  m_comm = pmesh ? pmesh->GetComm() : MPI_COMM_NULL;
505 #endif
506 
507  // hasBP is used to indicate if the data currently in the blueprint should be
508  // used to replace the data in the mesh.
509  bool hasBP = bp_grp->getNumViews() > 0 || bp_grp->getNumGroups() > 0;
510  bool has_bnd_elts = (new_mesh->GetNBE() > 0);
511 
512  createMeshBlueprintStubs(hasBP);
513  createMeshBlueprintState(hasBP);
514  createMeshBlueprintCoordset(hasBP);
515 
516  GridFunction *nodes = new_mesh->GetNodes();
517 
518  // register the "mesh" topology in the blueprint.
519  createMeshBlueprintTopologies(hasBP, "mesh");
520 
521  if (has_bnd_elts)
522  {
523  // Set the "boundary_topology" of "mesh" to "boundary".
524  bp_grp->createViewString("topologies/mesh/boundary_topology", "boundary");
525 
526  // register the "boundary" topology in the blueprint.
527  createMeshBlueprintTopologies(hasBP, "boundary");
528  }
529 
530  if (nodes)
531  {
532  // See the comment at the definition of 'hasBP' above.
533  if (hasBP)
534  {
535  // Get the bp mesh nodes name.
536  sidre::DataView *v_bp_nodes_name =
537  bp_grp->getView("topologies/mesh/grid_function");
538  std::string bp_nodes_name(v_bp_nodes_name->getString());
539 
540  // Check that the names match, e.g. when loading the collection.
541  MFEM_VERIFY(m_meshNodesGFName == bp_nodes_name,
542  "mismatch of requested and blueprint mesh nodes names");
543  // Support renaming bp_nodes_name --> m_meshNodesGFName ?
544  }
545 
546  if (m_owns_mesh_data)
547  {
548  // Make sure Sidre owns the data of the new_mesh's Nodes.
549  if (!GetNamedBuffer(m_meshNodesGFName))
550  {
551  int sz = new_mesh->GetNodalFESpace()->GetVSize();
552  double *gfData = AllocNamedBuffer(m_meshNodesGFName, sz)->getData();
553 
554  // See the comment at the definition of 'hasBP' above.
555  if (!hasBP)
556  {
557  MFEM_ASSERT(nodes->Size() == sz, "");
558  std::memcpy(gfData, nodes->GetData(), sizeof(double) * sz);
559  }
560  }
561  // Since the named buffer for m_meshNodesGFName exists, the call to
562  // RegisterField() below will replace the data of the nodes with the
563  // data from the named buffer.
564  }
565  else
566  {
567  // Make sure Sidre does not have a named buffer for m_meshNodesGFName.
568  MFEM_VERIFY(GetNamedBuffer(m_meshNodesGFName) == NULL, "");
569  }
570 
571  RegisterField(m_meshNodesGFName, nodes);
572 
573  if (own_data)
574  {
575  // Avoid double delete calls (for the nodes gf) when (own_data == true)
576  // and the new_mesh owns its Nodes --> take ownership from new_mesh.
577  // When new_mesh does not own its Nodes and (own_data == true), we can
578  // not take ownership --> verify that does not happen.
579  MFEM_VERIFY(new_mesh->OwnsNodes(), "mesh does not own its nodes, "
580  "can not take ownership");
581  new_mesh->SetNodesOwner(false);
582  }
583  }
584 }
585 
587 SetGroupPointers(asctoolkit::sidre::DataGroup *global_grp,
588  asctoolkit::sidre::DataGroup *domain_grp)
589 {
590  MFEM_VERIFY(domain_grp->hasGroup("blueprint"),
591  "Domain group does not contain a blueprint group.");
592  MFEM_VERIFY(global_grp->hasGroup("blueprint_index/" + name),
593  "Global group does not contain a blueprint index group.");
594 
595  bp_grp = domain_grp->getGroup("blueprint");
596  bp_index_grp = global_grp->getGroup("blueprint_index/" + name);
597  named_bufs_grp = domain_grp->getGroup("named_buffers");
598 }
599 
600 void SidreDataCollection::Load(const std::string& path,
601  const std::string& protocol)
602 {
604  // Reset DataStore?
605 
606 #ifdef MFEM_USE_MPI
607  if (m_comm != MPI_COMM_NULL)
608  {
609  asctoolkit::spio::IOManager reader(m_comm);
610  reader.read(bp_grp->getDataStore()->getRoot(), path);
611  }
612  else
613 #endif
614  {
615  bp_grp->load(path, protocol);
616  }
617 
618  // If the data collection created the datastore, it knows the layout of where
619  // the domain and global groups are, and can restore them after the Load().
620  //
621  // If the data collection did not create the datastore, the host code must
622  // reset these pointers after the load operation and also reset the state
623  // variables.
624  if (m_owns_datastore)
625  {
626  SetGroupPointers(m_datastore_ptr->getRoot()->getGroup(name + "_global"),
627  m_datastore_ptr->getRoot()->getGroup(name));
628 
630  }
631 }
632 
633 void SidreDataCollection::LoadExternalData(const std::string& path)
634 {
635 #ifdef MFEM_USE_MPI
636  if (m_comm != MPI_COMM_NULL)
637  {
638  asctoolkit::spio::IOManager reader(m_comm);
639  reader.loadExternalData(bp_grp->getDataStore()->getRoot(), path);
640  }
641  else
642 #endif
643  {
644  bp_grp->loadExternalData(path);
645  }
646 }
647 
649 {
650  SetTime( bp_grp->getView("state/time")->getData<double>() );
651  SetCycle( bp_grp->getView("state/cycle")->getData<int>() );
652  SetTimeStep( bp_grp->getView("state/time_step")->getData<double>() );
653 }
654 
656 {
657  bp_grp->getView("state/cycle")->setScalar(GetCycle());
658  bp_grp->getView("state/time")->setScalar(GetTime());
659  bp_grp->getView("state/time_step")->setScalar(GetTimeStep());
660 
661  if (myid == 0)
662  {
663  bp_index_grp->getView("state/cycle")->setScalar(GetCycle());
664  bp_index_grp->getView("state/time")->setScalar(time);
665  }
666 }
667 
669 {
670  verifyMeshBlueprint();
671  UpdateStateToDS();
672 }
673 
675 {
676  std::string filename = name;
677  std::string protocol = "sidre_hdf5";
678 
679  Save(filename, protocol);
680 }
681 
682 void SidreDataCollection::Save(const std::string& filename,
683  const std::string& protocol)
684 {
685  PrepareToSave();
686 
688 
689  std::string file_path = get_file_path(filename);
690 
691  sidre::DataGroup * blueprint_indicies_grp = bp_index_grp->getParent();
692 #ifdef MFEM_USE_MPI
693  if (m_comm != MPI_COMM_NULL)
694  {
695  asctoolkit::spio::IOManager writer(m_comm);
696  sidre::DataStore *datastore = bp_grp->getDataStore();
697  writer.write(datastore->getRoot(), num_procs, file_path, protocol);
698  if (myid == 0)
699  {
700  if (protocol == "sidre_hdf5")
701  {
702  writer.writeGroupToRootFile(blueprint_indicies_grp,
703  file_path + ".root");
704  }
705  // Root file support only available in hdf5.
706  else
707  {
708  writer.write(blueprint_indicies_grp, 1,
709  file_path + ".root", protocol);
710  }
711  }
712  }
713  else
714 #endif
715  {
716  // If serial, use sidre group writer.
717  bp_grp->save(file_path, protocol);
718 
719  blueprint_indicies_grp
720  ->save(file_path + ".root", protocol);
721  }
722 }
723 
724 // private method
725 void SidreDataCollection::
726 addScalarBasedGridFunction(const std::string &field_name, GridFunction *gf,
727  const std::string &buffer_name,
728  asctoolkit::sidre::SidreLength offset)
729 {
730  sidre::DataGroup* grp = bp_grp->getGroup("fields/" + field_name);
731  MFEM_ASSERT(grp != NULL, "field " << field_name << " does not exist");
732 
733  const int numDofs = gf->FESpace()->GetVSize();
734 
735  if (gf->GetData() == NULL)
736  {
737  AllocNamedBuffer(buffer_name, offset + numDofs);
738  // gf->data is set below.
739  }
740 
741  /*
742  * Mesh blueprint for a scalar-based grid function is of the form
743  * /fields/field_name/basis
744  * -- string value is GridFunction's FEC::Name
745  * /fields/field_name/values
746  * -- array of size numDofs
747  */
748 
749  // Make sure we have the DataView "values".
750  sidre::DataView *vv = alloc_view(grp, "values");
751 
752  // Describe and apply the "values" DataView.
753  // If the data store has buffer for field_name (e.g. AllocNamedBuffer was
754  // called, or it was loaded from file), use that buffer.
755  sidre::DataView *bv = named_buffers_grp()->getView(buffer_name);
756  if (bv)
757  {
758  MFEM_ASSERT(bv->hasBuffer() && bv->isDescribed(), "");
759 
760  // named buffers always have offset 0
761  MFEM_ASSERT(bv->getSchema().dtype().offset() == 0, "");
762  MFEM_ASSERT(bv->getNumElements() >= offset + numDofs, "");
763 
764  if (vv->isEmpty())
765  {
766  vv->attachBuffer(bv->getBuffer())
767  ->apply(sidre::DOUBLE_ID, numDofs, offset);
768  }
769 
770  gf->NewDataAndSize(vv->getData(), numDofs);
771  }
772  else
773  {
774  // If we are not managing the grid function's data,
775  // create a view with the external data
776  vv->setExternalDataPtr(sidre::DOUBLE_ID, numDofs, gf->GetData());
777  }
778  MFEM_ASSERT((numDofs > 0 && vv->isApplied()) ||
779  (numDofs == 0 && vv->isEmpty() && vv->isDescribed()),
780  "invlid DataView state");
781  MFEM_ASSERT(numDofs == 0 || vv->getData() == gf->GetData(),
782  "DataView data is different from GridFunction data");
783  MFEM_ASSERT(vv->getNumElements() == numDofs,
784  "DataView size is different from GridFunction size");
785 }
786 
787 // private method
788 void SidreDataCollection::
789 addVectorBasedGridFunction(const std::string& field_name, GridFunction *gf,
790  const std::string &buffer_name,
791  asctoolkit::sidre::SidreLength offset)
792 {
793  sidre::DataGroup* grp = bp_grp->getGroup("fields/" + field_name);
794  MFEM_ASSERT(grp != NULL, "field " << field_name << " does not exist");
795 
796  const int FLD_SZ = 20;
797  char fidxName[FLD_SZ];
798 
799  int vdim = gf->FESpace()->GetVDim();
800  int ndof = gf->FESpace()->GetNDofs();
801  Ordering::Type ordering = gf->FESpace()->GetOrdering();
802 
803  if (gf->GetData() == NULL)
804  {
805  AllocNamedBuffer(buffer_name, offset + vdim*ndof);
806  // gf->data is set below.
807  }
808 
809  /*
810  * Mesh blueprint for a vector-based grid function is of the form
811  * /fields/field_name/basis
812  * -- string value is GridFunction's FEC::Name
813  * /fields/field_name/values/x0
814  * /fields/field_name/values/x1
815  * ...
816  * /fields/field_name/values/xn
817  * -- each coordinate is an array of size ndof
818  */
819 
820  // Get/create the DataGroup "values".
821  sidre::DataGroup *vg = alloc_group(grp, "values");
822 
823  // Create the DataViews "x0", "x1", etc inside the "values" DataGroup, vg.
824  // If we have a named buffer for field_name, attach it to the DataViews;
825  // otherwise set the DataViews to use gf->GetData() as external data.
826  sidre::DataType dtype = sidre::DataType::c_double(ndof);
827  const int entry_stride = (ordering == Ordering::byNODES ? 1 : vdim);
828  const int vdim_stride = (ordering == Ordering::byNODES ? ndof : 1);
829  dtype.set_stride(dtype.stride()*entry_stride);
830 
831  sidre::DataView *bv = named_buffers_grp()->getView(buffer_name);
832  if (bv)
833  {
834  MFEM_ASSERT(bv->hasBuffer() && bv->isDescribed(), "");
835 
836  // named buffers always have offset 0
837  MFEM_ASSERT(bv->getSchema().dtype().offset() == 0, "");
838  dtype.set_offset(dtype.element_bytes()*offset);
839 
840  for (int d = 0; d < vdim; d++)
841  {
842  std::snprintf(fidxName, FLD_SZ, "x%d", d);
843  sidre::DataView *xv = alloc_view(vg, fidxName, dtype);
844  xv->attachBuffer(bv->getBuffer());
845  dtype.set_offset(dtype.offset() + dtype.element_bytes()*vdim_stride);
846  }
847 
848  gf->NewDataAndSize(bv->getData<double*>() + offset, vdim*ndof);
849  }
850  else
851  {
852  for (int d = 0; d < vdim; d++)
853  {
854  std::snprintf(fidxName, FLD_SZ, "x%d", d);
855  sidre::DataView *xv = alloc_view(vg, fidxName, dtype);
856  xv->setExternalDataPtr(gf->GetData());
857  dtype.set_offset(dtype.offset() + dtype.element_bytes()*vdim_stride);
858  }
859  }
860 
861 #ifdef MFEM_DEBUG
862  for (int d = 0; d < vdim; d++)
863  {
864  std::snprintf(fidxName, FLD_SZ, "x%d", d);
865  sidre::DataView *xv = vg->getView(fidxName);
866  MFEM_ASSERT((ndof > 0 && xv->isApplied()) ||
867  (ndof == 0 && xv->isEmpty() && xv->isDescribed()),
868  "invlid DataView state");
869  MFEM_ASSERT(ndof == 0 || xv->getData() == gf->GetData() + d*vdim_stride,
870  "DataView data is different from GridFunction data");
871  MFEM_ASSERT(xv->getNumElements() == ndof,
872  "DataView size is different from GridFunction size");
873  }
874 #endif
875 }
876 
877 // private method
878 // Should only be called on mpi rank 0 ( or if serial problem ).
879 void SidreDataCollection::
880 RegisterFieldInBPIndex(const std::string& field_name, GridFunction *gf)
881 {
882  sidre::DataGroup *bp_field_grp = bp_grp->getGroup("fields/" + field_name);
883  sidre::DataGroup *bp_index_field_grp =
884  bp_index_grp->createGroup("fields/" + field_name);
885 
886  bp_index_field_grp->createViewString( "path", bp_field_grp->getPathName() );
887  bp_index_field_grp->copyView( bp_field_grp->getView("topology") );
888  bp_index_field_grp->copyView( bp_field_grp->getView("basis") );
889 
890  // Note: The bp index requires GridFunction::VectorDim()
891  // since the GF might be scalar valued and have a vector basis
892  // (e.g. hdiv and hcurl spaces)
893  const int number_of_components = gf->VectorDim();
894  bp_index_field_grp->createViewScalar("number_of_components",
895  number_of_components);
896 }
897 
898 // private method
899 // Should only be called on mpi rank 0 ( or if serial problem ).
900 void SidreDataCollection::
901 DeregisterFieldInBPIndex(const std::string& field_name)
902 {
903  sidre::DataGroup * fields_grp = bp_index_grp->getGroup("fields");
904  MFEM_VERIFY(fields_grp->hasGroup(field_name),
905  "No field exists in blueprint index with name " << name);
906 
907  // Note: This will destroy all orphaned views or buffer classes under this
908  // group also. If sidre owns this field data, the memory will be deleted
909  // unless it's referenced somewhere else in sidre.
910  fields_grp->destroyGroup(field_name);
911 }
912 
913 void SidreDataCollection::RegisterField(const std::string &field_name,
914  GridFunction *gf,
915  const std::string &buffer_name,
916  asctoolkit::sidre::SidreLength offset)
917 {
918  if ( field_name.empty() || buffer_name.empty() ||
919  gf == NULL || gf->FESpace() == NULL )
920  {
921  return;
922  }
923 
924  // Register field_name in the blueprint group.
925  sidre::DataGroup* f = bp_grp->getGroup("fields");
926 
927  if (f->hasGroup( field_name ))
928  {
929  // There are two possibilities:
930  // 1. If HasField(field_name) is true - we are overwriting a field that
931  // was previously registered.
932  // 2. Otherwise, the field was loaded from a file, or defined outside of
933  // the data collection.
934  if (HasField(field_name))
935  {
936  MFEM_DEBUG_DO(
937  MFEM_WARNING("field with the name '" << field_name<< "' is already "
938  "registered, overwriting the old field"));
939  DeregisterField(field_name);
940  }
941  }
942 
943  sidre::DataGroup* grp = f->createGroup( field_name );
944 
945  // Set the "basis" string using the gf's finite element space, overwrite if
946  // necessary.
947  sidre::DataView *v = alloc_view(grp, "basis");
948  v->setString(gf->FESpace()->FEColl()->Name());
949 
950  // Set the topology of the GridFunction.
951  // This is always 'mesh' except for a special case with the boundary material
952  // attributes field.
953  v = alloc_view(grp, "topology")->setString("mesh");
954 
955  MFEM_ASSERT(gf->Size() == gf->FESpace()->GetVSize(),
956  "GridFunction size does not match FiniteElementSpace size");
957 
958  // Set the data views of the grid function
959  // e.g. the number of coefficients per DoF -- either scalar-valued or
960  // vector-valued
961  bool const isScalarValued = (gf->FESpace()->GetVDim() == 1);
962  if (isScalarValued)
963  {
964  // Set the DataView "<bp_grp>/fields/<field_name>/values"
965  addScalarBasedGridFunction(field_name, gf, buffer_name, offset);
966  }
967  else // vector valued
968  {
969  // Set the DataGroup "<bp_grp>/fields/<field_name>/values"
970  addVectorBasedGridFunction(field_name, gf, buffer_name, offset);
971  }
972 
973  // Register field_name in the blueprint_index group.
974  if (myid == 0)
975  {
976  RegisterFieldInBPIndex(field_name, gf);
977  }
978 
979  // Register field_name + gf in field_map.
980  DataCollection::RegisterField(field_name, gf);
981 }
982 
983 void SidreDataCollection::DeregisterField(const std::string& field_name)
984 {
985  // Deregister field_name from field_map.
987 
988  sidre::DataGroup * fields_grp = bp_grp->getGroup("fields");
989  MFEM_VERIFY(fields_grp->hasGroup(field_name),
990  "No field exists in blueprint with name " << field_name);
991 
992  // Delete field_name from the blueprint group.
993 
994  // Note: This will destroy all orphaned views or buffer classes under this
995  // group also. If sidre owns this field data, the memory will be deleted
996  // unless it's referenced somewhere else in sidre.
997  fields_grp->destroyGroup(field_name);
998 
999  // Delete field_name from the blueprint_index group.
1000  if (myid == 0)
1001  {
1002  DeregisterFieldInBPIndex(field_name);
1003  }
1004 
1005  // Delete field_name from the named_buffers group, if allocated.
1006  FreeNamedBuffer(field_name);
1007 }
1008 
1009 // private method
1010 std::string SidreDataCollection::getElementName(Element::Type elementEnum)
1011 {
1012  // Note -- the mapping from Element::Type to string is based on
1013  // enum Element::Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
1014  // TETRAHEDRON, HEXAHEDRON };
1015  // Note: -- the string names are from conduit's blueprint
1016 
1017  switch (elementEnum)
1018  {
1019  case Element::POINT: return "point";
1020  case Element::SEGMENT: return "line";
1021  case Element::TRIANGLE: return "tri";
1022  case Element::QUADRILATERAL: return "quad";
1023  case Element::TETRAHEDRON: return "tet";
1024  case Element::HEXAHEDRON: return "hex";
1025  }
1026 
1027  return "unknown";
1028 }
1029 
1030 } // end namespace mfem
1031 
1032 #endif
int GetVSize() const
Definition: fespace.hpp:163
asctoolkit::sidre::DataView * AllocNamedBuffer(const std::string &buffer_name, asctoolkit::sidre::SidreLength sz, asctoolkit::sidre::TypeID type=asctoolkit::sidre::DOUBLE_ID)
Return newly allocated or existing named buffer for buffer_name.
asctoolkit::sidre::DataGroup * alloc_group(asctoolkit::sidre::DataGroup *grp, const std::string &group_name)
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:617
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:94
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
void Load(const std::string &path, const std::string &protocol)
Load the Sidre DataStore from file.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
double GetTime() const
Get physical time (for time-dependent simulations)
bool appendRankToFileName
Append rank to any output file names.
virtual void SetMesh(Mesh *new_mesh)
Set/change the mesh associated with the collection.
void DeleteAll()
Delete data owned by the DataCollection including field information.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Register a GridFunction in the Sidre DataStore.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:587
const Geometry::Type geom
asctoolkit::sidre::DataView * GetNamedBuffer(const std::string &buffer_name) const
Get a pointer to the sidre::DataView holding the named buffer for buffer_name.
asctoolkit::sidre::DataView * alloc_view(asctoolkit::sidre::DataGroup *grp, const std::string &view_name)
int Size() const
Returns the size of the vector.
Definition: vector.hpp:106
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:584
bool HasField(const std::string &name) const
Check if a grid function is part of the collection.
virtual void DeregisterField(const std::string &field_name)
Remove a grid function from the collection.
SidreDataCollection(const std::string &collection_name, Mesh *the_mesh=NULL, bool owns_mesh_data=false)
Constructor that allocates and initializes a Sidre DataStore.
void LoadExternalData(const std::string &path)
Load external data after registering externally owned fields.
Data type for vertex.
Definition: vertex.hpp:21
virtual void DeregisterField(const std::string &field_name)
De-register field_name from the SidreDataCollection.
bool own_data
Should the collection delete its mesh and fields.
asctoolkit::sidre::DataGroup * named_buffers_grp() const
double * GetData() const
Definition: vector.hpp:114
double GetTimeStep() const
Get the simulation time step (for time-dependent simulations)
double time
Physical time (for time-dependent simulations)
Mesh * mesh
The (common) mesh for the collected fields.
int dim
Definition: ex3.cpp:47
bool serial
Serial or parallel run?
int GetCycle() const
Get time cycle (for time-dependent simulations)
Type
Constants for the classes derived from Element.
Definition: element.hpp:37
static int create_directory(const std::string &dir_name, const Mesh *mesh, int myid)
void SetNodesOwner(bool nodes_owner)
Set the mesh nodes ownership flag.
Definition: mesh.hpp:881
const Element * GetElement(int i) const
Definition: mesh.hpp:642
void SetGroupPointers(asctoolkit::sidre::DataGroup *global_grp, asctoolkit::sidre::DataGroup *domain_grp)
Reset the domain and global datastore group pointers.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:153
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:293
void UpdateStateToDS()
Updates the data store&#39;s cycle, time, and time-step variables with the values from the SidreDataColle...
void SetTime(double t)
Set physical time (for time-dependent simulations)
int SpaceDimension() const
Definition: mesh.hpp:612
void UpdateStateFromDS()
Updates the DataCollection&#39;s cycle, time, and time-step variables with the values from the data store...
MPI_Comm GetComm() const
Definition: pmesh.hpp:116
virtual ~SidreDataCollection()
Delete all owned data.
int myid
MPI rank (in parallel)
std::string get_file_path(const std::string &filename) const
virtual void PrepareToSave()
Prepare the DataStore for writing.
int GetElementBaseGeometry(int i=0) const
Definition: mesh.hpp:654
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:581
void SetTimeStep(double ts)
Set the simulation time step (for time-dependent simulations)
virtual const char * Name() const
Definition: fe_coll.hpp:117
bool OwnsNodes() const
Return the mesh nodes ownership flag.
Definition: mesh.hpp:879
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:3273
void GetBdrElementData(int geom, Array< int > &bdr_elem_vtx, Array< int > &bdr_attr) const
Definition: mesh.hpp:627
virtual void Save()
Save the collection to file.
virtual int GetNVertices() const =0
void FreeNamedBuffer(const std::string &buffer_name)
Deallocate the named buffer buffer_name.
int pad_digits
Number of digits used for the cycle and MPI rank in filenames.
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:177
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5368
void SetComm(MPI_Comm comm)
Associate an MPI communicator with the collection.
virtual void SetMesh(Mesh *new_mesh)
Set/change the mesh associated with the collection.
void ChangeVertexDataOwnership(double *vertices, int len_vertices, bool zerocopy=false)
Set the internal Vertex array to point to the given vertices array without assuming ownership of the ...
Definition: mesh.cpp:2392
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
std::string name
Name of the collection, used as a directory name when saving.
int num_procs
Number of MPI ranks (in parallel)
Class for parallel meshes.
Definition: pmesh.hpp:28
int GetBdrElementBaseGeometry(int i=0) const
Definition: mesh.hpp:657
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:646
void GetElementData(const Array< Element * > &elem_array, int geom, Array< int > &elem_vtx, Array< int > &attr) const
Definition: mesh.cpp:6218
virtual int GetType() const =0
Returns element&#39;s type.