44 #include "VtkMeshReader.hpp"
47 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
57 mBoundaryFacesRead(0),
58 mBoundaryFacesSkipped(0),
59 mCableElementsRead(0),
60 mNumElementAttributes(0),
61 mNumFaceAttributes(0),
62 mNumCableElementAttributes(0),
65 mVtkCellType(VTK_TETRA)
67 vtkXMLUnstructuredGridReader* vtk_xml_unstructured_grid_reader;
73 EXCEPTION(
"Could not open VTU file: " + pathBaseName);
78 vtk_xml_unstructured_grid_reader = vtkXMLUnstructuredGridReader::New();
79 vtk_xml_unstructured_grid_reader->SetFileName( pathBaseName.c_str() );
80 vtk_xml_unstructured_grid_reader->Update();
83 vtk_xml_unstructured_grid_reader->Delete();
86 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
90 mNumNodes = mpVtkUnstructuredGrid->GetNumberOfPoints();
91 unsigned num_cells = mpVtkUnstructuredGrid->GetNumberOfCells();
93 if (ELEMENT_DIM == 2u)
96 mVtkCellType = VTK_TRIANGLE;
98 else if (ELEMENT_DIM == 1u)
100 mNodesPerElement = 2;
101 mVtkCellType = VTK_LINE;
105 vtkCellTypes* cell_types = vtkCellTypes::New();
106 mpVtkUnstructuredGrid->GetCellTypes(cell_types);
108 if (cell_types->GetNumberOfTypes() > 1)
110 mNumCableElementAttributes = 1;
111 for (
unsigned cell_id = 0; cell_id < num_cells; ++cell_id)
113 if (mpVtkUnstructuredGrid->GetCellType(cell_id) == mVtkCellType)
116 assert(mNumCableElements == 0);
118 else if (mpVtkUnstructuredGrid->GetCellType(cell_id) == VTK_LINE)
131 mNumElements = num_cells;
134 cell_types->Delete();
138 if (ELEMENT_DIM == 2u)
140 vtkDataSetSurfaceFilter* p_surface = vtkDataSetSurfaceFilter::New();
141 mpVtkFilterEdges = vtkFeatureEdges::New();
142 #if VTK_MAJOR_VERSION >= 6
143 p_surface->SetInputData(mpVtkUnstructuredGrid);
144 mpVtkFilterEdges->SetInputConnection(p_surface->GetOutputPort());
146 p_surface->SetInput(mpVtkUnstructuredGrid);
147 mpVtkFilterEdges->SetInput(p_surface->GetOutput());
149 mpVtkFilterEdges->Update();
150 mNumFaces = mpVtkFilterEdges->GetOutput()->GetNumberOfCells();
153 else if (ELEMENT_DIM == 3u)
155 mpVtkGeometryFilter = vtkGeometryFilter::New();
156 #if VTK_MAJOR_VERSION >= 6
157 mpVtkGeometryFilter->SetInputData(mpVtkUnstructuredGrid);
159 mpVtkGeometryFilter->SetInput(mpVtkUnstructuredGrid);
161 mpVtkGeometryFilter->Update();
163 mNumFaces = mpVtkGeometryFilter->GetOutput()->GetNumberOfCells();
164 if (mNumCableElements > 0)
167 unsigned num_all_cells = mNumFaces;
168 for (
unsigned i=0; i<num_all_cells; i++)
170 if (mpVtkGeometryFilter->GetOutput()->GetCellType(i) == VTK_LINE)
177 else if (ELEMENT_DIM == 1u)
180 vtkSmartPointer<vtkIdList> enclosing_cells = vtkSmartPointer<vtkIdList>::New();
181 assert(mNumNodes == (
unsigned) mpVtkUnstructuredGrid->GetNumberOfPoints());
182 for (
unsigned point_index = 0; point_index < mNumNodes; ++point_index)
184 mpVtkUnstructuredGrid->GetPointCells(point_index, enclosing_cells);
186 if (enclosing_cells->GetNumberOfIds() == 1u)
198 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
200 mIndexFromZero(true),
204 mNumCableElements(0),
208 mBoundaryFacesRead(0),
209 mBoundaryFacesSkipped(0),
210 mCableElementsRead(0),
211 mNumElementAttributes(0),
212 mNumFaceAttributes(0),
213 mNumCableElementAttributes(0),
216 mVtkCellType(VTK_TETRA)
222 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
225 if (ELEMENT_DIM == 3u)
227 mpVtkGeometryFilter->Delete();
229 else if (ELEMENT_DIM == 2u)
231 mpVtkFilterEdges->Delete();
235 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
241 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
244 return mNumCableElements;
247 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
253 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
259 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
265 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
268 return mNumElementAttributes;
271 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
274 return mNumCableElementAttributes;
277 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
280 return mNumFaceAttributes;
283 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
289 mBoundaryFacesRead = 0;
290 mBoundaryFacesSkipped = 0;
291 mCableElementsRead = 0;
294 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
297 if (mNodesRead >= mNumNodes)
299 EXCEPTION(
"Trying to read data for a node that doesn't exist" );
302 std::vector<double> next_node;
304 for (
unsigned i = 0; i < 3; i++)
306 next_node.push_back( mpVtkUnstructuredGrid->GetPoint(mNodesRead)[i] );
313 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
316 if (mElementsRead >= mNumElements)
318 EXCEPTION(
"Trying to read data for an element that doesn't exist" );
321 if (mpVtkUnstructuredGrid->GetCellType(mElementsRead)!= mVtkCellType)
323 EXCEPTION(
"Element is not of expected type (vtkTetra/vtkTriangle)");
328 for (
unsigned i = 0; i < mNodesPerElement; i++)
330 next_element_data.
NodeIndices.push_back(mpVtkUnstructuredGrid->GetCell(mElementsRead)->GetPointId(i));
337 return next_element_data;
340 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
343 if (mCableElementsRead >= mNumCableElements)
345 EXCEPTION(
"Trying to read data for a cable element that doesn't exist" );
348 unsigned next_index = mNumElements + mCableElementsRead;
349 assert(mpVtkUnstructuredGrid->GetCellType(next_index)==VTK_LINE);
353 for (
unsigned i = 0; i < 2; i++)
355 next_element_data.
NodeIndices.push_back(mpVtkUnstructuredGrid->GetCell(next_index)->GetPointId(i));
358 vtkDataArray *p_scalars = mpVtkUnstructuredGrid->GetCellData()->GetArray(
"Cable radius" );
359 next_element_data.AttributeValue = p_scalars->GetTuple(next_index)[0];
361 mCableElementsRead++;
362 return next_element_data;
366 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
369 if (mBoundaryFacesRead >= mNumFaces)
371 EXCEPTION(
"Trying to read data for a boundary element that doesn't exist");
376 if (ELEMENT_DIM == 3u)
378 while (mpVtkGeometryFilter->GetOutput()->GetCellType(mBoundaryFacesRead + mBoundaryFacesSkipped) == VTK_LINE)
380 mBoundaryFacesSkipped++;
382 for (
unsigned i = 0; i < (mNodesPerElement-1); i++)
384 next_face_data.
NodeIndices.push_back(mpVtkGeometryFilter->GetOutput()->GetCell(mBoundaryFacesRead + mBoundaryFacesSkipped)->GetPointId(i));
387 else if (ELEMENT_DIM == 2u)
389 for (
unsigned i = 0; i < (mNodesPerElement-1); i++)
391 next_face_data.
NodeIndices.push_back(mpVtkFilterEdges->GetOutput()->GetCell(mBoundaryFacesRead)->GetPointId(i));
394 else if (ELEMENT_DIM == 1u)
396 vtkSmartPointer<vtkIdList> enclosing_cells = vtkSmartPointer<vtkIdList>::New();
397 assert(mNumNodes == (
unsigned) mpVtkUnstructuredGrid->GetNumberOfPoints());
398 for (
unsigned point_index = mBoundaryFacesSkipped; point_index < mNumNodes; ++point_index)
400 mBoundaryFacesSkipped++;
401 mpVtkUnstructuredGrid->GetPointCells(point_index, enclosing_cells);
403 if (enclosing_cells->GetNumberOfIds() == 1u)
415 mBoundaryFacesRead++;
416 return next_face_data;
420 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
423 vtkCellData *p_cell_data = mpVtkUnstructuredGrid->GetCellData();
425 if (!p_cell_data->HasArray(dataName.c_str()))
427 EXCEPTION(
"No cell data '" + dataName +
"'");
430 vtkDataArray *p_scalars = p_cell_data->GetArray( dataName.c_str() );
431 if (p_scalars->GetNumberOfComponents() != 1)
433 EXCEPTION(
"The cell data '" + dataName +
"' is not scalar data.");
437 for (
unsigned i = 0; i < mNumElements; i++)
439 dataPayload.push_back( p_scalars->GetTuple(i)[0] );
443 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
446 vtkCellData *p_cell_data = mpVtkUnstructuredGrid->GetCellData();
448 if (!p_cell_data->HasArray(dataName.c_str()))
450 EXCEPTION(
"No cell data '" + dataName +
"'");
453 vtkDataArray *p_scalars = p_cell_data->GetArray( dataName.c_str() );
454 if (p_scalars->GetNumberOfComponents() != 3)
456 EXCEPTION(
"The cell data '" + dataName +
"' is not 3-vector data.");
460 for (
unsigned i = 0; i < mNumElements; i++)
462 c_vector <double, SPACE_DIM> data;
463 for (
unsigned j = 0; j < SPACE_DIM; j++)
465 data[j]= p_scalars->GetTuple(i)[j];
467 dataPayload.push_back(data);
472 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
475 vtkPointData *p_point_data = mpVtkUnstructuredGrid->GetPointData();
477 if (!p_point_data->HasArray(dataName.c_str()))
479 EXCEPTION(
"No point data '" + dataName +
"'");
482 vtkDataArray *p_scalars = p_point_data->GetArray( dataName.c_str() );
483 if (p_scalars->GetNumberOfComponents() != 1)
485 EXCEPTION(
"The point data '" + dataName +
"' is not scalar data.");
490 for (
unsigned i = 0; i < mNumNodes; i++)
492 dataPayload.push_back( p_scalars->GetTuple(i)[0] );
496 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
499 vtkPointData *p_point_data = mpVtkUnstructuredGrid->GetPointData();
501 if (!p_point_data->HasArray(dataName.c_str()))
503 EXCEPTION(
"No point data '" + dataName +
"'");
506 vtkDataArray *p_scalars = p_point_data->GetArray( dataName.c_str() );
508 if (p_scalars->GetNumberOfComponents() != 3)
510 EXCEPTION(
"The point data '" + dataName +
"' is not 3-vector data.");
513 for (
unsigned i = 0; i < mNumNodes; i++)
515 c_vector<double, SPACE_DIM> data;
516 for (
unsigned j = 0; j< SPACE_DIM; j++)
518 data[j] = p_scalars->GetTuple(i)[j];
520 dataPayload.push_back( data );
525 template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
528 return mpVtkUnstructuredGrid;
unsigned GetNumCableElementAttributes() const
#define EXCEPTION(message)
unsigned GetNumEdges() const
void GetPointData(std::string dataName, std::vector< double > &dataPayload)
unsigned GetNumFaceAttributes() const
ElementData GetNextFaceData()
unsigned GetNumElementAttributes() const
ElementData GetNextElementData()
std::vector< unsigned > NodeIndices
unsigned GetNumElements() const
ElementData GetNextCableElementData()
unsigned GetNumCableElements() const
void GetCellData(std::string dataName, std::vector< double > &dataPayload)
unsigned GetNumFaces() const
unsigned GetNumNodes() const
vtkUnstructuredGrid * OutputMeshAsVtkUnstructuredGrid()
std::vector< double > GetNextNode()
vtkSmartPointer< vtkUnstructuredGrid > mpVtkUnstructuredGrid
VtkMeshReader(std::string pathBaseName)