1 #include <vtkPolyDataReader.h>
2 #include <vtkPolyData.h>
3 #include <vtkPointData.h>
5 #include <vtkCellArray.h>
6 #include <vtkCellDataToPointData.h>
7 #include <vtkContourFilter.h>
9 #include <boost/unordered_set.hpp>
10 #include <boost/program_options.hpp>
11 namespace po = boost::program_options;
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)),
48 return (x == v.
x && y == v.
y && z == v.
z);
56 return v1->operator==(*v2);
63 struct VertexHash : std::unary_function<VertexSharedPtr, std::size_t>
68 boost::hash_combine(seed, p -> x);
69 boost::hash_combine(seed, p -> y);
70 boost::hash_combine(seed, p -> z);
76 typedef boost::unordered_set<VertexSharedPtr, VertexHash>
VertexSet;
84 int main(
int argc,
char* argv[])
87 po::options_description desc(
"Available 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.");
97 po::options_description hidden(
"Hidden options");
99 (
"file", po::value<vector<string> >(),
"Input filename");
101 po::options_description cmdline_options;
102 cmdline_options.add(desc).add(hidden);
104 po::positional_options_description
p;
107 po::variables_map vm;
112 po::store(po::command_line_parser(argc, argv).
113 options(cmdline_options).positional(p).
run(), vm);
116 catch (
const std::exception& e)
118 cerr << e.what() << endl;
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]"
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>();
139 std::vector<std::string> vFilenames;
144 vFilenames.push_back(vFiles[0]);
145 vSession = LibUtilities::SessionReader::CreateInstance(2, argv, vFilenames);
163 int coordim = Exp->GetCoordim(0);
164 int nq = Exp->GetNpoints();
173 Exp->GetCoords(xc0,xc1);
176 Exp->GetCoords(xc0,xc1,xc2);
179 ASSERTL0(
false,
"Coordim not valid");
184 vtkPolyDataReader *vtkMeshReader = vtkPolyDataReader::New();
185 vtkMeshReader->SetFileName(infile.c_str());
186 vtkMeshReader->Update();
188 vtkPolyData *vtkMesh = vtkMeshReader->GetOutput();
189 vtkCellDataToPointData* c2p = vtkCellDataToPointData::New();
190 #if VTK_MAJOR_VERSION <= 5
191 c2p->SetInput(vtkMesh);
193 c2p->SetInputData(vtkMesh);
195 c2p->PassCellDataOn();
197 vtkPolyData *vtkDataAtPoints = c2p->GetPolyDataOutput();
199 vtkPoints *vtkPoints = vtkMesh->GetPoints();
200 ASSERTL0(vtkPoints,
"ERROR: cannot get points from mesh.");
202 vtkCellArray *vtkPolys = vtkMesh->GetPolys();
203 ASSERTL0(vtkPolys,
"ERROR: cannot get polygons from mesh.");
205 vtkPointData *vtkPData = vtkDataAtPoints->GetPointData();
206 ASSERTL0(vtkPolys,
"ERROR: cannot get point data from file.");
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)
225 << vtkDataAtPoints->GetPointData()->GetArray(i)->GetName()
234 n = vtkPoints->GetNumberOfPoints();
235 for (i = 0; i < n; ++i)
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));
245 for (i = 0; i < Exp->GetNumElmts(); ++i)
248 for (j = 0; j < e->GetNverts(); ++j)
251 coeff_idx = Exp->GetCoeff_Offset(i) + e->GetVertexMap(j);
259 boost::shared_ptr<Vertex> v(
new Vertex(x,y,z,0.0,factor));
260 vIter = points.find(v);
265 if (vIter == points.end())
267 cerr <<
"Vertex " << i <<
" not found. Looking for ("
268 << x <<
", " << y <<
", " << z <<
")" << endl;
272 Exp->UpdateCoeffs()[coeff_idx] = (*vIter)->scalar;
276 Exp->SetPhysState(
false);
280 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
281 = Exp->GetFieldDefinitions();
282 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
284 for(i = 0; i < FieldDef.size(); ++i)
286 FieldDef[i]->m_fields.push_back(outname);
287 Exp->AppendFieldData(FieldDef[i], FieldData[i]);
291 LibUtilities::FieldIO::CreateDefault(vSession);
292 vFld->Write(outfile, FieldDef, FieldData);
296 cout <<
"An error occurred." << endl;
std::size_t operator()(VertexSharedPtr const &p) const
#define ASSERTL0(condition, msg)
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.
int main(int argc, char *argv[])
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
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 ...
boost::shared_ptr< Vertex > VertexSharedPtr
Vertex(double pX, double pY, double pZ, double pScalar, double factor)
boost::shared_ptr< FieldIO > FieldIOSharedPtr
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Represents a vertex in the mesh.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
bool operator==(const Vertex &v)
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Hash function for the Vertex struct used for defining sets.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
boost::shared_ptr< PointGeom > PointGeomSharedPtr