vtk_io.c

一个用来实现偏微分方程中网格的计算库

// $Id: vtk_io.C 2971 2008-08-11 04:53:37Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library 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
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


// C++ includes
#include <fstream>

// Local includes
#include "vtk_io.h"
#include "mesh_base.h"
#include "mesh.h"
#include "equation_systems.h"
#include "cell_tet4.h"
#include "cell_tet10.h"
#include "cell_prism6.h"
#include "cell_pyramid5.h"
#include "cell_hex8.h"
#include "cell_hex20.h"
#include "numeric_vector.h"
#include "system.h"
//#include "cell_inf.h"
//#include "cell_inf_prism12.h"
//#include "cell_prism6.h"
//#include "cell_tet4.h"

//#include "cell_inf_hex16.h"
//#include "cell_inf_prism6.h"
//#include "cell_prism.h"
//#include "cell_tet.h"

//#include "cell_hex27.h"
//#include "cell_inf_hex18.h"
//#include "cell_inf_prism.h"
//#include "cell_pyramid5.h"

//#include "cell_hex8.h"
//#include "cell_inf_hex8.h"
//#include "cell_prism15.h"
//#include "cell_pyramid.h"

//#include "cell_hex.h"
//#include "cell_inf_hex.h"
//#include "cell_prism18.h"
//#include "cell_tet10.h"



#include "mesh_data.h"

#ifdef HAVE_VTK

#include "vtkXMLUnstructuredGridReader.h"
#include "vtkXMLUnstructuredGridWriter.h"
#include "vtkXMLPUnstructuredGridWriter.h"
#include "vtkUnstructuredGrid.h"
#include "vtkGenericGeometryFilter.h"
#include "vtkCellArray.h"
#include "vtkConfigure.h"
/*
#include "vtkTetra.h"
#include "vtkTriangle.h"
#include "vtkWedge.h"
#include "vtkPyramid.h"
#include "vtkQuad.h"
#include "vtkHexahedron.h"
#include "vtkQuadraticTriangle.h"
#include "vtkQuadraticQuad.h"
#include "vtkQuadraticHexahedron.h"
#include "vtkQuadraticWedge.h"
#include "vtkQuadraticTetra.h"
#include "vtkBiQuadraticQuad.h"
*/
#include "vtkFloatArray.h"
#include "vtkDoubleArray.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#endif //HAVE_VTK

//static unsigned int vtk_tet4_mapping[4]= {0,1,2,3};
//static unsigned int vtk_pyramid5_mapping[5]= {0,3,2,1,4};
//static unsigned int vtk_wedge6_mapping[6]= {0,2,1,3,5,4};
//static unsigned int vtk_hex8_mapping[8]= {0,1,2,3,4,5,6,7};

#ifdef HAVE_VTK // private functions
vtkPoints* VTKIO::nodes_to_vtk(const MeshBase& mesh){
	if (libMesh::processor_id() == 0)
		{

			vtkPoints* points = vtkPoints::New();
			vtkDoubleArray* pcoords = vtkDoubleArray::New();
			pcoords->SetNumberOfComponents(3);
			pcoords->SetNumberOfTuples(mesh.n_nodes());

			// FIXME change to local iterators when I've figured out how to write in parallel
			//  MeshBase::const_node_iterator nd = mesh.local_nodes_begin();
			//  MeshBase::const_node_iterator nd_end = mesh.local_nodes_end();
			//  points->SetNumberOfPoints(mesh.n_local_nodes()); // it seems that it needs this to prevent a segfault
			MeshBase::const_node_iterator nd = mesh.nodes_begin();
			MeshBase::const_node_iterator nd_end = mesh.nodes_end();
			// write out nodal data
			for (;nd!=nd_end;++nd)
			{
				Node* node = (*nd);
				float* pnt = new float[3];
				for(unsigned int i=0;i<3;++i){
					pnt[i] = (*node)(i);
				} 
				//     pcoords->InsertNextTuple(pnt);
				pcoords->SetTuple(node->id(),pnt);
				//     std::cout<<"pnt "<<pnt[0]<<" "<<pnt[1]<<" "<<pcoords-><<std::endl;  
				//     points->InsertPoint(node->id(),pnt);
				//     std::cout<<"point "<<node->id()<<" "<<(*node)(0)<<" "<<(*node)(1)<<std::endl;  
			}
			points->SetData(pcoords);
			return points;
	}
	return NULL;
}

