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