Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
XmlToVtk.cpp File Reference
#include <cstdlib>
#include <iomanip>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
Include dependency graph for XmlToVtk.cpp:

Go to the source code of this file.

Functions

Array< OneD, NekDouble > GetQ (SpatialDomains::DerivStorage &deriv, LocalRegions::ExpansionSharedPtr exp)
 
int main (int argc, char *argv[])
 

Function Documentation

Array< OneD, NekDouble > GetQ ( SpatialDomains::DerivStorage &  deriv,
LocalRegions::ExpansionSharedPtr  exp 
)

Definition at line 256 of file XmlToVtk.cpp.

References Vmath::Fill(), Vmath::Vvtvp(), and Vmath::Vvtvvtm().

Referenced by main().

259 {
260  StdRegions::StdExpansionSharedPtr chi = exp->GetGeom()->GetXmap();
261  const int pts = deriv[0][0].num_elements();
262  const int nq = chi->GetTotPoints();
263  const int expDim = chi->GetNumBases();
264 
265  Array<OneD, NekDouble> jac(pts, 0.0);
266  Array<OneD, NekDouble> eta(nq);
267 
268  switch (expDim)
269  {
270  case 2:
271  {
272  Vmath::Vvtvvtm(pts, &deriv[0][0][0], 1, &deriv[1][1][0], 1,
273  &deriv[1][0][0], 1, &deriv[0][1][0], 1,
274  &jac[0], 1);
275  break;
276  }
277  case 3:
278  {
279  Array<OneD, NekDouble> tmp(pts, 0.0);
280  Vmath::Vvtvvtm(pts, &deriv[1][1][0], 1, &deriv[2][2][0], 1,
281  &deriv[2][1][0], 1, &deriv[1][2][0], 1,
282  &tmp[0], 1);
283  Vmath::Vvtvp (pts, &deriv[0][0][0], 1, &tmp[0], 1,
284  &jac[0], 1, &jac[0], 1);
285  Vmath::Vvtvvtm(pts, &deriv[2][1][0], 1, &deriv[0][2][0], 1,
286  &deriv[0][1][0], 1, &deriv[2][2][0], 1,
287  &tmp[0], 1);
288  Vmath::Vvtvp (pts, &deriv[1][0][0], 1, &tmp[0], 1,
289  &jac[0], 1, &jac[0], 1);
290  Vmath::Vvtvvtm(pts, &deriv[0][1][0], 1, &deriv[1][2][0], 1,
291  &deriv[1][1][0], 1, &deriv[0][2][0], 1,
292  &tmp[0], 1);
293  Vmath::Vvtvp (pts, &deriv[2][0][0], 1, &tmp[0], 1,
294  &jac[0], 1, &jac[0], 1);
295  break;
296  }
297  }
298 
299  for (int k = 0; k < pts; ++k)
300  {
301  NekDouble frob = 0.0;
302 
303  for (int i = 0; i < deriv.num_elements(); ++i)
304  {
305  for (int j = 0; j < deriv[i].num_elements(); ++j)
306  {
307  frob += deriv[i][j][k]*deriv[i][j][k];
308  }
309  }
310 
311  NekDouble delta = 0.0; // fixme
312  NekDouble sigma = 0.5*(jac[k] + sqrt(jac[k]*jac[k] + 4*delta*delta));
313 
314  eta[k] = expDim * pow(sigma,2.0/expDim) / frob;
315  if(eta[k] > 1 || eta[k] < 0)
316  cout << eta[k] << endl;
317  }
318 
319  if (pts == 1)
320  {
321  Vmath::Fill(nq-1, eta[0], &eta[1], 1);
322  }
323 
324  return eta;
325 }
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Vmath::Vvtvp
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
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Vvtvvtm
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
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1873
int main ( int  argc,
char *  argv[] 
)

Definition at line 52 of file XmlToVtk.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::eFourier, Nektar::LibUtilities::ePolyEvenlySpaced, Nektar::LibUtilities::BasisKey::GetBasisType(), Nektar::LibUtilities::BasisKey::GetNumModes(), Nektar::LibUtilities::BasisKey::GetPointsKey(), GetQ(), Vmath::Imin(), npts, Nektar::SpatialDomains::MeshGraph::Read(), Nektar::LibUtilities::SessionReader::RegisterCmdLineFlag(), Vmath::Vcopy(), and Vmath::Vmin().

