16 int main(
int argc,
char *argv[])
19 int surfID, nfiles, nStart;
23 fprintf(stderr,
"Usage: FldAddWSS meshfile nfiles infld BoundaryID\n");
27 surfID = boost::lexical_cast<
int>(argv[argc - 1]);
28 nfiles = boost::lexical_cast<
int>(argv[2]);
30 vector<string> infiles(nfiles), outfiles(nfiles);
40 string basename = argv[3];
41 basename = basename.substr(0, basename.find_last_of(
"."));
42 stringstream filename2;
43 filename2 << basename <<
"_wss.fld";
44 filename2 >> outfiles[0];
48 string basename = argv[3];
49 basename = basename.substr(basename.find_last_of(
"t")+1, basename.find_last_of(
".")-basename.find_last_of(
"t"));
50 stringstream filename3;
51 filename3 << basename;
54 for (i = 0; i< nfiles; ++i)
57 string extension =
".fld";
58 basename = basename.substr(0, basename.find_first_of(
"_"));
59 stringstream filename, filename2;
60 filename << basename <<
"_t" << i + nStart << extension;
61 filename >> infiles[i];
62 filename2 << basename <<
"_t" << i +nStart <<
"_wss.fld";
63 filename2 >> outfiles[i];
67 argv[argc - 1] = argv[argc - 2];
68 argv[argc - 3] = argv[argc - 2];
71 = LibUtilities::SessionReader::CreateInstance(argc, argv);
75 string meshfile(argv[argc-4]);
82 m_kinvis = vSession->GetParameter(
"Kinvis");
87 string fieldfile(argv[3]);
88 vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
89 vector<vector<NekDouble> > fielddata;
95 int expdim = graphShPt->GetMeshDimension();
96 int nfields = fielddef[0]->m_fields.size()-1;
97 int addfields = (nfields == 3)? 4:3;
98 int nstress = (nfields == 3)? 6:3;
106 ASSERTL0(
false,
"Expansion dimension not recognised");
112 ASSERTL0(
false,
"Not implemented in 2D");
120 for(i = 1; i < nfields; ++i)
132 vSession->GetVariable(0));
134 m_locToGlobalMap = firstfield->GetLocalToGlobalMap();
137 for(i = 1; i < nfields; ++i)
141 vSession->GetVariable(i));
144 for(i = 0; i < addfields; ++i)
148 vSession->GetVariable(0));
153 ASSERTL0(
false,
"Expansion dimension not recognised");
160 int n, cnt, elmtid, nq, offset, nt, boundary, nfq;
172 for (
int fileNo = 0; fileNo < nfiles ; ++fileNo)
183 for(j = 0; j < nfields; ++j)
185 for(i = 0; i < fielddata.size(); ++i)
187 vel[j]->ExtractDataToCoeffs(fielddef[i],
189 fielddef[i]->m_fields[j],
190 vel[j]->UpdateCoeffs());
193 vel[j]->BwdTrans(vel[j]->GetCoeffs(),vel[j]->UpdatePhys());
196 for(j = 0; j < addfields; ++j)
198 for(i = 0; i < fielddata.size(); ++i)
200 shear[j]->ExtractDataToCoeffs(fielddef[i],
202 fielddef[i]->m_fields[0],
203 shear[j]->UpdateCoeffs());
207 shear[j]->BwdTrans(shear[j]->GetCoeffs(),shear[j]->UpdatePhys());
215 nt = shear[0]->GetNpoints();
218 for(j = 0; j < addfields; ++j)
224 shear[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
227 for(j = 0; j < addfields; ++j)
229 BndExp[j] = shear[j]->GetBndCondExpansions();
233 for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
238 for(i = 0; i < BndExp[0][n]->GetExpSize(); ++i, cnt++)
241 elmtid = BoundarytoElmtID[cnt];
242 elmt = shear[0]->GetExp(elmtid);
243 nq = elmt->GetTotPoints();
244 offset = shear[0]->GetPhys_Offset(elmtid);
248 for(j = 0; j < nfields*nfields; ++j)
253 for(j = 0; j < nstress; ++j)
270 boundary = BoundarytoTraceID[cnt];
276 = elmt->GetFaceNormal(boundary);
279 for(j = 0; j < nstress; ++j)
286 U = vel[0]->GetPhys() + offset;
287 V = vel[1]->GetPhys() + offset;
288 W = vel[2]->GetPhys() + offset;
291 elmt->PhysDeriv(U,grad[0],grad[1],grad[2]);
292 elmt->PhysDeriv(V,grad[3],grad[4],grad[5]);
293 elmt->PhysDeriv(W,grad[6],grad[7],grad[8]);
297 Vmath::Smul (nq,(2*m_kinvis),grad[0],1,stress[0],1);
299 Vmath::Smul (nq,(2*m_kinvis),grad[4],1,stress[1],1);
301 Vmath::Smul (nq,(2*m_kinvis),grad[8],1,stress[2],1);
313 for(j = 0; j < nstress; ++j)
315 elmt->GetFacePhysVals(boundary,bc,stress[j],fstress[j]);
319 for (j = 0; j< addfields; j++)
321 values[j] = BndExp[j][n]->UpdateCoeffs() + BndExp[j][n]->GetCoeff_Offset(i);
328 normals[1],1,fstress[3],1,Sx,1);
333 normals[1],1,fstress[1],1,Sy,1);
338 normals[1],1,fstress[5],1,Sz,1);
345 normals[1][0],fstress[3],1,Sx,1);
350 normals[1][0],fstress[1],1,Sy,1);
355 normals[1][0],fstress[5],1,Sz,1);
363 normals[1],1, Sy,1,values2,1);
364 Vmath::Vvtvp (nfq,normals[2],1, Sz,1,values2,1,values2,1);
371 bc->FwdTrans(Sx, values[0]);
375 bc->FwdTrans(Sy, values[1]);
379 bc->FwdTrans(Sz, values[2]);
385 normals[1][0],Sy,1,values2,1);
386 Vmath::Svtvp(nfq,normals[2][0],Sz,1,values2,1,values2,1);
392 bc->FwdTrans(Sx, values[0]);
396 bc->FwdTrans(Sy, values[1]);
400 bc->FwdTrans(Sz, values[2]);
405 Vmath::Vvtvvtp(nfq, Sx, 1, Sx, 1, Sy, 1, Sy, 1, S, 1);
408 bc->FwdTrans(S, values[3]);
414 cnt += BndExp[0][n]->GetExpSize();
419 for(j = 0; j < addfields; ++j)
421 int ncoeffs = shear[j]->GetNcoeffs();
424 output=shear[j]->UpdateCoeffs();
426 int nGlobal=m_locToGlobalMap->GetNumGlobalCoeffs();
434 for(i = 0; i < BndExp[j].num_elements(); ++i)
439 for(k = 0; k < (BndExp[j][i])->GetNcoeffs(); ++k)
441 sign = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
442 outarray[map[bndcnt++]] = sign * coeffs[k];
447 bndcnt += BndExp[j][i]->GetNcoeffs();
450 m_locToGlobalMap->GlobalToLocal(outarray,output);
456 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
457 = shear[0]->GetFieldDefinitions();
458 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
460 vector<string > outname;
468 outname.push_back(
"Tx");
469 outname.push_back(
"Ty");
470 outname.push_back(
"Tz");
471 outname.push_back(
"Tw");
474 for(j = 0; j < nfields + addfields; ++j)
476 for(i = 0; i < FieldDef.size(); ++i)
480 FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
481 vel[j]->AppendFieldData(FieldDef[i], FieldData[i]);
485 FieldDef[i]->m_fields.push_back(outname[j-nfields]);
486 shear[j-nfields]->AppendFieldData(FieldDef[i], FieldData[i]);
#define ASSERTL0(condition, msg)
#define sign(a, b)
return the sign(b)*a
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
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< ContField2D > ContField2DSharedPtr
int main(int argc, char *argv[])
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
boost::shared_ptr< StdExpansion2D > StdExpansion2DSharedPtr
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap, const bool backup)
This function allows for data to be written to an FLD file when a session and/or communicator is not ...
boost::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
void Vvtvvtp(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)
vvtvvtp (vector times vector plus vector times vector):
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
boost::shared_ptr< ContField3D > ContField3DSharedPtr
boost::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Geometry is curved or has non-constant factors.
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.