vtk_io.c

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

					const unsigned int dof_nr = mesh.node(k).dof_number(i,j,0);

#ifdef USE_COMPLEX_NUMBERS
					data->SetValue(k,sys.current_solution(dof_nr).real());
#else
					data->SetValue(k,sys.current_solution(dof_nr));
#endif

				} 
				grid->GetPointData()->AddArray(data);				
			} 
		}
	}
}

/*
 * FIXME now this is known to write nonsense on AMR meshes
 * and it strips the imaginary parts of complex Numbers
 */
void VTKIO::system_vectors_to_vtk(const EquationSystems& es,vtkUnstructuredGrid*& grid){
	if (libMesh::processor_id() == 0){

		std::map<std::string,std::vector<Number> > vecs; 
		for(unsigned int i=0;i<es.n_systems();++i){
			const System& sys = es.get_system(i);
			System::const_vectors_iterator v_end = sys.vectors_end();
			System::const_vectors_iterator it = sys.vectors_begin();
			for(;it!= v_end;++it){ // for all vectors on this system
				std::vector<Number> values; 	
				std::cout<<"it "<<it->first<<std::endl;  

				it->second->localize_to_one(values,0);
				std::cout<<"finish localize"<<std::endl;  
				vecs[it->first] = values;
			}
		}


		std::map<std::string,std::vector<Number> >::iterator it = vecs.begin(); 

		for(;it!=vecs.end();++it){

			vtkFloatArray *data = vtkFloatArray::New(); 

			data->SetName(it->first.c_str());

			libmesh_assert(it->second.size()==es.get_mesh().n_nodes());

			data->SetNumberOfValues(it->second.size());

			for(unsigned int i=0;i<it->second.size();++i){

#ifdef USE_COMPLEX_NUMBERS
				data->SetValue(i,it->second[i].real());
#else
				data->SetValue(i,it->second[i]);
#endif

			}

			grid->GetPointData()->AddArray(data);

		} 

	}

/*
	// write out the vectors added to the systems
	if (libMesh::processor_id() == 0)
		{
			std::cout<<"write system vectors"<<std::endl;  
		const MeshBase& mesh = (MeshBase&)es.get_mesh();
		const unsigned int n_nodes = mesh.n_nodes();
		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);
			std::cout<<"i "<<i<<std::endl;  
			System::const_vectors_iterator v_end = sys.vectors_end();
			System::const_vectors_iterator it = sys.vectors_begin();
			for(;it!= v_end;++it){ // for all vectors on this system
				std::cout<<"it- "<<it->second->size()<<" "<<n_nodes<<" "<<it->first<<std::endl;  
				vtkFloatArray *data = vtkFloatArray::New(); 
				data->SetName(it->first.c_str());
				std::vector<Number> values; 	
				it->second->localize(values);
				data->SetNumberOfValues(n_nodes);


	//         MeshBase::const_node_iterator it = mesh.active_nodes_begin();
	//         const MeshBase::const_node_iterator n_end = mesh.nodes_end();
	//         for(unsigned int count=0;it!=n_end;++it,++count){			
				for(unsigned int j=0;j<n_nodes;++j){
					std::cout<<"j "<<j<<" "<<(*it->second).size()<<std::endl;  
//               const unsigned int dof_nr = mesh.node(j).dof_number(i,0,0);
					data->SetValue(j,values[j]);
	//            it++;
				} 
				grid->GetPointData()->AddArray(data);		
			} 
		} 
	}*/
}

/*
// write out mesh data to the VTK file, this might come in handy to display
// boundary conditions and material data
inline void meshdata_to_vtk(const MeshData& meshdata,
										vtkUnstructuredGrid* grid){
	vtkPointData* pointdata = vtkPointData::New();
	
	const unsigned int n_vn = meshdata.n_val_per_node();
	const unsigned int n_dat = meshdata.n_node_data();

	pointdata->SetNumberOfTuples(n_dat);
}*/

#endif


// ------------------------------------------------------------
// vtkIO class members
//
void VTKIO::read (const std::string& name)
{
  // This is a serial-only process for now;
  // the Mesh should be read on processor 0 and
  // broadcast later
  libmesh_assert(libMesh::processor_id() == 0);

#ifndef HAVE_VTK
  std::cerr << "Cannot read VTK file: " << name
	    << "\nYou must have VTK installed and correctly configured to read VTK meshes."
	    << std::endl;
  libmesh_error();

#else
//  std::cout<<"read "<<name <<std::endl;  
  vtkXMLUnstructuredGridReader *reader = vtkXMLUnstructuredGridReader::New();
  reader->SetFileName( name.c_str() );
  //std::cout<<"force read"<<std::endl;  
  // Force reading
  reader->Update();

  // read in the grid   
//  vtkUnstructuredGrid *grid = reader->GetOutput();
  _vtk_grid = reader->GetOutput();
  _vtk_grid->Update();

  // Get a reference to the mesh
  MeshBase& mesh = MeshInput<MeshBase>::mesh();

  // Clear out any pre-existing data from the Mesh
  mesh.clear();
	
  // always numbered nicely??, so we can loop like this
  // I'm pretty sure it is numbered nicely
  for (unsigned int i=0; i < static_cast<unsigned int>(_vtk_grid->GetNumberOfPoints()); ++i)
    {
      // add to the id map
      // and add the actual point
      double * pnt = _vtk_grid->GetPoint(static_cast<vtkIdType>(i));
      Point xyz(pnt[0],pnt[1],pnt[2]);
      Node* newnode = mesh.add_point(xyz,i);

      // Add node to the nodes vector &
      // tell the MeshData object the foreign node id.
      if (this->_mesh_data != NULL)
	this->_mesh_data->add_foreign_node_id (newnode, i);
    }
	
  for (unsigned int i=0; i < static_cast<unsigned int>(_vtk_grid->GetNumberOfCells()); ++i)
    {
      vtkCell* cell = _vtk_grid->GetCell(i);
      Elem* elem = NULL;  // Initialize to avoid compiler warning
      switch(cell->GetCellType())
	{
	case VTK_TETRA:
	  elem = new Tet4();
	  break;
	case VTK_WEDGE:
	  elem = new Prism6();
	  break;
	case VTK_HEXAHEDRON:
	  elem = new Hex8();
	  break;
	case VTK_PYRAMID:	
	  elem = new Pyramid5();					
	  break;
	case VTK_QUADRATIC_HEXAHEDRON:
  	  elem = new Hex20();
	  break;
	case VTK_QUADRATIC_TETRA:
	  elem = new Tet10();
	  break;
	default:
	  std::cerr << "element type not implemented in vtkinterface " << cell->GetCellType() << std::endl;
	  libmesh_error();
	}

  // get the straightforward numbering from the VTK cells
  for(unsigned int j=0;j<elem->n_nodes();++j){
	  elem->set_node(j) = mesh.node_ptr(cell->GetPointId(j));
  } 
  // then get the connectivity 
  std::vector<unsigned int> conn;
  elem->connectivity(0,VTK,conn);
  // then reshuffle the nodes according to the connectivity, this
  // two-time-assign would evade the definition of the vtk_mapping
  for(unsigned int j=0;j<conn.size();++j){
	  elem->set_node(j) = mesh.node_ptr(conn[j]);
  } 
  elem->set_id(i);
  mesh.add_elem(elem);
  } // end loop over VTK cells
  
#endif // HAVE_VTK
}

