Nektar++
XmlToVtk.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 //
3 // File: XmlToVtk.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Output VTK file of XML mesh, optionally with Jacobian.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <cstdlib>
37 #include <iomanip>
38 
39 #include <MultiRegions/ExpList.h>
40 #include <MultiRegions/ExpList1D.h>
41 #include <MultiRegions/ExpList2D.h>
43 #include <MultiRegions/ExpList3D.h>
45 
46 using namespace Nektar;
47 
48 int main(int argc, char *argv[])
49 {
50  if(argc < 2)
51  {
52  cerr << "Usage: XmlToVtk meshfile" << endl;
53  exit(1);
54  }
55 
57  "jacobian", "j", "Output Jacobian as scalar field");
59  "quality", "q", "Output distribution of scaled Jacobians");
60 
63 
64  bool jac = vSession->DefinesCmdLineArgument("jacobian");
65  bool quality = vSession->DefinesCmdLineArgument("quality");
66 
67  jac = quality ? true : jac;
68 
69  // Read in mesh from input file
70  string meshfile(argv[argc-1]);
72  SpatialDomains::MeshGraph::Read(vSession); //meshfile);
73 
74  // Set up Expansion information
75  SpatialDomains::ExpansionMap emap = graphShPt->GetExpansions();
77 
78  for (it = emap.begin(); it != emap.end(); ++it)
79  {
80  for (int i = 0; i < it->second->m_basisKeyVector.size(); ++i)
81  {
82  LibUtilities::BasisKey tmp1 = it->second->m_basisKeyVector[i];
84  it->second->m_basisKeyVector[i] = LibUtilities::BasisKey(
85  tmp1.GetBasisType(), tmp1.GetNumModes(),
88  }
89  }
90 
91  // Define Expansion
92  int expdim = graphShPt->GetMeshDimension();
94 
95  switch(expdim)
96  {
97  case 1:
98  {
99  if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
100  {
101  std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
103 
104  ASSERTL0(
105  HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
106  HomoStr == "1D" || HomoStr == "Homo1D",
107  "Only 3DH1D supported for XML output currently.");
108 
109  int nplanes;
110  vSession->LoadParameter("HomModesZ", nplanes);
111 
112  // choose points to be at evenly spaced points at nplanes + 1
113  // points
114  const LibUtilities::PointsKey Pkey(
115  nplanes + 1, LibUtilities::ePolyEvenlySpaced);
116  const LibUtilities::BasisKey Bkey(
117  LibUtilities::eFourier, nplanes, Pkey);
118  NekDouble lz = vSession->GetParameter("LZ");
119 
122  vSession, Bkey, lz, false, false, graphShPt);
123  Exp[0] = Exp2DH1;
124  }
125  else
126  {
129  ::AllocateSharedPtr(vSession,graphShPt);
130  Exp[0] = Exp1D;
131  }
132 
133  break;
134  }
135  case 2:
136  {
137  if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
138  {
139  std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
141 
142  ASSERTL0(
143  HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
144  HomoStr == "1D" || HomoStr == "Homo1D",
145  "Only 3DH1D supported for XML output currently.");
146 
147  int nplanes;
148  vSession->LoadParameter("HomModesZ", nplanes);
149 
150  // choose points to be at evenly spaced points at nplanes + 1
151  // points
152  const LibUtilities::PointsKey Pkey(
153  nplanes + 1, LibUtilities::ePolyEvenlySpaced);
154  const LibUtilities::BasisKey Bkey(
155  LibUtilities::eFourier, nplanes, Pkey);
156  NekDouble lz = vSession->GetParameter("LZ");
157 
160  vSession, Bkey, lz, false, false, graphShPt);
161  Exp[0] = Exp3DH1;
162  }
163  else
164  {
167  ::AllocateSharedPtr(vSession,graphShPt);
168  Exp[0] = Exp2D;
169  }
170  break;
171  }
172  case 3:
173  {
176  ::AllocateSharedPtr(vSession,graphShPt);
177  Exp[0] = Exp3D;
178  break;
179  }
180  default:
181  {
182  ASSERTL0(false,"Expansion dimension not recognised");
183  break;
184  }
185  }
186 
187  // Write out VTK file.
188  string outname(strtok(argv[argc-1],"."));
189  outname += ".vtu";
190  ofstream outfile(outname.c_str());
191 
192  Exp[0]->WriteVtkHeader(outfile);
193 
194  if (jac)
195  {
196  // Find minimum Jacobian.
198  Array<OneD, NekDouble> x0 (Exp[0]->GetNpoints());
199  Array<OneD, NekDouble> x1 (Exp[0]->GetNpoints());
200  Array<OneD, NekDouble> x2 (Exp[0]->GetNpoints());
201  Exp[0]->GetCoords(x0, x1, x2);
202 
203  vector<NekDouble> jacDist;
204 
205  if (quality)
206  {
207  jacDist.resize(Exp[0]->GetExpSize());
208  }
209 
210  // Write out field containing Jacobian.
211  for(int i = 0; i < Exp[0]->GetExpSize(); ++i)
212  {
213  LocalRegions::ExpansionSharedPtr e = Exp[0]->GetExp(i);
214  SpatialDomains::GeomFactorsSharedPtr g = e->GetMetricInfo();
215  LibUtilities::PointsKeyVector ptsKeys = e->GetPointsKeys();
216  unsigned int npts = e->GetTotPoints();
217  NekDouble scaledJac = 1.0;
218 
219  if (g->GetGtype() == SpatialDomains::eDeformed)
220  {
221  const Array<OneD, const NekDouble> &jacobian
222  = g->GetJac(ptsKeys);
223  if (!quality)
224  {
225  Vmath::Vcopy(npts, jacobian, 1,
226  tmp = Exp[0]->UpdatePhys()
227  + Exp[0]->GetPhys_Offset(i), 1);
228  }
229  else
230  {
231  scaledJac = Vmath::Vmin(npts, jacobian, 1) /
232  Vmath::Vmax(npts, jacobian, 1);
233 
234  Vmath::Fill(npts, scaledJac,
235  tmp = Exp[0]->UpdatePhys()
236  + Exp[0]->GetPhys_Offset(i), 1);
237  }
238  }
239  else
240  {
241  Vmath::Fill (npts, g->GetJac(ptsKeys)[0],
242  tmp = Exp[0]->UpdatePhys()
243  + Exp[0]->GetPhys_Offset(i), 1);
244  }
245 
246  if (quality)
247  {
248  jacDist[i] = scaledJac;
249  }
250 
251  Exp[0]->WriteVtkPieceHeader(outfile, i);
252  Exp[0]->WriteVtkPieceData (outfile, i, "Jac");
253  Exp[0]->WriteVtkPieceFooter(outfile, i);
254  }
255 
256  unsigned int n
257  = Vmath::Imin(Exp[0]->GetNpoints(), Exp[0]->GetPhys(), 1);
258  cout << "- Minimum Jacobian: "
259  << Vmath::Vmin(Exp[0]->GetNpoints(), Exp[0]->GetPhys(), 1)
260  << " at coords (" << x0[n] << ", " << x1[n] << ", " << x2[n] << ")"
261  << endl;
262 
263  if (quality)
264  {
265  string distName = vSession->GetSessionName() + ".jac";
266  ofstream dist(distName.c_str());
267  dist.setf (ios::scientific, ios::floatfield);
268 
269  for (int i = 0; i < Exp[0]->GetExpSize(); ++i)
270  {
271  dist << setw(10) << i << " "
272  << setw(20) << setprecision(15) << jacDist[i] << endl;
273  }
274 
275  dist.close();
276 
277  cout << "- Minimum/maximum scaled Jacobian: "
278  << Vmath::Vmin(Exp[0]->GetExpSize(), &jacDist[0], 1) << " "
279  << Vmath::Vmax(Exp[0]->GetExpSize(), &jacDist[0], 1)
280  << endl;
281  }
282  }
283  else
284  {
285  // For each field write header and footer, since there is no field data.
286  for(int i = 0; i < Exp[0]->GetExpSize(); ++i)
287  {
288  Exp[0]->WriteVtkPieceHeader(outfile, i);
289  Exp[0]->WriteVtkPieceFooter(outfile, i);
290  }
291  }
292 
293  Exp[0]->WriteVtkFooter(outfile);
294 
295  return 0;
296 }
297 
static std::string RegisterCmdLineFlag(const std::string &pName, const std::string &pShortName, const std::string &pDescription)
Registers a command-line flag with the session reader.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:119
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:756
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:848
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:139
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:824
int main(int argc, char *argv[])
Definition: XmlToVtk.cpp:48
Fourier Expansion .
Definition: BasisType.h:52
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:63
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
static std::string npts
Definition: InputFld.cpp:43
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
boost::shared_ptr< ExpList2DHomogeneous1D > ExpList2DHomogeneous1DSharedPtr
Shared pointer to an ExpList2DHomogeneous1D object.
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
Definition: ExpList3D.h:110
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:145
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Geometry is curved or has non-constant factors.
Describes the specification for a Basis.
Definition: Basis.h:50
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:171
std::map< int, ExpansionShPtr >::iterator ExpansionMapIter
Definition: MeshGraph.h:172