OGS
NodePartitionedMeshReader.cpp
Go to the documentation of this file.
1 
17 
18 #include "BaseLib/Logging.h"
19 
20 #ifdef USE_PETSC
21 #include <mpi.h>
22 #endif
23 
24 #include "BaseLib/FileTools.h"
25 #include "BaseLib/RunTime.h"
27 #include "MeshLib/MeshEnums.h"
28 #include "MeshLib/Properties.h"
29 
30 // Check if the value can by converted to given type without overflow.
31 template <typename VALUE, typename TYPE>
32 bool is_safely_convertable(VALUE const& value)
33 {
34  bool const result = value <= std::numeric_limits<TYPE>::max();
35  if (!result)
36  {
37  ERR("The value {:d} is too large for conversion.", value);
38  ERR("Maximum available size is {:d}.",
39  std::numeric_limits<TYPE>::max());
40  }
41  return result;
42 }
43 
44 namespace MeshLib
45 {
46 namespace IO
47 {
49  : _mpi_comm(comm)
50 {
51  MPI_Comm_size(_mpi_comm, &_mpi_comm_size);
52  MPI_Comm_rank(_mpi_comm, &_mpi_rank);
53 
55 }
56 
58 {
59  MPI_Type_free(&_mpi_node_type);
60 }
61 
63 {
64  int const count = 2;
65  int blocks[count] = {1, 3};
66  MPI_Datatype types[count] = {MPI_UNSIGNED_LONG, MPI_DOUBLE};
67  MPI_Aint displacements[count] = {0, sizeof(NodeData::index)};
68 
69  MPI_Type_create_struct(count, blocks, displacements, types,
71  MPI_Type_commit(&_mpi_node_type);
72 }
73 
75  const std::string& file_name_base)
76 {
77  BaseLib::RunTime timer;
78  timer.start();
79 
80  MeshLib::NodePartitionedMesh* mesh = nullptr;
81 
82  // Always try binary file first
83  std::string const fname_new = file_name_base + "_partitioned_msh_cfg" +
84  std::to_string(_mpi_comm_size) + ".bin";
85 
86  if (!BaseLib::IsFileExisting(fname_new)) // binary file does not exist.
87  {
88  OGS_FATAL(
89  "Required binary file {:s} does not exist.\n"
90  "Reading of ASCII mesh file is not supported since OGS version "
91  "6.3.3.",
92  fname_new);
93  }
94 
95  INFO("Reading corresponding part of mesh data from binary file {:s} ...",
96  file_name_base);
97  mesh = readMesh(file_name_base);
98 
99  INFO("[time] Reading the mesh took {:f} s.", timer.elapsed());
100 
101  MPI_Barrier(_mpi_comm);
102 
103  return mesh;
104 }
105 
106 template <typename DATA>
107 bool NodePartitionedMeshReader::readDataFromFile(std::string const& filename,
108  MPI_Offset offset,
109  MPI_Datatype type,
110  DATA& data) const
111 {
112  // Check container size
113  if (!is_safely_convertable<std::size_t, int>(data.size()))
114  {
115  ERR("The container size is too large for MPI_File_read() call.");
116  return false;
117  }
118 
119  // Open file
120  MPI_File file;
121 
122  char* filename_char = const_cast<char*>(filename.data());
123  int const file_status = MPI_File_open(
124  _mpi_comm, filename_char, MPI_MODE_RDONLY, MPI_INFO_NULL, &file);
125 
126  if (file_status != 0)
127  {
128  ERR("Error opening file {:s}. MPI error code {:d}", filename,
129  file_status);
130  return false;
131  }
132 
133  // Read data
134  char file_mode[] = "native";
135  MPI_File_set_view(file, offset, type, type, file_mode, MPI_INFO_NULL);
136  // The static cast is checked above.
137  MPI_File_read(file, data.data(), static_cast<int>(data.size()), type,
138  MPI_STATUS_IGNORE);
139  MPI_File_close(&file);
140 
141  return true;
142 }
143 
145  const std::string& file_name_base)
146 {
147  //----------------------------------------------------------------------------------
148  // Read headers
149  const std::string fname_header = file_name_base + "_partitioned_msh_";
150  const std::string fname_num_p_ext = std::to_string(_mpi_comm_size) + ".bin";
151 
152  if (!readDataFromFile(
153  fname_header + "cfg" + fname_num_p_ext,
154  static_cast<MPI_Offset>(static_cast<unsigned>(_mpi_rank) *
155  sizeof(_mesh_info)),
156  MPI_LONG, _mesh_info))
157  return nullptr;
158 
159  //----------------------------------------------------------------------------------
160  // Read Nodes
161  std::vector<NodeData> nodes(_mesh_info.nodes);
162 
163  if (!readDataFromFile(fname_header + "nod" + fname_num_p_ext,
164  static_cast<MPI_Offset>(_mesh_info.offset[2]),
165  _mpi_node_type, nodes))
166  return nullptr;
167 
168  std::vector<MeshLib::Node*> mesh_nodes;
169  std::vector<unsigned long> glb_node_ids;
170  setNodes(nodes, mesh_nodes, glb_node_ids);
171 
172  //----------------------------------------------------------------------------------
173  // Read non-ghost elements
174  std::vector<unsigned long> elem_data(_mesh_info.regular_elements +
175  _mesh_info.offset[0]);
176  if (!readDataFromFile(fname_header + "ele" + fname_num_p_ext,
177  static_cast<MPI_Offset>(_mesh_info.offset[3]),
178  MPI_LONG, elem_data))
179  return nullptr;
180 
181  std::vector<MeshLib::Element*> mesh_elems(_mesh_info.regular_elements +
183  setElements(mesh_nodes, elem_data, mesh_elems);
184 
185  //----------------------------------------------------------------------------------
186  // Read ghost element
187  std::vector<unsigned long> ghost_elem_data(_mesh_info.ghost_elements +
188  _mesh_info.offset[1]);
189 
190  if (!readDataFromFile(fname_header + "ele_g" + fname_num_p_ext,
191  static_cast<MPI_Offset>(_mesh_info.offset[4]),
192  MPI_LONG, ghost_elem_data))
193  return nullptr;
194 
195  const bool process_ghost = true;
196  setElements(mesh_nodes, ghost_elem_data, mesh_elems, process_ghost);
197 
198  //----------------------------------------------------------------------------------
199  // read the properties
200  MeshLib::Properties p(readProperties(file_name_base));
201 
202  return newMesh(BaseLib::extractBaseName(file_name_base), mesh_nodes,
203  glb_node_ids, mesh_elems, p);
204 }
205 
207  const std::string& file_name_base) const
208 {
210  readProperties(file_name_base, MeshLib::MeshItemType::Node, p);
211  readProperties(file_name_base, MeshLib::MeshItemType::Cell, p);
212  return p;
213 }
214 
216  const std::string& file_name_base, MeshLib::MeshItemType t,
217  MeshLib::Properties& p) const
218 {
219  std::string const item_type =
220  t == MeshLib::MeshItemType::Node ? "node" : "cell";
221  const std::string fname_cfg = file_name_base + "_partitioned_" + item_type +
222  "_properties_cfg" +
223  std::to_string(_mpi_comm_size) + ".bin";
224  std::ifstream is(fname_cfg.c_str(), std::ios::binary | std::ios::in);
225  if (!is)
226  {
227  WARN(
228  "Could not open file '{:s}'.\n"
229  "\tYou can ignore this warning if the mesh does not contain {:s}-"
230  "wise property data.",
231  fname_cfg, item_type.data());
232  return;
233  }
234  std::size_t number_of_properties = 0;
235  is.read(reinterpret_cast<char*>(&number_of_properties),
236  sizeof(std::size_t));
237  std::vector<std::optional<MeshLib::IO::PropertyVectorMetaData>> vec_pvmd(
239  for (std::size_t i(0); i < number_of_properties; ++i)
240  {
241  vec_pvmd[i] = MeshLib::IO::readPropertyVectorMetaData(is);
242  if (!vec_pvmd[i])
243  {
244  OGS_FATAL(
245  "Error in NodePartitionedMeshReader::readProperties: "
246  "Could not read the meta data for the PropertyVector {:d}",
247  i);
248  }
249  }
250  for (std::size_t i(0); i < number_of_properties; ++i)
251  {
253  }
254  auto pos = is.tellg();
255  auto offset =
256  static_cast<long>(pos) +
257  static_cast<long>(_mpi_rank *
259  is.seekg(offset);
260  std::optional<MeshLib::IO::PropertyVectorPartitionMetaData> pvpmd(
262  bool pvpmd_read_ok = static_cast<bool>(pvpmd);
263  bool all_pvpmd_read_ok;
264  MPI_Allreduce(&pvpmd_read_ok, &all_pvpmd_read_ok, 1, MPI_C_BOOL, MPI_LOR,
265  _mpi_comm);
266  if (!all_pvpmd_read_ok)
267  {
268  OGS_FATAL(
269  "Error in NodePartitionedMeshReader::readProperties: "
270  "Could not read the partition meta data for the mpi process {:d}",
271  _mpi_rank);
272  }
273  DBUG("[{:d}] offset in the PropertyVector: {:d}", _mpi_rank, pvpmd->offset);
274  DBUG("[{:d}] {:d} tuples in partition.", _mpi_rank,
275  pvpmd->number_of_tuples);
276  is.close();
277 
278  const std::string fname_val = file_name_base + "_partitioned_" + item_type +
279  "_properties_val" +
280  std::to_string(_mpi_comm_size) + ".bin";
281  is.open(fname_val.c_str(), std::ios::binary | std::ios::in);
282  if (!is)
283  {
284  ERR("Could not open file '{:s}'\n."
285  "\tYou can ignore this warning if the mesh does not contain {:s}-"
286  "wise property data.",
287  fname_val, item_type.data());
288  }
289 
290  readDomainSpecificPartOfPropertyVectors(vec_pvmd, *pvpmd, t, is, p);
291 }
292 
294  std::vector<std::optional<MeshLib::IO::PropertyVectorMetaData>> const&
295  vec_pvmd,
298  std::istream& is,
299  MeshLib::Properties& p) const
300 {
301  unsigned long global_offset = 0;
302  std::size_t const number_of_properties = vec_pvmd.size();
303  for (std::size_t i(0); i < number_of_properties; ++i)
304  {
305  DBUG(
306  "[{:d}] global offset: {:d}, offset within the PropertyVector: "
307  "{:d}.",
308  _mpi_rank, global_offset,
309  global_offset + pvpmd.offset * vec_pvmd[i]->number_of_components *
310  vec_pvmd[i]->data_type_size_in_bytes);
311  if (vec_pvmd[i]->is_int_type)
312  {
313  if (vec_pvmd[i]->is_data_type_signed)
314  {
315  if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(char))
316  createPropertyVectorPart<char>(is, *vec_pvmd[i], pvpmd, t,
317  global_offset, p);
318  else if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(int))
319  createPropertyVectorPart<int>(is, *vec_pvmd[i], pvpmd, t,
320  global_offset, p);
321  else if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(long))
322  createPropertyVectorPart<long>(is, *vec_pvmd[i], pvpmd, t,
323  global_offset, p);
324  else
325  {
326  WARN(
327  "Implementation for reading signed integer property "
328  "vector '{:s}' is not available.",
329  vec_pvmd[i]->property_name);
330  }
331  }
332  else
333  {
334  if (vec_pvmd[i]->data_type_size_in_bytes ==
335  sizeof(unsigned char))
336  createPropertyVectorPart<unsigned char>(
337  is, *vec_pvmd[i], pvpmd, t, global_offset, p);
338  else if (vec_pvmd[i]->data_type_size_in_bytes ==
339  sizeof(unsigned int))
340  createPropertyVectorPart<unsigned int>(
341  is, *vec_pvmd[i], pvpmd, t, global_offset, p);
342  else if (vec_pvmd[i]->data_type_size_in_bytes ==
343  sizeof(unsigned long))
344  createPropertyVectorPart<unsigned long>(
345  is, *vec_pvmd[i], pvpmd, t, global_offset, p);
346  else
347  {
348  WARN(
349  "Implementation for reading unsigned property vector "
350  "'{:s}' is not available.",
351  vec_pvmd[i]->property_name);
352  }
353  }
354  }
355  else
356  {
357  if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(float))
358  createPropertyVectorPart<float>(is, *vec_pvmd[i], pvpmd, t,
359  global_offset, p);
360  else if (vec_pvmd[i]->data_type_size_in_bytes == sizeof(double))
361  createPropertyVectorPart<double>(is, *vec_pvmd[i], pvpmd, t,
362  global_offset, p);
363  else
364  {
365  WARN(
366  "Implementation for reading floating point property vector "
367  "'{:s}' is not available.",
368  vec_pvmd[i]->property_name);
369  }
370  }
371  global_offset += vec_pvmd[i]->data_type_size_in_bytes *
372  vec_pvmd[i]->number_of_tuples *
373  vec_pvmd[i]->number_of_components;
374  }
375 }
376 
378  std::string const& mesh_name,
379  std::vector<MeshLib::Node*> const& mesh_nodes,
380  std::vector<unsigned long> const& glb_node_ids,
381  std::vector<MeshLib::Element*> const& mesh_elems,
382  MeshLib::Properties const& properties) const
383 {
384  return new MeshLib::NodePartitionedMesh(
385  mesh_name, mesh_nodes, glb_node_ids, mesh_elems, properties,
388 }
389 
391  const std::vector<NodeData>& node_data,
392  std::vector<MeshLib::Node*>& mesh_node,
393  std::vector<unsigned long>& glb_node_ids) const
394 {
395  mesh_node.resize(_mesh_info.nodes);
396  glb_node_ids.resize(_mesh_info.nodes);
397 
398  for (std::size_t i = 0; i < mesh_node.size(); i++)
399  {
400  NodeData const& nd = node_data[i];
401  glb_node_ids[i] = nd.index;
402  mesh_node[i] = new MeshLib::Node(nd.x, nd.y, nd.z, i);
403  }
404 }
405 
407  const std::vector<MeshLib::Node*>& mesh_nodes,
408  const std::vector<unsigned long>& elem_data,
409  std::vector<MeshLib::Element*>& mesh_elems, const bool ghost) const
410 {
411  // Number of elements, either ghost or regular
412  unsigned long const ne =
414  unsigned long const id_offset_ghost =
415  ghost ? _mesh_info.regular_elements : 0;
416 
417  for (unsigned long i = 0; i < ne; i++)
418  {
419  unsigned long id_offset_elem = elem_data[i];
420 
421  // Unused for now, keep for elem_data documentation purpose here.
422  {
423  const unsigned mat_idx =
424  static_cast<unsigned>(elem_data[id_offset_elem++]);
425  (void)mat_idx;
426  }
427  const unsigned long e_type = elem_data[id_offset_elem++];
428  unsigned long const nnodes = elem_data[id_offset_elem++];
429 
430  MeshLib::Node** elem_nodes = new MeshLib::Node*[nnodes];
431  for (unsigned long k = 0; k < nnodes; k++)
432  elem_nodes[k] = mesh_nodes[elem_data[id_offset_elem++]];
433 
434  // The element types below are defined by the MeshLib::CellType.
435  switch (static_cast<CellType>(e_type))
436  {
437  case CellType::POINT1:
438  mesh_elems[i + id_offset_ghost] =
439  new MeshLib::Point(elem_nodes);
440  break;
441  case CellType::LINE2:
442  mesh_elems[i + id_offset_ghost] = new MeshLib::Line(elem_nodes);
443  break;
444  case CellType::LINE3:
445  mesh_elems[i + id_offset_ghost] =
446  new MeshLib::Line3(elem_nodes);
447  break;
448  case CellType::QUAD4:
449  mesh_elems[i + id_offset_ghost] = new MeshLib::Quad(elem_nodes);
450  break;
451  case CellType::QUAD8:
452  mesh_elems[i + id_offset_ghost] =
453  new MeshLib::Quad8(elem_nodes);
454  break;
455  case CellType::QUAD9:
456  mesh_elems[i + id_offset_ghost] =
457  new MeshLib::Quad9(elem_nodes);
458  break;
459  case CellType::HEX8:
460  mesh_elems[i + id_offset_ghost] = new MeshLib::Hex(elem_nodes);
461  break;
462  case CellType::HEX20:
463  mesh_elems[i + id_offset_ghost] =
464  new MeshLib::Hex20(elem_nodes);
465  break;
466  case CellType::HEX27:
467  OGS_FATAL(
468  "NodePartitionedMeshReader: construction of HEX27 element "
469  "with id {:d} is not implemented.",
470  i);
471  break;
472  case CellType::TRI3:
473  mesh_elems[i + id_offset_ghost] = new MeshLib::Tri(elem_nodes);
474  break;
475  case CellType::TRI6:
476  mesh_elems[i + id_offset_ghost] = new MeshLib::Tri6(elem_nodes);
477  break;
478  case CellType::TET4:
479  mesh_elems[i + id_offset_ghost] = new MeshLib::Tet(elem_nodes);
480  break;
481  case CellType::TET10:
482  mesh_elems[i + id_offset_ghost] =
483  new MeshLib::Tet10(elem_nodes);
484  break;
485  case CellType::PRISM6:
486  mesh_elems[i + id_offset_ghost] =
487  new MeshLib::Prism(elem_nodes);
488  break;
489  case CellType::PRISM15:
490  mesh_elems[i + id_offset_ghost] =
491  new MeshLib::Prism15(elem_nodes);
492  break;
493  case CellType::PYRAMID5:
494  mesh_elems[i + id_offset_ghost] =
495  new MeshLib::Pyramid(elem_nodes);
496  break;
497  case CellType::PYRAMID13:
498  mesh_elems[i + id_offset_ghost] =
499  new MeshLib::Pyramid13(elem_nodes);
500  break;
501  case CellType::INVALID:
502  OGS_FATAL(
503  "NodePartitionedMeshReader: construction of INVALID "
504  "element type with id {:d} is not possible.",
505  i);
506  break;
507  default:
508  OGS_FATAL(
509  "NodePartitionedMeshReader: construction of element type "
510  "{:d} is not implemented.",
511  e_type);
512  }
513  }
514 }
515 } // namespace IO
516 } // namespace MeshLib
#define OGS_FATAL(...)
Definition: Error.h:26
Filename manipulation routines.
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
Definition of mesh-related Enumerations.
Definition of the class Properties that implements a container of properties.
bool is_safely_convertable(VALUE const &value)
Declare a class to read node-wise partitioned mesh with MPI functions.
Definition of the RunTime class.
Count the running time.
Definition: RunTime.h:29
double elapsed() const
Get the elapsed time in seconds.
Definition: RunTime.h:42
void start()
Start the timer.
Definition: RunTime.h:32
void setElements(const std::vector< MeshLib::Node * > &mesh_nodes, const std::vector< unsigned long > &elem_data, std::vector< MeshLib::Element * > &mesh_elems, const bool ghost=false) const
Set mesh elements from a temporary array containing node data read from file.
bool readDataFromFile(std::string const &filename, MPI_Offset offset, MPI_Datatype type, DATA &data) const
void readDomainSpecificPartOfPropertyVectors(std::vector< std::optional< MeshLib::IO::PropertyVectorMetaData >> const &vec_pvmd, MeshLib::IO::PropertyVectorPartitionMetaData const &pvpmd, MeshLib::MeshItemType t, std::istream &is, MeshLib::Properties &p) const
MeshLib::NodePartitionedMesh * readMesh(const std::string &file_name_base)
Create a NodePartitionedMesh object, read binary mesh data in the manner of parallel,...
MeshLib::NodePartitionedMesh * read(const std::string &file_name_base)
Create a NodePartitionedMesh object, read data to it, and return a pointer to it. Data files are in b...
void setNodes(const std::vector< NodeData > &node_data, std::vector< MeshLib::Node * > &mesh_node, std::vector< unsigned long > &glb_node_ids) const
Set mesh nodes from a temporary array containing node data read from file.
int _mpi_comm_size
Number of processes in the communicator: _mpi_comm.
MPI_Comm _mpi_comm
Pointer to MPI communicator.
MPI_Datatype _mpi_node_type
MPI data type for struct NodeData.
void registerNodeDataMpiType()
Define MPI data type for NodeData struct.
struct MeshLib::IO::NodePartitionedMeshReader::PartitionedMeshInfo _mesh_info
MeshLib::Properties readProperties(const std::string &file_name_base) const
MeshLib::NodePartitionedMesh * newMesh(std::string const &mesh_name, std::vector< MeshLib::Node * > const &mesh_nodes, std::vector< unsigned long > const &glb_node_ids, std::vector< MeshLib::Element * > const &mesh_elems, MeshLib::Properties const &properties) const
Create a new mesh of NodePartitionedMesh after reading and processing the data.
Property manager on mesh items. Class Properties manages scalar, vector or matrix properties....
Definition: Properties.h:36
bool IsFileExisting(const std::string &strFilename)
Returns true if given file exists.
Definition: FileTools.cpp:43
std::string extractBaseName(std::string const &pathname)
Definition: FileTools.cpp:175
std::optional< PropertyVectorPartitionMetaData > readPropertyVectorPartitionMetaData(std::istream &is)
void writePropertyVectorMetaData(std::ostream &os, PropertyVectorMetaData const &pvmd)
std::optional< PropertyVectorMetaData > readPropertyVectorMetaData(std::istream &is)
TemplateElement< MeshLib::PyramidRule13 > Pyramid13
Definition: Pyramid.h:26
CellType
Types of mesh elements supported by OpenGeoSys.
Definition: MeshEnums.h:43
TemplateElement< MeshLib::QuadRule8 > Quad8
Definition: Quad.h:29
TemplateElement< MeshLib::TetRule4 > Tet
Definition: Tet.h:25
TemplateElement< MeshLib::LineRule2 > Line
Definition: Line.h:25
MeshItemType
Definition: Location.h:21
TemplateElement< MeshLib::LineRule3 > Line3
Definition: Line.h:26
TemplateElement< MeshLib::TriRule3 > Tri
Definition: Tri.h:26
TemplateElement< MeshLib::QuadRule9 > Quad9
Definition: Quad.h:30
TemplateElement< MeshLib::PrismRule6 > Prism
Definition: Prism.h:25
TemplateElement< MeshLib::PyramidRule5 > Pyramid
Definition: Pyramid.h:25
TemplateElement< PointRule1 > Point
Definition: Point.h:20
TemplateElement< MeshLib::HexRule20 > Hex20
Definition: Hex.h:26
TemplateElement< MeshLib::TriRule6 > Tri6
Definition: Tri.h:27
TemplateElement< MeshLib::TetRule10 > Tet10
Definition: Tet.h:26
TemplateElement< MeshLib::PrismRule15 > Prism15
Definition: Prism.h:26
TemplateElement< MeshLib::HexRule8 > Hex
Definition: Hex.h:25
TemplateElement< MeshLib::QuadRule4 > Quad
Definition: Quad.h:28
unsigned long active_nodes
5: Number of all active nodes a partition,
unsigned long ghost_elements
3: Number of ghost element of a partition,
unsigned long nodes
0: Number of all nodes of a partition,
unsigned long global_nodes
7: Number of all nodes of global mesh,