Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
VtkToFld.cpp
Go to the documentation of this file.
6 
7 #include <vtkPolyDataReader.h>
8 #include <vtkPolyData.h>
9 #include <vtkPointData.h>
10 #include <vtkPoints.h>
11 #include <vtkCellArray.h>
12 #include <vtkCellDataToPointData.h>
13 #include <vtkContourFilter.h>
14 
15 #include <boost/unordered_set.hpp>
16 #include <boost/program_options.hpp>
17 
18 namespace po = boost::program_options;
19 
20 /**
21  * @brief Represents a vertex in the mesh.
22  *
23  * Each vertex has a 3-component coordinate and a scalar value associated with
24  * it. The factor provides a mechanism to specify precision of the coordinate
25  * comparison. Although coordinates are provided as floating-point numbers, they
26  * are stored as integers. The integer value is computed as
27  * floor (x * factor)
28  * Therefore a factor of 100 would ensure a precision of 0.01 in the comparison.
29  */
30 struct Vertex
31 {
32  Vertex(double pX, double pY, double pZ, double pScalar, double factor)
33  : x((int)floor(pX*factor)),
34  y((int)floor(pY*factor)),
35  z((int)floor(pZ*factor)),
36  scalar(pScalar) {}
37 
38  int x;
39  int y;
40  int z;
41  double scalar;
42 
43  bool operator==(const Vertex& v)
44  {
45  return (x == v.x && y == v.y && z == v.z);
46  }
47 };
48 typedef boost::shared_ptr<Vertex> VertexSharedPtr;
49 
50 /// Define comparison operator for the vertex struct
51 bool operator==(const VertexSharedPtr& v1, const VertexSharedPtr& v2)
52 {
53  return v1->operator==(*v2);
54 }
55 
56 
57 /**
58  * @brief Hash function for the Vertex struct used for defining sets.
59  */
60 struct VertexHash : std::unary_function<VertexSharedPtr, std::size_t>
61 {
62  std::size_t operator()(VertexSharedPtr const& p) const
63  {
64  std::size_t seed = 0;
65  boost::hash_combine(seed, p -> x);
66  boost::hash_combine(seed, p -> y);
67  boost::hash_combine(seed, p -> z);
68  return seed;
69  }
70 };
71 
72 /// Define a set of Vertices with associated hash function
73 typedef boost::unordered_set<VertexSharedPtr, VertexHash> VertexSet;
74 
75 
76 /**
77  * Main function.
78  *
79  * Usage: VtkToFld session.xml input.vtk output.fld [options]
80  */
81 int main(int argc, char* argv[])
82 {
83  // Set up available options
84  po::options_description desc("Available options");
85  desc.add_options()
86  ("help,h", "Produce this help message.")
87  ("name,n", po::value<string>()->default_value("Intensity"),
88  "Name of field in VTK file to use for intensity.")
89  ("outname,m", po::value<string>()->default_value("intensity"),
90  "Name of field in output FLD file.")
91  ("precision,p", po::value<double>()->default_value(1),
92  "Precision of vertex matching.");
93 
94  po::options_description hidden("Hidden options");
95  hidden.add_options()
96  ("file", po::value<vector<string> >(), "Input filename");
97 
98  po::options_description cmdline_options;
99  cmdline_options.add(desc).add(hidden);
100 
101  po::positional_options_description p;
102  p.add("file", -1);
103 
104  po::variables_map vm;
105 
106  // Parse command-line options
107  try
108  {
109  po::store(po::command_line_parser(argc, argv).
110  options(cmdline_options).positional(p).run(), vm);
111  po::notify(vm);
112  }
113  catch (const std::exception& e)
114  {
115  cerr << e.what() << endl;
116  cerr << desc;
117  return 1;
118  }
119 
120  if ( vm.count("help") || vm.count("file") == 0 ||
121  vm["file"].as<vector<string> >().size() != 3) {
122  cerr << "Usage: VtkToFld session.xml intensity.vtk output.fld [options]"
123  << endl;
124  cerr << desc;
125  return 1;
126  }
127 
128  // Extract command-line argument values
129  std::vector<std::string> vFiles = vm["file"].as<vector<string> >();
130  const string infile = vFiles[1];
131  const string outfile = vFiles[2];
132  const double factor = vm["precision"].as<double>();
133  const string name = vm["name"].as<string>();
134  const string outname = vm["outname"].as<string>();
135 
136  std::vector<std::string> vFilenames;
140 
141  vFilenames.push_back(vFiles[0]);
142  vSession = LibUtilities::SessionReader::CreateInstance(2, argv, vFilenames);
143 
144  try
145  {
146  //----------------------------------------------
147  // Read in mesh from input file
149  AllocateSharedPtr(vSession);
150  //----------------------------------------------
151 
152  //----------------------------------------------
153  // Define Expansion
155  AllocateSharedPtr(vSession,graph2D);
156  //----------------------------------------------
157 
158  //----------------------------------------------
159  // Set up coordinates of mesh
160  int coordim = Exp->GetCoordim(0);
161  int nq = Exp->GetNpoints();
162 
163  Array<OneD, NekDouble> xc0(nq,0.0);
164  Array<OneD, NekDouble> xc1(nq,0.0);
165  Array<OneD, NekDouble> xc2(nq,0.0);
166 
167  switch(coordim)
168  {
169  case 2:
170  Exp->GetCoords(xc0,xc1);
171  break;
172  case 3:
173  Exp->GetCoords(xc0,xc1,xc2);
174  break;
175  default:
176  ASSERTL0(false,"Coordim not valid");
177  break;
178  }
179  //----------------------------------------------
180 
181  vtkPolyDataReader *vtkMeshReader = vtkPolyDataReader::New();
182  vtkMeshReader->SetFileName(infile.c_str());
183  vtkMeshReader->Update();
184 
185  vtkPolyData *vtkMesh = vtkMeshReader->GetOutput();
186  vtkCellDataToPointData* c2p = vtkCellDataToPointData::New();
187 #if VTK_MAJOR_VERSION <= 5
188  c2p->SetInput(vtkMesh);
189 #else
190  c2p->SetInputData(vtkMesh);
191 #endif
192  c2p->PassCellDataOn();
193  c2p->Update();
194  vtkPolyData *vtkDataAtPoints = c2p->GetPolyDataOutput();
195 
196  vtkPoints *vtkPoints = vtkMesh->GetPoints();
197  ASSERTL0(vtkPoints, "ERROR: cannot get points from mesh.");
198 
199  vtkCellArray *vtkPolys = vtkMesh->GetPolys();
200  ASSERTL0(vtkPolys, "ERROR: cannot get polygons from mesh.");
201 
202  vtkPointData *vtkPData = vtkDataAtPoints->GetPointData();
203  ASSERTL0(vtkPolys, "ERROR: cannot get point data from file.");
204 
205  VertexSet points;
206  VertexSet::iterator vIter;
207  double p[3];
208  double val;
209  double x, y, z;
210  int coeff_idx;
211  int i,j,n;
212 
213  if (!vtkDataAtPoints->GetPointData()->HasArray(name.c_str())) {
214  n = vtkDataAtPoints->GetPointData()->GetNumberOfArrays();
215  cerr << "Input file '" << infile
216  << "' does not have a field named '"
217  << name << "'" << endl;
218  cerr << "There are " << n << " arrays in this file." << endl;
219  for (int i = 0; i < n; ++i)
220  {
221  cerr << " "
222  << vtkDataAtPoints->GetPointData()->GetArray(i)->GetName()
223  << endl;
224  }
225  return 1;
226  }
227 
228  // Build up an unordered set of vertices from the VTK file. For each
229  // vertex a hashed value of the coordinates is generated to within a
230  // given tolerance.
231  n = vtkPoints->GetNumberOfPoints();
232  for (i = 0; i < n; ++i)
233  {
234  vtkPoints->GetPoint(i,p);
235  val = vtkPData->GetScalars(name.c_str())->GetTuple1(i);
236  boost::shared_ptr<Vertex> v(new Vertex(p[0],p[1],p[2],val,factor));
237  points.insert(v);
238  }
239 
240  // Now process each vertex of each element in the mesh
242  for (i = 0; i < Exp->GetNumElmts(); ++i)
243  {
244  StdRegions::StdExpansionSharedPtr e = Exp->GetExp(i);
245  for (j = 0; j < e->GetNverts(); ++j)
246  {
247  // Get the index of the coefficient corresponding to this vertex
248  coeff_idx = Exp->GetCoeff_Offset(i) + e->GetVertexMap(j);
249 
250  // Get the coordinates of the vertex
251  vert = e->as<LocalRegions::Expansion2D>()->GetGeom2D()
252  ->GetVertex(j);
253  vert->GetCoords(x,y,z);
254 
255  // Look up the vertex in the VertexSet
256  boost::shared_ptr<Vertex> v(new Vertex(x,y,z,0.0,factor));
257  vIter = points.find(v);
258 
259  // If not found, maybe the tolerance should be reduced?
260  // If found, record the scalar value from the VTK file in the
261  // corresponding coefficient.
262  if (vIter == points.end())
263  {
264  cerr << "Vertex " << i << " not found. Looking for ("
265  << x << ", " << y << ", " << z << ")" << endl;
266  }
267  else
268  {
269  Exp->UpdateCoeffs()[coeff_idx] = (*vIter)->scalar;
270  }
271  }
272  }
273  Exp->SetPhysState(false);
274 
275  //-----------------------------------------------
276  // Write solution to file
277  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
278  = Exp->GetFieldDefinitions();
279  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
280 
281  for(i = 0; i < FieldDef.size(); ++i)
282  {
283  FieldDef[i]->m_fields.push_back(outname);
284  Exp->AppendFieldData(FieldDef[i], FieldData[i]);
285  }
286 
287  LibUtilities::FieldIO vFld(vSession->GetComm());
288  vFld.Write(outfile, FieldDef, FieldData);
289  //-----------------------------------------------
290  }
291  catch (...) {
292  cout << "An error occurred." << endl;
293  }
294 }