DOLFIN-X
DOLFIN-X C++ interface
xdmf_read.h
1 // Copyright (C) 2012-2018 Chris N. Richardson and Garth N. Wells
2 //
3 // This file is part of DOLFINX (https://www.fenicsproject.org)
4 //
5 // SPDX-License-Identifier: LGPL-3.0-or-later
6 
7 #pragma once
8 
9 #include "HDF5File.h"
10 #include "pugixml.hpp"
11 #include "xdmf_utils.h"
12 #include <Eigen/Dense>
13 #include <boost/algorithm/string.hpp>
14 #include <boost/lexical_cast.hpp>
15 #include <dolfinx/common/IndexMap.h>
16 #include <dolfinx/mesh/Geometry.h>
17 #include <dolfinx/mesh/Mesh.h>
18 #include <dolfinx/mesh/Topology.h>
19 
22 {
23 
25 template <typename T>
26 std::vector<T> get_dataset(MPI_Comm comm, const pugi::xml_node& dataset_node,
27  const hid_t h5_id,
28  std::array<std::int64_t, 2> range = {{0, 0}})
29 {
30  // FIXME: Need to sort out datasset dimensions - can't depend on HDF5
31  // shape, and a Topology data item is not required to have a
32  // 'Dimensions' attribute since the dimensions can be determined from
33  // the number of cells and the cell type (for topology, one must
34  // supply cell type + (number of cells or dimensions).
35  //
36  // A geometry data item must have 'Dimensions' attribute.
37 
38  assert(dataset_node);
39  pugi::xml_attribute format_attr = dataset_node.attribute("Format");
40  assert(format_attr);
41 
42  // Get data set shape from 'Dimensions' attribute (empty if not
43  // available)
44  const std::vector<std::int64_t> shape_xml
45  = xdmf_utils::get_dataset_shape(dataset_node);
46 
47  const std::string format = format_attr.as_string();
48  std::vector<T> data_vector;
49  // Only read ASCII on process 0
50  const int rank = dolfinx::MPI::rank(comm);
51  if (format == "XML")
52  {
53  if (rank == 0)
54  {
55  // Read data and trim any leading/trailing whitespace
56  pugi::xml_node data_node = dataset_node.first_child();
57  assert(data_node);
58  std::string data_str = data_node.value();
59 
60  // Split data based on spaces and line breaks
61  std::vector<boost::iterator_range<std::string::iterator>> data_vector_str;
62  boost::split(data_vector_str, data_str, boost::is_any_of(" \n"));
63 
64  // Add data to numerical vector
65  data_vector.reserve(data_vector_str.size());
66  for (auto& v : data_vector_str)
67  {
68  if (v.begin() != v.end())
69  data_vector.push_back(
70  boost::lexical_cast<T>(boost::copy_range<std::string>(v)));
71  }
72  }
73  }
74  else if (format == "HDF")
75  {
76  // Get file and data path
77  auto paths = xdmf_utils::get_hdf5_paths(dataset_node);
78 
79  // Get data shape from HDF5 file
80  const std::vector<std::int64_t> shape_hdf5
81  = HDF5Interface::get_dataset_shape(h5_id, paths[1]);
82 
83  // FIXME: should we support empty data sets?
84  // Check that data set is not empty
85  assert(!shape_hdf5.empty());
86  assert(shape_hdf5[0] != 0);
87 
88  // Determine range of data to read from HDF5 file. This is
89  // complicated by the XML Dimension attribute and the HDF5 storage
90  // possibly having different shapes, e.g. the HDF5 storage may be a
91  // flat array.
92 
93  // If range = {0, 0} then no range is supplied and we must determine
94  // the range
95  if (range[0] == 0 and range[1] == 0)
96  {
97  if (shape_xml == shape_hdf5)
98  {
99  range = dolfinx::MPI::local_range(rank, shape_hdf5[0],
100  dolfinx::MPI::size(comm));
101  }
102  else if (!shape_xml.empty() and shape_hdf5.size() == 1)
103  {
104  // Size of dims > 0
105  std::int64_t d = std::accumulate(shape_xml.begin(), shape_xml.end(), 1,
106  std::multiplies<std::int64_t>());
107 
108  // Check for data size consistency
109  if (d * shape_xml[0] != shape_hdf5[0])
110  {
111  throw std::runtime_error("Data size in XDMF/XML and size of HDF5 "
112  "dataset are inconsistent");
113  }
114 
115  // Compute data range to read
116  range = dolfinx::MPI::local_range(rank, shape_xml[0],
117  dolfinx::MPI::rank(comm));
118  range[0] *= d;
119  range[1] *= d;
120  }
121  else
122  {
123  throw std::runtime_error(
124  "This combination of array shapes in XDMF and HDF5 "
125  "is not supported");
126  }
127  }
128 
129  // Retrieve data
130  data_vector = HDF5Interface::read_dataset<T>(h5_id, paths[1], range);
131  }
132  else
133  throw std::runtime_error("Storage format \"" + format + "\" is unknown");
134 
135  // Get dimensions for consistency (if available in DataItem node)
136  if (shape_xml.empty())
137  {
138  std::int64_t size = 1;
139  for (auto dim : shape_xml)
140  size *= dim;
141 
142  std::int64_t size_global = 0;
143  const std::int64_t size_local = data_vector.size();
144  MPI_Allreduce(&size_local, &size_global, 1, MPI_INT64_T, MPI_SUM, comm);
145  if (size != size_global)
146  {
147  throw std::runtime_error(
148  "Data sizes in attribute and size of data read are inconsistent");
149  }
150  }
151 
152  return data_vector;
153 }
154 //----------------------------------------------------------------------------
155 
156 } // namespace dolfinx::io::xdmf_read
dolfinx::io::xdmf_read
Low-level methods for reading XDMF files.
Definition: xdmf_read.h:21
dolfinx::MPI::rank
static int rank(MPI_Comm comm)
Return process rank for the communicator.
Definition: MPI.cpp:75
dolfinx::io::HDF5Interface::get_dataset_shape
static std::vector< std::int64_t > get_dataset_shape(const hid_t hdf5_file_handle, const std::string dataset_path)
Get dataset shape (size of each dimension)
Definition: HDF5Interface.cpp:364
dolfinx::MPI::size
static int size(MPI_Comm comm)
Return size of the group (number of processes) associated with the communicator.
Definition: MPI.cpp:83
dolfinx::io::xdmf_read::get_dataset
std::vector< T > get_dataset(MPI_Comm comm, const pugi::xml_node &dataset_node, const hid_t h5_id, std::array< std::int64_t, 2 > range={{0, 0}})
Return data associated with a data set node.
Definition: xdmf_read.h:26
dolfinx::MPI::local_range
static std::array< std::int64_t, 2 > local_range(int process, std::int64_t N, int size)
Return local range for given process, splitting [0, N - 1] into size() portions of almost equal size.
Definition: MPI.cpp:101