MEPP2 Project
VtkFileReader.h
Go to the documentation of this file.
1 // Copyright (c) 2012-2019 University of Lyon and CNRS (France).
2 // All rights reserved.
3 //
4 // This file is part of MEPP2; you can redistribute it and/or modify
5 // it under the terms of the GNU Lesser General Public License as
6 // published by the Free Software Foundation; either version 3 of
7 // the License, or (at your option) any later version.
8 //
9 // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
10 // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
11 #pragma once
12 
13 #include <iostream>
14 #include <exception>
15 
16 #include "FEVV/Tools/IO/FileUtilities.hpp"
17 #include "FEVV/Tools/IO/StringUtilities.hpp"
18 
19 #if defined(__APPLE__) && defined(__clang__)
20 // In order to silence out warnings about VTK deprecated code invocation
21 #define VTK_LEGACY_SILENT
22 // We also want to silence out vtk complains about member function not being
23 // marked override:
24 #pragma clang diagnostic push
25 #pragma clang diagnostic ignored "-Winconsistent-missing-override"
26 #endif
27 #include <vtkVersion.h>
28 
29 #include <vtkSmartPointer.h>
30 #include <vtkPointData.h>
31 #include <vtkCellData.h>
32 #include <vtkCellArray.h>
33 #include <vtkPolyData.h>
34 #include <vtkUnstructuredGrid.h>
35 
36 #include <vtkPolyDataReader.h>
37 #include <vtkXMLPolyDataReader.h>
38 #include <vtkUnstructuredGridReader.h>
39 #include <vtkXMLUnstructuredGridReader.h>
40 #include <vtkDataSetSurfaceFilter.h>
41 #if defined(__APPLE__) && defined(__clang__)
42 #pragma clang diagnostic pop
43 #endif
44 
45 #ifdef _MSC_VER
46 #define IO_TOOLS_EXPORT //"DataStructures/IO_Tools/IO_Tools_export.h"
47 #else
48 #define IO_TOOLS_EXPORT
49 #endif
50 
51 namespace FEVV {
52 namespace IO {
53 
54 using namespace StrUtils;
55 using namespace FileUtils;
56 
57 inline std::string
58 get_line_containing_dataset(std::string file_path)
59 {
60  std::ifstream myfile(file_path);
61  if(myfile.is_open())
62  {
63  std::string line;
64  while(myfile.good())
65  {
66  getline(myfile, line);
67  if(line.find("DATASET") != std::string::npos)
68  {
69  myfile.close();
70  return line;
71  }
72  }
73  myfile.close();
74  }
75  else
76  {
77  throw std::runtime_error(
78  "Reader::getLineContainingDATASET -> unable to open file");
79  }
80  return std::string();
81 }
82 
83 template< typename CoordType,
84  typename CoordNType,
85  typename CoordCType,
86  typename IndexType >
87 void
88 read_vtk_poly_data(
89  const vtkSmartPointer< vtkPolyData > &poly_data,
90  std::vector< std::vector< CoordType > > &points_coords,
91  std::vector< std::vector< CoordNType > >
92  &normals_coords,
93  std::vector< std::vector< CoordCType > > &vertex_color_coords,
94  std::vector< std::vector< IndexType > >
95  &line_indices,
96  std::vector< std::vector< CoordCType > > &lines_color_coords,
97  std::vector< std::vector< IndexType > >
98  &face_indices,
99  std::vector< std::vector< CoordCType > > &face_color_coords,
100  std::vector< std::vector< std::vector< double > > >
101  &field_attributes, // array of (2D) arrays: first point data, then cell
102  // data
103  std::vector< std::string > &field_names)
104 {
106  // clearing vectors...
107  points_coords.clear();
108  normals_coords.clear();
109  vertex_color_coords.clear();
110  line_indices.clear();
111  lines_color_coords.clear();
112  face_indices.clear();
113  face_color_coords.clear();
114  field_attributes.clear();
115  field_names.clear();
118  // some checks:
119  // we do not manage triangle strips independently from cells yet
120  if(poly_data->GetNumberOfStrips() > 0)
121  {
122  std::cerr << "Reader::read_vtkPolyData -> Warning: triangle strips not "
123  "taken into account yet.\n";
124  }
125 
126  // each point there is at maximum one normal (some points may have no normal):
127  assert((poly_data->GetPointData()->GetNormals() == NULL) ||
128  (poly_data->GetNumberOfPoints() >=
129  poly_data->GetPointData()->GetNormals()->GetNumberOfTuples()));
131  // Get points and point normals
132  if(poly_data->GetNumberOfPoints() > 0)
133  {
134  points_coords.reserve(static_cast< size_t >(
135  poly_data->GetNumberOfPoints())); // use exactly the nb of points
136 
137  vtkSmartPointer< vtkPoints > ptr_points = poly_data->GetPoints();
138  vtkSmartPointer< vtkDataArray > ptr_normals =
139  poly_data->GetPointData()->GetNormals();
140  bool use_norm =
141  (ptr_normals != NULL &&
142  ptr_points->GetNumberOfPoints() ==
143  ptr_normals->GetNumberOfTuples()); // use normals only if one
144  // normal per point/vertex
145 
146  if(use_norm)
147  normals_coords.reserve(static_cast< size_t >(
148  poly_data->GetPointData()
149  ->GetNormals()
150  ->GetNumberOfTuples())); // use exactly the nb of normals
151 
152  for(vtkIdType id = 0; id < ptr_points->GetNumberOfPoints(); id++)
153  {
154  std::vector< double > vd(ptr_points->GetPoint(id),
155  ptr_points->GetPoint(id) +
156  3); // ** 3 coordinates of type double
157  // hard-coded in vtkPoints class ** ;
158  std::vector< CoordType > final_vc;
159  final_vc.reserve(vd.size());
160  std::vector< double >::iterator it(vd.begin()), ite(vd.end());
161  for(; it != ite; ++it)
162  final_vc.push_back(static_cast< CoordType >(*it));
163 
164  points_coords.push_back(final_vc);
165 
166  if(use_norm)
167  {
168  std::vector< double > vdn(ptr_normals->GetTuple(id),
169  ptr_normals->GetTuple(id) +
170  ptr_normals->GetNumberOfComponents());
171  std::vector< CoordNType > final_vcn;
172  final_vcn.reserve(vdn.size());
173  it = vdn.begin();
174  ite = vdn.end();
175  for(; it != ite; ++it)
176  final_vcn.push_back(static_cast< CoordNType >(*it));
177 
178  normals_coords.push_back(final_vcn);
179  }
180  }
181  }
183  // Get facet indices
184  if(poly_data->GetNumberOfPolys() > 0)
185  {
186  face_indices.reserve(static_cast< size_t >(
187  poly_data->GetNumberOfPolys())); // use exactly the nb of faces
188 
189  vtkSmartPointer< vtkCellArray > ptr_cell_polygons = poly_data->GetPolys();
190  vtkIdType npts;
191 #if (VTK_MAJOR_VERSION >= 9) // not tested with VTK 8...
192  const vtkIdType *pts_poly = new vtkIdType[ptr_cell_polygons->GetMaxCellSize()];
193 #else
194  vtkIdType *pts_poly = new vtkIdType[ptr_cell_polygons->GetMaxCellSize()];
195 #endif
196 
197  ptr_cell_polygons->InitTraversal();
198  while(ptr_cell_polygons->GetNextCell(
199  npts, pts_poly)) // 0 (i.e. false) is encountered at the list' end
200  {
201  std::vector< IndexType > facet_indices(npts);
202  for(vtkIdType i = 0; i < npts; i++)
203  facet_indices[i] = static_cast< IndexType >(pts_poly[i]);
204 
205  face_indices.push_back(facet_indices);
206  }
207 
208  delete[] pts_poly;
209  }
211  // Get line indices
212  if(poly_data->GetNumberOfLines() > 0)
213  {
214  line_indices.reserve(static_cast< size_t >(
215  poly_data->GetNumberOfLines())); // use exactly the nb of line_indices
216 
217  vtkSmartPointer< vtkCellArray > ptr_cell_lines = poly_data->GetLines();
218  vtkIdType npts;
219 #if (VTK_MAJOR_VERSION >= 9) // not tested with VTK 8...
220  const vtkIdType *pts_line = new vtkIdType[ptr_cell_lines->GetMaxCellSize()];
221 #else
222  vtkIdType *pts_line = new vtkIdType[ptr_cell_lines->GetMaxCellSize()];
223 #endif
224 
225  ptr_cell_lines->InitTraversal();
226  while(ptr_cell_lines->GetNextCell(
227  npts, pts_line)) // 0 (i.e. false) is encountered at the list' end
228  {
229  std::vector< IndexType > line_indices_tmp(npts);
230  for(vtkIdType i = 0; i < npts; i++)
231  line_indices_tmp[i] = static_cast< IndexType >(pts_line[i]);
232 
233  line_indices.push_back(line_indices_tmp);
234  }
235 
236  delete[] pts_line;
237  }
239  // Get 2D arrays of FIELD data when present
240  vtkFieldData *fd = poly_data->GetFieldData();
241  if(fd != NULL)
242  {
243  field_attributes.resize(static_cast< size_t >(
244  fd->GetNumberOfArrays())); // use exactly the nb of 2D arrays for FIELD
245  // data
246  field_names.resize(static_cast< size_t >(fd->GetNumberOfArrays()));
247 
248  for(int id_array = 0; id_array < fd->GetNumberOfArrays(); id_array++)
249  {
250  field_names[id_array] = fd->GetArrayName(id_array);
251  // for each FIELD data (arrays)
252  vtkSmartPointer< vtkDataArray > ptr_data = fd->GetArray(
253  id_array); // poly_data->GetFieldData()->GetArrayName(idArray) ) ;
254 
255  field_attributes[id_array].resize(
256  static_cast< size_t >(ptr_data->GetNumberOfTuples()));
257  // std::cout << "ptrData->GetNumberOfTuples() = " <<
258  // ptrData->GetNumberOfTuples() << std::endl ;
259  for(vtkIdType id_tuple_in_array = 0;
260  id_tuple_in_array < ptr_data->GetNumberOfTuples();
261  id_tuple_in_array++)
262  {
263  // field_attributes[idArray] is a 2D array
264  field_attributes[id_array][id_tuple_in_array].resize(
265  static_cast< size_t >(ptr_data->GetNumberOfComponents()));
266  std::copy(ptr_data->GetTuple(id_tuple_in_array),
267  ptr_data->GetTuple(id_tuple_in_array) +
268  ptr_data->GetNumberOfComponents(),
269 
270  field_attributes[id_array][id_tuple_in_array].begin());
271  }
272  }
273  }
275  // Get point data when present
276  vtkPointData *pd = poly_data->GetPointData();
277  if(pd != NULL)
278  {
279  size_t old_nb_e = field_attributes.size();
280  field_attributes.resize(static_cast< size_t >(
281  old_nb_e + pd->GetNumberOfArrays())); // use exactly the nb of 2D arrays
282  // for FIELD data
283  field_names.resize(
284  static_cast< size_t >(old_nb_e + pd->GetNumberOfArrays()));
285  for(int id_array = 0; id_array < pd->GetNumberOfArrays(); id_array++)
286  {
287  field_names[id_array + old_nb_e] =
288  std::string("POINT_DATA_") + pd->GetArrayName(id_array);
289  // std::cout << "POINT_DATA: array name is " << pd->GetArrayName(idArray)
290  // << std::endl;
291  // for each FIELD data (arrays)
292  vtkSmartPointer< vtkDataArray > ptr_data = pd->GetArray(id_array);
293 
294  field_attributes[id_array + old_nb_e].resize(
295  static_cast< size_t >(ptr_data->GetNumberOfTuples()));
296  // std::cout << "ptrData->GetNumberOfTuples() = " <<
297  // ptrData->GetNumberOfTuples() << std::endl ;
298  for(vtkIdType id_tuple_in_array = 0;
299  id_tuple_in_array < ptr_data->GetNumberOfTuples();
300  id_tuple_in_array++)
301  {
302  // field_attributes[idArray] is a 2D array
303  field_attributes[id_array + old_nb_e][id_tuple_in_array].resize(
304  static_cast< size_t >(ptr_data->GetNumberOfComponents()));
305  std::copy(
306  ptr_data->GetTuple(id_tuple_in_array),
307  ptr_data->GetTuple(id_tuple_in_array) +
308  ptr_data->GetNumberOfComponents(),
309 
310  field_attributes[id_array + old_nb_e][id_tuple_in_array].begin());
311  }
312  }
313  }
315  // Get cell data when present
316  vtkCellData *cd = poly_data->GetCellData();
317  if(cd != NULL)
318  {
320  // arrays (SCALARS are taken into account there)
321  size_t old_nb_e = field_attributes.size();
322  field_attributes.resize(static_cast< size_t >(
323  old_nb_e + cd->GetNumberOfArrays())); // use exactly the nb of 2D arrays
324  // for FIELD data
325  field_names.resize(
326  static_cast< size_t >(old_nb_e + cd->GetNumberOfArrays()));
327  for(int id_array = 0; id_array < cd->GetNumberOfArrays(); id_array++)
328  {
329  field_names[old_nb_e + id_array] =
330  std::string("CELL_DATA_") + cd->GetArrayName(id_array);
331  // std::cout << "CELL_DATA: array name is " << cd->GetArrayName(idArray)
332  // << std::endl;
333  // for each FIELD data (arrays)
334  vtkSmartPointer< vtkDataArray > ptr_data = cd->GetArray(id_array);
335 
336  field_attributes[old_nb_e + id_array].resize(
337  static_cast< size_t >(ptr_data->GetNumberOfTuples()));
338  // std::cout << "ptrData->GetNumberOfTuples() = " <<
339  // ptrData->GetNumberOfTuples() << std::endl ;
340  for(vtkIdType id_tuple_in_array = 0;
341  id_tuple_in_array < ptr_data->GetNumberOfTuples();
342  id_tuple_in_array++)
343  {
344  // field_attributes[oldNbE+idArray] is a 2D array
345  field_attributes[old_nb_e + id_array][id_tuple_in_array].resize(
346  static_cast< size_t >(ptr_data->GetNumberOfComponents()));
347  std::copy(
348  ptr_data->GetTuple(id_tuple_in_array),
349  ptr_data->GetTuple(id_tuple_in_array) +
350  ptr_data->GetNumberOfComponents(),
351 
352  field_attributes[old_nb_e + id_array][id_tuple_in_array].begin());
353  }
354  }
355  }
356  if((line_indices.size() > 0) && // line_indices are present
357  (line_indices[0].size() >
358  2) && // but at least 3 vertex for an edge (strange)
359  (face_indices.size() == 0) // and no face_indices provided...
360  )
361  { // line_indices in 2D can sometimes be misused as polygons by some users of
362  // MEPP (but for vtk line_indices refer to edges)
363  std::swap(line_indices, face_indices);
364  }
366  if(field_names.size() == field_attributes.size())
367  {
368  size_t cpt = 0, nb_e_glob = field_names.size();
369  for(; cpt < nb_e_glob; ++cpt)
370  {
371  if(field_names[cpt].find("colo") != std::string::npos)
372  {
373  if((vertex_color_coords.size() == 0) &&
374  ((field_names[cpt].find("vertex") != std::string::npos) ||
375  (field_names[cpt].find("point") != std::string::npos)))
376  {
377  size_t nb_e = points_coords.size();
378  vertex_color_coords.resize(nb_e);
379  for(size_t i = 0; i < nb_e; ++i)
380  {
381  size_t nb_c = field_attributes[cpt][i].size();
382  vertex_color_coords[i].resize(nb_c);
383  for(size_t j = 0; j < nb_c; ++j)
384  vertex_color_coords[i][j] =
385  static_cast< CoordCType >(field_attributes[cpt][i][j]);
386  // std::copy(field_attributes[cpt][i].begin(),
387  // field_attributes[cpt][i].end(), vertex_color_coords[i].begin()) ;
388  // std::copy(field_attributes[cpt][i].begin(),
389  // field_attributes[cpt][i].end(), std::back_inserter(
390  // vertex_color_coords[i] )) ;
391  }
392  }
393  else if((face_color_coords.size() == 0) &&
394  ((field_names[cpt].find("face") != std::string::npos) ||
395  (field_names[cpt].find("polygon") != std::string::npos) ||
396  (line_indices.size() == 0)))
397  {
398  size_t nb_e = face_indices.size();
399  face_color_coords.resize(nb_e);
400  for(size_t i = 0; i < nb_e; ++i)
401  {
402  size_t nb_c = field_attributes[cpt][i].size();
403  face_color_coords[i].resize(nb_c);
404  for(size_t j = 0; j < nb_c; ++j)
405  face_color_coords[i][j] =
406  static_cast< CoordCType >(field_attributes[cpt][i][j]);
407  // std::copy(field_attributes[cpt][i].begin(),
408  // field_attributes[cpt][i].end(), face_color_coords[i].begin()) ;
409  // std::copy(field_attributes[cpt][i].begin(),
410  // field_attributes[cpt][i].end(), std::back_inserter(
411  // face_color_coords[i] )) ;
412  }
413  }
414  else if((lines_color_coords.size() == 0) && (line_indices.size() > 0))
415  {
416  size_t nb_e = line_indices.size();
417  lines_color_coords.resize(nb_e);
418  for(size_t i = 0; i < nb_e; ++i)
419  {
420  size_t nb_c = field_attributes[cpt][i].size();
421  lines_color_coords[i].resize(nb_c);
422  for(size_t j = 0; j < nb_c; ++j)
423  lines_color_coords[i][j] =
424  static_cast< CoordCType >(field_attributes[cpt][i][j]);
425  // std::copy(field_attributes[cpt][i].begin(),
426  // field_attributes[cpt][i].end(), lines_color_coords[i].begin()) ;
427  // std::copy(field_attributes[cpt][i].begin(),
428  // field_attributes[cpt][i].end(), std::back_inserter(
429  // lines_color_coords[i] )) ;
430  }
431  }
432  field_attributes[cpt].clear(); // to not take into account twice colors!
433 
434  size_t last_index = nb_e_glob - 1;
435  if(cpt < last_index)
436  {
437  std::swap(field_attributes[cpt], field_attributes[last_index]);
438  std::swap(field_names[cpt], field_names[last_index]);
439  }
440  field_attributes.resize(last_index);
441  field_names.resize(last_index);
442  --nb_e_glob; // update end iterator
443  --cpt;
444  }
445  else if(field_names[cpt].find("norm") != std::string::npos)
446  { // we do not take into account normals twice!
447  field_attributes[cpt].clear(); // to not take into account twice colors!
448 
449  size_t last_index = nb_e_glob - 1;
450  if(cpt < last_index)
451  {
452  std::swap(field_attributes[cpt], field_attributes[last_index]);
453  std::swap(field_names[cpt], field_names[last_index]);
454  }
455  field_attributes.resize(last_index);
456  field_names.resize(last_index);
457  --nb_e_glob; // update end iterator
458  --cpt;
459  }
460  }
461  }
462 }
463 
464 template< typename CoordType,
465  typename CoordNType,
466  typename CoordCType,
467  typename IndexType >
468 void
469 read_vtk_unstructured_grid(
470  const vtkSmartPointer< vtkUnstructuredGrid > &unstructured_grid,
471  std::vector< std::vector< CoordType > > &points_coords,
472  std::vector< std::vector< CoordNType > >
473  &normals_coords,
474  std::vector< std::vector< CoordCType > > &vertex_color_coords,
476  std::vector< std::vector< IndexType > >
477  &line_indices,
478  std::vector< std::vector< CoordCType > > &lines_color_coords,
479  std::vector< std::vector< IndexType > >
480  &face_indices,
481  std::vector< std::vector< CoordCType > > &face_color_coords,
482  std::vector< std::vector< std::vector< double > > >
483  &field_attributes, // array of (2D) arrays: first point data, then cell
484  // data
485  std::vector< std::string > &field_names)
486 {
488  // clearing vectors...
489  points_coords.clear();
490  normals_coords.clear();
491  vertex_color_coords.clear();
492  line_indices.clear();
493  lines_color_coords.clear();
494  face_indices.clear();
495  face_color_coords.clear();
496  field_attributes.clear();
498 
499  // each point there is at maximum one normal (some points may have no normal):
500  assert(
501  (unstructured_grid->GetPointData()->GetNormals() == NULL) ||
502  (unstructured_grid->GetNumberOfPoints() >=
503  unstructured_grid->GetPointData()->GetNormals()->GetNumberOfTuples()));
505  // Get points and point normals
506  if(unstructured_grid->GetNumberOfPoints() > 0)
507  {
508  points_coords.reserve(static_cast< size_t >(
509  unstructured_grid
510  ->GetNumberOfPoints())); // use exactly the nb of points
511 
512  vtkSmartPointer< vtkPoints > ptr_points = unstructured_grid->GetPoints();
513  vtkSmartPointer< vtkDataArray > ptr_normals =
514  unstructured_grid->GetPointData()->GetNormals();
515  bool use_norm =
516  (ptr_normals != NULL &&
517  ptr_points->GetNumberOfPoints() ==
518  ptr_normals->GetNumberOfTuples()); // use normals only if one
519  // normal per point/vertex
520 
521  if(use_norm)
522  normals_coords.reserve(static_cast< size_t >(
523  unstructured_grid->GetPointData()
524  ->GetNormals()
525  ->GetNumberOfTuples())); // use exactly the nb of normals
526 
527  for(vtkIdType id = 0; id < ptr_points->GetNumberOfPoints(); id++)
528  {
529  std::vector< double > vd(ptr_points->GetPoint(id),
530  ptr_points->GetPoint(id) +
531  3); // ** 3 coordinates of type double
532  // hard-coded in vtkPoints class ** ;
533  std::vector< CoordType > final_vc;
534  final_vc.reserve(vd.size());
535  std::vector< double >::iterator it(vd.begin()), ite(vd.end());
536  for(; it != ite; ++it)
537  final_vc.push_back(static_cast< CoordType >(*it));
538 
539  points_coords.push_back(final_vc);
540 
541  if(use_norm)
542  {
543  std::vector< double > vdn(ptr_normals->GetTuple(id),
544  ptr_normals->GetTuple(id) +
545  ptr_normals->GetNumberOfComponents());
546  std::vector< CoordNType > final_vcn;
547  final_vcn.reserve(vdn.size());
548  it = vdn.begin();
549  ite = vdn.end();
550  for(; it != ite; ++it)
551  final_vcn.push_back(static_cast< CoordNType >(*it));
552 
553  normals_coords.push_back(final_vcn);
554  }
555  }
556  }
558  // Get facet indices
559  if(unstructured_grid->GetNumberOfCells() > 0)
560  {
561  face_indices.reserve(static_cast< size_t >(
562  unstructured_grid->GetNumberOfCells())); // use exactly the nb of faces
563  line_indices.reserve(
564  static_cast< size_t >(unstructured_grid->GetNumberOfCells()));
565 
566  size_t nb_face_cell = 0, nb_vol_cell = 0;
567 
568  vtkSmartPointer< vtkCellArray > ptr_cell_polygons =
569  unstructured_grid->GetCells();
570  vtkIdType npts;
571 #if (VTK_MAJOR_VERSION >= 9) // not tested with VTK 8...
572  const vtkIdType *pts_poly = new vtkIdType[ptr_cell_polygons->GetMaxCellSize()];
573 #else
574  vtkIdType *pts_poly = new vtkIdType[ptr_cell_polygons->GetMaxCellSize()];
575 #endif
576 
577  ptr_cell_polygons->InitTraversal();
578  vtkIdType cell_id = 0;
579  while(ptr_cell_polygons->GetNextCell(
580  npts, pts_poly)) // 0 (i.e. false) is encountered at the list' end
581  {
582  std::vector< IndexType > facet_indices(npts);
583  for(vtkIdType i = 0; i < npts; i++)
584  {
585  facet_indices[i] = static_cast< IndexType >(pts_poly[i]);
586  // std::cout << static_cast<Index>( ptsPoly[i] ) << std::endl ;
587  }
588 
589  switch(unstructured_grid->GetCellType(cell_id))
590  {
591  case VTK_VERTEX: // vertex (nothing to do)
592  throw std::runtime_error("Reader::read_vtkUnstructuredGrid -> vertex "
593  "cells are not managed.");
594  // break;
595  case VTK_LINE: // line/edge (nothing to do, because edge alone are not of
596  // interest for numerical simulation)
597  throw std::runtime_error(
598  "Reader::read_vtkUnstructuredGrid -> edge cells are not managed.");
599  // break;
600  case VTK_TRIANGLE: // triangle
601  face_indices.push_back(facet_indices);
602  nb_face_cell++;
603  break;
604  case VTK_POLYGON: // polygon
605  face_indices.push_back(facet_indices);
606  nb_face_cell++;
607  break;
608  case VTK_QUAD: // quad
609  face_indices.push_back(facet_indices);
610  nb_face_cell++;
611  break;
612  case VTK_TETRA: // tetra
613  {
614  line_indices.push_back(facet_indices);
615  nb_vol_cell++;
616  std::vector< IndexType > face1, face2, face3, face4;
617  face1.push_back(facet_indices[0]);
618  face1.push_back(facet_indices[1]);
619  face1.push_back(facet_indices[3]);
620  face_indices.push_back(face1);
621  nb_face_cell++; // ok
622 
623  face2.push_back(facet_indices[1]);
624  face2.push_back(facet_indices[2]);
625  face2.push_back(facet_indices[3]);
626  face_indices.push_back(face2);
627  nb_face_cell++; // ok
628 
629  face3.push_back(facet_indices[3]);
630  face3.push_back(facet_indices[2]);
631  face3.push_back(facet_indices[0]);
632  face_indices.push_back(face3);
633  nb_face_cell++; // ok
634 
635  face4.push_back(facet_indices[2]);
636  face4.push_back(facet_indices[1]);
637  face4.push_back(facet_indices[0]);
638  face_indices.push_back(face4);
639  nb_face_cell++; // ok
640  }
641  break;
642  case VTK_HEXAHEDRON: // hexa
643  {
644  line_indices.push_back(facet_indices);
645  nb_vol_cell++;
646  std::vector< IndexType > face1, face2, face3, face4, face5, face6;
647  face1.push_back(facet_indices[0]);
648  face1.push_back(facet_indices[1]);
649  face1.push_back(facet_indices[5]);
650  face1.push_back(facet_indices[4]);
651  face_indices.push_back(face1);
652  nb_face_cell++; // ok
653 
654  face2.push_back(facet_indices[1]);
655  face2.push_back(facet_indices[2]);
656  face2.push_back(facet_indices[6]);
657  face2.push_back(facet_indices[5]);
658  face_indices.push_back(face2);
659  nb_face_cell++; // ok
660 
661  face3.push_back(facet_indices[2]);
662  face3.push_back(facet_indices[3]);
663  face3.push_back(facet_indices[7]);
664  face3.push_back(facet_indices[6]);
665  face_indices.push_back(face3);
666  nb_face_cell++; // ok
667 
668  face4.push_back(facet_indices[3]);
669  face4.push_back(facet_indices[0]);
670  face4.push_back(facet_indices[4]);
671  face4.push_back(facet_indices[7]);
672  face_indices.push_back(face4);
673  nb_face_cell++; // ok
674 
675  face5.push_back(facet_indices[4]);
676  face5.push_back(facet_indices[5]);
677  face5.push_back(facet_indices[6]);
678  face5.push_back(facet_indices[7]);
679  face_indices.push_back(face5);
680  nb_face_cell++; // ok
681 
682  face6.push_back(facet_indices[3]);
683  face6.push_back(facet_indices[2]);
684  face6.push_back(facet_indices[1]);
685  face6.push_back(facet_indices[0]);
686  face_indices.push_back(face6);
687  nb_face_cell++; // ok
688  }
689  break;
690  /*
691  case VTK_WEDGE: //wedge/triangular prism
692  line_indices.push_back( facetIndices ) ; nbVolCell++; break ;
693  case VTK_PYRAMID: //pyramid
694  line_indices.push_back( facetIndices ) ; nbVolCell++; break ;*/
695  default:
696  throw std::runtime_error(
697  "Reader::read_vtkUnstructuredGrid -> cell type is not known.");
698  }
699 
700  cell_id++;
701  }
702  face_indices.resize(nb_face_cell);
703  line_indices.resize(nb_vol_cell);
704 
705  delete[] pts_poly;
706  }
708  // Get 2D arrays of FIELD data when present
709  vtkFieldData *fd = unstructured_grid->GetFieldData();
710  if(fd != NULL)
711  {
712  field_attributes.resize(static_cast< size_t >(
713  fd->GetNumberOfArrays())); // use exactly the nb of 2D arrays for FIELD
714  // data
715  field_names.resize(static_cast< size_t >(fd->GetNumberOfArrays()));
716 
717  for(int id_array = 0; id_array < fd->GetNumberOfArrays(); id_array++)
718  {
719  field_names[id_array] = fd->GetArrayName(id_array);
720  // for each FIELD data (arrays)
721  vtkSmartPointer< vtkDataArray > ptr_data = fd->GetArray(
722  id_array); // poly_data->GetFieldData()->GetArrayName(idArray) ) ;
723 
724  field_attributes[id_array].resize(
725  static_cast< size_t >(ptr_data->GetNumberOfTuples()));
726  // std::cout << "ptrData->GetNumberOfTuples() = " <<
727  // ptrData->GetNumberOfTuples() << std::endl ;
728  for(vtkIdType id_tuple_in_array = 0;
729  id_tuple_in_array < ptr_data->GetNumberOfTuples();
730  id_tuple_in_array++)
731  {
732  // field_attributes[idArray] is a 2D array
733  field_attributes[id_array][id_tuple_in_array].resize(
734  static_cast< size_t >(ptr_data->GetNumberOfComponents()));
735  std::copy(ptr_data->GetTuple(id_tuple_in_array),
736  ptr_data->GetTuple(id_tuple_in_array) +
737  ptr_data->GetNumberOfComponents(),
738 
739  field_attributes[id_array][id_tuple_in_array].begin());
740  }
741  }
742  }
744  // Get point data when present
745  vtkPointData *pd = unstructured_grid->GetPointData();
746  if(pd != NULL)
747  {
748  size_t old_nb_e = field_attributes.size();
749  field_attributes.resize(static_cast< size_t >(
750  old_nb_e + pd->GetNumberOfArrays())); // use exactly the nb of 2D arrays
751  // for FIELD data
752  field_names.resize(
753  static_cast< size_t >(old_nb_e + pd->GetNumberOfArrays()));
754  for(int id_array = 0; id_array < pd->GetNumberOfArrays(); id_array++)
755  {
756  field_names[id_array + old_nb_e] =
757  std::string("POINT_DATA_") + pd->GetArrayName(id_array);
758  // std::cout << "POINT_DATA: array name is " << pd->GetArrayName(idArray)
759  // << std::endl;
760  // for each FIELD data (arrays)
761  vtkSmartPointer< vtkDataArray > ptr_data = pd->GetArray(id_array);
762 
763  field_attributes[id_array + old_nb_e].resize(
764  static_cast< size_t >(ptr_data->GetNumberOfTuples()));
765  // std::cout << "ptrData->GetNumberOfTuples() = " <<
766  // ptrData->GetNumberOfTuples() << std::endl ;
767  for(vtkIdType id_tuple_in_array = 0;
768  id_tuple_in_array < ptr_data->GetNumberOfTuples();
769  id_tuple_in_array++)
770  {
771  // field_attributes[idArray] is a 2D array
772  field_attributes[id_array + old_nb_e][id_tuple_in_array].resize(
773  static_cast< size_t >(ptr_data->GetNumberOfComponents()));
774  std::copy(
775  ptr_data->GetTuple(id_tuple_in_array),
776  ptr_data->GetTuple(id_tuple_in_array) +
777  ptr_data->GetNumberOfComponents(),
778 
779  field_attributes[id_array + old_nb_e][id_tuple_in_array].begin());
780  }
781  }
782  }
784  // Get cell data when present
785  vtkCellData *cd = unstructured_grid->GetCellData();
786  if(cd != NULL)
787  {
789  // arrays (SCALARS are taken into account there)
790  size_t old_nb_e = field_attributes.size();
791  field_attributes.resize(static_cast< size_t >(
792  old_nb_e + cd->GetNumberOfArrays())); // use exactly the nb of 2D arrays
793  // for FIELD data
794  field_names.resize(
795  static_cast< size_t >(old_nb_e + cd->GetNumberOfArrays()));
796  for(int id_array = 0; id_array < cd->GetNumberOfArrays(); id_array++)
797  {
798  field_names[old_nb_e + id_array] =
799  std::string("CELL_DATA_") + cd->GetArrayName(id_array);
800  // std::cout << "CELL_DATA: array name is " << cd->GetArrayName(idArray)
801  // << std::endl;
802  // for each FIELD data (arrays)
803  vtkSmartPointer< vtkDataArray > ptr_data = cd->GetArray(id_array);
804 
805  field_attributes[old_nb_e + id_array].resize(
806  static_cast< size_t >(ptr_data->GetNumberOfTuples()));
807  // std::cout << "ptrData->GetNumberOfTuples() = " <<
808  // ptrData->GetNumberOfTuples() << std::endl ;
809  for(vtkIdType id_tuple_in_array = 0;
810  id_tuple_in_array < ptr_data->GetNumberOfTuples();
811  id_tuple_in_array++)
812  {
813  // field_attributes[oldNbE+idArray] is a 2D array
814  field_attributes[old_nb_e + id_array][id_tuple_in_array].resize(
815  static_cast< size_t >(ptr_data->GetNumberOfComponents()));
816  std::copy(
817  ptr_data->GetTuple(id_tuple_in_array),
818  ptr_data->GetTuple(id_tuple_in_array) +
819  ptr_data->GetNumberOfComponents(),
820 
821  field_attributes[old_nb_e + id_array][id_tuple_in_array].begin());
822  }
823  }
824  }
826  if(field_names.size() == field_attributes.size())
827  {
828  size_t cpt = 0, nb_e_glob = field_names.size();
829  for(; cpt < nb_e_glob; ++cpt)
830  {
831  if(field_names[cpt].find("colo") != std::string::npos)
832  {
833  if((vertex_color_coords.size() == 0) &&
834  ((field_names[cpt].find("vertex") != std::string::npos) ||
835  (field_names[cpt].find("point") != std::string::npos)))
836  {
837  size_t nb_e = points_coords.size();
838  vertex_color_coords.resize(nb_e);
839  for(size_t i = 0; i < nb_e; ++i)
840  {
841  size_t nb_c = field_attributes[cpt][i].size();
842  vertex_color_coords[i].resize(nb_c);
843  for(size_t j = 0; j < nb_c; ++j)
844  vertex_color_coords[i][j] =
845  static_cast< CoordCType >(field_attributes[cpt][i][j]);
846  // std::copy(field_attributes[cpt][i].begin(),
847  // field_attributes[cpt][i].end(), vertex_color_coords[i].begin()) ;
848  // std::copy(field_attributes[cpt][i].begin(),
849  // field_attributes[cpt][i].end(), std::back_inserter(
850  // vertex_color_coords[i] )) ;
851  }
852  }
853  else if((face_color_coords.size() == 0) &&
854  ((field_names[cpt].find("face") != std::string::npos) ||
855  (field_names[cpt].find("polygon") != std::string::npos) ||
856  (line_indices.size() == 0)))
857  {
858  size_t nb_e = face_indices.size();
859  face_color_coords.resize(nb_e);
860  for(size_t i = 0; i < nb_e; ++i)
861  {
862  size_t nb_c = field_attributes[cpt][i].size();
863  face_color_coords[i].resize(nb_c);
864  for(size_t j = 0; j < nb_c; ++j)
865  face_color_coords[i][j] =
866  static_cast< CoordCType >(field_attributes[cpt][i][j]);
867  // std::copy(field_attributes[cpt][i].begin(),
868  // field_attributes[cpt][i].end(), face_color_coords[i].begin()) ;
869  // std::copy(field_attributes[cpt][i].begin(),
870  // field_attributes[cpt][i].end(), std::back_inserter(
871  // face_color_coords[i] )) ;
872  }
873  }
874  else if((lines_color_coords.size() == 0) && (line_indices.size() > 0))
875  {
876  size_t nb_e = line_indices.size();
877  lines_color_coords.resize(nb_e);
878  for(size_t i = 0; i < nb_e; ++i)
879  {
880  size_t nb_c = field_attributes[cpt][i].size();
881  lines_color_coords[i].resize(nb_c);
882  for(size_t j = 0; j < nb_c; ++j)
883  lines_color_coords[i][j] =
884  static_cast< CoordCType >(field_attributes[cpt][i][j]);
885  // std::copy(field_attributes[cpt][i].begin(),
886  // field_attributes[cpt][i].end(), lines_color_coords[i].begin()) ;
887  // std::copy(field_attributes[cpt][i].begin(),
888  // field_attributes[cpt][i].end(), std::back_inserter(
889  // lines_color_coords[i] )) ;
890  }
891  }
892  field_attributes[cpt].clear(); // to not take into account twice colors!
893 
894  size_t last_index = nb_e_glob - 1;
895  if(cpt < last_index)
896  {
897  std::swap(field_attributes[cpt], field_attributes[last_index]);
898  std::swap(field_names[cpt], field_names[last_index]);
899  }
900  field_attributes.resize(last_index);
901  field_names.resize(last_index);
902  --nb_e_glob; // update end iterator
903  --cpt;
904  }
905  else if(field_names[cpt].find("norm") != std::string::npos)
906  { // we do not take into account normals twice!
907  field_attributes[cpt].clear(); // to not take into account twice colors!
908 
909  size_t last_index = nb_e_glob - 1;
910  if(cpt < last_index)
911  {
912  std::swap(field_attributes[cpt], field_attributes[last_index]);
913  std::swap(field_names[cpt], field_names[last_index]);
914  }
915  field_attributes.resize(last_index);
916  field_names.resize(last_index);
917  --nb_e_glob; // update end iterator
918  --cpt;
919  }
920  }
921  }
922 }
923 
924 template< typename CoordType,
925  typename CoordNType,
926  typename CoordCType,
927  typename IndexType >
928 void
929 read_vtk_or_vtp_or_vtu_file(
930  std::string file_path,
931  std::vector< std::vector< CoordType > > &points_coords,
932  std::vector< std::vector< CoordNType > >
933  &normals_coords,
934  std::vector< std::vector< CoordCType > > &vertex_color_coords,
936  std::vector< std::vector< IndexType > >
937  &line_indices,
938  std::vector< std::vector< CoordCType > > &lines_color_coords,
939  std::vector< std::vector< IndexType > >
940  &face_indices,
941  std::vector< std::vector< CoordCType > > &face_color_coords,
942  std::vector< std::vector< std::vector< double > > >
943  &field_attributes, // array of 2D arrays: first point data, then cell
944  // data
945  std::vector< std::string > &field_names,
946  bool if_unstructured_grid_then_take_surface = false)
947 {
949  // reading and getting polydata/unstructured mesh
950  bool is_a_poly_data = true;
951  vtkSmartPointer< vtkPolyData > poly_data;
952  vtkSmartPointer< vtkUnstructuredGrid > unstructured_grid;
953 
954  std::string dataset_line = get_line_containing_dataset(file_path);
955 
956  if(FileUtils::has_extension(file_path, ".vtp"))
957  {
958  vtkSmartPointer< vtkXMLPolyDataReader > reader =
959  vtkSmartPointer< vtkXMLPolyDataReader >::New();
960  reader->SetFileName(file_path.c_str());
961  reader->Update();
962 
963  poly_data = vtkSmartPointer< vtkPolyData >::New();
964  poly_data->DeepCopy(reader->GetOutput());
965  }
966  else if(FileUtils::has_extension(file_path, ".vtu"))
967  {
968  vtkSmartPointer< vtkXMLUnstructuredGridReader > reader =
969  vtkSmartPointer< vtkXMLUnstructuredGridReader >::New();
970  reader->SetFileName(file_path.c_str());
971  reader->Update();
972 
973  if(if_unstructured_grid_then_take_surface)
974  {
975  vtkSmartPointer< vtkDataSetSurfaceFilter > surface_filter =
976  vtkSmartPointer< vtkDataSetSurfaceFilter >::New();
977 #if VTK_MAJOR_VERSION <= 5
978  surfaceFilter->SetInput(reader->GetOutput());
979 #else
980  surface_filter->SetInputData(reader->GetOutput());
981 #endif
982  surface_filter->Update();
983 
984  poly_data = vtkSmartPointer< vtkPolyData >::New();
985  poly_data->DeepCopy(surface_filter->GetOutput());
986  }
987  else
988  {
989  is_a_poly_data = false;
990  unstructured_grid = vtkSmartPointer< vtkUnstructuredGrid >::New();
991  unstructured_grid->DeepCopy(reader->GetOutput());
992  }
993  }
994  else
995  {
996  if(FileUtils::has_extension(file_path, ".vtk"))
997  {
998  if(dataset_line.find("POLYDATA") != std::string::npos)
999  {
1000  vtkSmartPointer< vtkPolyDataReader > reader =
1001  vtkSmartPointer< vtkPolyDataReader >::New();
1002  reader->SetFileName(file_path.c_str());
1003  reader->Update();
1004 
1005  poly_data = vtkSmartPointer< vtkPolyData >::New();
1006  poly_data->DeepCopy(reader->GetOutput());
1007  }
1008  else
1009  {
1010  if(dataset_line.find("UNSTRUCTURED_GRID") != std::string::npos)
1011  {
1012  vtkSmartPointer< vtkUnstructuredGridReader > reader =
1013  vtkSmartPointer< vtkUnstructuredGridReader >::New();
1014  reader->SetFileName(file_path.c_str());
1015  reader->Update();
1016 
1017  if(if_unstructured_grid_then_take_surface)
1018  {
1019  vtkSmartPointer< vtkDataSetSurfaceFilter > surface_filter =
1020  vtkSmartPointer< vtkDataSetSurfaceFilter >::New();
1021 #if VTK_MAJOR_VERSION <= 5
1022  surfaceFilter->SetInput(reader->GetOutput());
1023 #else
1024  surface_filter->SetInputData(reader->GetOutput());
1025 #endif
1026  surface_filter->Update();
1027 
1028  poly_data = vtkSmartPointer< vtkPolyData >::New();
1029  poly_data->DeepCopy(surface_filter->GetOutput());
1030  }
1031  else
1032  {
1033  is_a_poly_data = false;
1034  unstructured_grid = vtkSmartPointer< vtkUnstructuredGrid >::New();
1035  unstructured_grid->DeepCopy(reader->GetOutput());
1036  }
1037  }
1038  else
1039  {
1040  throw std::runtime_error(
1041  "Reader::read_vtk_or_vtp_or_vtu_file -> file dataset type cannot "
1042  "be processed (neither POLYDATA nor UNSTRUCTURED_GRID)");
1043  }
1044  }
1045  }
1046  else
1047  {
1048  throw std::runtime_error(
1049  "Reader::read_vtk_or_vtp_or_vtu_file -> file extension is "
1050  "inappropriate (neither .vtk nor .vtp)");
1051  }
1052  }
1054  if(is_a_poly_data)
1055  return read_vtk_poly_data< CoordType,
1056  CoordNType,
1057  /*coordT_type,*/ CoordCType,
1058  IndexType >(poly_data,
1059  points_coords,
1060  normals_coords,
1061  vertex_color_coords,
1062  line_indices,
1063  lines_color_coords,
1064  face_indices,
1065  face_color_coords,
1066  field_attributes,
1067  field_names);
1068  else
1069  return read_vtk_unstructured_grid< CoordType,
1070  CoordNType,
1071  /*coordT_type,*/ CoordCType,
1072  IndexType >(unstructured_grid,
1073  points_coords,
1074  normals_coords,
1075  vertex_color_coords,
1076  line_indices,
1077  lines_color_coords,
1078  face_indices,
1079  face_color_coords,
1080  field_attributes,
1081  field_names);
1083 }
1084 
1085 
1086 } // namespace IO
1087 } // namespace FEVV
1088 
FEVV::IO::read_vtk_or_vtp_or_vtu_file
void read_vtk_or_vtp_or_vtu_file(std::string file_path, std::vector< std::vector< CoordType > > &points_coords, std::vector< std::vector< CoordNType > > &normals_coords, std::vector< std::vector< CoordCType > > &vertex_color_coords, std::vector< std::vector< IndexType > > &line_indices, std::vector< std::vector< CoordCType > > &lines_color_coords, std::vector< std::vector< IndexType > > &face_indices, std::vector< std::vector< CoordCType > > &face_color_coords, std::vector< std::vector< std::vector< double > > > &field_attributes, std::vector< std::string > &field_names, bool if_unstructured_grid_then_take_surface=false)
Definition: VtkFileReader.h:929
FEVV::IO::read_vtk_unstructured_grid
void read_vtk_unstructured_grid(const vtkSmartPointer< vtkUnstructuredGrid > &unstructured_grid, std::vector< std::vector< CoordType > > &points_coords, std::vector< std::vector< CoordNType > > &normals_coords, std::vector< std::vector< CoordCType > > &vertex_color_coords, std::vector< std::vector< IndexType > > &line_indices, std::vector< std::vector< CoordCType > > &lines_color_coords, std::vector< std::vector< IndexType > > &face_indices, std::vector< std::vector< CoordCType > > &face_color_coords, std::vector< std::vector< std::vector< double > > > &field_attributes, std::vector< std::string > &field_names)
Definition: VtkFileReader.h:469
FEVV
Interfaces for plugins These interfaces will be used for different plugins.
Definition: Assert.h:16
FEVV::IO::get_line_containing_dataset
std::string get_line_containing_dataset(std::string file_path)
Definition: VtkFileReader.h:58
StringUtilities.hpp
FEVV::FileUtils::has_extension
bool has_extension(const std::string &file_name)
Definition: FileUtilities.hpp:58
FEVV::IO::read_vtk_poly_data
void read_vtk_poly_data(const vtkSmartPointer< vtkPolyData > &poly_data, std::vector< std::vector< CoordType > > &points_coords, std::vector< std::vector< CoordNType > > &normals_coords, std::vector< std::vector< CoordCType > > &vertex_color_coords, std::vector< std::vector< IndexType > > &line_indices, std::vector< std::vector< CoordCType > > &lines_color_coords, std::vector< std::vector< IndexType > > &face_indices, std::vector< std::vector< CoordCType > > &face_color_coords, std::vector< std::vector< std::vector< double > > > &field_attributes, std::vector< std::string > &field_names)
Definition: VtkFileReader.h:88
FileUtilities.hpp
Generated by doxygen 1.8.20