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