38 #include <boost/algorithm/string.hpp>
63 ParamMap::const_iterator it;
66 it = pParams.find(
"OutputFile");
67 if (it == pParams.end())
73 ASSERTL0(it->second.length() > 0,
"Missing parameter 'OutputFile'.");
83 it = pParams.find(
"OutputFrequency");
84 if (it == pParams.end())
95 it = pParams.find(
"StartTime");
96 if (it == pParams.end())
107 m_session->MatchSolverInfo(
"Homogeneous",
"1D",
111 it = pParams.find(
"OutputAllPlanes");
112 if (it == pParams.end())
118 std::string sOption =
121 ( boost::iequals(sOption,
"yes"));
130 it = pParams.find(
"Boundary");
131 ASSERTL0(it != pParams.end(),
"Missing parameter 'Boundary");
132 ASSERTL0(it->second.length() > 0,
"Empty parameter 'Boundary'.");
142 for (
int i = 0; i < 3; ++i)
147 std::stringstream directionStream;
148 std::string directionString;
150 for (
int i = 0; i < 3; ++i)
152 std::stringstream tmp;
154 std::string dir =
"Direction" + tmp.str();
155 it = pParams.find(dir);
156 if ( it != pParams.end() )
159 "Missing parameter '"+dir+
"'.");
160 directionStream.str(it->second);
163 directionStream.clear();
166 for (
int j = 0; j < 3; j++)
168 directionStream >> directionString;
169 if (!directionString.empty())
177 for(
int j = 0; j < 3; j++)
207 std::string::size_type FirstInd =
209 std::string::size_type LastInd =
213 (std::string(
"Error reading boundary region definition:") +
216 std::string IndString =
221 (std::string(
"Unable to read boundary regions index "
222 "range for FilterAeroForces: ") + IndString).c_str());
226 unsigned int numBoundaryRegions =
227 pFields[0]->GetBndConditions().num_elements();
229 numBoundaryRegions, 0);
232 pFields[0]->GetGraph());
235 SpatialDomains::BoundaryRegionCollection::const_iterator it;
237 for (cnt = 0, it = bregions.begin(); it != bregions.end();
251 pFields[0]->GetPlane(0)->GetBoundaryToElmtMap
265 m_nPlanes = pFields[0]->GetHomogeneousBasis()->
266 GetZ().num_elements();
283 for(j = 0; j < IDs.num_elements(); ++j)
293 if(j != IDs.num_elements())
307 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
308 if (vComm->GetRank() == 0)
313 for(
int i = 0; i < expdim; i++ )
324 m_outputStream <<
"# Boundary regions: " << IndString.c_str() << endl;
328 for(
int i = 1; i <= expdim; i++ )
342 if (
m_session->DefinesSolverInfo(
"HomoStrip"))
345 "Output forces on all planes not compatible with HomoStrips");
374 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
378 for(
int i = 0; i < expdim; i++)
382 Ft[i] = Fp[i] + Fv[i];
389 if (vComm->GetRank() == 0)
394 for(
int plane = 0; plane <
m_nPlanes; plane++)
400 for(
int i = 0; i < expdim; i++ )
420 for(
int i = 0; i < expdim; i++)
435 if(
m_session->DefinesSolverInfo(
"HomoStrip"))
444 m_session->LoadParameter(
"Strip_Z", nstrips);
445 for(
int i = 1; i<nstrips; i++)
447 vComm->GetColumnComm()->Recv(i, Fp);
448 vComm->GetColumnComm()->Recv(i, Fv);
449 vComm->GetColumnComm()->Recv(i, Ft);
453 for(
int j = 0; j < expdim; j++)
475 if (
m_session->DefinesSolverInfo(
"HomoStrip") &&
476 (vComm->GetRowComm()->GetRank() == 0) )
478 vComm->GetColumnComm()->Send(0, Fp);
479 vComm->GetColumnComm()->Send(0, Fv);
480 vComm->GetColumnComm()->Send(0, Ft);
494 if (pFields[0]->GetComm()->GetRank() == 0)
525 int local_planes = ZIDs.num_elements();
526 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
531 for(
int plane = 0 ; plane < local_planes; plane++)
533 for(
int dir =0; dir < expdim; dir++)
535 Aeroforces[plane + dir*local_planes] =
542 for(
int plane = 0 ; plane < local_planes; plane++)
544 for(
int dir =0; dir < expdim; dir++)
546 Aeroforces[plane + dir*local_planes] =
573 int i, j, k, n, cnt, elmtid, nq, offset, boundary, plane;
575 int physTot = pFields[0]->GetNpoints();
576 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
587 fields( pFields.num_elements() );
604 m_session->DefinesSolverInfo(
"HomoStrip") ?
605 vComm->GetColumnComm()->GetColumnComm():
606 vComm->GetColumnComm();
612 for( i = 0; i < expdim; i++)
631 for(i = 0; i < pFields.num_elements(); ++i)
635 pFields[i]->SetWaveSpace(
false);
637 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
638 pFields[i]->UpdatePhys());
639 pFields[i]->SetPhysState(
true);
647 BndConds = pFields[0]->GetPlane(0)->GetBndConditions();
648 BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
652 BndConds = pFields[0]->GetBndConditions();
653 BndExp = pFields[0]->GetBndCondExpansions();
658 for(plane = 0; plane <
m_nPlanes; plane++ )
667 for(n = 0; n < pFields.num_elements(); n++)
669 fields[n] = pFields[n]->GetPlane(
m_planesID[plane]);
674 for(n = 0; n < pFields.num_elements(); n++)
676 fields[n] = pFields[n];
681 for( cnt = n = 0; n < BndConds.num_elements(); n++)
685 for (i=0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
688 elmt = fields[0]->GetExp(elmtid);
689 nq = elmt->GetTotPoints();
690 offset = fields[0]->GetPhys_Offset(elmtid);
693 for( j=0; j<expdim; j++)
695 velocity[j] = fields[j]->GetPhys() + offset;
697 P = fields[nVel]->GetPhys() + offset;
700 for (j=0; j<expdim; j++)
702 for (k=0; k<expdim; k++)
706 elmt->PhysDeriv(k,velocity[j],
724 bc = BndExp[n]->GetExp(i)->
725 as<LocalRegions::Expansion1D> ();
728 nbc = bc->GetTotPoints();
731 normals = elmt->GetEdgeNormal(boundary);
735 elmt->GetEdgePhysVals(boundary,bc,P,Pb);
736 for(
int j = 0; j < expdim*expdim; ++j)
740 elmt->GetEdgePhysVals(boundary,
741 bc,grad[j],gradb[j]);
750 bc = BndExp[n]->GetExp(i)->
751 as<LocalRegions::Expansion2D> ();
754 nbc = bc->GetTotPoints();
757 normals = elmt->GetFaceNormal(boundary);
761 elmt->GetFacePhysVals(boundary,bc,P,Pb);
762 for(
int j = 0; j < expdim*expdim; ++j)
766 elmt->GetFacePhysVals(boundary,
767 bc,grad[j],gradb[j]);
774 "Expansion not supported by FilterForces");
781 for ( j = 0; j < expdim; j++)
791 for ( j = 0; j < expdim; j++ )
794 for ( k = 0; k < expdim; k++ )
809 for ( j = 0; j < expdim; j++)
811 m_Fpplane[j][plane] += BndExp[n]->GetExp(i)->
813 m_Fvplane[j][plane] += BndExp[n]->GetExp(i)->
820 cnt += BndExp[n]->GetExpSize();
829 for( i = 0; i < expdim; i++)
839 for(plane = 0; plane <
m_nPlanes; plane++)
843 for( i = 0; i < expdim; i++)
845 for( j = 0; j < expdim; j++ )
848 tmpV[i] += m_Fvplane[j][plane]*m_directions[i][j];
852 for( i = 0; i < expdim; i++)
855 m_Fvplane[i][plane] = tmpV[i];
860 for(plane = 0; plane <
m_nPlanes; plane++)
862 for( i = 0; i < expdim; i++)
864 m_Ftplane[i][plane] =
m_Fpplane[i][plane] + m_Fvplane[i][plane];
871 for (i = 0; i < pFields.num_elements(); ++i)
873 pFields[i]->SetWaveSpace(
true);
874 pFields[i]->HomogeneousFwdTrans(pFields[i]->GetPhys(),
875 pFields[i]->UpdatePhys());
888 int i, j, k, n, cnt, elmtid, offset, boundary, plane;
890 int physTot = pFields[0]->GetNpoints();
891 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
928 m_session->DefinesSolverInfo(
"HomoStrip") ?
929 vComm->GetColumnComm()->GetColumnComm():
930 vComm->GetColumnComm();
936 for( i = 0; i < expdim; i++)
955 for(i = 0; i < pFields.num_elements(); ++i)
959 pFields[i]->SetWaveSpace(
false);
961 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
962 pFields[i]->UpdatePhys());
963 pFields[i]->SetPhysState(
true);
971 BndConds = pFields[0]->GetPlane(0)->GetBndConditions();
972 BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
976 BndConds = pFields[0]->GetBndConditions();
977 BndExp = pFields[0]->GetBndCondExpansions();
992 Exp3DH1 = boost::dynamic_pointer_cast
995 for(i = 0; i < nVel*nVel; i++)
1012 Exp2D = boost::dynamic_pointer_cast
1015 for(i = 0; i < nVel*nVel; i++)
1034 Exp3D = boost::dynamic_pointer_cast
1037 for(i = 0; i < nVel*nVel; i++)
1054 ASSERTL0(
false,
"Expansion dimension not supported by FilterAeroForces");
1064 for (i = 0; i < nVel*nVel; i++)
1067 pFields[nVel]->GetPhys(), 1,
1068 P[i]->UpdatePhys(), 1);
1073 for (i=0; i<nVel; i++)
1079 m_mapping->ApplyChristoffelContravar(wk, tmp);
1080 for (i=0; i< nVel; i++)
1082 for (k=0; k< nVel; k++)
1085 wk[i], grad[i*nVel+k]->UpdatePhys());
1088 grad[i*nVel+k]->UpdatePhys(), 1,
1089 grad[i*nVel+k]->UpdatePhys(), 1);
1093 for (i=0; i< nVel; i++)
1095 for (k=0; k< nVel; k++)
1101 for (j=0; j<nVel; j++)
1104 grad[i*nVel+j]->UpdatePhys(), 1);
1109 m_mapping->GetJacobian( Jac->UpdatePhys());
1112 for (i=0; i< nVel; i++)
1115 for (k=0; k< nVel; k++)
1124 for (k=0; k< nVel; k++)
1127 C[k*nVel+i]->UpdatePhys(), 1);
1137 for(plane = 0; plane <
m_nPlanes; plane++ )
1146 for(n = 0; n < nVel*nVel; n++)
1148 PPlane[n] = P[n]->GetPlane(
m_planesID[plane]);
1149 gradPlane[n] = grad[n]->GetPlane(
m_planesID[plane]);
1150 CPlane[n] = C[n]->GetPlane(
m_planesID[plane]);
1156 for(n = 0; n < nVel*nVel; n++)
1159 gradPlane[n] = grad[n];
1166 for( cnt = n = 0; n < BndConds.num_elements(); n++)
1170 for (i=0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
1173 elmt = PPlane[0]->GetExp(elmtid);
1174 offset = PPlane[0]->GetPhys_Offset(elmtid);
1177 for( j=0; j<nVel*nVel; j++)
1179 PElmt[j] = PPlane[j]->GetPhys()
1181 gradElmt[j] = gradPlane[j]->GetPhys()
1183 CElmt[j] = CPlane[j]->GetPhys()
1186 JacElmt = JacPlane->GetPhys() + offset;
1201 bc = BndExp[n]->GetExp(i)->
1202 as<LocalRegions::Expansion1D> ();
1205 nbc = bc->GetTotPoints();
1208 normals = elmt->GetEdgeNormal(boundary);
1211 for(
int j = 0; j < nVel*nVel; ++j)
1215 elmt->GetEdgePhysVals(boundary,
1221 elmt->GetEdgePhysVals(boundary,
1226 elmt->GetEdgePhysVals(boundary,
1232 elmt->GetEdgePhysVals(boundary,
1242 bc = BndExp[n]->GetExp(i)->
1243 as<LocalRegions::Expansion2D> ();
1246 nbc = bc->GetTotPoints();
1249 normals = elmt->GetFaceNormal(boundary);
1252 for(
int j = 0; j < nVel*nVel; ++j)
1256 elmt->GetFacePhysVals(boundary,
1262 elmt->GetFacePhysVals(boundary,
1268 elmt->GetFacePhysVals(boundary,
1274 elmt->GetFacePhysVals(boundary,
1282 "Expansion not supported by FilterForces");
1289 for ( j = 0; j < nVel; j++)
1293 for ( k = 0; k < expdim; k++)
1304 for ( j = 0; j < nVel; j++ )
1307 for ( k = 0; k < expdim; k++ )
1322 for ( k = 0; k < nVel; k++ )
1331 for ( k = 0; k < nVel; k++ )
1334 for ( j = 0; j < nVel; j++ )
1342 for ( k = 0; k < nVel; k++ )
1347 for ( k = 0; k < nVel; k++ )
1350 for ( j = 0; j < nVel; j++ )
1358 for ( k = 0; k < nVel; k++ )
1364 for ( j = 0; j < expdim; j++)
1366 m_Fpplane[j][plane] += BndExp[n]->GetExp(i)->
1368 m_Fvplane[j][plane] += BndExp[n]->GetExp(i)->
1375 cnt += BndExp[n]->GetExpSize();
1384 for( i = 0; i < expdim; i++)
1394 for(plane = 0; plane <
m_nPlanes; plane++)
1398 for( i = 0; i < expdim; i++)
1400 for( j = 0; j < expdim; j++ )
1403 tmpV[i] += m_Fvplane[j][plane]*m_directions[i][j];
1407 for( i = 0; i < expdim; i++)
1410 m_Fvplane[i][plane] = tmpV[i];
1415 for(plane = 0; plane <
m_nPlanes; plane++)
1417 for( i = 0; i < expdim; i++)
1419 m_Ftplane[i][plane] =
m_Fpplane[i][plane] + m_Fvplane[i][plane];
1426 for (i = 0; i < pFields.num_elements(); ++i)
1428 pFields[i]->SetWaveSpace(
true);
1429 pFields[i]->HomogeneousFwdTrans(pFields[i]->GetPhys(),
1430 pFields[i]->UpdatePhys());
std::ofstream m_outputStream
#define ASSERTL0(condition, msg)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual bool v_IsTimeDependent()
Array< OneD, int > m_BCtoTraceID
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.
Array< OneD, int > m_BCtoElmtID
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.
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.
GlobalMapping::MappingSharedPtr m_mapping
Array< OneD, Array< OneD, NekDouble > > m_directions
std::string m_BoundaryString
SOLVER_UTILS_EXPORT void GetForces(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Aeroforces, const NekDouble &time)
NekDouble Evaluate() const
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
void CalculateForces(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
std::map< std::string, std::string > ParamMap
LibUtilities::SessionReaderSharedPtr m_session
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
int m_nPlanes
number of planes for homogeneous1D expansion
boost::shared_ptr< ExpList3D > ExpList3DSharedPtr
Shared pointer to an ExpList3D object.
MultiRegions::Direction const DirCartesianMap[]
Array< OneD, Array< OneD, NekDouble > > m_Ftplane
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
unsigned int m_outputFrequency
Array< OneD, Array< OneD, NekDouble > > m_Fpplane
Abstraction of a three-dimensional multi-elemental expansion which is merely a collection of local ex...
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
FilterFactory & GetFilterFactory()
const BoundaryRegionCollection & GetBoundaryRegions(void) const
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Array< OneD, Array< OneD, NekDouble > > m_Fvplane
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
bool m_outputAllPlanes
if using a homogeneous1D expansion, determine if should output all planes or just the average ...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
void CalculateForcesMapping(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
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.
Array< OneD, int > m_planesID