54 const std::map<std::string, std::string> &pParams) :
57 if (pParams.find(
"OutputFile") == pParams.end())
63 ASSERTL0(!(pParams.find(
"OutputFile")->second.empty()),
64 "Missing parameter 'OutputFile'.");
73 if (pParams.find(
"OutputFrequency") == pParams.end())
80 atoi(pParams.find(
"OutputFrequency")->second.c_str());
84 m_session->MatchSolverInfo(
"Homogeneous",
"1D",
89 if (pParams.find(
"OutputPlane") == pParams.end())
96 atoi(pParams.find(
"OutputPlane")->second.c_str());
101 if (pParams.find(
"Boundary") == pParams.end())
103 ASSERTL0(
false,
"Missing parameter 'Boundary'.");
107 ASSERTL0(!(pParams.find(
"Boundary")->second.empty()),
108 "Missing parameter 'Boundary'.");
131 std::string::size_type FirstInd =
133 std::string::size_type LastInd =
137 (std::string(
"Error reading boundary region definition:") +
140 std::string IndString =
145 (std::string(
"Unable to read boundary regions index "
146 "range for FilterAeroForces: ") + IndString).c_str());
150 unsigned int numBoundaryRegions =
151 pFields[0]->GetBndConditions().num_elements();
153 numBoundaryRegions, 0);
156 pFields[0]->GetGraph());
159 SpatialDomains::BoundaryRegionCollection::const_iterator it;
161 for (cnt = 0, it = bregions.begin(); it != bregions.end();
174 if (vComm->GetRank() == 0)
219 int n, cnt, elmtid, nq, offset, nt, boundary;
220 nt = pFields[0]->GetNpoints();
221 int dim = pFields.num_elements()-1;
244 NekDouble Fx,Fy,Fz,Fxp,Fxv,Fyp,Fyv,Fzp,Fzv;
263 for(
int i = 0; i < pFields.num_elements(); ++i)
265 pFields[i]->SetWaveSpace(
false);
266 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
267 pFields[i]->UpdatePhys());
268 pFields[i]->SetPhysState(
true);
275 if(vComm->GetColumnComm()->GetRank() == 0)
277 pFields[0]->GetPlane(0)->GetBoundaryToElmtMap(
278 BoundarytoElmtID,BoundarytoTraceID);
279 BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
283 for(cnt = n = 0; n < BndExp.num_elements(); ++n)
287 for(
int i = 0; i < BndExp[n]->GetExpSize();
291 elmtid = BoundarytoElmtID[cnt];
292 elmt = pFields[0]->GetPlane(0)->GetExp(elmtid);
293 nq = elmt->GetTotPoints();
294 offset = pFields[0]->GetPlane(0)->GetPhys_Offset(elmtid);
299 for(
int j = 0; j < dim; ++j)
307 boundary = BoundarytoTraceID[cnt];
310 U = pFields[0]->GetPlane(0)->GetPhys() + offset;
311 V = pFields[1]->GetPlane(0)->GetPhys() + offset;
312 P = pFields[3]->GetPlane(0)->GetPhys() + offset;
315 elmt->PhysDeriv(U,gradU[0],gradU[1]);
316 elmt->PhysDeriv(V,gradV[0],gradV[1]);
319 bc = BndExp[n]->GetExp(i)
328 for(
int j = 0; j < dim; ++j)
342 boundary = BoundarytoTraceID[cnt];
346 elmt->GetEdgePhysVals(boundary,bc,P,Pb);
348 for(
int j = 0; j < dim; ++j)
350 elmt->GetEdgePhysVals(boundary,bc,gradU[j],
352 elmt->GetEdgePhysVals(boundary,bc,gradV[j],
358 = elmt->GetEdgeNormal(boundary);
415 Fxv += bc->Integral(drag_t);
416 Fyv += bc->Integral(lift_t);
418 Fxp += bc->Integral(drag_p);
419 Fyp += bc->Integral(lift_p);
424 cnt += BndExp[n]->GetExpSize();
429 for(
int i = 0; i < pFields.num_elements(); ++i)
431 pFields[i]->SetWaveSpace(
true);
432 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
433 pFields[i]->UpdatePhys());
434 pFields[i]->SetPhysState(
false);
440 pFields[0]->GetBoundaryToElmtMap(BoundarytoElmtID,
442 BndExp = pFields[0]->GetBndCondExpansions();
446 for(cnt = n = 0; n < BndExp.num_elements(); ++n)
450 for(
int i = 0; i < BndExp[n]->GetExpSize(); ++i, cnt++)
453 elmtid = BoundarytoElmtID[cnt];
454 elmt = pFields[0]->GetExp(elmtid);
455 nq = elmt->GetTotPoints();
456 offset = pFields[0]->GetPhys_Offset(elmtid);
461 for(
int j = 0; j < dim; ++j)
469 boundary = BoundarytoTraceID[cnt];
472 U = pFields[0]->GetPhys() + offset;
473 V = pFields[1]->GetPhys() + offset;
474 W = pFields[2]->GetPhys() + offset;
475 P = pFields[3]->GetPhys() + offset;
478 elmt->PhysDeriv(U,gradU[0],gradU[1],gradU[2]);
479 elmt->PhysDeriv(V,gradV[0],gradV[1],gradV[2]);
480 elmt->PhysDeriv(W,gradW[0],gradW[1],gradW[2]);
483 bc = BndExp[n]->GetExp(i)
492 for(
int j = 0; j < dim; ++j)
510 boundary = BoundarytoTraceID[cnt];
514 elmt->GetFacePhysVals(boundary,bc,P,Pb);
516 for(
int j = 0; j < dim; ++j)
518 elmt->GetFacePhysVals(boundary,bc,gradU[j],
520 elmt->GetFacePhysVals(boundary,bc,gradV[j],
522 elmt->GetFacePhysVals(boundary,bc,gradW[j],
528 = elmt->GetFaceNormal(boundary);
621 Fxv += bc->Expansion::Integral(drag_t);
622 Fyv += bc->Expansion::Integral(lift_t);
623 Fzv += bc->Expansion::Integral(side_t);
625 Fxp += bc->Expansion::Integral(drag_p);
626 Fyp += bc->Expansion::Integral(lift_p);
627 Fzp += bc->Expansion::Integral(side_p);
632 cnt += BndExp[n]->GetExpSize();
639 pFields[0]->GetBoundaryToElmtMap(BoundarytoElmtID,
641 BndExp = pFields[0]->GetBndCondExpansions();
645 for(cnt = n = 0; n < BndExp.num_elements(); ++n)
649 for(
int i = 0; i < BndExp[n]->GetExpSize(); ++i, cnt++)
652 elmtid = BoundarytoElmtID[cnt];
653 elmt = pFields[0]->GetExp(elmtid);
654 nq = elmt->GetTotPoints();
655 offset = pFields[0]->GetPhys_Offset(elmtid);
657 for(
int j = 0; j < dim; ++j)
663 boundary = BoundarytoTraceID[cnt];
665 U = pFields[0]->GetPhys() + offset;
666 V = pFields[1]->GetPhys() + offset;
667 P = pFields[2]->GetPhys() + offset;
669 elmt->PhysDeriv(U,gradU[0],gradU[1]);
670 elmt->PhysDeriv(V,gradV[0],gradV[1]);
672 bc = BndExp[n]->GetExp(i)
684 boundary = BoundarytoTraceID[cnt];
686 elmt->GetEdgePhysVals(boundary,bc,P,Pb);
688 for(
int j = 0; j < dim; ++j)
695 for(
int j = 0; j < dim; ++j)
697 elmt->GetEdgePhysVals(boundary,bc,gradU[j],
699 elmt->GetEdgePhysVals(boundary,bc,gradV[j],
704 = elmt->GetEdgeNormal(boundary);
728 Fxp += bc->Integral(drag_p);
729 Fyp += bc->Integral(lift_p);
731 Fxv += bc->Integral(drag_t);
732 Fyv += bc->Integral(lift_t);
737 cnt += BndExp[n]->GetExpSize();
757 if (vComm->GetRank() == 0)
795 if (pFields[0]->GetComm()->GetRank() == 0)
std::ofstream m_outputStream
#define ASSERTL0(condition, msg)
virtual bool v_IsTimeDependent()
SOLVER_UTILS_EXPORT FilterAeroForces(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
virtual SOLVER_UTILS_EXPORT ~FilterAeroForces()
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
virtual void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
virtual void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
vector< unsigned int > m_boundaryRegionsIdList
ID's of boundary regions where we want the forces.
unsigned int m_outputPlane
plane to take history point from if using a homogeneous1D expansion
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
static std::string className
Name of the class.
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
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.
std::string m_BoundaryString
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
boost::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
void Neg(int n, T *x, const int incx)
Negate x = -x.
LibUtilities::SessionReaderSharedPtr m_session
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
unsigned int m_outputFrequency
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
FilterFactory & GetFilterFactory()
const BoundaryRegionCollection & GetBoundaryRegions(void) const
void Zero(int n, T *x, const int incx)
Zero vector.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
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.
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.