Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Classes | Typedefs | Functions
VtkToFld.cpp File Reference
#include <vtkPolyDataReader.h>
#include <vtkPolyData.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkCellArray.h>
#include <vtkCellDataToPointData.h>
#include <vtkContourFilter.h>
#include <boost/unordered_set.hpp>
#include <boost/program_options.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <LibUtilities/Communication/Comm.h>
#include <SpatialDomains/MeshGraph2D.h>
#include <MultiRegions/ExpList2D.h>
#include <LocalRegions/Expansion2D.h>
Include dependency graph for VtkToFld.cpp:

Go to the source code of this file.

Classes

struct  Vertex
 Represents a vertex in the mesh. More...
 
struct  VertexHash
 Hash function for the Vertex struct used for defining sets. More...
 

Typedefs

typedef boost::shared_ptr< VertexVertexSharedPtr
 
typedef boost::unordered_set
< VertexSharedPtr, VertexHash
VertexSet
 Define a set of Vertices with associated hash function. More...
 

Functions

bool operator== (const VertexSharedPtr &v1, const VertexSharedPtr &v2)
 Define comparison operator for the vertex struct. More...
 
int main (int argc, char *argv[])
 

Typedef Documentation

typedef boost::unordered_set<VertexSharedPtr, VertexHash> VertexSet

Define a set of Vertices with associated hash function.

Definition at line 76 of file VtkToFld.cpp.

typedef boost::shared_ptr<Vertex> VertexSharedPtr

Definition at line 51 of file VtkToFld.cpp.

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Main function.

Usage: VtkToFld session.xml input.vtk output.fld [options]

Definition at line 84 of file VtkToFld.cpp.

References ASSERTL0, Nektar::StdRegions::StdExpansion::GetCoords(), Nektar::iterator, CellMLToNektar.cellml_metadata::p, and CellMLToNektar.translators::run().

