Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 std;
47 using namespace Nektar;
48 
52 
53 int main(int argc, char *argv[])
54 {
55  if(argc < 2)
56  {
57  cerr << "Usage: XmlToVtk meshfile" << endl;
58  exit(1);
59  }
60 
61  LibUtilities::SessionReader::RegisterCmdLineFlag(
62  "jacobian", "j", "Output Jacobian as scalar field");
63  LibUtilities::SessionReader::RegisterCmdLineFlag(
64  "quality", "q", "Output distribution of scaled Jacobians");
65 
67  = LibUtilities::SessionReader::CreateInstance(argc, argv);
68 
69  bool jac = vSession->DefinesCmdLineArgument("jacobian");
70  bool quality = vSession->DefinesCmdLineArgument("quality");
71 
72  jac = quality ? true : jac;
73 
74  // Read in mesh from input file
75  string meshfile(argv[argc-1]);
77  SpatialDomains::MeshGraph::Read(vSession); //meshfile);
78 
79  // Set up Expansion information
80  SpatialDomains::ExpansionMap emap = graphShPt->GetExpansions();
82 
83  for (it = emap.begin(); it != emap.end(); ++it)
84  {
85  for (int i = 0; i < it->second->m_basisKeyVector.size(); ++i)
86  {
87  LibUtilities::BasisKey tmp1 = it->second->m_basisKeyVector[i];
89  it->second->m_basisKeyVector[i] = LibUtilities::BasisKey(
90  tmp1.GetBasisType(), tmp1.GetNumModes(),
93  }
94  }
95 
96  // Define Expansion
97  int expdim = graphShPt->GetMeshDimension();
99 
100  switch(expdim)
101  {
102  case 1:
103  {
104  if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
105  {
106  std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
108 
109  ASSERTL0(
110  HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
111  HomoStr == "1D" || HomoStr == "Homo1D",
112  "Only 3DH1D supported for XML output currently.");
113 
114  int nplanes;
115  vSession->LoadParameter("HomModesZ", nplanes);
116 
117  // choose points to be at evenly spaced points at nplanes + 1
118  // points
119  const LibUtilities::PointsKey Pkey(
120  nplanes + 1, LibUtilities::ePolyEvenlySpaced);
121  const LibUtilities::BasisKey Bkey(
122  LibUtilities::eFourier, nplanes, Pkey);
123  NekDouble lz = vSession->GetParameter("LZ");
124 
127  vSession, Bkey, lz, false, false, graphShPt);
128  Exp[0] = Exp2DH1;
129  }
130  else
131  {
134  ::AllocateSharedPtr(vSession,graphShPt);
135  Exp[0] = Exp1D;
136  }
137 
138  break;
139  }
140  case 2:
141  {
142  if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
143  {
144  std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
146 
147  ASSERTL0(
148  HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
149  HomoStr == "1D" || HomoStr == "Homo1D",
150  "Only 3DH1D supported for XML output currently.");
151 
152  int nplanes;
153  vSession->LoadParameter("HomModesZ", nplanes);
154 
155  // choose points to be at evenly spaced points at nplanes + 1
156  // points
157  const LibUtilities::PointsKey Pkey(
158  nplanes + 1, LibUtilities::ePolyEvenlySpaced);
159  const LibUtilities::BasisKey Bkey(
160  LibUtilities::eFourier, nplanes, Pkey);
161  NekDouble lz = vSession->GetParameter("LZ");
162 
165  vSession, Bkey, lz, false, false, graphShPt);
166  Exp[0] = Exp3DH1;
167  }
168  else
169  {
172  ::AllocateSharedPtr(vSession,graphShPt);
173  Exp[0] = Exp2D;
174  }
175  break;
176  }
177  case 3:
178  {
181  ::AllocateSharedPtr(vSession,graphShPt);
182  Exp[0] = Exp3D;
183  break;
184  }
185  default:
186  {
187  ASSERTL0(false,"Expansion dimension not recognised");
188  break;
189  }
190  }
191 
192  // Write out VTK file.
193  string outname(strtok(argv[argc-1],"."));
194  outname += ".vtu";
195  ofstream outfile(outname.c_str());
196 
197  Exp[0]->WriteVtkHeader(outfile);
198 
199  if (jac)
200  {
201  // Find minimum Jacobian.
203  Array<OneD, NekDouble> x0 (Exp[0]->GetNpoints());
204  Array<OneD, NekDouble> x1 (Exp[0]->GetNpoints());
205  Array<OneD, NekDouble> x2 (Exp[0]->GetNpoints());
206  Exp[0]->GetCoords(x0, x1, x2);
207 
208  // Write out field containing Jacobian.
209  for(int i = 0; i < Exp[0]->GetExpSize(); ++i)
210  {
211  LocalRegions::ExpansionSharedPtr e = Exp[0]->GetExp(i);
212  SpatialDomains::GeometrySharedPtr geom = e->GetGeom();
213  StdRegions::StdExpansionSharedPtr chi = geom->GetXmap();
214  LibUtilities::PointsKeyVector p = chi->GetPointsKeys();
215  SpatialDomains::GeomFactorsSharedPtr gfac = geom->GetGeomFactors();
216 
217  // Grab Jacobian entries (on collapsed quadrature points)
218  SpatialDomains::DerivStorage deriv = gfac->GetDeriv(p);
219  const int npts = chi->GetTotPoints();
220 
221 
222  Array<OneD, const NekDouble> q = GetQ(deriv, e);
223  Vmath::Vcopy(npts, q, 1,
224  tmp = Exp[0]->UpdatePhys()
225  + Exp[0]->GetPhys_Offset(i), 1);
226 
227 
228 
229  Exp[0]->WriteVtkPieceHeader(outfile, i);
230  Exp[0]->WriteVtkPieceData (outfile, i, "Q");
231  Exp[0]->WriteVtkPieceFooter(outfile, i);
232  }
233 
234  unsigned int n
235  = Vmath::Imin(Exp[0]->GetNpoints(), Exp[0]->GetPhys(), 1);
236  cout << "- Minimum Jacobian: "
237  << Vmath::Vmin(Exp[0]->GetNpoints(), Exp[0]->GetPhys(), 1)
238  << " at coords (" << x0[n] << ", " << x1[n] << ", " << x2[n] << ")"
239  << endl;
240 
241  }
242  else
243  {
244  // For each field write header and footer, since there is no field data.
245  for(int i = 0; i < Exp[0]->GetExpSize(); ++i)
246  {
247  Exp[0]->WriteVtkPieceHeader(outfile, i);
248  Exp[0]->WriteVtkPieceFooter(outfile, i);
249  }
250  }
251 
252  Exp[0]->WriteVtkFooter(outfile);
253 
254  return 0;
255 }
256 
260 {
261  StdRegions::StdExpansionSharedPtr chi = exp->GetGeom()->GetXmap();
262  const int pts = deriv[0][0].num_elements();
263  const int nq = chi->GetTotPoints();
264  const int expDim = chi->GetNumBases();
265 
266  Array<OneD, NekDouble> jac(pts, 0.0);
267  Array<OneD, NekDouble> eta(nq);
268 
269  switch (expDim)
270  {
271  case 2:
272  {
273  Vmath::Vvtvvtm(pts, &deriv[0][0][0], 1, &deriv[1][1][0], 1,
274  &deriv[1][0][0], 1, &deriv[0][1][0], 1,
275  &jac[0], 1);
276  break;
277  }
278  case 3:
279  {
280  Array<OneD, NekDouble> tmp(pts, 0.0);
281  Vmath::Vvtvvtm(pts, &deriv[1][1][0], 1, &deriv[2][2][0], 1,
282  &deriv[2][1][0], 1, &deriv[1][2][0], 1,
283  &tmp[0], 1);
284  Vmath::Vvtvp (pts, &deriv[0][0][0], 1, &tmp[0], 1,
285  &jac[0], 1, &jac[0], 1);
286  Vmath::Vvtvvtm(pts, &deriv[2][1][0], 1, &deriv[0][2][0], 1,
287  &deriv[0][1][0], 1, &deriv[2][2][0], 1,
288  &tmp[0], 1);
289  Vmath::Vvtvp (pts, &deriv[1][0][0], 1, &tmp[0], 1,
290  &jac[0], 1, &jac[0], 1);
291  Vmath::Vvtvvtm(pts, &deriv[0][1][0], 1, &deriv[1][2][0], 1,
292  &deriv[1][1][0], 1, &deriv[0][2][0], 1,
293  &tmp[0], 1);
294  Vmath::Vvtvp (pts, &deriv[2][0][0], 1, &tmp[0], 1,
295  &jac[0], 1, &jac[0], 1);
296  break;
297  }
298  }
299 
300  for (int k = 0; k < pts; ++k)
301  {
302  NekDouble frob = 0.0;
303 
304  for (int i = 0; i < deriv.num_elements(); ++i)
305  {
306  for (int j = 0; j < deriv[i].num_elements(); ++j)
307  {
308  frob += deriv[i][j][k]*deriv[i][j][k];
309  }
310  }
311 
312  NekDouble delta = 0.0; // fixme
313  NekDouble sigma = 0.5*(jac[k] + sqrt(jac[k]*jac[k] + 4*delta*delta));
314 
315  eta[k] = expDim * pow(sigma,2.0/expDim) / frob;
316  if(eta[k] > 1 || eta[k] < 0)
317  cout << eta[k] << endl;
318  }
319 
320  if (pts == 1)
321  {
322  Vmath::Fill(nq-1, eta[0], &eta[1], 1);
323  }
324 
325  return eta;
326 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
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:857
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:833
int main(int argc, char *argv[])
Definition: XmlToVtk.cpp:53
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
STL namespace.
Fourier Expansion .
Definition: BasisType.h:52
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:63
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< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void Vvtvvtm(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtm (vector times vector minus vector times vector):
Definition: Vmath.cpp:550
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
Array< OneD, NekDouble > GetQ(SpatialDomains::DerivStorage &deriv, LocalRegions::ExpansionSharedPtr exp)
Definition: XmlToVtk.cpp:257
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:114
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< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
Describes the specification for a Basis.
Definition: Basis.h:50
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:174
std::map< int, ExpansionShPtr >::iterator ExpansionMapIter
Definition: MeshGraph.h:175