49 RegisterCreatorFunction(
"MovingBody",
50 FilterMovingBody::create,
51 "Moving Body Filter");
55 FilterMovingBody::FilterMovingBody(
61 ParamMap::const_iterator it;
64 it = pParams.find(
"OutputFile");
65 if (it == pParams.end())
72 ASSERTL0(it->second.length() > 0,
"Missing parameter 'OutputFile'.");
89 it = pParams.find(
"OutputFrequency");
90 if (it == pParams.end())
102 "Homogeneous 1D discetisations.");
105 it = pParams.find(
"Boundary");
106 ASSERTL0(it != pParams.end(),
"Missing parameter 'Boundary'.");
107 ASSERTL0(it->second.length() > 0,
"Missing parameter 'Boundary'.");
136 (std::string(
"Error reading boundary region definition:") +
140 LastInd - FirstInd + 1);
146 (std::string(
"Unable to read boundary regions index "
147 "range for FilterAeroForces: ") + IndString).c_str());
152 unsigned int numBoundaryRegions
153 = pFields[0]->GetBndConditions().num_elements();
156 numBoundaryRegions, 0);
163 SpatialDomains::BoundaryRegionCollection::const_iterator it;
165 for (cnt = 0, it = bregions.begin(); it != bregions.end();
178 if (vComm->GetRank() == 0)
234 int n, cnt, elmtid, nq, offset, boundary;
235 int nt = pFields[0]->GetNpoints();
236 int dim = pFields.num_elements()-1;
264 int Num_z_pos = pFields[0]->GetHomogeneousBasis()->GetNumModes();
273 NekDouble rho = (pSession->DefinesParameter(
"rho"))
274 ? (pSession->GetParameter(
"rho"))
276 NekDouble mu = rho*pSession->GetParameter(
"Kinvis");
278 for(
int i = 0; i < pFields.num_elements(); ++i)
280 pFields[i]->SetWaveSpace(
false);
281 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
282 pFields[i]->UpdatePhys());
283 pFields[i]->SetPhysState(
true);
289 ZIDs = pFields[0]->GetZIDs();
290 int local_planes = ZIDs.num_elements();
294 for(
int plane = 0 ; plane < local_planes; plane++)
296 pFields[0]->GetPlane(plane)->GetBoundaryToElmtMap(BoundarytoElmtID,
298 BndExp = pFields[0]->GetPlane(plane)->GetBndCondExpansions();
302 for(cnt = n = 0; n < BndExp.num_elements(); ++n)
306 for(
int i = 0; i < BndExp[n]->GetExpSize(); ++i, cnt++)
309 elmtid = BoundarytoElmtID[cnt];
310 elmt = pFields[0]->GetPlane(plane)->GetExp(elmtid);
311 nq = elmt->GetTotPoints();
312 offset = pFields[0]->GetPlane(plane)
313 ->GetPhys_Offset(elmtid);
318 for(
int j = 0; j < dim; ++j)
326 boundary = BoundarytoTraceID[cnt];
329 U = pFields[0]->GetPlane(plane)->GetPhys() + offset;
330 V = pFields[1]->GetPlane(plane)->GetPhys() + offset;
331 P = pFields[3]->GetPlane(plane)->GetPhys() + offset;
334 elmt->PhysDeriv(U,gradU[0],gradU[1]);
335 elmt->PhysDeriv(V,gradV[0],gradV[1]);
346 for(
int j = 0; j < dim; ++j)
360 boundary = BoundarytoTraceID[cnt];
364 elmt->GetEdgePhysVals(boundary,bc,P,Pb);
366 for(
int j = 0; j < dim; ++j)
368 elmt->GetEdgePhysVals(boundary,bc,gradU[j],fgradU[j]);
369 elmt->GetEdgePhysVals(boundary,bc,gradV[j],fgradV[j]);
374 = elmt->GetEdgeNormal(boundary);
393 Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, drag_t, 1);
394 Vmath::Vmul(nbc, drag_t, 1, normals[1], 1, drag_t, 1);
397 Vmath::Vmul(nbc, fgradU[0], 1, normals[0], 1, temp2, 1);
411 Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, lift_t, 1);
412 Vmath::Vmul(nbc, lift_t, 1, normals[0], 1, lift_t, 1);
415 Vmath::Vmul(nbc, fgradV[1], 1, normals[1], 1, temp2, 1);
416 Vmath::Smul(nbc, -0.5, fgradV[1], 1, fgradV[1], 1);
431 Fxv[ZIDs[plane]] += bc->Integral(drag_t);
432 Fyv[ZIDs[plane]] += bc->Integral(lift_t);
434 Fxp[ZIDs[plane]] += bc->Integral(drag_p);
435 Fyp[ZIDs[plane]] += bc->Integral(lift_p);
440 cnt += BndExp[n]->GetExpSize();
445 for(
int i = 0; i < pFields.num_elements(); ++i)
447 pFields[i]->SetWaveSpace(
true);
448 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
449 pFields[i]->UpdatePhys());
450 pFields[i]->SetPhysState(
false);
459 if(vComm->GetRowComm()->GetSize() > 0)
466 for(
int plane = 0 ; plane < local_planes; plane++)
481 if(!pSession->DefinesSolverInfo(
"HomoStrip"))
483 if(vComm->GetRowComm()->GetRank() == 0)
485 for(
int z = 0 ; z < Num_z_pos; z++)
496 if(vComm->GetRowComm()->GetRank() == 0)
498 for(
int z = 0 ; z < Num_z_pos; z++)
508 if(!pSession->DefinesSolverInfo(
"HomoStrip"))
511 for(
int plane = 0 ; plane < local_planes; plane++)
513 Aeroforces[plane] = Fxp[ZIDs[plane]]
515 Aeroforces[plane + local_planes] = Fyp[ZIDs[plane]]
532 = pFields[0]->GetHomogeneousBasis()->GetZ();
535 pSession->LoadParameter(
"LZ", LZ);
537 Vmath::Sadd(Num_z_pos,LZ/2.0,z_coords,1,z_coords,1);
538 if (vComm->GetRank() == 0)
544 for(
int i = 0 ; i < Num_z_pos; i++)
573 for(
int plane = 0 ; plane < local_planes; plane++)
575 fces[0] += Fxp[ZIDs[plane]] + Fxv[ZIDs[plane]];
576 fces[1] += Fyp[ZIDs[plane]] + Fyv[ZIDs[plane]];
577 fces[2] += Fxp[ZIDs[plane]] ;
578 fces[3] += Fyp[ZIDs[plane]] ;
579 fces[4] += Fxv[ZIDs[plane]] ;
580 fces[5] += Fyv[ZIDs[plane]] ;
583 fces[0] = fces[0]/local_planes;
584 fces[1] = fces[1]/local_planes;
585 fces[2] = fces[2]/local_planes;
586 fces[3] = fces[3]/local_planes;
587 fces[4] = fces[4]/local_planes;
588 fces[5] = fces[5]/local_planes;
598 int npts = vComm->GetColumnComm()->GetColumnComm()->GetSize();
600 fces[0] = fces[0]/
npts;
601 fces[1] = fces[1]/
npts;
602 fces[2] = fces[2]/
npts;
603 fces[3] = fces[3]/
npts;
604 fces[4] = fces[4]/
npts;
605 fces[5] = fces[5]/
npts;
607 for(
int plane = 0 ; plane < local_planes; plane++)
609 Aeroforces[plane] = fces[0];
610 Aeroforces[plane+local_planes] = fces[1];
619 int colrank = vColComm->GetRank();
624 pSession->LoadParameter(
"Strip_Z", nstrips);
625 pSession->LoadParameter(
"DistStrip", DistStrip);
628 for(
int i = 0; i < nstrips; i++)
630 z_coords[i] = i * DistStrip;
656 for(
int i = 1; i < nstrips; i++)
658 vColComm->Recv(i, fces);
684 for(
int i = 1; i < nstrips; i++)
688 vColComm->Send(0, fces);
715 if(!pSession->DefinesSolverInfo(
"HomoStrip"))
717 pSession->LoadParameter(
"LZ", Length);
718 npts =
m_session->GetParameter(
"HomModesZ");
722 pSession->LoadParameter(
"LC", Length);
723 npts =
m_session->GetParameter(
"HomStructModesZ");
727 for(
int n = 0; n <
npts; n++)
729 z_coords = Length/npts*n;
758 if (pFields[0]->GetComm()->GetRank() == 0)
#define ASSERTL0(condition, msg)
std::string m_BoundaryString
Array< OneD, std::ofstream > m_outputStream
virtual bool v_IsTimeDependent()
std::vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
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
std::string m_outputFile_mot
LibUtilities::SessionReaderSharedPtr m_session
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
virtual void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
void UpdateMotion(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &MotionVars, const NekDouble &time)
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
unsigned int m_outputFrequency
NekDouble Evaluate() const
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
std::string m_outputFile_fce
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
std::map< std::string, std::string > ParamMap
void UpdateForce(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Aeroforces, const NekDouble &time)
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
FilterFactory & GetFilterFactory()
const BoundaryRegionCollection & GetBoundaryRegions(void) const
std::vector< unsigned int > m_boundaryRegionsIdList
ID's of boundary regions where we want the forces.
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
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.