59 mBoundaryFacesRead(0),
60 mBoundaryFacesSkipped(0),
61 mCableElementsRead(0),
62 mNumElementAttributes(0),
63 mNumFaceAttributes(0),
64 mNumCableElementAttributes(0),
67 mVtkCellType(VTK_TETRA)
69 vtkXMLUnstructuredGridReader* vtk_xml_unstructured_grid_reader;
75 EXCEPTION(
"Could not open VTU file: " + pathBaseName);
80 vtk_xml_unstructured_grid_reader = vtkXMLUnstructuredGridReader::New();
81 vtk_xml_unstructured_grid_reader->SetFileName( pathBaseName.c_str() );
82 vtk_xml_unstructured_grid_reader->Update();
85 vtk_xml_unstructured_grid_reader->Delete();
92 mNumNodes = mpVtkUnstructuredGrid->GetNumberOfPoints();
93 unsigned num_cells = mpVtkUnstructuredGrid->GetNumberOfCells();
95 if (ELEMENT_DIM == 2u)
98 mVtkCellType = VTK_TRIANGLE;
100 else if (ELEMENT_DIM == 1u)
102 mNodesPerElement = 2;
103 mVtkCellType = VTK_LINE;
106#if (VTK_MAJOR_VERSION >= 9 && VTK_MINOR_VERSION >= 2)
107 const auto num_distinct_cell_types =
static_cast<unsigned>(mpVtkUnstructuredGrid->GetDistinctCellTypesArray()->GetNumberOfTuples());
109 vtkCellTypes* cell_types = vtkCellTypes::New();
110 mpVtkUnstructuredGrid->GetCellTypes(cell_types);
111 const auto num_distinct_cell_types =
static_cast<unsigned>(cell_types->GetNumberOfTypes());
112 cell_types->Delete();
115 if (num_distinct_cell_types > 1u)
117 mNumCableElementAttributes = 1;
118 for (
unsigned cell_id = 0; cell_id < num_cells; ++cell_id)
120 if (mpVtkUnstructuredGrid->GetCellType(cell_id) == mVtkCellType)
123 assert(mNumCableElements == 0);
125 else if (mpVtkUnstructuredGrid->GetCellType(cell_id) == VTK_LINE)
138 mNumElements = num_cells;
142 if (ELEMENT_DIM == 2u)
144 vtkDataSetSurfaceFilter* p_surface = vtkDataSetSurfaceFilter::New();
145 mpVtkFilterEdges = vtkFeatureEdges::New();
146#if VTK_MAJOR_VERSION >= 6
147 p_surface->SetInputData(mpVtkUnstructuredGrid);
148 mpVtkFilterEdges->SetInputConnection(p_surface->GetOutputPort());
150 p_surface->SetInput(mpVtkUnstructuredGrid);
151 mpVtkFilterEdges->SetInput(p_surface->GetOutput());
153 mpVtkFilterEdges->Update();
154 mNumFaces = mpVtkFilterEdges->GetOutput()->GetNumberOfCells();
157 else if (ELEMENT_DIM == 3u)
159 mpVtkGeometryFilter = vtkGeometryFilter::New();
160#if VTK_MAJOR_VERSION >= 6
161 mpVtkGeometryFilter->SetInputData(mpVtkUnstructuredGrid);
163 mpVtkGeometryFilter->SetInput(mpVtkUnstructuredGrid);
165 mpVtkGeometryFilter->Update();
167 mNumFaces = mpVtkGeometryFilter->GetOutput()->GetNumberOfCells();
168 if (mNumCableElements > 0)
171 unsigned num_all_cells = mNumFaces;
172 for (
unsigned i=0; i<num_all_cells; i++)
174 if (mpVtkGeometryFilter->GetOutput()->GetCellType(i) == VTK_LINE)
181 else if (ELEMENT_DIM == 1u)
184 vtkSmartPointer<vtkIdList> enclosing_cells = vtkSmartPointer<vtkIdList>::New();
185 assert(mNumNodes == (
unsigned) mpVtkUnstructuredGrid->GetNumberOfPoints());
186 for (
unsigned point_index = 0; point_index < mNumNodes; ++point_index)
188 mpVtkUnstructuredGrid->GetPointCells(point_index, enclosing_cells);
190 if (enclosing_cells->GetNumberOfIds() == 1u)
347 if (mCableElementsRead >= mNumCableElements)
349 EXCEPTION(
"Trying to read data for a cable element that doesn't exist" );
352 unsigned next_index = mNumElements + mCableElementsRead;
353 assert(mpVtkUnstructuredGrid->GetCellType(next_index)==VTK_LINE);
357 for (
unsigned i = 0; i < 2; i++)
359 next_element_data.
NodeIndices.push_back(mpVtkUnstructuredGrid->GetCell(next_index)->GetPointId(i));
362 vtkDataArray *p_scalars = mpVtkUnstructuredGrid->GetCellData()->GetArray(
"Cable radius" );
363 next_element_data.
AttributeValue = p_scalars->GetTuple(next_index)[0];
365 mCableElementsRead++;
366 return next_element_data;
373 if (mBoundaryFacesRead >= mNumFaces)
375 EXCEPTION(
"Trying to read data for a boundary element that doesn't exist");
380 if (ELEMENT_DIM == 3u)
382 while (mpVtkGeometryFilter->GetOutput()->GetCellType(mBoundaryFacesRead + mBoundaryFacesSkipped) == VTK_LINE)
384 mBoundaryFacesSkipped++;
386 for (
unsigned i = 0; i < (mNodesPerElement-1); i++)
388 next_face_data.
NodeIndices.push_back(mpVtkGeometryFilter->GetOutput()->GetCell(mBoundaryFacesRead + mBoundaryFacesSkipped)->GetPointId(i));
391 else if (ELEMENT_DIM == 2u)
393 for (
unsigned i = 0; i < (mNodesPerElement-1); i++)
395 next_face_data.
NodeIndices.push_back(mpVtkFilterEdges->GetOutput()->GetCell(mBoundaryFacesRead)->GetPointId(i));
398 else if (ELEMENT_DIM == 1u)
400 vtkSmartPointer<vtkIdList> enclosing_cells = vtkSmartPointer<vtkIdList>::New();
401 assert(mNumNodes == (
unsigned) mpVtkUnstructuredGrid->GetNumberOfPoints());
402 for (
unsigned point_index = mBoundaryFacesSkipped; point_index < mNumNodes; ++point_index)
404 mBoundaryFacesSkipped++;
405 mpVtkUnstructuredGrid->GetPointCells(point_index, enclosing_cells);
407 if (enclosing_cells->GetNumberOfIds() == 1u)
419 mBoundaryFacesRead++;
420 return next_face_data;
427 vtkCellData *p_cell_data = mpVtkUnstructuredGrid->GetCellData();
429 if (!p_cell_data->HasArray(dataName.c_str()))
431 EXCEPTION(
"No cell data '" + dataName +
"'");
434 vtkDataArray *p_scalars = p_cell_data->GetArray( dataName.c_str() );
435 if (p_scalars->GetNumberOfComponents() != 1)
437 EXCEPTION(
"The cell data '" + dataName +
"' is not scalar data.");
441 for (
unsigned i = 0; i < mNumElements; i++)
443 dataPayload.push_back( p_scalars->GetTuple(i)[0] );
450 vtkCellData *p_cell_data = mpVtkUnstructuredGrid->GetCellData();
452 if (!p_cell_data->HasArray(dataName.c_str()))
454 EXCEPTION(
"No cell data '" + dataName +
"'");
457 vtkDataArray *p_scalars = p_cell_data->GetArray( dataName.c_str() );
458 if (p_scalars->GetNumberOfComponents() != 3)
460 EXCEPTION(
"The cell data '" + dataName +
"' is not 3-vector data.");
464 for (
unsigned i = 0; i < mNumElements; i++)
466 c_vector <double, SPACE_DIM> data;
467 for (
unsigned j = 0; j < SPACE_DIM; j++)
469 data[j]= p_scalars->GetTuple(i)[j];
471 dataPayload.push_back(data);
479 vtkPointData *p_point_data = mpVtkUnstructuredGrid->GetPointData();
481 if (!p_point_data->HasArray(dataName.c_str()))
483 EXCEPTION(
"No point data '" + dataName +
"'");
486 vtkDataArray *p_scalars = p_point_data->GetArray( dataName.c_str() );
487 if (p_scalars->GetNumberOfComponents() != 1)
489 EXCEPTION(
"The point data '" + dataName +
"' is not scalar data.");
494 for (
unsigned i = 0; i < mNumNodes; i++)
496 dataPayload.push_back( p_scalars->GetTuple(i)[0] );
503 vtkPointData *p_point_data = mpVtkUnstructuredGrid->GetPointData();
505 if (!p_point_data->HasArray(dataName.c_str()))
507 EXCEPTION(
"No point data '" + dataName +
"'");
510 vtkDataArray *p_scalars = p_point_data->GetArray( dataName.c_str() );
512 if (p_scalars->GetNumberOfComponents() != 3)
514 EXCEPTION(
"The point data '" + dataName +
"' is not 3-vector data.");
517 for (
unsigned i = 0; i < mNumNodes; i++)
519 c_vector<double, SPACE_DIM> data;
520 for (
unsigned j = 0; j< SPACE_DIM; j++)
522 data[j] = p_scalars->GetTuple(i)[j];
524 dataPayload.push_back( data );