53 {
54  if(argc < 2)
55  {
56  cerr << "Usage: XmlToVtk meshfile" << endl;
57  exit(1);
58  }
59 
60  LibUtilities::SessionReader::RegisterCmdLineFlag(
61  "jacobian", "j", "Output Jacobian as scalar field");
62  LibUtilities::SessionReader::RegisterCmdLineFlag(
63  "quality", "q", "Output distribution of scaled Jacobians");
64 
65  LibUtilities::SessionReaderSharedPtr vSession
66  = LibUtilities::SessionReader::CreateInstance(argc, argv);
67 
68  bool jac = vSession->DefinesCmdLineArgument("jacobian");
69  bool quality = vSession->DefinesCmdLineArgument("quality");
70 
71  jac = quality ? true : jac;
72 
73  // Read in mesh from input file
74  string meshfile(argv[argc-1]);
75  SpatialDomains::MeshGraphSharedPtr graphShPt =
76  SpatialDomains::MeshGraph::Read(vSession); //meshfile);
77 
78  // Set up Expansion information
79  SpatialDomains::ExpansionMap emap = graphShPt->GetExpansions();
80  SpatialDomains::ExpansionMapIter it;
81 
82  for (it = emap.begin(); it != emap.end(); ++it)
83  {
84  for (int i = 0; i < it->second->m_basisKeyVector.size(); ++i)
85  {
86  LibUtilities::BasisKey tmp1 = it->second->m_basisKeyVector[i];
87  LibUtilities::PointsKey tmp2 = tmp1.GetPointsKey();
88  it->second->m_basisKeyVector[i] = LibUtilities::BasisKey(
89  tmp1.GetBasisType(), tmp1.GetNumModes(),
90  LibUtilities::PointsKey(tmp1.GetNumModes(),
91  LibUtilities::ePolyEvenlySpaced));
92  }
93  }
94 
95  // Define Expansion
96  int expdim = graphShPt->GetMeshDimension();
97  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(1);
98 
99  switch(expdim)
100  {
101  case 1:
102  {
103  if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
104  {
105  std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
106  MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
107 
108  ASSERTL0(
109  HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
110  HomoStr == "1D" || HomoStr == "Homo1D",
111  "Only 3DH1D supported for XML output currently.");
112 
113  int nplanes;
114  vSession->LoadParameter("HomModesZ", nplanes);
115 
116  // choose points to be at evenly spaced points at nplanes + 1
117  // points
118  const LibUtilities::PointsKey Pkey(
119  nplanes + 1, LibUtilities::ePolyEvenlySpaced);
120  const LibUtilities::BasisKey Bkey(
121  LibUtilities::eFourier, nplanes, Pkey);
122  NekDouble lz = vSession->GetParameter("LZ");
123 
124  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>
125  ::AllocateSharedPtr(
126  vSession, Bkey, lz, false, false, graphShPt);
127  Exp[0] = Exp2DH1;
128  }
129  else
130  {
131  MultiRegions::ExpList1DSharedPtr Exp1D;
132  Exp1D = MemoryManager<MultiRegions::ExpList1D>
133  ::AllocateSharedPtr(vSession,graphShPt);
134  Exp[0] = Exp1D;
135  }
136 
137  break;
138  }
139  case 2:
140  {
141  if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
142  {
143  std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
144  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
145 
146  ASSERTL0(
147  HomoStr == "HOMOGENEOUS1D" || HomoStr == "Homogeneous1D" ||
148  HomoStr == "1D" || HomoStr == "Homo1D",
149  "Only 3DH1D supported for XML output currently.");
150 
151  int nplanes;
152  vSession->LoadParameter("HomModesZ", nplanes);
153 
154  // choose points to be at evenly spaced points at nplanes + 1
155  // points
156  const LibUtilities::PointsKey Pkey(
157  nplanes + 1, LibUtilities::ePolyEvenlySpaced);
158  const LibUtilities::BasisKey Bkey(
159  LibUtilities::eFourier, nplanes, Pkey);
160  NekDouble lz = vSession->GetParameter("LZ");
161 
162  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>
163  ::AllocateSharedPtr(
164  vSession, Bkey, lz, false, false, graphShPt);
165  Exp[0] = Exp3DH1;
166  }
167  else
168  {
169  MultiRegions::ExpList2DSharedPtr Exp2D;
170  Exp2D = MemoryManager<MultiRegions::ExpList2D>
171  ::AllocateSharedPtr(vSession,graphShPt);
172  Exp[0] = Exp2D;
173  }
174  break;
175  }
176  case 3:
177  {
178  MultiRegions::ExpList3DSharedPtr Exp3D;
179  Exp3D = MemoryManager<MultiRegions::ExpList3D>
180  ::AllocateSharedPtr(vSession,graphShPt);
181  Exp[0] = Exp3D;
182  break;
183  }
184  default:
185  {
186  ASSERTL0(false,"Expansion dimension not recognised");
187  break;
188  }
189  }
190 
191  // Write out VTK file.
192  string outname(strtok(argv[argc-1],"."));
193  outname += ".vtu";
194  ofstream outfile(outname.c_str());
195 
196  Exp[0]->WriteVtkHeader(outfile);
197 
198  if (jac)
199  {
200  // Find minimum Jacobian.
201  Array<OneD, NekDouble> tmp;
202  Array<OneD, NekDouble> x0 (Exp[0]->GetNpoints());
203  Array<OneD, NekDouble> x1 (Exp[0]->GetNpoints());
204  Array<OneD, NekDouble> x2 (Exp[0]->GetNpoints());
205  Exp[0]->GetCoords(x0, x1, x2);
206 
207  // Write out field containing Jacobian.
208  for(int i = 0; i < Exp[0]->GetExpSize(); ++i)
209  {
210  LocalRegions::ExpansionSharedPtr e = Exp[0]->GetExp(i);
211  SpatialDomains::GeometrySharedPtr geom = e->GetGeom();
212  StdRegions::StdExpansionSharedPtr chi = geom->GetXmap();
213  LibUtilities::PointsKeyVector p = chi->GetPointsKeys();
214  SpatialDomains::GeomFactorsSharedPtr gfac = geom->GetGeomFactors();
215 
216  // Grab Jacobian entries (on collapsed quadrature points)
217  SpatialDomains::DerivStorage deriv = gfac->GetDeriv(p);
218  const int npts = chi->GetTotPoints();
219 
220 
221  Array<OneD, const NekDouble> q = GetQ(deriv, e);
222  Vmath::Vcopy(npts, q, 1,
223  tmp = Exp[0]->UpdatePhys()
224  + Exp[0]->GetPhys_Offset(i), 1);
225 
226 
227 
228  Exp[0]->WriteVtkPieceHeader(outfile, i);
229  Exp[0]->WriteVtkPieceData (outfile, i, "Q");
230  Exp[0]->WriteVtkPieceFooter(outfile, i);
231  }
232 
233  unsigned int n
234  = Vmath::Imin(Exp[0]->GetNpoints(), Exp[0]->GetPhys(), 1);
235  cout << "- Minimum Jacobian: "
236  << Vmath::Vmin(Exp[0]->GetNpoints(), Exp[0]->GetPhys(), 1)
237  << " at coords (" << x0[n] << ", " << x1[n] << ", " << x2[n] << ")"
238  << endl;
239 
240  }
241  else
242  {
243  // For each field write header and footer, since there is no field data.
244  for(int i = 0; i < Exp[0]->GetExpSize(); ++i)
245  {
246  Exp[0]->WriteVtkPieceHeader(outfile, i);
247  Exp[0]->WriteVtkPieceFooter(outfile, i);
248  }
249  }
250 
251  Exp[0]->WriteVtkFooter(outfile);
252 
253  return 0;
254 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Vmath::Vmin
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
Nektar::LibUtilities::BasisKey::GetBasisType
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:139
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:833
Nektar::LibUtilities::eFourier
Fourier Expansion .
Definition: BasisType.h:52
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::LibUtilities::ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:63
Nektar::MultiRegions::ExpList1DSharedPtr
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::LocalRegions::ExpansionSharedPtr
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Nektar::MultiRegions::ExpList2DSharedPtr
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
GetQ
Array< OneD, NekDouble > GetQ(SpatialDomains::DerivStorage &deriv, LocalRegions::ExpansionSharedPtr exp)
Definition: XmlToVtk.cpp:256
Nektar::MultiRegions::ExpList2DHomogeneous1DSharedPtr
boost::shared_ptr< ExpList2DHomogeneous1D > ExpList2DHomogeneous1DSharedPtr
Shared pointer to an ExpList2DHomogeneous1D object.
Definition: ExpList2DHomogeneous1D.h:50
Nektar::MultiRegions::ExpList3DSharedPtr
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
Definition: ExpList3D.h:114
Nektar::LibUtilities::BasisKey::GetPointsKey
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:145
Nektar::SpatialDomains::GeomFactorsSharedPtr
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
Nektar::MultiRegions::ExpList3DHomogeneous1DSharedPtr
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
Definition: ExpList3DHomogeneous1D.h:50
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1873
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Nektar::SpatialDomains::GeometrySharedPtr
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::SpatialDomains::ExpansionMap
std::map< int, ExpansionShPtr > ExpansionMap
Definition: MeshGraph.h:174
Nektar::SpatialDomains::ExpansionMapIter
std::map< int, ExpansionShPtr >::iterator ExpansionMapIter
Definition: MeshGraph.h:175
Generated on Wed Mar 23 2016 20:14:27 for Nektar++ by   doxygen 1.8.8