vtkCellArray* VTKIO::cells_to_vtk(const MeshBase& mesh, int*& types){
	//, vtkUnstructuredGrid*& grid){
	//  MeshBase::const_element_iterator       it  = mesh.active_local_elements_begin();
	//  const MeshBase::const_element_iterator end = mesh.active_local_elements_end(); 
	if (libMesh::processor_id() == 0)
	{

		vtkCellArray* cells = vtkCellArray::New();
		//     cells->Allocate(100000);
		vtkIdList *pts = vtkIdList::New();

		// for some reason the iterator variant of this code throws an error
		//     MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
		//     MeshBase::const_element_iterator end = mesh.active_elements_end(); 

		//     for ( ; it != end; ++it)
		//     {
		for(unsigned int el_nr =0; el_nr< mesh.n_active_elem(); ++el_nr){
			Elem *elem  = mesh.elem(el_nr);
			//         if(elem->active()){ 
			//       std::cout<<"sub elem "<<elem->has_children()<<std::endl;  
			//       (*it);
			vtkIdType celltype = VTK_EMPTY_CELL; // initialize to something to avoid compiler warning

			switch(elem->type())
			{
				case EDGE2:				
					celltype = VTK_LINE;
					break;
				case EDGE3:      
					celltype = VTK_QUADRATIC_EDGE;
					break;// 1
				case TRI3:       
					celltype = VTK_TRIANGLE;
					break;// 3
				case TRI6:       
					celltype = VTK_QUADRATIC_TRIANGLE;
					break;// 4
				case QUAD4:      
					celltype = VTK_QUAD;
					break;// 5
				case QUAD8:      
					celltype = VTK_QUADRATIC_QUAD;
					break;// 6
				case TET4:      
					celltype = VTK_TETRA;
					break;// 8
				case TET10:      
					celltype = VTK_QUADRATIC_TETRA;
					break;// 9
				case HEX8:    
					celltype = VTK_HEXAHEDRON;
					break;// 10
				case HEX20:      
					celltype = VTK_QUADRATIC_HEXAHEDRON;
					break;// 12
				case PRISM6:     
					celltype = VTK_WEDGE;
					break;// 13
				case PRISM15:   
					celltype = VTK_HIGHER_ORDER_WEDGE;
					break;// 14
				case PRISM18:    
					break;// 15
				case PYRAMID5:
					celltype = VTK_PYRAMID;
					break;// 16
#if VTK_MAJOR_VERSION > 5 || (VTK_MAJOR_VERSION == 5 && VTK_MINOR_VERSION > 0)
				case QUAD9:      
					celltype = VTK_BIQUADRATIC_QUAD;
					break;
#else
				case QUAD9:
#endif 
				case EDGE4:      
				case HEX27:      
				case INFEDGE2:   
				case INFQUAD4:   
				case INFQUAD6:   
				case INFHEX8:    
				case INFHEX16:   
				case INFHEX18:   
				case INFPRISM6:  
				case INFPRISM12: 
				case NODEELEM:   
				case INVALID_ELEM:
				default:
					{
						std::cerr<<"element type "<<elem->type()<<" not implemented"<<std::endl;
						libmesh_error();
					}
			}
			//        std::cout<<"elem "<<elem->n_nodes()<<" "<<elem->n_sub_elem()<<std::endl;  
			//        std::cout<<"type "<<celltype<<" "<<VTK_BIQUADRATIC_QUAD<<std::endl;  
			types[elem->id()] = celltype;
			pts->SetNumberOfIds(elem->n_nodes());
			// get the connectivity for this element
			std::vector<unsigned int> conn;
			elem->connectivity(0,VTK,conn);
			//        std::cout<<"conn size "<<conn.size()<<std::endl;  
			for(unsigned int i=0;i<conn.size();++i)
			{
				//           std::cout<<" sz "<<pts->GetNumberOfIds()<<std::endl;  
				pts->InsertId(i,conn[i]);
			} 
			//        std::cout<<"set cell"<<std::endl;  
			cells->InsertNextCell(pts);
			//         }
			//         std::cout<<"finish set cell"<<std::endl;  
		} // end loop over active elements
		pts->Delete();

		return cells;
	}
	return NULL;
}

/*
 * single processor implementation, this can be done in parallel, but requires a
 * bit more thinking on how to deal with overlapping local solution vectors
 */
/*
void VTKIO::solution_to_vtk(const EquationSystems& es,vtkUnstructuredGrid*& grid){
//   if (libMesh::processor_id() == 0)
//      {
//   std::vector<Number> soln;
//   std::cout<<"build"<<std::endl;  
//   es.build_solution_vector(soln);
//   std::cout<<"add solution "<<soln.size()<<std::endl;  

	if(libMesh::processor_id()==0){
		const MeshBase& mesh = (MeshBase&)es.get_mesh();
		const unsigned int n_nodes = es.get_mesh().n_nodes();
		const unsigned int n_vars = es.n_vars();
			// write the solutions belonging to the system
			for(unsigned int i=0;i<es.n_systems();++i){ // for all systems, regardless of whether they are active or not
				for(unsigned int j=0;j<n_vars;++j){
					const System& sys = es.get_system(i);
					const std::string& varname = sys.variable_name(j); 
					vtkFloatArray *data = vtkFloatArray::New(); 
					data->SetName(varname.c_str());
					data->SetNumberOfValues(n_nodes);
					MeshBase::const_node_iterator it = mesh.nodes_begin();
					const MeshBase::const_node_iterator n_end = mesh.nodes_end();
					for(;it!=n_end;++it){
						if((*it)->n_comp(i,j)>0){
							const unsigned int nd_id = (*it)->id();
	//                  const unsigned int dof_nr = (*it)->dof_number(i,j,0);
	//                  data->InsertValue((*it)->id(),sys.current_solution(count++));
							data->InsertValue(nd_id,libmesh_real(13.0));
//                     soln[nd_id*n_vars+j]));
	//               data->InsertValue((*it)->id(),sys.current_solution(dof_nr));
						}else{
							data->InsertValue((*it)->id(),0);
						}
					} 
					grid->GetPointData()->AddArray(data);
				} 
			} 
   }
}*/

void VTKIO::solution_to_vtk(const EquationSystems& es, vtkUnstructuredGrid*& grid){
	// write only single processor
	if(libMesh::processor_id()==0){

		const MeshBase& mesh = (MeshBase&)es.get_mesh();

		const unsigned int n_nodes = es.get_mesh().n_nodes();
		// loop over the systems
		for(unsigned int i=0;i<es.n_systems();++i){ // for all systems, regardless of whether they are active or not
			
			const System& sys = es.get_system(i);

			const unsigned int n_vars = sys.n_vars();
			// loop over variables
			for(unsigned int j=0;j<n_vars;++j){

				std::string name = sys.variable_name(j);

				vtkFloatArray *data = vtkFloatArray::New(); 

				data->SetName(name.c_str());

				data->SetNumberOfValues(sys.solution->size());

				for(unsigned int k=0;k<n_nodes;++k){