Imported Upstream version 0.0.34

1 files changed, 340 insertions, 0 deletions

diff --git a/openEMS/tools/vtk_file_writer.cpp b/openEMS/tools/vtk_file_writer.cpp
new file mode 100644
index 0000000..79c40f3
--- /dev/null
+++ b/openEMS/tools/vtk_file_writer.cpp
@@ -0,0 +1,340 @@
+/*
+*	Copyright (C) 2011,2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+using namespace std;
+
+#include "vtk_file_writer.h"
+
+#include <vtkRectilinearGrid.h>
+#include <vtkRectilinearGridWriter.h>
+#include <vtkXMLRectilinearGridWriter.h>
+#include <vtkStructuredGrid.h>
+#include <vtkStructuredGridWriter.h>
+#include <vtkXMLStructuredGridWriter.h>
+#include <vtkZLibDataCompressor.h>
+#include <vtkFloatArray.h>
+#include <vtkDoubleArray.h>
+#include <vtkFieldData.h>
+#include <vtkPointData.h>
+
+#include <sstream>
+#include <iomanip>
+
+
+VTK_File_Writer::VTK_File_Writer(string filename, int meshType)
+{
+	SetFilename(filename);
+	m_MeshType = meshType;
+	m_NativeDump = false;
+	m_Binary = true;
+	m_Compress = true;
+	m_AppendMode = false;
+	m_ActiveTS = false;
+
+	if (m_MeshType==0) //cartesian mesh
+		m_GridData = vtkRectilinearGrid::New();
+	else if (m_MeshType==1) //cylindrical mesh
+		m_GridData = vtkStructuredGrid::New();
+	else
+	{
+		cerr << "VTK_File_Writer::VTK_File_Writer: Error, unknown mesh type: " << m_MeshType << endl;
+		m_GridData=NULL;
+	}
+}
+
+VTK_File_Writer::~VTK_File_Writer()
+{
+	if (m_GridData)
+		m_GridData->Delete();
+	m_GridData = NULL;
+}
+
+void VTK_File_Writer::SetMeshLines(double const* const* lines, unsigned int const* count, double scaling)
+{
+	if (m_MeshType==0) //cartesian mesh
+	{
+		vtkRectilinearGrid* RectGrid = dynamic_cast<vtkRectilinearGrid*>(m_GridData);
+		if (RectGrid==NULL)
+		{
+			cerr << "VTK_File_Writer::SetMeshLines: Error, grid invalid, this should not have happend! " << endl;
+			exit(1);
+		}
+		RectGrid->SetDimensions(count[0],count[1],count[2]);
+		vtkDoubleArray *Coords[3];
+		for (int n=0;n<3;++n)
+		{
+			m_MeshLines[n].clear();
+			m_MeshLines[n].reserve(count[n]);
+			Coords[n] = vtkDoubleArray::New();
+			for (unsigned int i=0; i<count[n]; i++)
+			{
+				Coords[n]->InsertNextValue(lines[n][i]*scaling);
+				m_MeshLines[n].push_back(lines[n][i]*scaling);
+			}
+		}
+		RectGrid->SetXCoordinates(Coords[0]);
+		RectGrid->SetYCoordinates(Coords[1]);
+		RectGrid->SetZCoordinates(Coords[2]);
+		for (int n=0;n<3;++n)
+			Coords[n]->Delete();
+	}
+	else if (m_MeshType==1) //cylindrical mesh
+	{
+		vtkStructuredGrid* StructGrid = dynamic_cast<vtkStructuredGrid*>(m_GridData);
+		if (StructGrid==NULL)
+		{
+			cerr << "VTK_File_Writer::SetMeshLines: Error, grid invalid, this should not have happend! " << endl;
+			exit(1);
+		}
+
+		for (int n=0;n<3;++n)
+		{
+			m_MeshLines[n].clear();
+			m_MeshLines[n].reserve(count[n]);
+			double scale=1;
+			if (n!=1)
+				scale*=scaling;
+			for (unsigned int i=0; i<count[n]; i++)
+				m_MeshLines[n].push_back(lines[n][i]*scale);
+		}
+
+		StructGrid->SetDimensions(count[0],count[1],count[2]);
+		vtkPoints *points = vtkPoints::New();
+		points->SetNumberOfPoints(count[0]*count[1]*count[2]);
+		double r[3];
+		int id=0;
+		for (unsigned int k=0; k<count[2]; ++k)
+			for (unsigned int j=0; j<count[1]; ++j)
+				for (unsigned int i=0; i<count[0]; ++i)
+				{
+					r[0] = lines[0][i] * cos(lines[1][j]) * scaling;
+					r[1] = lines[0][i] * sin(lines[1][j]) * scaling;
+					r[2] = lines[2][k] * scaling;
+					points->SetPoint(id++,r);
+				}
+		StructGrid->SetPoints(points);
+		points->Delete();
+	}
+	else
+	{
+		cerr << "VTK_File_Writer::SetMeshLines: Error, unknown mesh type: " << m_MeshType << endl;
+	}
+}
+
+void VTK_File_Writer::AddScalarField(string fieldname, double const* const* const* field)
+{
+	vtkDoubleArray* array = vtkDoubleArray::New();
+	array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
+	array->SetName(fieldname.c_str());
+	int id=0;
+	for (unsigned int k=0;k<m_MeshLines[2].size();++k)
+	{
+		for (unsigned int j=0;j<m_MeshLines[1].size();++j)
+		{
+			for (unsigned int i=0;i<m_MeshLines[0].size();++i)
+			{
+				array->SetTuple1(id++,field[i][j][k]);
+			}
+		}
+	}
+	m_GridData->GetPointData()->AddArray(array);
+	array->Delete();
+}
+
+void VTK_File_Writer::AddScalarField(string fieldname, float const* const* const* field)
+{
+	vtkFloatArray* array = vtkFloatArray::New();
+	array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
+	array->SetName(fieldname.c_str());
+	int id=0;
+	for (unsigned int k=0;k<m_MeshLines[2].size();++k)
+	{
+		for (unsigned int j=0;j<m_MeshLines[1].size();++j)
+		{
+			for (unsigned int i=0;i<m_MeshLines[0].size();++i)
+			{
+				array->SetTuple1(id++,field[i][j][k]);
+			}
+		}
+	}
+	m_GridData->GetPointData()->AddArray(array);
+	array->Delete();
+}
+
+void VTK_File_Writer::AddVectorField(string fieldname, double const* const* const* const* field)
+{
+	vtkDoubleArray* array = vtkDoubleArray::New();
+	array->SetNumberOfComponents(3);
+	array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
+	array->SetName(fieldname.c_str());
+	int id=0;
+	double out[3];
+	for (unsigned int k=0;k<m_MeshLines[2].size();++k)
+	{
+		for (unsigned int j=0;j<m_MeshLines[1].size();++j)
+		{
+			double cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
+			double sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
+			for (unsigned int i=0;i<m_MeshLines[0].size();++i)
+			{
+				if ((m_MeshType==0) || (m_NativeDump))
+					array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]);
+				else
+				{
+					out[0] = field[0][i][j][k] * cos_a - field[1][i][j][k] * sin_a;
+					out[1] = field[0][i][j][k] * sin_a + field[1][i][j][k] * cos_a;
+					out[2] = field[2][i][j][k];
+					array->SetTuple3(id++,out[0],out[1],out[2]);
+				}
+			}
+		}
+	}
+	m_GridData->GetPointData()->AddArray(array);
+	array->Delete();
+}
+
+void VTK_File_Writer::AddVectorField(string fieldname, float const* const* const* const* field)
+{
+	vtkFloatArray* array = vtkFloatArray::New();
+	array->SetNumberOfComponents(3);
+	array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
+	array->SetName(fieldname.c_str());
+	int id=0;
+	float out[3];
+	for (unsigned int k=0;k<m_MeshLines[2].size();++k)
+	{
+		for (unsigned int j=0;j<m_MeshLines[1].size();++j)
+		{
+			float cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
+			float sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
+			for (unsigned int i=0;i<m_MeshLines[0].size();++i)
+			{
+				if ((m_MeshType==0) || (m_NativeDump))
+					array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]);
+				else
+				{
+					out[0] = field[0][i][j][k] * cos_a - field[1][i][j][k] * sin_a;
+					out[1] = field[0][i][j][k] * sin_a + field[1][i][j][k] * cos_a;
+					out[2] = field[2][i][j][k];
+					array->SetTuple3(id++,out[0],out[1],out[2]);
+				}
+			}
+		}
+	}
+	m_GridData->GetPointData()->AddArray(array);
+	array->Delete();
+}
+
+
+int VTK_File_Writer::GetNumberOfFields() const
+{
+	return m_GridData->GetPointData()->GetNumberOfArrays();
+}
+
+void VTK_File_Writer::ClearAllFields()
+{
+	while (m_GridData->GetPointData()->GetNumberOfArrays()>0)
+	{
+		const char* name = m_GridData->GetPointData()->GetArrayName(0);
+		m_GridData->GetPointData()->RemoveArray(name);
+	}
+}
+
+bool VTK_File_Writer::Write()
+{
+	return WriteXML();
+}
+
+string VTK_File_Writer::GetTimestepFilename(int pad_length) const
+{
+	if (m_ActiveTS==false)
+		return m_filename;
+
+	stringstream ss;
+	ss << m_filename << "_" << std::setw( pad_length ) << std::setfill( '0' ) << m_timestep;
+
+	return ss.str();
+}
+
+
+bool VTK_File_Writer::WriteASCII()
+{
+	vtkDataWriter* writer = NULL;
+	if (m_MeshType==0) //cartesian mesh
+		writer = vtkRectilinearGridWriter::New();
+	else if (m_MeshType==1) //cylindrical mesh
+		writer = vtkStructuredGridWriter::New();
+	else
+	{
+		cerr << "VTK_File_Writer::WriteASCII: Error, unknown mesh type: " << m_MeshType << endl;
+		return false;
+	}
+
+	writer->SetHeader(m_header.c_str());
+#if VTK_MAJOR_VERSION>=6
+	writer->SetInputData(m_GridData);
+#else
+	writer->SetInput(m_GridData);
+#endif
+
+	string filename = GetTimestepFilename() + ".vtk";
+	writer->SetFileName(filename.c_str());
+	if (m_Binary)
+		writer->SetFileTypeToBinary();
+	else
+		writer->SetFileTypeToASCII();
+
+	writer->Write();
+	writer->Delete();
+	return true;
+}
+
+bool VTK_File_Writer::WriteXML()
+{
+	vtkXMLStructuredDataWriter* writer = NULL;
+	if (m_MeshType==0) //cartesian mesh
+		writer = vtkXMLRectilinearGridWriter::New();
+	else if (m_MeshType==1) //cylindrical mesh
+		writer = vtkXMLStructuredGridWriter::New();
+	else
+	{
+		cerr << "VTK_File_Writer::WriteXML: Error, unknown mesh type: " << m_MeshType << endl;
+		return false;
+	}
+
+#if VTK_MAJOR_VERSION>=6
+	writer->SetInputData(m_GridData);
+#else
+	writer->SetInput(m_GridData);
+#endif
+
+	string filename = GetTimestepFilename() + "." + writer->GetDefaultFileExtension();
+	writer->SetFileName(filename.c_str());
+	if (m_Compress)
+		writer->SetCompressor(vtkZLibDataCompressor::New());
+	else
+		writer->SetCompressor(NULL);
+
+	if (m_Binary)
+		writer->SetDataModeToBinary();
+	else
+		writer->SetDataModeToAscii();
+
+	writer->Write();
+	writer->Delete();
+	return true;
+}