Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
VtkMeshReader.cpp
1/*
2
3Copyright (C) Fujitsu Laboratories of Europe, 2009
4
5*/
6
7/*
8
9Copyright (c) 2005-2024, University of Oxford.
10All rights reserved.
11
12University of Oxford means the Chancellor, Masters and Scholars of the
13University of Oxford, having an administrative office at Wellington
14Square, Oxford OX1 2JD, UK.
15
16This file is part of Chaste.
17
18Redistribution and use in source and binary forms, with or without
19modification, are permitted provided that the following conditions are met:
20 * Redistributions of source code must retain the above copyright notice,
21 this list of conditions and the following disclaimer.
22 * Redistributions in binary form must reproduce the above copyright notice,
23 this list of conditions and the following disclaimer in the documentation
24 and/or other materials provided with the distribution.
25 * Neither the name of the University of Oxford nor the names of its
26 contributors may be used to endorse or promote products derived from this
27 software without specific prior written permission.
28
29THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
30AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
31IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
32ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
33LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
34CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
35GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
36HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
37LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
38OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
39
40*/
41
42#ifdef CHASTE_VTK
43
44#include <vtkCellTypes.h>
45
46#include "VtkMeshReader.hpp"
47#include "Exception.hpp"
48
49template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
51 mIndexFromZero(true),
52 mNumNodes(0),
53 mNumElements(0),
54 mNumFaces(0),
55 mNumCableElements(0),
56 mNodesRead(0),
57 mElementsRead(0),
58 mFacesRead(0),
59 mBoundaryFacesRead(0),
60 mBoundaryFacesSkipped(0),
61 mCableElementsRead(0),
62 mNumElementAttributes(0),
63 mNumFaceAttributes(0),
64 mNumCableElementAttributes(0),
65 mOrderOfElements(1),
66 mNodesPerElement(4),
67 mVtkCellType(VTK_TETRA)
68{
69 vtkXMLUnstructuredGridReader* vtk_xml_unstructured_grid_reader;
70
71 // Check file exists
72 mVtuFile.open(pathBaseName.c_str());
73 if (!mVtuFile.is_open())
74 {
75 EXCEPTION("Could not open VTU file: " + pathBaseName);
76 }
77 mVtuFile.close();
78
79 // Load the mesh geometry and data from a file
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();
83 mpVtkUnstructuredGrid = vtk_xml_unstructured_grid_reader->GetOutput();
85 vtk_xml_unstructured_grid_reader->Delete();
86}
87
88template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
90{
91
92 mNumNodes = mpVtkUnstructuredGrid->GetNumberOfPoints();
93 unsigned num_cells = mpVtkUnstructuredGrid->GetNumberOfCells();
94
95 if (ELEMENT_DIM == 2u)
96 {
97 mNodesPerElement = 3;
98 mVtkCellType = VTK_TRIANGLE;
99 }
100 else if (ELEMENT_DIM == 1u)
101 {
102 mNodesPerElement = 2;
103 mVtkCellType = VTK_LINE;
104 }
105
106#if (VTK_MAJOR_VERSION >= 9 && VTK_MINOR_VERSION >= 2)
107 const auto num_distinct_cell_types = static_cast<unsigned>(mpVtkUnstructuredGrid->GetDistinctCellTypesArray()->GetNumberOfTuples());
108#else // VTK older than 9.2
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();
113#endif
114
115 if (num_distinct_cell_types > 1u)
116 {
117 mNumCableElementAttributes = 1;
118 for (unsigned cell_id = 0; cell_id < num_cells; ++cell_id)
119 {
120 if (mpVtkUnstructuredGrid->GetCellType(cell_id) == mVtkCellType)
121 {
122 ++mNumElements;
123 assert(mNumCableElements == 0); //We expect all the simplices first and then the cables at the end of the array
124 }
125 else if (mpVtkUnstructuredGrid->GetCellType(cell_id) == VTK_LINE)
126 {
127 ++mNumCableElements;
128 }
129 else
130 {
132 }
133 }
134 }
135 else
136 {
137 //There is only 1 cell type, so all cells are elements
138 mNumElements = num_cells;
139 }
140
141 // Extract the surface faces
142 if (ELEMENT_DIM == 2u)
143 {
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());
149#else
150 p_surface->SetInput(mpVtkUnstructuredGrid);
151 mpVtkFilterEdges->SetInput(p_surface->GetOutput());
152#endif
153 mpVtkFilterEdges->Update();
154 mNumFaces = mpVtkFilterEdges->GetOutput()->GetNumberOfCells();
155 p_surface->Delete();
156 }
157 else if (ELEMENT_DIM == 3u)
158 {
159 mpVtkGeometryFilter = vtkGeometryFilter::New();
160#if VTK_MAJOR_VERSION >= 6
161 mpVtkGeometryFilter->SetInputData(mpVtkUnstructuredGrid);
162#else
163 mpVtkGeometryFilter->SetInput(mpVtkUnstructuredGrid);
164#endif
165 mpVtkGeometryFilter->Update();
166
167 mNumFaces = mpVtkGeometryFilter->GetOutput()->GetNumberOfCells();
168 if (mNumCableElements > 0)
169 {
170 //The boundary face filter includes the cable elements - get rid of them
171 unsigned num_all_cells = mNumFaces;
172 for (unsigned i=0; i<num_all_cells; i++)
173 {
174 if (mpVtkGeometryFilter->GetOutput()->GetCellType(i) == VTK_LINE)
175 {
176 mNumFaces--;
177 }
178 }
179 }
180 }
181 else if (ELEMENT_DIM == 1u)
182 {
183 mNumFaces = 0;
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)
187 {
188 mpVtkUnstructuredGrid->GetPointCells(point_index, enclosing_cells);
189
190 if (enclosing_cells->GetNumberOfIds() == 1u)
191 {
192 mNumFaces++;
193 }
194 }
195 }
196 else
197 {
199 }
200}
201
202template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
203VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::VtkMeshReader(vtkUnstructuredGrid* p_vtkUnstructuredGrid) :
204 mIndexFromZero(true),
205 mNumNodes(0),
206 mNumElements(0),
207 mNumFaces(0),
208 mNumCableElements(0),
209 mNodesRead(0),
210 mElementsRead(0),
211 mFacesRead(0),
212 mBoundaryFacesRead(0),
213 mBoundaryFacesSkipped(0),
214 mCableElementsRead(0),
215 mNumElementAttributes(0),
216 mNumFaceAttributes(0),
217 mNumCableElementAttributes(0),
218 mOrderOfElements(1),
219 mNodesPerElement(4),
220 mVtkCellType(VTK_TETRA)
221{
222 mpVtkUnstructuredGrid = p_vtkUnstructuredGrid;
224}
225
226template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
228{
229 if (ELEMENT_DIM == 3u)
230 {
231 mpVtkGeometryFilter->Delete();
232 }
233 else if (ELEMENT_DIM == 2u)
234 {
235 mpVtkFilterEdges->Delete();
236 }
237}
238
239template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
241{
242 return mNumElements;
243}
244
245template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
247{
248 return mNumCableElements;
249}
250
251template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
253{
254 return mNumNodes;
255}
256
257template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
259{
260 return mNumFaces;
261}
262
263template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
265{
266 return mNumFaces;
267}
268
269template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
271{
272 return mNumElementAttributes;
273}
274
275template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
277{
278 return mNumCableElementAttributes;
279}
280
281template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
283{
284 return mNumFaceAttributes;
285}
286
287template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
289{
290 mNodesRead = 0;
291 mElementsRead = 0;
292 mFacesRead = 0;
293 mBoundaryFacesRead = 0;
294 mBoundaryFacesSkipped = 0;
295 mCableElementsRead = 0;
296}
297
298template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
300{
301 if (mNodesRead >= mNumNodes)
302 {
303 EXCEPTION( "Trying to read data for a node that doesn't exist" );
304 }
305
306 std::vector<double> next_node;
307
308 for (unsigned i = 0; i < 3; i++)
309 {
310 next_node.push_back( mpVtkUnstructuredGrid->GetPoint(mNodesRead)[i] );
311 }
312
313 mNodesRead++;
314 return next_node;
315}
316
317template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
319{
320 if (mElementsRead >= mNumElements)
321 {
322 EXCEPTION( "Trying to read data for an element that doesn't exist" );
323 }
324
325 if (mpVtkUnstructuredGrid->GetCellType(mElementsRead)!= mVtkCellType)
326 {
327 EXCEPTION("Element is not of expected type (vtkTetra/vtkTriangle)");
328 }
329
330 ElementData next_element_data;
331
332 for (unsigned i = 0; i < mNodesPerElement; i++)
333 {
334 next_element_data.NodeIndices.push_back(mpVtkUnstructuredGrid->GetCell(mElementsRead)->GetPointId(i));
335 }
336
337 // \todo implement method to read element data properly (currently returns zero always...)
338 next_element_data.AttributeValue = 0;
339
340 mElementsRead++;
341 return next_element_data;
342}
343
344template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
346{
347 if (mCableElementsRead >= mNumCableElements)
348 {
349 EXCEPTION( "Trying to read data for a cable element that doesn't exist" );
350 }
351
352 unsigned next_index = mNumElements + mCableElementsRead;
353 assert(mpVtkUnstructuredGrid->GetCellType(next_index)==VTK_LINE);
354
355 ElementData next_element_data;
356
357 for (unsigned i = 0; i < 2; i++)
358 {
359 next_element_data.NodeIndices.push_back(mpVtkUnstructuredGrid->GetCell(next_index)->GetPointId(i));
360 }
361
362 vtkDataArray *p_scalars = mpVtkUnstructuredGrid->GetCellData()->GetArray( "Cable radius" );
363 next_element_data.AttributeValue = p_scalars->GetTuple(next_index)[0];
364
365 mCableElementsRead++;
366 return next_element_data;
367}
368
369
370template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
372{
373 if (mBoundaryFacesRead >= mNumFaces)
374 {
375 EXCEPTION( "Trying to read data for a boundary element that doesn't exist");
376 }
377
378 ElementData next_face_data;
379
380 if (ELEMENT_DIM == 3u)
381 {
382 while (mpVtkGeometryFilter->GetOutput()->GetCellType(mBoundaryFacesRead + mBoundaryFacesSkipped) == VTK_LINE)
383 {
384 mBoundaryFacesSkipped++;
385 }
386 for (unsigned i = 0; i < (mNodesPerElement-1); i++)
387 {
388 next_face_data.NodeIndices.push_back(mpVtkGeometryFilter->GetOutput()->GetCell(mBoundaryFacesRead + mBoundaryFacesSkipped)->GetPointId(i));
389 }
390 }
391 else if (ELEMENT_DIM == 2u)
392 {
393 for (unsigned i = 0; i < (mNodesPerElement-1); i++)
394 {
395 next_face_data.NodeIndices.push_back(mpVtkFilterEdges->GetOutput()->GetCell(mBoundaryFacesRead)->GetPointId(i));
396 }
397 }
398 else if (ELEMENT_DIM == 1u)
399 {
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)
403 {
404 mBoundaryFacesSkipped++;
405 mpVtkUnstructuredGrid->GetPointCells(point_index, enclosing_cells);
406
407 if (enclosing_cells->GetNumberOfIds() == 1u)
408 {
409 next_face_data.NodeIndices.push_back(point_index);
410 break;
411 }
412 }
413 }
414 else
415 {
417 }
418
419 mBoundaryFacesRead++;
420 return next_face_data;
421}
422
423
424template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
425void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetCellData(std::string dataName, std::vector<double>& dataPayload)
426{
427 vtkCellData *p_cell_data = mpVtkUnstructuredGrid->GetCellData();
428
429 if (!p_cell_data->HasArray(dataName.c_str()))
430 {
431 EXCEPTION("No cell data '" + dataName + "'");
432 }
433
434 vtkDataArray *p_scalars = p_cell_data->GetArray( dataName.c_str() );
435 if (p_scalars->GetNumberOfComponents() != 1)
436 {
437 EXCEPTION("The cell data '" + dataName + "' is not scalar data.");
438 }
439
440 dataPayload.clear();
441 for (unsigned i = 0; i < mNumElements; i++)
442 {
443 dataPayload.push_back( p_scalars->GetTuple(i)[0] );
444 }
445}
446
447template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
448void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetCellData(std::string dataName, std::vector<c_vector<double,SPACE_DIM> >& dataPayload)
449{
450 vtkCellData *p_cell_data = mpVtkUnstructuredGrid->GetCellData();
451
452 if (!p_cell_data->HasArray(dataName.c_str()))
453 {
454 EXCEPTION("No cell data '" + dataName + "'");
455 }
456
457 vtkDataArray *p_scalars = p_cell_data->GetArray( dataName.c_str() );
458 if (p_scalars->GetNumberOfComponents() != 3)
459 {
460 EXCEPTION("The cell data '" + dataName + "' is not 3-vector data.");
461 }
462
463 dataPayload.clear();
464 for (unsigned i = 0; i < mNumElements; i++)
465 {
466 c_vector <double, SPACE_DIM> data;
467 for (unsigned j = 0; j < SPACE_DIM; j++)
468 {
469 data[j]= p_scalars->GetTuple(i)[j];
470 }
471 dataPayload.push_back(data);
472 }
473}
474
475
476template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
477void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetPointData(std::string dataName, std::vector<double>& dataPayload)
478{
479 vtkPointData *p_point_data = mpVtkUnstructuredGrid->GetPointData();
480
481 if (!p_point_data->HasArray(dataName.c_str()))
482 {
483 EXCEPTION("No point data '" + dataName + "'");
484 }
485
486 vtkDataArray *p_scalars = p_point_data->GetArray( dataName.c_str() );
487 if (p_scalars->GetNumberOfComponents() != 1)
488 {
489 EXCEPTION("The point data '" + dataName + "' is not scalar data.");
490 }
491
492 dataPayload.clear();
493
494 for (unsigned i = 0; i < mNumNodes; i++)
495 {
496 dataPayload.push_back( p_scalars->GetTuple(i)[0] );
497 }
498}
499
500template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
501void VtkMeshReader<ELEMENT_DIM,SPACE_DIM>::GetPointData(std::string dataName, std::vector<c_vector<double,SPACE_DIM> >& dataPayload)
502{
503 vtkPointData *p_point_data = mpVtkUnstructuredGrid->GetPointData();
504
505 if (!p_point_data->HasArray(dataName.c_str()))
506 {
507 EXCEPTION("No point data '" + dataName + "'");
508 }
509
510 vtkDataArray *p_scalars = p_point_data->GetArray( dataName.c_str() );
511
512 if (p_scalars->GetNumberOfComponents() != 3)
513 {
514 EXCEPTION("The point data '" + dataName + "' is not 3-vector data.");
515 }
516 dataPayload.clear();
517 for (unsigned i = 0; i < mNumNodes; i++)
518 {
519 c_vector<double, SPACE_DIM> data;
520 for (unsigned j = 0; j< SPACE_DIM; j++)
521 {
522 data[j] = p_scalars->GetTuple(i)[j];
523 }
524 dataPayload.push_back( data );
525 }
526}
527
528
529template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
531{
532 return mpVtkUnstructuredGrid;
533}
534
535// Explicit instantiation
536template class VtkMeshReader<1,1>;
537template class VtkMeshReader<1,2>;
538template class VtkMeshReader<1,3>;
539template class VtkMeshReader<2,2>;
540template class VtkMeshReader<2,3>;
541template class VtkMeshReader<3,3>;
542#endif // CHASTE_VTK
#define EXCEPTION(message)
#define NEVER_REACHED
ElementData GetNextCableElementData()
unsigned GetNumCableElementAttributes() const
vtkSmartPointer< vtkUnstructuredGrid > mpVtkUnstructuredGrid
void GetPointData(std::string dataName, std::vector< double > &dataPayload)
unsigned GetNumElementAttributes() const
unsigned GetNumCableElements() const
void GetCellData(std::string dataName, std::vector< double > &dataPayload)
VtkMeshReader(std::string pathBaseName)
void CommonConstructor()
unsigned GetNumFaceAttributes() const
std::ifstream mVtuFile
std::vector< double > GetNextNode()
virtual ~VtkMeshReader()
vtkUnstructuredGrid * OutputMeshAsVtkUnstructuredGrid()
unsigned GetNumNodes() const
unsigned GetNumEdges() const
ElementData GetNextFaceData()
unsigned GetNumFaces() const
unsigned GetNumElements() const
ElementData GetNextElementData()
std::vector< unsigned > NodeIndices