Skip to content

Instantly share code, notes, and snippets.

@ChristosT

ChristosT/CMakeLists.txt

Last active May 28, 2022 06:33
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save ChristosT/1c355b7ce4e39523c4e2774d16496fb0 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save ChristosT/1c355b7ce4e39523c4e2774d16496fb0 to your computer and use it in GitHub Desktop.
Download ZIP
Read triangular/tetrahedral mesh saved as vtk unstructured grid in python/ C++
Raw
CMakeLists.txt
cmake_minimum_required(VERSION 3.0)
project(ReadMesh)
find_package(VTK)
if(NOT VTK_FOUND)
message(FATAL "VTK not found")
endif()
include(${VTK_USE_FILE})
add_executable(read_2d_mesh read_2d_mesh.cpp)
target_link_libraries(read_2d_mesh ${VTK_LIBRARIES})
Raw
read_2d_mesh.cpp
#include <vtkUnstructuredGridReader.h>
#include <vtkUnstructuredGrid.h>
#include <vtkCellIterator.h>
#include <vtkPoints.h>
#include <vtkDataArray.h>
#include <vtkPointData.h>
#include <vtkNew.h>
// 2019.12.05 Christos Tsolakis
int main(int argc, char** argv)
{
if(argc == 1)
return 1;
const char* filename = argv[1];
vtkNew<vtkUnstructuredGridReader> reader;
reader->SetFileName(filename);
reader->Update();
vtkUnstructuredGrid* mesh = reader->GetOutput();
std::cout << "Number of elements " << mesh->GetNumberOfCells() << "\n";
std::cout << "Number of points " << mesh->GetNumberOfPoints() << "\n";
vtkIdType npoints = mesh->GetNumberOfPoints();
vtkPoints* points = mesh->GetPoints();
for(vtkIdType id = 0 ; id < npoints ; id++)
{
// both 2D and 3D meshes use 3D points, for 2D meshes z = 0
double xyz[3];
points->GetPoint(id,xyz);
//std::cout << xyz[0] << " " << xyz[1] << " " << xyz[2] << "\n";
}
// cell iterator
vtkCellIterator* iter = mesh->NewCellIterator();
for(iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextCell())
{
vtkIdList* point_ids = iter->GetPointIds();
size_t a,b,c;
a = point_ids->GetId(0);
b = point_ids->GetId(1);
c = point_ids->GetId(2);
//std::cout << a << " " << b << " " << c << "\n";
}
iter->Delete();
// get point weights
vtkDataArray* weights = mesh->GetPointData()->GetArray("weights");
for(vtkIdType id = 0 ; id < npoints ; id++)
{
float w = weights->GetTuple1(id);
//std::cout << w << "\n";
}
return 0;
}
Raw
read_2d_mesh.py
import sys
from vtk import vtkUnstructuredGridReader
# 2019.12.05 Christos Tsolakis
reader = vtkUnstructuredGridReader()
reader.SetFileName(sys.argv[1])
reader.Update()
mesh = reader.GetOutput()
print("Number of elements {}".format(mesh.GetNumberOfCells()))
print("Number of points {}".format(mesh.GetNumberOfPoints()))
# this is the array of weights
print("Number of points arrays {}".format(mesh.GetPointData().GetNumberOfArrays()))
# One way is to iterate the mesh
npoints = mesh.GetNumberOfPoints()
ncells = mesh.GetNumberOfPoints()
for i in range(npoints):
x,y,z = mesh.GetPoint(i) # for 2D meshes z is always zero
# all ids are 0-based
for i in range(ncells):
cell = mesh.GetCell(i)
a = cell.GetPointId(0)
b = cell.GetPointId(1)
c = cell.GetPointId(2)
# d = cell.GetPointId(3) for 3D
# get point weights
weights = mesh.GetPointData().GetArray("weights")
for i in range(npoints):
w = weights.GetValue(i)
#--------------------------------------
# OR use numpy to get the data as arrays
import numpy as np
# a 2D matrix of size npoints x 3
coordinates = np.array(mesh.GetPoints().GetData())
print(coordinates.shape)
print(coordinates[0])
# a flat array of size (3 + 1) * ncells ((4 +1 )*ncells for 3D mesh)
# holding : [ numpoints for element, point_id_0, point_id_1, point_id_2, ...]
connectivity = np.array(mesh.GetCells().GetData()).reshape(-1,4)
print(connectivity.shape)
print(connectivity[0])
weights = np.array(mesh.GetPointData().GetArray("weights"))
print(weights.shape)
print(weights[0])
#--------------------------------------
# OR use wrap object
# source https://blog.kitware.com/improved-vtk-numpy-integration-part-2/
from vtk.numpy_interface import dataset_adapter as dsa
data_object = dsa.WrapDataObject(mesh)
points = data_object.Points
cells = data_object.Cells
weights = data_object.PointData['weights']
print(points.shape)
print(points[0])
# turn flat array to matrix
cells = cells.reshape(-1,4)
print(cells.shape)
print(cells[0])
print(weights.shape)
print(weights[0])
Raw
read_ascii_vtk.cpp
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>
#include <cassert>
// Christos Tsolakis 2019.12.05
// Read a triangular VTK mesh saved as Legacy ASCII unstructured grid
int main(int argc, char** argv )
{
if(argc == 1)
return 1;
std::ifstream file(argv[1]);
size_t npoints, ncells;
// skip first four lines
std::string line,tmp;
for(int i = 0 ; i < 4; i++)
{
std::getline(file,line);
}
// Read points
file >> line ; // POINTS;
file >> npoints;
file >> line; // point type
std::vector<double> points(3 * npoints);
for(size_t i = 0 ; i < npoints; i++)
{
file >> points[3*i];
file >> points[3*i + 1];
file >> points[3*i + 2];
}
std::cout << npoints <<"\n";
// first
std::cout << points[0] << " " << points[1] << " " << points[2] <<"\n";
// last
//std::cout << points[3*(npoints - 1) ] << " " << points[3*(npoints - 1) +1 ] << " " << points[3*(npoints -1) + 2] <<"\n";
// Read Connectivity
// read until you find the CELLS line
// this is needed because meshes exported from Paraview have sometimes
// a METADATA section
while(line.find("CELLS") == std::string::npos)
{
std::getline(file,line);
}
std::stringstream ss(line);
std::getline(ss,tmp,' '); // skip CELLS
ss >> ncells;
const int cell_npoints = 3; // 4 for 3D tetrahedral mesh
std::vector<size_t> cells( cell_npoints * ncells);
for(size_t i = 0 ; i < ncells; i++)
{
size_t points_per_cell;
file >> points_per_cell;
assert(points_per_cell == cell_npoints);
file >> cells[cell_npoints*i ];
file >> cells[cell_npoints*i + 1];
file >> cells[cell_npoints*i + 2];
//file >> cells[cell_npoints*i + 3]; for 3D
}
std::cout << ncells <<"\n";
std::cout << cells[0] << " " << cells[1] << " " << cells[2] <<"\n";
//std::cout << cells[3*(ncells - 1) ] << " " << cells[3*(ncells - 1) +1 ] << " " << cells[3*(ncells -1) + 2] <<"\n";
// Read weights
// Read until you find the POINTDATA line
while(line.find("POINT_DATA") == std::string::npos)
{
std::getline(file,line);
}
std::stringstream ss1(line);
size_t nweights;
std::getline(ss1,tmp,' '); // skip POINT_DATA
ss1 >> nweights;
assert(nweights == npoints);
// skip next two lines
std::getline(file,line);
std::getline(file,line);
std::vector<float> weights(npoints);
for(size_t i = 0 ; i < npoints; i++)
{
file >> weights[i];
}
std::cout << weights[0] << "\n";
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment