16 int main(
int argc,
char *argv[])
19 int nfiles, nStart, surfID;
22 fprintf(stderr,
"Usage: FldAddMultiShear meshfile nfiles FirstInFile BoundaryID\n");
26 nfiles = boost::lexical_cast<
int>(argv[2]);
27 surfID = boost::lexical_cast<
int>(argv[4]);
29 vector<string> infiles(nfiles);
33 stringstream filename2;
34 filename2 <<
"multishear.fld";
38 string basename = argv[3];
39 basename = basename.substr(basename.find_last_of(
"t")+1, basename.find_last_of(
".")-basename.find_last_of(
"t"));
40 stringstream filename3;
41 filename3 << basename;
44 for (i = 0; i< nfiles; ++i)
47 string extension =
".fld";
48 basename = basename.substr(0, basename.find_first_of(
"_"));
49 stringstream filename;
50 filename << basename <<
"_t" << i+nStart <<
"_wss.fld";
51 filename >> infiles[i];
52 cout << infiles[i]<<endl;
59 = LibUtilities::SessionReader::CreateInstance(argc, argv);
64 string meshfile(argv[1]);
71 vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
72 vector<vector<NekDouble> > fielddata;
80 int expdim = graphShPt->GetMeshDimension();
81 int nfields = fielddef[0]->m_fields.size();
83 int sfields = nfields - expdim;
91 ASSERTL0(
false,
"Expansion dimension not recognised");
96 ASSERTL0(
false,
"Expansion dimension not recognised");
105 vSession->GetVariable(i));
107 m_locToGlobalMap = firstfield->GetLocalToGlobalMap();
110 for(i = 1; i < expdim; ++i)
114 vSession->GetVariable(i));
117 for(i = 0; i < sfields; ++i)
121 vSession->GetVariable(i));
124 for(i = 0; i < addfields; ++i)
128 vSession->GetVariable(0));
134 ASSERTL0(
false,
"Expansion dimension not recognised");
151 for(i = 0; i < fielddata.size(); ++i)
153 Exp[0]->ExtractDataToCoeffs(fielddef[i],
155 fielddef[i]->m_fields[0],
156 Exp[0]->UpdateCoeffs());
158 Exp[0]->BwdTrans(Exp[0]->GetCoeffs(),Exp[0]->UpdatePhys());
160 Exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
161 BndExp[0] = Exp[0]->GetBndCondExpansions();
167 int nelem = BndExp[0][surfID]->GetExpSize();
170 int nt = Exp[0]->GetNpoints();
175 int n, cnt, elmtid, offset, boundary, bndOffset;
180 for (i = 0; i < sfields; i++)
194 for (
int fileNo = 0; fileNo < nfiles ; ++fileNo)
202 for(j = 0; j < nfields; ++j)
204 for(
int i = 0; i < fielddata.size(); ++i)
206 Exp[j]->ExtractDataToCoeffs(fielddef[i],
208 fielddef[i]->m_fields[j],
209 Exp[j]->UpdateCoeffs());
211 Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
214 Exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
215 BndExp[0] = Exp[0]->GetBndCondExpansions();
217 for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
221 for(i = 0; i < nelem; ++i, cnt++)
224 elmtid = BoundarytoElmtID[cnt];
225 elmt = Exp[0]->GetExp(elmtid);
226 offset = Exp[0]->GetPhys_Offset(elmtid);
239 boundary = BoundarytoTraceID[cnt];
242 for (
int t = 0; t< expdim; t++)
244 elmt->GetFacePhysVals(boundary,bc, Exp[t+expdim]->GetPhys() + offset, temp[t]);
247 for (
int m = 0; m < nfq; m++)
249 Sxr[bndOffset + m] = temp[0][m];
250 Syr[bndOffset + m] = temp[1][m];
251 Szr[bndOffset + m] = temp[2][m];
267 cnt += BndExp[0][n]->GetExpSize();
279 Vmath::Vvtvvtp(nbq, Sx, 1, Sx, 1, Sy, 1, Sy, 1, Save, 1);
299 for (
int fileNo = 0; fileNo < nfiles ; ++fileNo)
311 for(j = 0; j < nfields; ++j)
313 for(
int i = 0; i < fielddata.size(); ++i)
315 Exp[j]->ExtractDataToCoeffs(fielddef[i],
317 fielddef[i]->m_fields[j],
318 Exp[j]->UpdateCoeffs());
320 Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
323 Exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
324 BndExp[0] = Exp[0]->GetBndCondExpansions();
326 for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
330 for(i = 0; i < nelem; ++i, cnt++)
333 elmtid = BoundarytoElmtID[cnt];
334 elmt = Exp[0]->GetExp(elmtid);
335 offset = Exp[0]->GetPhys_Offset(elmtid);
347 boundary = BoundarytoTraceID[cnt];
350 for (
int t = 0; t< sfields; t++)
352 elmt->GetFacePhysVals(boundary,bc,Exp[t+expdim]->GetPhys() + offset, temp[t]);
355 for (
int m = 0; m < nfq; m++)
357 Sx[bndOffset + m] = temp[0][m];
358 Sy[bndOffset + m] = temp[1][m];
359 Sz[bndOffset + m] = temp[2][m];
365 Vmath::Vvtvvtp(nbq, Sx, 1, Sx, 1, Sy, 1, Sy, 1, S, 1);
371 Vmath::Vvtvvtp(nbq, Sx, 1, Sxr, 1, Sy, 1, Syr, 1, temp2, 1);
375 for (
int m = 0; m < nbq; m++)
379 trs[m] = trs[m] + sqrt(temp2[m]);
394 cnt += BndExp[0][n]->GetExpSize();
402 Vmath::Smul (nbq, (1.0/nfiles), TAwss, 1, TAwss, 1);
414 for(j = 0; j < nfields; ++j)
416 for(
int i = 0; i < fielddata.size(); ++i)
418 Exp[j]->ExtractDataToCoeffs(fielddef[i],
420 fielddef[i]->m_fields[j],
421 Exp[j]->UpdateCoeffs());
423 Exp[j]->BwdTrans(Exp[j]->GetCoeffs(),Exp[j]->UpdatePhys());
426 for(j = 0; j < addfields; ++j)
428 for(
int i = 0; i < fielddata.size(); ++i)
430 Exp[j+nfields]->ExtractDataToCoeffs(fielddef[i],
432 fielddef[i]->m_fields[j],
433 Exp[j+nfields]->UpdateCoeffs());
435 Exp[j+nfields]->BwdTrans(Exp[j+nfields]->GetCoeffs(),Exp[j+nfields]->UpdatePhys());
438 Exp[0]->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
441 for(j = 0; j < nfields+addfields; ++j)
443 BndExp[j] = Exp[j]->GetBndCondExpansions();
446 for(cnt = n = 0; n < BndExp[0].num_elements(); ++n)
450 for(i = 0; i < nelem; ++i, cnt++)
453 elmtid = BoundarytoElmtID[cnt];
454 offset = Exp[0]->GetPhys_Offset(elmtid);
466 boundary = BoundarytoTraceID[cnt];
469 for (j = 0; j < nfq; j++)
471 temp[0][j] = trs[bndOffset + j];
472 temp[1][j] = TAwss[bndOffset + j];
473 temp[2][j] = osi[bndOffset + j];
476 for (j = 0; j < 3; j++)
478 values[j] = BndExp[j+nfields][n]->UpdateCoeffs() + BndExp[j+nfields][n]->GetCoeff_Offset(i);
479 bc->FwdTrans(temp[j], values[j]);
484 for (j = 0; j < nfq; j++)
486 temp[0][j] = Sxr[bndOffset + j];
487 temp[1][j] = Syr[bndOffset + j];
488 temp[2][j] = Szr[bndOffset + j];
491 for (j = 0; j < 3; j++)
493 values[j] = BndExp[j+nfields+3][n]->UpdateCoeffs() + BndExp[j+nfields+3][n]->GetCoeff_Offset(i);
494 bc->FwdTrans(temp[j], values[j]);
500 for (j = 0; j < nfq; j++)
502 temp[0][j] = Save[bndOffset + j];
505 values[0] = BndExp[nfields+addfields-1][n]->UpdateCoeffs() + BndExp[nfields+addfields-1][n]->GetCoeff_Offset(i);
506 bc->FwdTrans(temp[0], values[0]);
513 cnt += BndExp[0][n]->GetExpSize();
517 for(j = 0; j < nfields + addfields; ++j)
519 int ncoeffs = Exp[j]->GetNcoeffs();
522 output=Exp[j]->UpdateCoeffs();
524 int nGlobal=m_locToGlobalMap->GetNumGlobalCoeffs();
532 for(i = 0; i < BndExp[j].num_elements(); ++i)
537 for(
int k = 0; k < (BndExp[j][i])->GetNcoeffs(); ++k)
539 sign = m_locToGlobalMap->GetBndCondCoeffsToGlobalCoeffsSign(bndcnt);
540 outarray[map[bndcnt++]] = sign * coeffs[k];
545 bndcnt += BndExp[j][i]->GetNcoeffs();
548 m_locToGlobalMap->GlobalToLocal(outarray,output);
555 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
556 = Exp[0]->GetFieldDefinitions();
557 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
559 vector<string > outname;
561 outname.push_back(
"TransWSS");
562 outname.push_back(
"TAWSS");
563 outname.push_back(
"OSI");
564 outname.push_back(
"norm_mean_x");
565 outname.push_back(
"norm_mean_y");
566 outname.push_back(
"norm_mean_z");
567 outname.push_back(
"mean_mag");
569 for(j = 0; j < nfields+addfields; ++j)
571 for(i = 0; i < FieldDef.size(); ++i)
575 FieldDef[i]->m_fields.push_back(outname[j-nfields]);
579 FieldDef[i]->m_fields.push_back(fielddef[i]->m_fields[j]);
581 Exp[j]->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 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
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
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 Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
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 ...
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 Vvtvm(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)
vvtvm (vector times vector plus vector): z = w*x - y
boost::shared_ptr< ContField3D > ContField3DSharedPtr
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
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.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
int main(int argc, char *argv[])