85 {
86  // Set up available options
87  po::options_description desc("Available options");
88  desc.add_options()
89  ("help,h", "Produce this help message.")
90  ("name,n", po::value<string>()->default_value("Intensity"),
91  "Name of field in VTK file to use for intensity.")
92  ("outname,m", po::value<string>()->default_value("intensity"),
93  "Name of field in output FLD file.")
94  ("precision,p", po::value<double>()->default_value(1),
95  "Precision of vertex matching.");
96 
97  po::options_description hidden("Hidden options");
98  hidden.add_options()
99  ("file", po::value<vector<string> >(), "Input filename");
100 
101  po::options_description cmdline_options;
102  cmdline_options.add(desc).add(hidden);
103 
104  po::positional_options_description p;
105  p.add("file", -1);
106 
107  po::variables_map vm;
108 
109  // Parse command-line options
110  try
111  {
112  po::store(po::command_line_parser(argc, argv).
113  options(cmdline_options).positional(p).run(), vm);
114  po::notify(vm);
115  }
116  catch (const std::exception& e)
117  {
118  cerr << e.what() << endl;
119  cerr << desc;
120  return 1;
121  }
122 
123  if ( vm.count("help") || vm.count("file") == 0 ||
124  vm["file"].as<vector<string> >().size() != 3) {
125  cerr << "Usage: VtkToFld session.xml intensity.vtk output.fld [options]"
126  << endl;
127  cerr << desc;
128  return 1;
129  }
130 
131  // Extract command-line argument values
132  std::vector<std::string> vFiles = vm["file"].as<vector<string> >();
133  const string infile = vFiles[1];
134  const string outfile = vFiles[2];
135  const double factor = vm["precision"].as<double>();
136  const string name = vm["name"].as<string>();
137  const string outname = vm["outname"].as<string>();
138 
139  std::vector<std::string> vFilenames;
143 
144  vFilenames.push_back(vFiles[0]);
145  vSession = LibUtilities::SessionReader::CreateInstance(2, argv, vFilenames);
146 
147  try
148  {
149  //----------------------------------------------
150  // Read in mesh from input file
152  AllocateSharedPtr(vSession);
153  //----------------------------------------------
154 
155  //----------------------------------------------
156  // Define Expansion
158  AllocateSharedPtr(vSession,graph2D);
159  //----------------------------------------------
160 
161  //----------------------------------------------
162  // Set up coordinates of mesh
163  int coordim = Exp->GetCoordim(0);
164  int nq = Exp->GetNpoints();
165 
166  Array<OneD, NekDouble> xc0(nq,0.0);
167  Array<OneD, NekDouble> xc1(nq,0.0);
168  Array<OneD, NekDouble> xc2(nq,0.0);
169 
170  switch(coordim)
171  {
172  case 2:
173  Exp->GetCoords(xc0,xc1);
174  break;
175  case 3:
176  Exp->GetCoords(xc0,xc1,xc2);
177  break;
178  default:
179  ASSERTL0(false,"Coordim not valid");
180  break;
181  }
182  //----------------------------------------------
183 
184  vtkPolyDataReader *vtkMeshReader = vtkPolyDataReader::New();
185  vtkMeshReader->SetFileName(infile.c_str());
186  vtkMeshReader->Update();
187 
188  vtkPolyData *vtkMesh = vtkMeshReader->GetOutput();
189  vtkCellDataToPointData* c2p = vtkCellDataToPointData::New();
190 #if VTK_MAJOR_VERSION <= 5
191  c2p->SetInput(vtkMesh);
192 #else
193  c2p->SetInputData(vtkMesh);
194 #endif
195  c2p->PassCellDataOn();
196  c2p->Update();
197  vtkPolyData *vtkDataAtPoints = c2p->GetPolyDataOutput();
198 
199  vtkPoints *vtkPoints = vtkMesh->GetPoints();
200  ASSERTL0(vtkPoints, "ERROR: cannot get points from mesh.");
201 
202  vtkCellArray *vtkPolys = vtkMesh->GetPolys();
203  ASSERTL0(vtkPolys, "ERROR: cannot get polygons from mesh.");
204 
205  vtkPointData *vtkPData = vtkDataAtPoints->GetPointData();
206  ASSERTL0(vtkPolys, "ERROR: cannot get point data from file.");
207 
208  VertexSet points;
209  VertexSet::iterator vIter;
210  double p[3];
211  double val;
212  double x, y, z;
213  int coeff_idx;
214  int i,j,n;
215 
216  if (!vtkDataAtPoints->GetPointData()->HasArray(name.c_str())) {
217  n = vtkDataAtPoints->GetPointData()->GetNumberOfArrays();
218  cerr << "Input file '" << infile
219  << "' does not have a field named '"
220  << name << "'" << endl;
221  cerr << "There are " << n << " arrays in this file." << endl;
222  for (int i = 0; i < n; ++i)
223  {
224  cerr << " "
225  << vtkDataAtPoints->GetPointData()->GetArray(i)->GetName()
226  << endl;
227  }
228  return 1;
229  }
230 
231  // Build up an unordered set of vertices from the VTK file. For each
232  // vertex a hashed value of the coordinates is generated to within a
233  // given tolerance.
234  n = vtkPoints->GetNumberOfPoints();
235  for (i = 0; i < n; ++i)
236  {
237  vtkPoints->GetPoint(i,p);
238  val = vtkPData->GetScalars(name.c_str())->GetTuple1(i);
239  boost::shared_ptr<Vertex> v(new Vertex(p[0],p[1],p[2],val,factor));
240  points.insert(v);
241  }
242 
243  // Now process each vertex of each element in the mesh
245  for (i = 0; i < Exp->GetNumElmts(); ++i)
246  {
247  StdRegions::StdExpansionSharedPtr e = Exp->GetExp(i);
248  for (j = 0; j < e->GetNverts(); ++j)
249  {
250  // Get the index of the coefficient corresponding to this vertex
251  coeff_idx = Exp->GetCoeff_Offset(i) + e->GetVertexMap(j);
252 
253  // Get the coordinates of the vertex
254  vert = e->as<LocalRegions::Expansion2D>()->GetGeom2D()
255  ->GetVertex(j);
256  vert->GetCoords(x,y,z);
257 
258  // Look up the vertex in the VertexSet
259  boost::shared_ptr<Vertex> v(new Vertex(x,y,z,0.0,factor));
260  vIter = points.find(v);
261 
262  // If not found, maybe the tolerance should be reduced?
263  // If found, record the scalar value from the VTK file in the
264  // corresponding coefficient.
265  if (vIter == points.end())
266  {
267  cerr << "Vertex " << i << " not found. Looking for ("
268  << x << ", " << y << ", " << z << ")" << endl;
269  }
270  else
271  {
272  Exp->UpdateCoeffs()[coeff_idx] = (*vIter)->scalar;
273  }
274  }
275  }
276  Exp->SetPhysState(false);
277 
278  //-----------------------------------------------
279  // Write solution to file
280  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
281  = Exp->GetFieldDefinitions();
282  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
283 
284  for(i = 0; i < FieldDef.size(); ++i)
285  {
286  FieldDef[i]->m_fields.push_back(outname);
287  Exp->AppendFieldData(FieldDef[i], FieldData[i]);
288  }
289 
291  LibUtilities::FieldIO::CreateDefault(vSession);
292  vFld->Write(outfile, FieldDef, FieldData);
293  //-----------------------------------------------
294  }
295  catch (...) {
296  cout << "An error occurred." << endl;
297  }
298 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
boost::unordered_set< VertexSharedPtr, VertexHash > VertexSet
Define a set of Vertices with associated hash function.
Definition: VtkToFld.cpp:76
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
void GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2=NullNekDouble1DArray, Array< OneD, NekDouble > &coords_3=NullNekDouble1DArray)
this function returns the physical coordinates of the quadrature points of the expansion ...
Definition: StdExpansion.h:682
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:309
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
Represents a vertex in the mesh.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
bool operator== ( const VertexSharedPtr v1,
const VertexSharedPtr v2 
)

Define comparison operator for the vertex struct.

Definition at line 54 of file VtkToFld.cpp.

55 {
56  return v1->operator==(*v2);
57 }