Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
VtkToFld.cpp
Go to the documentation of this file.
1 #include <vtkPolyDataReader.h>
2 #include <vtkPolyData.h>
3 #include <vtkPointData.h>
4 #include <vtkPoints.h>
5 #include <vtkCellArray.h>
6 #include <vtkCellDataToPointData.h>
7 #include <vtkContourFilter.h>
8 
9 #include <boost/unordered_set.hpp>
10 #include <boost/program_options.hpp>
11 namespace po = boost::program_options;
12 
17 #include <MultiRegions/ExpList2D.h>
19 
20 using namespace std;
21 using namespace Nektar;
22 
23 /**
24  * @brief Represents a vertex in the mesh.
25  *
26  * Each vertex has a 3-component coordinate and a scalar value associated with
27  * it. The factor provides a mechanism to specify precision of the coordinate
28  * comparison. Although coordinates are provided as floating-point numbers, they
29  * are stored as integers. The integer value is computed as
30  * floor (x * factor)
31  * Therefore a factor of 100 would ensure a precision of 0.01 in the comparison.
32  */
33 struct Vertex
34 {
35  Vertex(double pX, double pY, double pZ, double pScalar, double factor)
36  : x((int)floor(pX*factor)),
37  y((int)floor(pY*factor)),
38  z((int)floor(pZ*factor)),
39  scalar(pScalar) {}
40 
41  int x;
42  int y;
43  int z;
44  double scalar;
45 
46  bool operator==(const Vertex& v)
47  {
48  return (x == v.x && y == v.y && z == v.z);
49  }
50 };
51 typedef boost::shared_ptr<Vertex> VertexSharedPtr;
52 
53 /// Define comparison operator for the vertex struct
54 bool operator==(const VertexSharedPtr& v1, const VertexSharedPtr& v2)
55 {
56  return v1->operator==(*v2);
57 }
58 
59 
60 /**
61  * @brief Hash function for the Vertex struct used for defining sets.
62  */
63 struct VertexHash : std::unary_function<VertexSharedPtr, std::size_t>
64 {
65  std::size_t operator()(VertexSharedPtr const& p) const
66  {
67  std::size_t seed = 0;
68  boost::hash_combine(seed, p -> x);
69  boost::hash_combine(seed, p -> y);
70  boost::hash_combine(seed, p -> z);
71  return seed;
72  }
73 };
74 
75 /// Define a set of Vertices with associated hash function
76 typedef boost::unordered_set<VertexSharedPtr, VertexHash> VertexSet;
77 
78 
79 /**
80  * Main function.
81  *
82  * Usage: VtkToFld session.xml input.vtk output.fld [options]
83  */
84 int main(int argc, char* argv[])
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 }
std::size_t operator()(VertexSharedPtr const &p) const
Definition: VtkToFld.cpp:65
#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
STL namespace.
int main(int argc, char *argv[])
Definition: VtkToFld.cpp:84
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
int z
Definition: VtkToFld.cpp:43
int y
Definition: VtkToFld.cpp:42
boost::shared_ptr< Vertex > VertexSharedPtr
Definition: VtkToFld.cpp:51
Vertex(double pX, double pY, double pZ, double pScalar, double factor)
Definition: VtkToFld.cpp:35
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
double scalar
Definition: VtkToFld.cpp:44
bool operator==(const Vertex &v)
Definition: VtkToFld.cpp:46
int x
Definition: VtkToFld.cpp:41
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Hash function for the Vertex struct used for defining sets.
Definition: VtkToFld.cpp:63
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
Definition: VtkToFld.cpp:54
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60