MEPP2 Project
VtkFileWriter.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 <vtkSmartPointer.h>
28 #include <vtkPointData.h>
29 #include <vtkCellData.h>
30 #include <vtkCellArray.h>
31 #include <vtkPolyData.h>
32 #include <vtkUnstructuredGrid.h>
33 #include <vtkDoubleArray.h>
34 #include <vtkFieldData.h>
35 
36 #include <vtkPolyDataWriter.h> // for vtk-files
37 #include <vtkXMLPolyDataWriter.h> // for vtp-files
38 #include <vtkUnstructuredGridWriter.h> // for vtu-files
39 #include <vtkXMLUnstructuredGridWriter.h>
40 #if defined(__APPLE__) && defined(__clang__)
41 #pragma clang diagnostic pop
42 #endif
43 
44 namespace FEVV {
45 namespace IO {
46 template< typename CoordType,
47  typename CoordNType,
48  /*typename coordT_type,*/ typename CoordCType,
49  typename IndexType >
50 void
51 load_vtk_poly_data(
52  const std::vector< std::vector< CoordType > > &points_coords,
53  const std::vector< std::vector< CoordNType > > &normals_coords,
54  const std::vector< std::vector< CoordCType > > &vertex_color_coords,
55  const std::vector< std::vector< IndexType > >
56  &line_indices,
57  const std::vector< std::vector< CoordCType > > &lines_color_coords,
59  const std::vector< std::vector< IndexType > >
60  &face_indices,
61  const std::vector< std::vector< CoordCType > > &face_color_coords,
62  const std::vector< std::vector< std::vector< double > > >
63  &field_attributes, // array of 2D arrays
64  const std::vector< std::string > &field_names,
65  vtkSmartPointer< vtkPolyData > &poly_data)
66 {
68  size_t nb_points = points_coords.size(), nb_normals = normals_coords.size(),
69  nb_point_colors = vertex_color_coords.size(),
70  nb_lines = line_indices.size(),
71  nb_line_colors = lines_color_coords.size(),
72  nb_faces = face_indices.size(),
73  nb_face_colors = face_color_coords.size(),
74  nb_arrays_in_field = field_attributes.size();
75  bool has_one_name_for_each_field = (nb_arrays_in_field == field_names.size());
77  // some checks
78  assert((nb_normals == 0u) || (nb_points == nb_normals));
80  poly_data = vtkSmartPointer< vtkPolyData >::New();
82  // add points
83  // std::cout << " before adding points..." << std::endl ;
84  vtkSmartPointer< vtkPoints > ptr_points = vtkSmartPointer< vtkPoints >::New();
85  for(size_t i = 0; i < nb_points;
86  ++i) // vtk writer will align point coordinates using 9 columns (so 3 3D
87  // points for the first line etc.)
88  {
89  // ptrPoints->InsertPoint( static_cast<vtkIdType>(i), points_coords[i][0],
90  // points_coords[i][1], points_coords[i][2] );
91  ptr_points->InsertNextPoint(
92  points_coords[i][0], points_coords[i][1], points_coords[i][2]);
93  }
94  // ptrPoints->Modified();
95  poly_data->SetPoints(ptr_points);
96  // std::cout << " after adding points..." << std::endl ;
98  // add point normals
99  // std::cout << " before adding point normals..." << std::endl ;
100  if((nb_normals > 0u) &&
101  (nb_points ==
102  nb_normals)) // vtk writer will align normal coordinates using 9 columns
103  // (so 3 3D normals for the first line etc.)
104  {
105  vtkSmartPointer< vtkDoubleArray > ptr_normals =
106  vtkSmartPointer< vtkDoubleArray >::New();
107  ptr_normals->SetName("vertex_normals"); // set the point normal name (the
108  // default name is "normals")
109  ptr_normals->SetNumberOfComponents(3); // 3d normals (ie x,y,z)
110  ptr_normals->SetNumberOfTuples(poly_data->GetNumberOfPoints());
111 
112  for(size_t i = 0; i < nb_normals; ++i)
113  {
114  ptr_normals->SetTuple3(
115  i, normals_coords[i][0], normals_coords[i][1], normals_coords[i][2]);
116  // ptrNormals->InsertTuple3(i, normals_coords[i][0], normals_coords[i][1],
117  // normals_coords[i][2] ) ; ptrNormals->InsertNextTuple3(
118  // normals_coords[i][0], normals_coords[i][1], normals_coords[i][2] );
119  }
120 
121  poly_data->GetPointData()->SetNormals(ptr_normals);
122  // std::cout << " nb of points = " << poly_data->GetNumberOfPoints() ;
123  // std::cout << " ; nb of point normals = " <<
124  // poly_data->GetPointData()->GetNormals()->GetNumberOfTuples() ;
125  assert(poly_data->GetNumberOfPoints() ==
126  poly_data->GetPointData()->GetNormals()->GetNumberOfTuples());
127  }
128  // std::cout << " after adding point normals..." << std::endl ;
130  // add point colors
131  if((nb_point_colors > 0u) && (vertex_color_coords[0].size() > 2u) &&
132  (nb_points == nb_point_colors))
133  {
134  vtkSmartPointer< vtkDoubleArray > ptr_colors =
135  vtkSmartPointer< vtkDoubleArray >::New();
136  ptr_colors->SetName("vertex_colors"); // set the point color name
137  ptr_colors->SetNumberOfComponents(
138  static_cast< int >(vertex_color_coords[0].size())); // 3d or 4d
139  ptr_colors->SetNumberOfTuples(nb_point_colors);
140 
141  for(size_t i = 0; i < nb_point_colors; ++i)
142  {
143  if(vertex_color_coords[0].size() == 3)
144  ptr_colors->SetTuple3(i,
145  vertex_color_coords[i][0],
146  vertex_color_coords[i][1],
147  vertex_color_coords[i][2]);
148  else
149  ptr_colors->SetTuple4(i,
150  vertex_color_coords[i][0],
151  vertex_color_coords[i][1],
152  vertex_color_coords[i][2],
153  vertex_color_coords[i][3]);
154  }
155 
156  poly_data->GetPointData()->AddArray(ptr_colors);
157  }
158  // std::cout << " after adding point colors..." << std::endl ;
160  // add polygons
161  // std::cout << " before adding polygons..." << std::endl ;
162  vtkSmartPointer< vtkCellArray > polys =
163  vtkSmartPointer< vtkCellArray >::New();
164  for(size_t i = 0; i < nb_faces; ++i)
165  {
166  vtkIdType *pt_indices = new vtkIdType[face_indices[i].size()];
167  typename std::vector< IndexType >::const_iterator it_begin(
168  face_indices[i].begin()),
169  it_end(face_indices[i].end());
170  unsigned long cpt = 0;
171  while(it_begin != it_end)
172  {
173  pt_indices[cpt++] = static_cast< vtkIdType >(*it_begin);
174 
175  ++it_begin;
176  }
177  polys->InsertNextCell(static_cast< int >(face_indices[i].size()),
178  pt_indices);
179 
180  delete[] pt_indices;
181  }
182  poly_data->SetPolys(polys);
183  // std::cout << " after adding polygons..." << std::endl ;
185  // add polygon colors
186  if((nb_face_colors > 0u) && (face_color_coords[0].size() > 2u) &&
187  (nb_faces == nb_face_colors))
188  {
189  vtkSmartPointer< vtkDoubleArray > ptr_colors =
190  vtkSmartPointer< vtkDoubleArray >::New();
191  ptr_colors->SetName("face_colors"); // set the point color name
192  ptr_colors->SetNumberOfComponents(
193  static_cast< int >(face_color_coords[0].size())); // 3d or 4d
194  ptr_colors->SetNumberOfTuples(nb_face_colors);
195 
196  for(size_t i = 0; i < nb_face_colors; ++i)
197  {
198  if(face_color_coords[0].size() == 3)
199  ptr_colors->SetTuple3(i,
200  face_color_coords[i][0],
201  face_color_coords[i][1],
202  face_color_coords[i][2]);
203  else
204  ptr_colors->SetTuple4(i,
205  face_color_coords[i][0],
206  face_color_coords[i][1],
207  face_color_coords[i][2],
208  face_color_coords[i][3]);
209  }
210  if(nb_lines == 0)
211  poly_data->GetCellData()->AddArray(ptr_colors);
212  else
213  poly_data->GetFieldData()->AddArray(ptr_colors);
214  }
215  // std::cout << " after adding face colors..." << std::endl ;
217  // add lines
218  vtkSmartPointer< vtkCellArray > lines_of_volume_cell =
219  vtkSmartPointer< vtkCellArray >::New();
220  for(size_t i = 0; i < nb_lines; ++i)
221  {
222  vtkIdType *pt_indices = new vtkIdType[line_indices[i].size()];
223  typename std::vector< IndexType >::const_iterator it_begin(
224  line_indices[i].begin()),
225  it_end(line_indices[i].end());
226  unsigned long cpt = 0;
227  while(it_begin != it_end)
228  {
229  pt_indices[cpt++] = static_cast< vtkIdType >(*it_begin);
230 
231  ++it_begin;
232  }
233  lines_of_volume_cell->InsertNextCell(
234  static_cast< int >(line_indices[i].size()), pt_indices);
235 
236  delete[] pt_indices;
237  }
238  poly_data->SetLines(lines_of_volume_cell);
240  // add lines colors
241  if((nb_line_colors > 0u) && (lines_color_coords[0].size() > 2u) &&
242  (nb_lines == nb_line_colors))
243  {
244  vtkSmartPointer< vtkDoubleArray > ptr_colors =
245  vtkSmartPointer< vtkDoubleArray >::New();
246  if(line_indices[0].size() > 2)
247  ptr_colors->SetName("cell_colors");
248  else
249  ptr_colors->SetName("edge_colors");
250 
251  ptr_colors->SetNumberOfComponents(
252  static_cast< int >(lines_color_coords[0].size())); // 3d or 4d
253  ptr_colors->SetNumberOfTuples(nb_line_colors);
254 
255  for(size_t i = 0; i < nb_line_colors; ++i)
256  {
257  if(lines_color_coords[0].size() == 3)
258  ptr_colors->SetTuple3(i,
259  lines_color_coords[i][0],
260  lines_color_coords[i][1],
261  lines_color_coords[i][2]);
262  else
263  ptr_colors->SetTuple4(i,
264  lines_color_coords[i][0],
265  lines_color_coords[i][1],
266  lines_color_coords[i][2],
267  lines_color_coords[i][3]);
268  }
269  if(nb_faces == 0)
270  poly_data->GetCellData()->AddArray(ptr_colors);
271  else
272  poly_data->GetFieldData()->AddArray(ptr_colors);
273  }
275  // Set 2D arrays of FIELD data when present
276  // std::cout << " before adding data fields..." << std::endl ;
277  if(nb_arrays_in_field > 0u) // vtk writer will align arrays using 9 columns
278  {
279  int choice = 0;
280  // std::cout << nbArraysInField << " data fields to add..." << std::endl ;
281  for(size_t id_array = 0; id_array < nb_arrays_in_field; id_array++)
282  {
283  // std::cout << " data fields number " << id_array << std::endl ;
284  // std::cout << field_names[id_array] << std::endl;
285  vtkSmartPointer< vtkDoubleArray > ptr_data =
286  vtkSmartPointer< vtkDoubleArray >::New();
287  if(has_one_name_for_each_field)
288  {
289  if(field_names[id_array].find("POINT_DATA_") != std::string::npos)
290  {
291  ptr_data->SetName((field_names[id_array].substr(11)).c_str());
292  choice = 1;
293  }
294  else if(field_names[id_array].find("CELL_DATA_") != std::string::npos)
295  {
296  ptr_data->SetName((field_names[id_array].substr(10)).c_str());
297  choice = 2;
298  }
299  else if(field_names[id_array].find("ELEMENT_DATA_") != std::string::npos)
300  {
301  ptr_data->SetName((field_names[id_array].substr(13)).c_str());
302  choice = 2;
303  }
304  else
305  {
306  ptr_data->SetName(field_names[id_array].c_str());
307  choice = 0;
308  }
309  }
310  else
311  {
312  ptr_data->SetName(
313  (std::string("cell_law_") + StrUtils::convert(id_array))
314  .c_str()); // name of array id_array
315  choice = 0;
316  }
317  ptr_data->SetNumberOfComponents(
318  ((field_attributes[id_array].size() > 0u)
319  ? static_cast< int >(field_attributes[id_array][0].size())
320  : 0));
321  ptr_data->SetNumberOfTuples(field_attributes[id_array].size());
322 
323  // std::cout << " it has " << ptrData->GetNumberOfTuples() << " tuples and
324  // " << ptrData->GetNumberOfComponents() << " components" << std::endl ;
325  // Set 2D array values:
326  for(vtkIdType i = 0; i < ptr_data->GetNumberOfTuples();
327  ++i) // for each line
328  {
329  assert(i < static_cast< vtkIdType >(field_attributes[id_array].size()));
330  for(vtkIdType j = 0; j < ptr_data->GetNumberOfComponents();
331  ++j) // for each column
332  {
333  if((j >=
334  static_cast< vtkIdType >(field_attributes[id_array][i].size())) ||
335  (ptr_data->GetNumberOfComponents() !=
336  static_cast< vtkIdType >(field_attributes[id_array][i].size())))
337  {
338  std::cerr << " Trying to access index " << j << " but only "
339  << field_attributes[id_array][i].size()
340  << " components in the data field." << std::endl;
341  std::cerr << " We should find " << ptr_data->GetNumberOfComponents()
342  << " components in the provided data field. "
343  << std::endl;
344  assert(false);
345  }
346  // if( field_attributes[id_array][i][j]<0. ||
347  // field_attributes[id_array][i][j]>1.0 )
348  //{
349  // std::cout << "Warning for GenSim project: attributes must be
350  //within [0;1] interval for the time being. " << std::endl ;
351  //}
352  ptr_data->SetComponent(i, j, field_attributes[id_array][i][j]);
353  }
354  }
356  switch(choice)
357  {
358  case 1:
359  poly_data->GetPointData()->AddArray(ptr_data);
360  break;
361  case 2:
362  poly_data->GetCellData()->AddArray(ptr_data);
363  break;
364  default:
365  poly_data->GetFieldData()->AddArray(ptr_data);
366  }
367  }
368  }
369  // std::cout << " after adding data fields..." << std::endl ;
370 }
371 template< typename CoordType,
372  typename CoordNType,
373  /*typename coordT_type,*/ typename CoordCType,
374  typename IndexType >
375 void
376 load_vtk_unstructured_grid(
377  const std::vector< std::vector< CoordType > > &points_coords,
378  const std::vector< std::vector< CoordNType > > &normals_coords,
379  const std::vector< std::vector< CoordCType > > &vertex_color_coords,
380  const std::vector< std::vector< IndexType > >
381  &line_indices,
382  const std::vector< std::vector< CoordCType > > &lines_color_coords,
384  const std::vector< std::vector< IndexType > >
385  &face_indices,
386  const std::vector< std::vector< CoordCType > > &face_color_coords,
387  const std::vector< std::vector< std::vector< double > > >
388  &field_attributes, // array of 2D arrays
389  const std::vector< std::string > &field_names,
390  vtkSmartPointer< vtkUnstructuredGrid > &unstructured_grid)
391 {
393  size_t nb_points = points_coords.size(), nb_normals = normals_coords.size(),
394  nb_point_colors = vertex_color_coords.size(),
395  nb_lines = line_indices.size(),
396  nb_line_colors = lines_color_coords.size(),
397  nb_faces = face_indices.size(),
398  nb_face_colors = face_color_coords.size(),
399  nb_arrays_in_field = field_attributes.size();
400  bool has_one_name_for_each_field = (nb_arrays_in_field == field_names.size());
402  // some checks
403  assert((nb_normals == 0u) || (nb_points == nb_normals));
405  unstructured_grid = vtkSmartPointer< vtkUnstructuredGrid >::New();
407  // add points
408  // std::cout << " before adding points..." << std::endl ;
409  vtkSmartPointer< vtkPoints > ptr_points = vtkSmartPointer< vtkPoints >::New();
410  for(size_t i = 0; i < nb_points;
411  ++i) // vtk writer will align point coordinates using 9 columns (so 3 3D
412  // points for the first line etc.)
413  {
414  // ptrPoints->InsertPoint( static_cast<vtkIdType>(i), points_coords[i][0],
415  // points_coords[i][1], points_coords[i][2] );
416  ptr_points->InsertNextPoint(
417  points_coords[i][0], points_coords[i][1], points_coords[i][2]);
418  }
419  unstructured_grid->SetPoints(ptr_points);
420  // std::cout << " after adding points..." << std::endl ;
422  // add point normals
423  // std::cout << " before adding point normals..." << std::endl ;
424  if((nb_normals > 0u) &&
425  (nb_points ==
426  nb_normals)) // vtk writer will align normal coordinates using 9 columns
427  // (so 3 3D normals for the first line etc.)
428  {
429  vtkSmartPointer< vtkDoubleArray > ptr_normals =
430  vtkSmartPointer< vtkDoubleArray >::New();
431  ptr_normals->SetName("vertex_normals"); // set the point normal name (the
432  // default name is "normals")
433  ptr_normals->SetNumberOfComponents(3); // 3d normals (ie x,y,z)
434  ptr_normals->SetNumberOfTuples(unstructured_grid->GetNumberOfPoints());
435 
436  for(size_t i = 0; i < nb_normals; ++i)
437  {
438  ptr_normals->SetTuple3(
439  i, normals_coords[i][0], normals_coords[i][1], normals_coords[i][2]);
440  // ptrNormals->InsertTuple3(i, normals_coords[i][0], normals_coords[i][1],
441  // normals_coords[i][2] ) ; ptrNormals->InsertNextTuple3(
442  // normals_coords[i][0], normals_coords[i][1], normals_coords[i][2] );
443  }
444 
445  unstructured_grid->GetPointData()->SetNormals(ptr_normals);
446  // std::cout << " nb of points = " << poly_data->GetNumberOfPoints() ;
447  // std::cout << " ; nb of point normals = " <<
448  // poly_data->GetPointData()->GetNormals()->GetNumberOfTuples() ;
449  assert(
450  unstructured_grid->GetNumberOfPoints() ==
451  unstructured_grid->GetPointData()->GetNormals()->GetNumberOfTuples());
452  }
453  // std::cout << " after adding point normals..." << std::endl ;
455  // add point colors
456  if((nb_point_colors > 0u) && (vertex_color_coords[0].size() > 2u) &&
457  (nb_points == nb_point_colors))
458  {
459  vtkSmartPointer< vtkDoubleArray > ptr_colors =
460  vtkSmartPointer< vtkDoubleArray >::New();
461  ptr_colors->SetName("vertex_colors"); // set the point color name
462  ptr_colors->SetNumberOfComponents(
463  static_cast< int >(vertex_color_coords[0].size())); // 3d or 4d
464  ptr_colors->SetNumberOfTuples(nb_point_colors);
465 
466  for(size_t i = 0; i < nb_point_colors; ++i)
467  {
468  if(vertex_color_coords[0].size() == 3)
469  ptr_colors->SetTuple3(i,
470  vertex_color_coords[i][0],
471  vertex_color_coords[i][1],
472  vertex_color_coords[i][2]);
473  else
474  ptr_colors->SetTuple4(i,
475  vertex_color_coords[i][0],
476  vertex_color_coords[i][1],
477  vertex_color_coords[i][2],
478  vertex_color_coords[i][3]);
479  }
480 
481  unstructured_grid->GetPointData()->AddArray(ptr_colors);
482  }
483  // std::cout << " after adding point colors..." << std::endl ;
485  // add polygons
486  // std::cout << " before adding polygons..." << std::endl ;
487  int *type_array = NULL;
488  if(nb_faces > 0 || nb_lines > 0)
489  type_array = new(std::nothrow) int[nb_faces + nb_lines];
490  vtkSmartPointer< vtkCellArray > cells =
491  vtkSmartPointer< vtkCellArray >::New();
492  for(size_t i = 0; i < nb_faces; ++i)
493  {
494  vtkIdType *pt_indices = new vtkIdType[face_indices[i].size()];
495  typename std::vector< IndexType >::const_iterator it_begin(
496  face_indices[i].begin()),
497  it_end(face_indices[i].end());
498  unsigned long cpt = 0;
499  while(it_begin != it_end)
500  {
501  pt_indices[cpt++] = static_cast< vtkIdType >(*it_begin);
502 
503  ++it_begin;
504  }
505  cells->InsertNextCell(static_cast< int >(face_indices[i].size()),
506  pt_indices);
507 
508  switch(face_indices[i].size())
509  {
510  case 0:
511  case 1:
512  case 2:
513  throw std::runtime_error(
514  "Writer::load_vtkUnstructuredGrid: a face seems to be not valid.");
515  case 3:
516  type_array[i] = VTK_TRIANGLE;
517  break;
518  case 4:
519  type_array[i] = VTK_QUAD;
520  break;
521  default:
522  type_array[i] = VTK_POLYGON;
523  }
524  delete[] pt_indices;
525  }
527  // add polygon colors
528  if((nb_face_colors > 0u) && (face_color_coords[0].size() > 2u) &&
529  (nb_faces == nb_face_colors))
530  {
531  vtkSmartPointer< vtkDoubleArray > ptr_colors =
532  vtkSmartPointer< vtkDoubleArray >::New();
533  ptr_colors->SetName("face_colors"); // set the point color name
534  ptr_colors->SetNumberOfComponents(
535  static_cast< int >(face_color_coords[0].size())); // 3d or 4d
536  ptr_colors->SetNumberOfTuples(nb_face_colors);
537 
538  for(size_t i = 0; i < nb_face_colors; ++i)
539  {
540  if(face_color_coords[0].size() == 3)
541  ptr_colors->SetTuple3(i,
542  face_color_coords[i][0],
543  face_color_coords[i][1],
544  face_color_coords[i][2]);
545  else
546  ptr_colors->SetTuple4(i,
547  face_color_coords[i][0],
548  face_color_coords[i][1],
549  face_color_coords[i][2],
550  face_color_coords[i][3]);
551  }
552  if(nb_lines == 0)
553  unstructured_grid->GetCellData()->AddArray(ptr_colors);
554  else
555  unstructured_grid->GetFieldData()->AddArray(ptr_colors);
556  }
557  // std::cout << " after adding face colors..." << std::endl ;
559  for(size_t i = 0; i < nb_lines; ++i)
560  {
561  vtkIdType *pt_indices = new vtkIdType[line_indices[i].size()];
562  typename std::vector< IndexType >::const_iterator it_begin(
563  line_indices[i].begin()),
564  it_end(line_indices[i].end());
565  unsigned long cpt = 0;
566  while(it_begin != it_end)
567  {
568  pt_indices[cpt++] = static_cast< vtkIdType >(*it_begin);
569 
570  ++it_begin;
571  }
572  cells->InsertNextCell(static_cast< int >(line_indices[i].size()),
573  pt_indices);
574  switch(cpt)
575  {
576  case 4:
577  type_array[nb_faces + i] = VTK_TETRA;
578  break;
579  case 6:
580  type_array[nb_faces + i] = VTK_HEXAHEDRON;
581  break;
582  default:
583  throw std::runtime_error(
584  "Writer::load_vtkUnstructuredGrid: a face seems to be not valid.");
585  }
586  delete[] pt_indices;
587  }
589  unstructured_grid->SetCells(type_array, cells);
590  if(type_array != NULL)
591  delete type_array;
592  // std::cout << " after adding cells..." << std::endl ;
594  // add lines colors
595  if((nb_line_colors > 0u) && (lines_color_coords[0].size() > 2u) &&
596  (nb_lines == nb_line_colors))
597  {
598  vtkSmartPointer< vtkDoubleArray > ptr_colors =
599  vtkSmartPointer< vtkDoubleArray >::New();
600  ptr_colors->SetName("cell_colors");
601 
602  ptr_colors->SetNumberOfComponents(
603  static_cast< int >(lines_color_coords[0].size())); // 3d or 4d
604  ptr_colors->SetNumberOfTuples(nb_line_colors);
605 
606  for(size_t i = 0; i < nb_line_colors; ++i)
607  {
608  if(lines_color_coords[0].size() == 3)
609  ptr_colors->SetTuple3(i,
610  lines_color_coords[i][0],
611  lines_color_coords[i][1],
612  lines_color_coords[i][2]);
613  else
614  ptr_colors->SetTuple4(i,
615  lines_color_coords[i][0],
616  lines_color_coords[i][1],
617  lines_color_coords[i][2],
618  lines_color_coords[i][3]);
619  }
620  unstructured_grid->GetCellData()->AddArray(ptr_colors);
621  }
622  // std::cout << " after adding line colors..." << std::endl ;
624  // Set 2D arrays of FIELD data when present
625  // std::cout << " before adding data fields..." << std::endl ;
626  if(nb_arrays_in_field > 0u) // vtk writer will align arrays using 9 columns
627  {
628  int choice = 0;
629  // std::cout << nbArraysInField << " data fields to add..." << std::endl ;
630  for(size_t id_array = 0; id_array < nb_arrays_in_field; id_array++)
631  {
632  // std::cout << " data fields number " << id_array << std::endl ;
633  vtkSmartPointer< vtkDoubleArray > ptr_data =
634  vtkSmartPointer< vtkDoubleArray >::New();
635  if(has_one_name_for_each_field)
636  {
637  if(field_names[id_array].find("POINT_DATA_") != std::string::npos)
638  {
639  ptr_data->SetName((field_names[id_array].substr(11)).c_str());
640  choice = 1;
641  }
642  else if(field_names[id_array].find("CELL_DATA_") != std::string::npos)
643  {
644  ptr_data->SetName((field_names[id_array].substr(10)).c_str());
645  choice = 2;
646  }
647  else if(field_names[id_array].find("ELEMENT_DATA_") != std::string::npos)
648  {
649  ptr_data->SetName((field_names[id_array].substr(13)).c_str());
650  choice = 2;
651  }
652  else
653  {
654  ptr_data->SetName(field_names[id_array].c_str());
655  choice = 0;
656  }
657  }
658  else
659  {
660  ptr_data->SetName(
661  (std::string("cell_law_") + StrUtils::convert(id_array))
662  .c_str()); // name of array id_array
663  choice = 0;
664  }
665  ptr_data->SetNumberOfComponents(
666  ((field_attributes[id_array].size() > 0u)
667  ? static_cast< int >(field_attributes[id_array][0].size())
668  : 0));
669  ptr_data->SetNumberOfTuples(field_attributes[id_array].size());
670 
671  // std::cout << " it has " << ptrData->GetNumberOfTuples() << " tuples and
672  // " << ptrData->GetNumberOfComponents() << " components" << std::endl ;
673  // Set 2D array values:
674  for(vtkIdType i = 0; i < ptr_data->GetNumberOfTuples();
675  ++i) // for each line
676  {
677  assert(i < static_cast< vtkIdType >(field_attributes[id_array].size()));
678  for(vtkIdType j = 0; j < ptr_data->GetNumberOfComponents();
679  ++j) // for each column
680  {
681  if((j >=
682  static_cast< vtkIdType >(field_attributes[id_array][i].size())) ||
683  (ptr_data->GetNumberOfComponents() !=
684  static_cast< vtkIdType >(field_attributes[id_array][i].size())))
685  {
686  std::cout << " Trying to access index " << j << " but only "
687  << field_attributes[id_array][i].size()
688  << " components in the data field." << std::endl;
689  std::cout << " We should find " << ptr_data->GetNumberOfComponents()
690  << " components in the provided data field. "
691  << std::endl;
692  assert(false);
693  }
694  // if( field_attributes[id_array][i][j]<0. ||
695  // field_attributes[id_array][i][j]>1.0 )
696  //{
697  // std::cout << "Warning for GenSim project: attributes must be
698  //within [0;1] interval for the time being. " << std::endl ;
699  //}
700  ptr_data->SetComponent(i, j, field_attributes[id_array][i][j]);
701  }
702  }
704  switch(choice)
705  {
706  case 1:
707  unstructured_grid->GetPointData()->AddArray(ptr_data);
708  break;
709  case 2:
710  unstructured_grid->GetCellData()->AddArray(ptr_data);
711  // unstructured_grid->GetCellData()->SetScalars(ptrData); // only the
712  // last ptrData set as scalars will be kept!
713  // for the time being we only manage data as fields...
714  break;
715  default:
716  unstructured_grid->GetFieldData()->AddArray(ptr_data);
717  }
718  }
719  }
720 }
721 
722 template< typename CoordType,
723  typename CoordNType,
724  /*typename coordT_type,*/ typename CoordCType,
725  typename IndexType >
726 void
727 write_vtk_or_vtp_or_vtu_file(
728  std::string file_path,
729  const std::vector< std::vector< CoordType > > &points_coords,
730  const std::vector< std::vector< CoordNType > > &normals_coords,
731  const std::vector< std::vector< CoordCType > > &vertex_color_coords,
732  const std::vector< std::vector< IndexType > >
733  &line_indices,
734  const std::vector< std::vector< CoordCType > > &lines_color_coords,
736  const std::vector< std::vector< IndexType > >
737  &face_indices,
738  const std::vector< std::vector< CoordCType > > &face_color_coords,
739  const std::vector< std::vector< std::vector< double > > >
740  &field_attributes, // array of 2D arrays
741  const std::vector< std::string > &field_names)
742 {
743  bool is_a_poly_data =
744  ((line_indices.size() == 0) ||
745  ((line_indices.size() > 0) && (face_indices.size() == 0) &&
746  !FileUtils::has_extension(
747  file_path,
748  ".vtu"))); // isAPolyData is true if there is no line indices [thus
749  // only vertices+faces] or if there is line indices
750  // without any faces not considered as volume [thus needs
751  // to ensure extension is not .vtu]
752  vtkSmartPointer< vtkPolyData > poly_data;
753  vtkSmartPointer< vtkUnstructuredGrid > unstructured_grid;
755  std::vector<std::string> names;
756  if (field_names.size() != field_attributes.size())
757  {
758  names.resize(field_attributes.size());
759  // When data field names need to be set, we assume that nD (n>=2)
760  // fields associated with positions are displacement fields while others
761  // are understood as physical laws.
762  for (unsigned long i(0); i<field_attributes.size(); ++i){
763  if( (points_coords.size() == field_attributes[i].size()) &&
764  ((points_coords.size() != face_indices.size()) || (field_attributes[i][0].size()>1)))
765  names[i] = std::string("Shifting_") + convert(i);
766  else
767  names[i] = std::string("cell_law_")+ convert(i);
768  }
769  }
770  else
771  names = field_names;
773  if(is_a_poly_data)
774  {
775  load_vtk_poly_data(points_coords,
776  normals_coords,
777  vertex_color_coords,
778  line_indices,
779  lines_color_coords,
780  face_indices,
781  face_color_coords,
782  field_attributes,
783  names,
784  poly_data);
785  }
786  else
787  {
788  // there is no face and no face color for unstructured_grid meshes.
789  std::vector< std::vector< IndexType > > new_facets;
790  std::vector< std::vector< CoordCType > > new_facets_c;
792  load_vtk_unstructured_grid(points_coords,
793  normals_coords,
794  vertex_color_coords,
795  line_indices,
796  lines_color_coords,
797  new_facets,
798  new_facets_c,
799  field_attributes,
800  names,
801  unstructured_grid);
802  }
804  if(FileUtils::has_extension(file_path, ".vtp"))
805  {
806  if(is_a_poly_data)
807  {
808  vtkSmartPointer< vtkXMLPolyDataWriter > writer =
809  vtkSmartPointer< vtkXMLPolyDataWriter >::New();
810  writer->SetFileName(file_path.c_str());
811 #if VTK_MAJOR_VERSION <= 5
812  writer->SetInput(poly_data);
813  // writer->SetInputConnection(poly_data->GetProducerPort()) ;
814 #else
815  // VTK 6 code
816  writer->SetInputData(poly_data);
817 #endif
818  writer->Update();
819  writer->Write();
820  }
821  else
822  {
823  throw std::runtime_error("Writer::write_vtk_or_vtp_or_vtu_file -> file "
824  "extension .vtp is reserved for PolyData.");
825  }
826  }
827  else
828  {
829  if(FileUtils::has_extension(file_path, ".vtk"))
830  {
831  if(is_a_poly_data)
832  {
833  vtkSmartPointer< vtkPolyDataWriter > writer =
834  vtkSmartPointer< vtkPolyDataWriter >::New();
835  writer->SetFileName(file_path.c_str());
836 
837  if(face_indices.size() > 0)
838  writer->SetHeader(
839  "vtk file: 2D data (surface mesh) generated by MEPP software "
840  "available at https://github.com/MEPP-team/MEPP2");
841  else
842  {
843  if(line_indices.size() == 0)
844  writer->SetHeader(
845  "vtk file: 1D data (point cloud) generated by MEPP software "
846  "available at https://github.com/MEPP-team/MEPP2");
847  else
848  {
849  if(line_indices[0].size() == 2)
850  writer->SetHeader(
851  "vtk file: 1D data (polyline(s)) generated by MEPP software "
852  "available at https://github.com/MEPP-team/MEPP2");
853  else
854  writer->SetHeader(
855  "vtk file: 3D data (volume cells) generated by MEPP software "
856  "available at https://github.com/MEPP-team/MEPP2");
857  }
858  }
859 
860 #if VTK_MAJOR_VERSION <= 5
861  writer->SetInput(poly_data);
862  // writer->SetInputConnection(poly_data->GetProducerPort()) ;
863 #else
864  // VTK 6 code
865  writer->SetInputData(poly_data);
866 #endif
867  writer->Update();
868  writer->Write();
869  }
870  else
871  {
872  vtkSmartPointer< vtkUnstructuredGridWriter > writer =
873  vtkSmartPointer< vtkUnstructuredGridWriter >::New();
874  writer->SetFileName(file_path.c_str());
875  writer->SetHeader("vtk file: 3D data (volume mesh) generated by MEPP "
876  "software available at https://github.com/MEPP-team/MEPP2");
877 
878 #if VTK_MAJOR_VERSION <= 5
879  writer->SetInput(unstructured_grid);
880  // writer->SetInputConnection(unstructured_grid->GetProducerPort()) ;
881 #else
882  // VTK 6 code
883  writer->SetInputData(unstructured_grid);
884 #endif
885  writer->Update();
886  writer->Write();
887  }
888  }
889  else if(FileUtils::has_extension(file_path, ".vtu"))
890  {
891  if(!is_a_poly_data)
892  {
893  vtkSmartPointer< vtkXMLUnstructuredGridWriter > writer =
894  vtkSmartPointer< vtkXMLUnstructuredGridWriter >::New();
895  writer->SetFileName(file_path.c_str());
896 
897 #if VTK_MAJOR_VERSION <= 5
898  writer->SetInput(unstructured_grid);
899  // writer->SetInputConnection(unstructured_grid->GetProducerPort()) ;
900 #else
901  // VTK 6 code
902  writer->SetInputData(unstructured_grid);
903 #endif
904  writer->Update();
905  writer->Write();
906  }
907  else
908  {
909  throw std::runtime_error(
910  "Writer::write_vtk_or_vtp_or_vtu_file -> file extension .vtu is "
911  "reserved for UnstructuredGrid.");
912  }
913  }
914  else
915  {
916  throw std::runtime_error(
917  "Writer::write_vtk_or_vtp_or_vtu_file -> file extension is "
918  "inappropriate (neither .vtk, .vtp nor .vtu)");
919  }
920  }
921 }
922 } // namespace IO
923 } // namespace FEVV
924 
FEVV::StrUtils::convert
void convert(const std::string &str, ConvertType &elem)
Definition: StringUtilities.hpp:81
FEVV::IO::load_vtk_poly_data
void load_vtk_poly_data(const std::vector< std::vector< CoordType > > &points_coords, const std::vector< std::vector< CoordNType > > &normals_coords, const std::vector< std::vector< CoordCType > > &vertex_color_coords, const std::vector< std::vector< IndexType > > &line_indices, const std::vector< std::vector< CoordCType > > &lines_color_coords, const std::vector< std::vector< IndexType > > &face_indices, const std::vector< std::vector< CoordCType > > &face_color_coords, const std::vector< std::vector< std::vector< double > > > &field_attributes, const std::vector< std::string > &field_names, vtkSmartPointer< vtkPolyData > &poly_data)
Definition: VtkFileWriter.h:51
FEVV
Interfaces for plugins These interfaces will be used for different plugins.
Definition: Assert.h:16
FEVV::IO::write_vtk_or_vtp_or_vtu_file
void write_vtk_or_vtp_or_vtu_file(std::string file_path, const std::vector< std::vector< CoordType > > &points_coords, const std::vector< std::vector< CoordNType > > &normals_coords, const std::vector< std::vector< CoordCType > > &vertex_color_coords, const std::vector< std::vector< IndexType > > &line_indices, const std::vector< std::vector< CoordCType > > &lines_color_coords, const std::vector< std::vector< IndexType > > &face_indices, const std::vector< std::vector< CoordCType > > &face_color_coords, const std::vector< std::vector< std::vector< double > > > &field_attributes, const std::vector< std::string > &field_names)
Definition: VtkFileWriter.h:727
FEVV::IO::load_vtk_unstructured_grid
void load_vtk_unstructured_grid(const std::vector< std::vector< CoordType > > &points_coords, const std::vector< std::vector< CoordNType > > &normals_coords, const std::vector< std::vector< CoordCType > > &vertex_color_coords, const std::vector< std::vector< IndexType > > &line_indices, const std::vector< std::vector< CoordCType > > &lines_color_coords, const std::vector< std::vector< IndexType > > &face_indices, const std::vector< std::vector< CoordCType > > &face_color_coords, const std::vector< std::vector< std::vector< double > > > &field_attributes, const std::vector< std::string > &field_names, vtkSmartPointer< vtkUnstructuredGrid > &unstructured_grid)
Definition: VtkFileWriter.h:376
StringUtilities.hpp
FEVV::FileUtils::has_extension
bool has_extension(const std::string &file_name)
Definition: FileUtilities.hpp:58
FileUtilities.hpp
Generated by doxygen 1.8.20