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