Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 257 of file XmlToVtk.cpp.

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

Referenced by main().

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 }
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:442
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:564
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1878
int main ( int  argc,
char *  argv[] 
)

Definition at line 53 of file XmlToVtk.cpp.

References ASSERTL0, Nektar::LibUtilities::eFourier, Nektar::LibUtilities::ePolyEvenlySpaced, Nektar::LibUtilities::BasisKey::GetBasisType(), Nektar::LibUtilities::BasisKey::GetNumModes(), Nektar::LibUtilities::BasisKey::GetPointsKey(), GetQ(), Vmath::Imin(), npts, CellMLToNektar.cellml_metadata::p, Vmath::Vcopy(), and Vmath::Vmin().

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 
66  LibUtilities::SessionReaderSharedPtr vSession
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]);
76  SpatialDomains::MeshGraphSharedPtr graphShPt =
77  SpatialDomains::MeshGraph::Read(vSession); //meshfile);
78 
79  // Set up Expansion information
80  SpatialDomains::ExpansionMap emap = graphShPt->GetExpansions();
81  SpatialDomains::ExpansionMapIter it;
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];
88  LibUtilities::PointsKey tmp2 = tmp1.GetPointsKey();
89  it->second->m_basisKeyVector[i] = LibUtilities::BasisKey(
90  tmp1.GetBasisType(), tmp1.GetNumModes(),
91  LibUtilities::PointsKey(tmp1.GetNumModes(),
92  LibUtilities::ePolyEvenlySpaced));
93  }
94  }
95 
96  // Define Expansion
97  int expdim = graphShPt->GetMeshDimension();
98  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(1);
99 
100  switch(expdim)
101  {
102  case 1:
103  {
104  if(vSession->DefinesSolverInfo("HOMOGENEOUS"))
105  {
106  std::string HomoStr = vSession->GetSolverInfo("HOMOGENEOUS");
107  MultiRegions::ExpList2DHomogeneous1DSharedPtr Exp2DH1;
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 
125  Exp2DH1 = MemoryManager<MultiRegions::ExpList2DHomogeneous1D>
126  ::AllocateSharedPtr(
127  vSession, Bkey, lz, false, false, graphShPt);
128  Exp[0] = Exp2DH1;
129  }
130  else
131  {
132  MultiRegions::ExpList1DSharedPtr Exp1D;
133  Exp1D = MemoryManager<MultiRegions::ExpList1D>
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");
145  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
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 
163  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>
164  ::AllocateSharedPtr(
165  vSession, Bkey, lz, false, false, graphShPt);
166  Exp[0] = Exp3DH1;
167  }
168  else
169  {
170  MultiRegions::ExpList2DSharedPtr Exp2D;
171  Exp2D = MemoryManager<MultiRegions::ExpList2D>
172  ::AllocateSharedPtr(vSession,graphShPt);
173  Exp[0] = Exp2D;
174  }
175  break;
176  }
177  case 3:
178  {
179  MultiRegions::ExpList3DSharedPtr Exp3D;
180  Exp3D = MemoryManager<MultiRegions::ExpList3D>
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.
202  Array<OneD, NekDouble> tmp;
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 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
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:871
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:847
Nektar::LibUtilities::eFourier
Fourier Expansion .
Definition: BasisType.h:52
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
Nektar::LibUtilities::ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
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:257
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:112
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:1878
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:1061
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 Sat Nov 25 2017 16:18:36 for Nektar++ by   doxygen 1.8.8