53 int main(
int argc,
char *argv[])
57 cerr <<
"Usage: XmlToVtk meshfile" << endl;
61 LibUtilities::SessionReader::RegisterCmdLineFlag(
62 "jacobian",
"j",
"Output Jacobian as scalar field");
63 LibUtilities::SessionReader::RegisterCmdLineFlag(
64 "quality",
"q",
"Output distribution of scaled Jacobians");
67 = LibUtilities::SessionReader::CreateInstance(argc, argv);
69 bool jac = vSession->DefinesCmdLineArgument(
"jacobian");
70 bool quality = vSession->DefinesCmdLineArgument(
"quality");
72 jac = quality ?
true : jac;
75 string meshfile(argv[argc-1]);
77 SpatialDomains::MeshGraph::Read(vSession);
83 for (it = emap.begin(); it != emap.end(); ++it)
85 for (
int i = 0; i < it->second->m_basisKeyVector.size(); ++i)
97 int expdim = graphShPt->GetMeshDimension();
104 if(vSession->DefinesSolverInfo(
"HOMOGENEOUS"))
106 std::string HomoStr = vSession->GetSolverInfo(
"HOMOGENEOUS");
110 HomoStr ==
"HOMOGENEOUS1D" || HomoStr ==
"Homogeneous1D" ||
111 HomoStr ==
"1D" || HomoStr ==
"Homo1D",
112 "Only 3DH1D supported for XML output currently.");
115 vSession->LoadParameter(
"HomModesZ", nplanes);
123 NekDouble lz = vSession->GetParameter(
"LZ");
127 vSession, Bkey, lz,
false,
false, graphShPt);
142 if(vSession->DefinesSolverInfo(
"HOMOGENEOUS"))
144 std::string HomoStr = vSession->GetSolverInfo(
"HOMOGENEOUS");
148 HomoStr ==
"HOMOGENEOUS1D" || HomoStr ==
"Homogeneous1D" ||
149 HomoStr ==
"1D" || HomoStr ==
"Homo1D",
150 "Only 3DH1D supported for XML output currently.");
153 vSession->LoadParameter(
"HomModesZ", nplanes);
161 NekDouble lz = vSession->GetParameter(
"LZ");
165 vSession, Bkey, lz,
false,
false, graphShPt);
187 ASSERTL0(
false,
"Expansion dimension not recognised");
193 string outname(strtok(argv[argc-1],
"."));
195 ofstream outfile(outname.c_str());
197 Exp[0]->WriteVtkHeader(outfile);
206 Exp[0]->GetCoords(x0, x1, x2);
209 for(
int i = 0; i < Exp[0]->GetExpSize(); ++i)
219 const int npts = chi->GetTotPoints();
224 tmp = Exp[0]->UpdatePhys()
225 + Exp[0]->GetPhys_Offset(i), 1);
229 Exp[0]->WriteVtkPieceHeader(outfile, i);
230 Exp[0]->WriteVtkPieceData (outfile, i,
"Q");
231 Exp[0]->WriteVtkPieceFooter(outfile, i);
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] <<
")"
245 for(
int i = 0; i < Exp[0]->GetExpSize(); ++i)
247 Exp[0]->WriteVtkPieceHeader(outfile, i);
248 Exp[0]->WriteVtkPieceFooter(outfile, i);
252 Exp[0]->WriteVtkFooter(outfile);
262 const int pts = deriv[0][0].num_elements();
263 const int nq = chi->GetTotPoints();
264 const int expDim = chi->GetNumBases();
274 &deriv[1][0][0], 1, &deriv[0][1][0], 1,
282 &deriv[2][1][0], 1, &deriv[1][2][0], 1,
285 &jac[0], 1, &jac[0], 1);
287 &deriv[0][1][0], 1, &deriv[2][2][0], 1,
290 &jac[0], 1, &jac[0], 1);
292 &deriv[1][1][0], 1, &deriv[0][2][0], 1,
295 &jac[0], 1, &jac[0], 1);
300 for (
int k = 0; k < pts; ++k)
304 for (
int i = 0; i < deriv.num_elements(); ++i)
306 for (
int j = 0; j < deriv[i].num_elements(); ++j)
308 frob += deriv[i][j][k]*deriv[i][j][k];
313 NekDouble sigma = 0.5*(jac[k] + sqrt(jac[k]*jac[k] + 4*delta*delta));
315 eta[k] = expDim * pow(sigma,2.0/expDim) / frob;
316 if(eta[k] > 1 || eta[k] < 0)
317 cout << eta[k] << endl;
#define ASSERTL0(condition, msg)
std::vector< PointsKey > PointsKeyVector
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.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
BasisType GetBasisType() const
Return type of expansion basis.
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
int main(int argc, char *argv[])
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
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
1D Evenly-spaced points using Lagrange polynomial
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Defines a specification for a set of points.
boost::shared_ptr< Expansion > ExpansionSharedPtr
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):
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Array< OneD, NekDouble > GetQ(SpatialDomains::DerivStorage &deriv, LocalRegions::ExpansionSharedPtr exp)
boost::shared_ptr< ExpList2DHomogeneous1D > ExpList2DHomogeneous1DSharedPtr
Shared pointer to an ExpList2DHomogeneous1D object.
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
PointsKey GetPointsKey() const
Return distribution of points.
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
int GetNumModes() const
Returns the order of the basis.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
boost::shared_ptr< Geometry > GeometrySharedPtr
Describes the specification for a Basis.
std::map< int, ExpansionShPtr > ExpansionMap
std::map< int, ExpansionShPtr >::iterator ExpansionMapIter