2 #include <vtkPolyDataReader.h>
4 #if VTK_MAJOR_VERSION >= 9
9 #include <vtkCellArray.h>
10 #include <vtkCellDataToPointData.h>
11 #include <vtkContourFilter.h>
12 #include <vtkPointData.h>
13 #include <vtkPoints.h>
14 #include <vtkPolyData.h>
16 #include <boost/program_options.hpp>
17 namespace po = boost::program_options;
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)
55 return (x == v.
x && y == v.
y && z == v.
z);
63 return v1->operator==(*v2);
69 struct VertexHash : std::unary_function<VertexSharedPtr, std::size_t>
78 typedef std::unordered_set<VertexSharedPtr, VertexHash>
VertexSet;
85 int main(
int argc,
char *argv[])
88 po::options_description desc(
"Available 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.");
101 po::options_description hidden(
"Hidden options");
105 (
"file", po::value<vector<string> >(),
"Input filename");
108 po::options_description cmdline_options;
109 cmdline_options.add(desc).add(hidden);
111 po::positional_options_description
p;
114 po::variables_map vm;
119 po::store(po::command_line_parser(argc, argv)
120 .options(cmdline_options)
126 catch (
const std::exception &e)
128 cerr << e.what() << endl;
133 if (vm.count(
"help") || vm.count(
"file") == 0 ||
134 vm[
"file"].as<vector<string>>().size() != 3)
136 cerr <<
"Usage: VtkToFld session.xml intensity.vtk output.fld [options]"
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>();
150 std::vector<std::string> vFilenames;
155 vFilenames.push_back(vFiles[0]);
156 vSession = LibUtilities::SessionReader::CreateInstance(2, argv, vFilenames);
162 graph2D = SpatialDomains::MeshGraph::Read(vSession);
173 int coordim = Exp->GetCoordim(0);
174 int nq = Exp->GetNpoints();
183 Exp->GetCoords(xc0, xc1);
186 Exp->GetCoords(xc0, xc1, xc2);
189 ASSERTL0(
false,
"Coordim not valid");
194 vtkPolyDataReader *vtkMeshReader = vtkPolyDataReader::New();
195 vtkMeshReader->SetFileName(infile.c_str());
196 vtkMeshReader->Update();
198 vtkPolyData *vtkMesh = vtkMeshReader->GetOutput();
199 vtkCellDataToPointData *c2p = vtkCellDataToPointData::New();
200 #if VTK_MAJOR_VERSION <= 5
201 c2p->SetInput(vtkMesh);
203 c2p->SetInputData(vtkMesh);
205 c2p->PassCellDataOn();
207 vtkPolyData *vtkDataAtPoints = c2p->GetPolyDataOutput();
209 vtkPoints *vtkPoints = vtkMesh->GetPoints();
210 ASSERTL0(vtkPoints,
"ERROR: cannot get points from mesh.");
212 vtkCellArray *vtkPolys = vtkMesh->GetPolys();
213 ASSERTL0(vtkPolys,
"ERROR: cannot get polygons from mesh.");
215 vtkPointData *vtkPData = vtkDataAtPoints->GetPointData();
216 ASSERTL0(vtkPolys,
"ERROR: cannot get point data from file.");
219 VertexSet::iterator vIter;
226 if (!vtkDataAtPoints->GetPointData()->HasArray(
name.c_str()))
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)
235 << vtkDataAtPoints->GetPointData()->GetArray(i)->GetName()
244 n = vtkPoints->GetNumberOfPoints();
245 for (i = 0; i < n; ++i)
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));
256 for (i = 0; i < Exp->GetNumElmts(); ++i)
259 for (j = 0; j < e->GetNverts(); ++j)
262 coeff_idx = Exp->GetCoeff_Offset(i) + e->GetVertexMap(j);
271 std::shared_ptr<Vertex> v(
new Vertex(x, y, z, 0.0, factor));
272 vIter = points.find(v);
277 if (vIter == points.end())
279 cerr <<
"Vertex " << i <<
" not found. Looking for (" << x
280 <<
", " << y <<
", " << z <<
")" << endl;
284 Exp->UpdateCoeffs()[coeff_idx] = (*vIter)->scalar;
288 Exp->SetPhysState(
false);
292 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
293 Exp->GetFieldDefinitions();
294 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
296 for (i = 0; i < FieldDef.size(); ++i)
298 FieldDef[i]->m_fields.push_back(outname);
299 Exp->AppendFieldData(FieldDef[i], FieldData[i]);
303 LibUtilities::FieldIO::CreateDefault(vSession);
304 vFld->Write(outfile, FieldDef, FieldData);
309 cout <<
"An error occurred." << endl;
#define ASSERTL0(condition, msg)
int main(int argc, char *argv[])
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
std::unordered_set< VertexSharedPtr, VertexHash > VertexSet
Define a set of Vertices with associated hash function.
std::shared_ptr< Vertex > VertexSharedPtr
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
std::shared_ptr< FieldIO > FieldIOSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
The above copyright notice and this permission notice shall be included.
void hash_combine(std::size_t &seed)
Hash function for the Vertex struct used for defining sets.
std::size_t operator()(VertexSharedPtr const &p) const
Represents a vertex in the mesh.
Vertex(double pX, double pY, double pZ, double pScalar, double factor)
bool operator==(const Vertex &v)