46 #include <vtkPolyData.h>
47 #include <vtkPoints.h>
48 #include <vtkXMLPolyDataReader.h>
49 #include <vtkXMLPolyDataWriter.h>
50 #include <vtkCellArray.h>
51 #include <vtkDoubleArray.h>
52 #include <vtkUnsignedIntArray.h>
53 #include <vtkPointData.h>
54 #include <vtkCellData.h>
55 #include <vtkDataArray.h>
56 #include <vtkAbstractArray.h>
57 #include <vtkStringArray.h>
58 #include <vtkSmartPointer.h>
67 vtkSmartPointer<vtkXMLPolyDataReader> reader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
68 reader->SetFileName(filename.c_str());
71 vtkSmartPointer<vtkPolyData> vtkMesh = reader->GetOutput();
74 vtkSmartPointer<vtkUnsignedIntArray> v_indices = vtkUnsignedIntArray::SafeDownCast(vtkMesh->GetPointData()->GetArray(
"Indices", trash));
75 vtkSmartPointer<vtkUnsignedIntArray> c_indices = vtkUnsignedIntArray::SafeDownCast(vtkMesh->GetCellData()->GetArray(
"Indices", trash));
77 const unsigned npts = vtkMesh->GetNumberOfPoints();
80 for (
unsigned i = 0; i < npts; ++i) {
82 vertices_.push_back(
Vertex(vtkMesh->GetPoint(i)[0], vtkMesh->GetPoint(i)[1], vtkMesh->GetPoint(i)[2]));
84 vertices_.push_back(
Vertex(vtkMesh->GetPoint(i)[0], vtkMesh->GetPoint(i)[1], vtkMesh->GetPoint(i)[2], v_indices->GetValue(i)));
89 std::vector<std::string> meshes_name;
92 vtkSmartPointer<vtkStringArray> cell_id = vtkStringArray::SafeDownCast(vtkMesh->GetCellData()->GetAbstractArray(
"Names"));
94 const unsigned ntrgs = vtkMesh->GetNumberOfCells();
95 for (
unsigned i = 0; i < ntrgs; ++i) {
96 if (std::find(meshes_name.begin(), meshes_name.end(), cell_id->GetValue(i)) == meshes_name.end()) {
97 meshes_name.push_back(cell_id->GetValue(i));
102 for (std::vector<std::string>::const_iterator vit = meshes_name.begin(); vit != meshes_name.end(); ++vit) {
107 vtkSmartPointer<vtkIdList> l;
110 for (
unsigned i = 0; i < ntrgs; ++i) {
112 if (cell_id->GetValue(i) == mit->name()) {
113 if (vtkMesh->GetCellType(i) == VTK_TRIANGLE) {
114 l = vtkMesh->GetCell(i)->GetPointIds();
119 c_indices->GetValue(i)));
126 std::cerr <<
"This is not a triangulation" << std::endl;
132 mit->build_mesh_vertices();
135 const unsigned nc = vtkMesh->GetPointData()->GetNumberOfArrays()-1;
136 const unsigned nl = vtkMesh->GetNumberOfPoints()+vtkMesh->GetNumberOfCells();
139 for (
unsigned ic = 0; ic < nc; ++ic) {
140 std::ostringstream sic;
142 for (
unsigned iil = 0; iil < vtkMesh->GetPointData()->GetArray((
"Potentials-"+sic.str()).c_str(), trash)->GetSize(); ++iil) {
143 data(vtkUnsignedIntArray::SafeDownCast(vtkMesh->GetPointData()->GetArray(
"Indices", trash))->GetValue(iil), ic) =
144 vtkDoubleArray::SafeDownCast(vtkMesh->GetPointData()->GetArray((
"Potentials-"+sic.str()).c_str(), trash))->GetValue(iil);
146 for (
unsigned iil = 0; iil < vtkMesh->GetCellData()->GetArray((
"Currents-"+sic.str()).c_str(), trash)->GetSize(); ++iil) {
147 data(vtkUnsignedIntArray::SafeDownCast(vtkMesh->GetCellData()->GetArray(
"Indices", trash))->GetValue(iil), ic) =
148 vtkDoubleArray::SafeDownCast(vtkMesh->GetCellData()->GetArray((
"Currents-"+sic.str()).c_str(), trash))->GetValue(iil);
155 std::cerr <<
"Error: please specify USE_VTK to cmake" << std::endl;
164 vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
165 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
166 vtkSmartPointer<vtkDoubleArray> normals = vtkSmartPointer<vtkDoubleArray>::New();
167 vtkSmartPointer<vtkStringArray> cell_id = vtkSmartPointer<vtkStringArray>::New();
168 vtkSmartPointer<vtkUnsignedIntArray> cell_indices = vtkSmartPointer<vtkUnsignedIntArray>::New();
169 vtkSmartPointer<vtkUnsignedIntArray> point_indices = vtkSmartPointer<vtkUnsignedIntArray>::New();
170 std::vector<vtkSmartPointer<vtkDoubleArray> > potentials(data.
ncol());
171 std::vector<vtkSmartPointer<vtkDoubleArray> > currents(data.
ncol());
173 normals->SetNumberOfComponents(3);
174 normals->SetName(
"Normals");
175 cell_id->SetName(
"Names");
176 cell_indices->SetName(
"Indices");
177 point_indices->SetName(
"Indices");
179 if (data.
nlin() != 0) {
181 bool HAS_OUTERMOST =
false;
183 HAS_OUTERMOST =
true;
187 for (
unsigned j = 0; j < data.
ncol(); ++j) {
188 std::stringstream sdip;
190 potentials[j] = vtkSmartPointer<vtkDoubleArray>::New();
191 currents[j] = vtkSmartPointer<vtkDoubleArray>::New();
192 potentials[j]->SetName((
"Potentials-"+sdip.str()).c_str());
193 currents[j]->SetName((
"Currents-"+sdip.str()).c_str());
196 potentials[j]->InsertNextValue(data(vit->index(), j));
199 for (Meshes::const_iterator mit =
meshes_.begin(); mit !=
meshes_.end(); ++mit) {
200 for (Mesh::const_iterator tit = mit->begin(); tit != mit->end(); ++tit) {
201 currents[j]->InsertNextValue(((mit->outermost() && !HAS_OUTERMOST) ? 0.0 : data(tit->index(), j)));
207 std::map<const Vertex*, unsigned> map;
211 points->InsertNextPoint((*vit)(0), (*vit)(1), (*vit)(2));
212 point_indices->InsertNextValue(vit->index());
216 polydata->SetPoints(points);
218 vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New();
220 for (Meshes::const_iterator mit =
meshes_.begin(); mit !=
meshes_.end(); ++mit) {
221 for ( Mesh::const_iterator tit = mit->begin(); tit != mit->end(); ++tit) {
222 vtkIdType triangle[3] = { map[&(tit->s1())], map[&(tit->s2())], map[&(tit->s3())] };
223 polys->InsertNextCell(3, triangle);
224 cell_id->InsertNextValue(mit->name());
225 cell_indices->InsertNextValue(tit->index());
230 polydata->SetPolys(polys);
232 polydata->GetCellData()->AddArray(cell_id);
234 polydata->GetCellData()->AddArray(cell_indices);
236 polydata->GetPointData()->AddArray(point_indices);
238 if (data.
nlin() != 0) {
239 for (
unsigned j = 0; j < data.
ncol(); ++j) {
240 polydata->GetPointData()->AddArray(potentials[j]);
241 polydata->GetCellData()->AddArray(currents[j]);
245 vtkSmartPointer<vtkXMLPolyDataWriter> writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
246 writer->SetFileName(filename.c_str());
248 #if VTK_MAJOR_VERSION <= 5
249 writer->SetInput(polydata);
251 writer->SetInputData(polydata);
256 std::cerr <<
"Error: please specify USE_VTK to cmake" << std::endl;