/*
 * FIXME this operates on the mesh it "gets" from the ES only, this would
 * prevent passing in a mesh that doesn't belong to the ES
 */
/**
 * This method writes out the equationsystems to a .pvtu file (VTK parallel
 * unstructured grid). 
 */
void VTKIO::write_equation_systems(const std::string& fname, const EquationSystems& es)
{
#ifndef HAVE_VTK

  // Do something with the es to avoid a compiler warning.
  es.n_systems();
  
  std::cerr << "Cannot write VTK file: " << fname
	    << "\nYou must have VTK installed and correctly configured to read VTK meshes."
	    << std::endl;
  libmesh_error();
  
#else
  if (libMesh::processor_id() == 0)
		{

	  // check if the filename extension is pvtu
	  libmesh_assert(fname.substr(fname.rfind("."),fname.size())==".pvtu");
	  /*
		* we only use Unstructured grids
		*/
	  _vtk_grid = vtkUnstructuredGrid::New();
	  vtkXMLPUnstructuredGridWriter* writer= vtkXMLPUnstructuredGridWriter::New();
	  std::cout<<"get points "<<std::endl;  
	  vtkPoints* pnts = nodes_to_vtk((const MeshBase&)es.get_mesh());
	  _vtk_grid->SetPoints(pnts);
	  int * types = new int[es.get_mesh().n_active_elem()];
	  std::cout<<"get cells"<<std::endl;  
	  vtkCellArray* cells = cells_to_vtk((const MeshBase&)es.get_mesh(), types);

	  std::cout<<"set cells"<<std::endl;  
	  _vtk_grid->SetCells(types,cells);
	  
	  // I'd like to write out meshdata, but this requires some coding, in
	  // particular, non_const meshdata iterators
	  //   const MeshData& md = es.get_mesh_data();
	  //   if(es.has_mesh_data())
	  //      meshdata_to_vtk(md,_vtk_grid);
	  //   libmesh_assert (soln.size() ==mesh.n_nodes()*names.size());
	  std::cout<<"write solution"<<std::endl;  
	  solution_to_vtk(es,_vtk_grid);

//#ifdef DEBUG
//     if(true) // add some condition here, although maybe it is more sensible to give each vector a flag on whether it is to be written out or not
//       system_vectors_to_vtk(es,_vtk_grid);
//#endif
	  writer->SetInput(_vtk_grid);
	  writer->SetFileName(fname.c_str());
	  writer->SetDataModeToAscii();
	  writer->Write();

	  _vtk_grid->Delete();
	  writer->Delete();
	}
#endif
}

/**
 * This method implements writing to a .vtu (VTK Unstructured Grid) file. 
 * This is one of the new style XML dataformats. 
 */
void VTKIO::write (const std::string& name)
{	
#ifndef HAVE_VTK
  std::cerr << "Cannot write VTK file: " << name
	    << "\nYou must have VTK installed and correctly configured to write VTK meshes."
	    << std::endl;
  libmesh_error();

#else
  if (libMesh::processor_id() == 0)
  {

	  MeshBase& mesh = MeshInput<MeshBase>::mesh();
	  _vtk_grid = vtkUnstructuredGrid::New();
	  vtkXMLUnstructuredGridWriter* writer = vtkXMLUnstructuredGridWriter::New();
	  std::cout<<"write nodes "<<std::endl;  
	  vtkPoints* pnts = nodes_to_vtk(mesh);
	  _vtk_grid->SetPoints(pnts);

	  std::cout<<"write elements "<<std::endl;  
	  int * types = new int[mesh.n_active_elem()];
	  vtkCellArray* cells = cells_to_vtk(mesh,types);
	  _vtk_grid->SetCells(types,cells);
	  //  , _vtk_grid);
	  writer->SetInput(_vtk_grid);
	  writer->SetDataModeToAscii();
	  writer->SetFileName(name.c_str());
	  writer->Write();
  }
#endif // HAVE_VTK
}

//  vim: sw=3 ts=3