38 #include <boost/algorithm/string.hpp>
53 std::string FilterAeroForces::className =
55 "AeroForces", FilterAeroForces::create);
60 FilterAeroForces::FilterAeroForces(
65 ParamMap::const_iterator it;
68 it = pParams.find(
"OutputFile");
69 if (it == pParams.end())
75 ASSERTL0(it->second.length() > 0,
"Missing parameter 'OutputFile'.");
85 it = pParams.find(
"OutputFrequency");
86 if (it == pParams.end())
97 it = pParams.find(
"StartTime");
98 if (it == pParams.end())
109 m_session->MatchSolverInfo(
"Homogeneous",
"1D",
113 it = pParams.find(
"OutputAllPlanes");
114 if (it == pParams.end())
120 std::string sOption =
123 ( boost::iequals(sOption,
"yes"));
132 it = pParams.find(
"Boundary");
133 ASSERTL0(it != pParams.end(),
"Missing parameter 'Boundary");
134 ASSERTL0(it->second.length() > 0,
"Empty parameter 'Boundary'.");
144 for (
int i = 0; i < 3; ++i)
149 std::stringstream directionStream;
150 std::string directionString;
152 for (
int i = 0; i < 3; ++i)
154 std::stringstream tmp;
156 std::string dir =
"Direction" + tmp.str();
157 it = pParams.find(dir);
158 if ( it != pParams.end() )
161 "Missing parameter '"+dir+
"'.");
162 directionStream.str(it->second);
165 directionStream.clear();
168 for (
int j = 0; j < 3; j++)
170 directionStream >> directionString;
171 if (!directionString.empty())
179 for(
int j = 0; j < 3; j++)
209 std::string::size_type FirstInd =
211 std::string::size_type LastInd =
215 (std::string(
"Error reading boundary region definition:") +
218 std::string IndString =
223 (std::string(
"Unable to read boundary regions index "
224 "range for FilterAeroForces: ") + IndString).c_str());
228 unsigned int numBoundaryRegions =
229 pFields[0]->GetBndConditions().num_elements();
231 numBoundaryRegions, 0);
234 pFields[0]->GetGraph());
237 SpatialDomains::BoundaryRegionCollection::const_iterator it;
239 for (cnt = 0, it = bregions.begin(); it != bregions.end();
253 pFields[0]->GetPlane(0)->GetBoundaryToElmtMap
267 m_nPlanes = pFields[0]->GetHomogeneousBasis()->
268 GetZ().num_elements();
285 for(j = 0; j < IDs.num_elements(); ++j)
295 if(j != IDs.num_elements())
309 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
310 if (vComm->GetRank() == 0)
314 m_session->MatchSolverInfo(
"Driver",
"Adaptive",
325 for(
int i = 0; i < expdim; i++ )
336 m_outputStream <<
"# Boundary regions: " << IndString.c_str() << endl;
340 for(
int i = 1; i <= expdim; i++ )
354 if (
m_session->DefinesSolverInfo(
"HomoStrip"))
357 "Output forces on all planes not compatible with HomoStrips");
386 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
390 for(
int i = 0; i < expdim; i++)
394 Ft[i] = Fp[i] + Fv[i];
401 if (vComm->GetRank() == 0)
406 for(
int plane = 0; plane <
m_nPlanes; plane++)
412 for(
int i = 0; i < expdim; i++ )
432 for(
int i = 0; i < expdim; i++)
447 if(
m_session->DefinesSolverInfo(
"HomoStrip"))
456 m_session->LoadParameter(
"Strip_Z", nstrips);
457 for(
int i = 1; i<nstrips; i++)
459 vComm->GetColumnComm()->Recv(i, Fp);
460 vComm->GetColumnComm()->Recv(i, Fv);
461 vComm->GetColumnComm()->Recv(i, Ft);
465 for(
int j = 0; j < expdim; j++)
487 if (
m_session->DefinesSolverInfo(
"HomoStrip") &&
488 (vComm->GetRowComm()->GetRank() == 0) )
490 vComm->GetColumnComm()->Send(0, Fp);
491 vComm->GetColumnComm()->Send(0, Fv);
492 vComm->GetColumnComm()->Send(0, Ft);
506 if (pFields[0]->GetComm()->GetRank() == 0)
537 int local_planes = ZIDs.num_elements();
538 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
543 for(
int plane = 0 ; plane < local_planes; plane++)
545 for(
int dir =0; dir < expdim; dir++)
547 Aeroforces[plane + dir*local_planes] =
554 for(
int plane = 0 ; plane < local_planes; plane++)
556 for(
int dir =0; dir < expdim; dir++)
558 Aeroforces[plane + dir*local_planes] =
585 int i, j, k, n, cnt, elmtid, nq, offset, boundary, plane;
587 int physTot = pFields[0]->GetNpoints();
588 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
599 fields( pFields.num_elements() );
616 m_session->DefinesSolverInfo(
"HomoStrip") ?
617 vComm->GetColumnComm()->GetColumnComm():
618 vComm->GetColumnComm();
624 for( i = 0; i < expdim; i++)
643 for(i = 0; i < pFields.num_elements(); ++i)
647 pFields[i]->SetWaveSpace(
false);
649 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
650 pFields[i]->UpdatePhys());
651 pFields[i]->SetPhysState(
true);
659 BndConds = pFields[0]->GetPlane(0)->GetBndConditions();
660 BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
664 BndConds = pFields[0]->GetBndConditions();
665 BndExp = pFields[0]->GetBndCondExpansions();
670 for(plane = 0; plane <
m_nPlanes; plane++ )
679 for(n = 0; n < pFields.num_elements(); n++)
681 fields[n] = pFields[n]->GetPlane(
m_planesID[plane]);
686 for(n = 0; n < pFields.num_elements(); n++)
688 fields[n] = pFields[n];
693 for( cnt = n = 0; n < BndConds.num_elements(); n++)
697 for (i=0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
700 elmt = fields[0]->GetExp(elmtid);
701 nq = elmt->GetTotPoints();
702 offset = fields[0]->GetPhys_Offset(elmtid);
705 for( j=0; j<expdim; j++)
707 velocity[j] = fields[j]->GetPhys() + offset;
709 P = fields[nVel]->GetPhys() + offset;
712 for (j=0; j<expdim; j++)
714 for (k=0; k<expdim; k++)
718 elmt->PhysDeriv(k,velocity[j],
736 bc = BndExp[n]->GetExp(i)->
737 as<LocalRegions::Expansion1D> ();
740 nbc = bc->GetTotPoints();
743 normals = elmt->GetEdgeNormal(boundary);
747 elmt->GetEdgePhysVals(boundary,bc,P,Pb);
748 for(
int j = 0; j < expdim*expdim; ++j)
752 elmt->GetEdgePhysVals(boundary,
753 bc,grad[j],gradb[j]);
762 bc = BndExp[n]->GetExp(i)->
763 as<LocalRegions::Expansion2D> ();
766 nbc = bc->GetTotPoints();
769 normals = elmt->GetFaceNormal(boundary);
773 elmt->GetFacePhysVals(boundary,bc,P,Pb);
774 for(
int j = 0; j < expdim*expdim; ++j)
778 elmt->GetFacePhysVals(boundary,
779 bc,grad[j],gradb[j]);
786 "Expansion not supported by FilterForces");
793 for ( j = 0; j < expdim; j++)
803 for ( j = 0; j < expdim; j++ )
806 for ( k = 0; k < expdim; k++ )
821 for ( j = 0; j < expdim; j++)
823 m_Fpplane[j][plane] += BndExp[n]->GetExp(i)->
825 m_Fvplane[j][plane] += BndExp[n]->GetExp(i)->
832 cnt += BndExp[n]->GetExpSize();
841 for( i = 0; i < expdim; i++)
851 for(plane = 0; plane <
m_nPlanes; plane++)
855 for( i = 0; i < expdim; i++)
857 for( j = 0; j < expdim; j++ )
860 tmpV[i] += m_Fvplane[j][plane]*m_directions[i][j];
864 for( i = 0; i < expdim; i++)
867 m_Fvplane[i][plane] = tmpV[i];
872 for(plane = 0; plane <
m_nPlanes; plane++)
874 for( i = 0; i < expdim; i++)
876 m_Ftplane[i][plane] =
m_Fpplane[i][plane] + m_Fvplane[i][plane];
883 for (i = 0; i < pFields.num_elements(); ++i)
885 pFields[i]->SetWaveSpace(
true);
886 pFields[i]->HomogeneousFwdTrans(pFields[i]->GetPhys(),
887 pFields[i]->UpdatePhys());
900 int i, j, k, n, cnt, elmtid, offset, boundary, plane;
902 int physTot = pFields[0]->GetNpoints();
903 int expdim = pFields[0]->GetGraph()->GetMeshDimension();
940 m_session->DefinesSolverInfo(
"HomoStrip") ?
941 vComm->GetColumnComm()->GetColumnComm():
942 vComm->GetColumnComm();
948 for( i = 0; i < expdim; i++)
967 for(i = 0; i < pFields.num_elements(); ++i)
971 pFields[i]->SetWaveSpace(
false);
973 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
974 pFields[i]->UpdatePhys());
975 pFields[i]->SetPhysState(
true);
983 BndConds = pFields[0]->GetPlane(0)->GetBndConditions();
984 BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
988 BndConds = pFields[0]->GetBndConditions();
989 BndExp = pFields[0]->GetBndCondExpansions();
1004 Exp3DH1 = boost::dynamic_pointer_cast
1007 for(i = 0; i < nVel*nVel; i++)
1024 Exp2D = boost::dynamic_pointer_cast
1027 for(i = 0; i < nVel*nVel; i++)
1046 Exp3D = boost::dynamic_pointer_cast
1049 for(i = 0; i < nVel*nVel; i++)
1066 ASSERTL0(
false,
"Expansion dimension not supported by FilterAeroForces");
1076 for (i = 0; i < nVel*nVel; i++)
1079 pFields[nVel]->GetPhys(), 1,
1080 P[i]->UpdatePhys(), 1);
1085 for (i=0; i<nVel; i++)
1091 m_mapping->ApplyChristoffelContravar(wk, tmp);
1092 for (i=0; i< nVel; i++)
1094 for (k=0; k< nVel; k++)
1097 wk[i], grad[i*nVel+k]->UpdatePhys());
1100 grad[i*nVel+k]->UpdatePhys(), 1,
1101 grad[i*nVel+k]->UpdatePhys(), 1);
1105 for (i=0; i< nVel; i++)
1107 for (k=0; k< nVel; k++)
1113 for (j=0; j<nVel; j++)
1116 grad[i*nVel+j]->UpdatePhys(), 1);
1121 m_mapping->GetJacobian( Jac->UpdatePhys());
1124 for (i=0; i< nVel; i++)
1127 for (k=0; k< nVel; k++)
1136 for (k=0; k< nVel; k++)
1139 C[k*nVel+i]->UpdatePhys(), 1);
1149 for(plane = 0; plane <
m_nPlanes; plane++ )
1158 for(n = 0; n < nVel*nVel; n++)
1160 PPlane[n] = P[n]->GetPlane(
m_planesID[plane]);
1161 gradPlane[n] = grad[n]->GetPlane(
m_planesID[plane]);
1162 CPlane[n] = C[n]->GetPlane(
m_planesID[plane]);
1168 for(n = 0; n < nVel*nVel; n++)
1171 gradPlane[n] = grad[n];
1178 for( cnt = n = 0; n < BndConds.num_elements(); n++)
1182 for (i=0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
1185 elmt = PPlane[0]->GetExp(elmtid);
1186 offset = PPlane[0]->GetPhys_Offset(elmtid);
1189 for( j=0; j<nVel*nVel; j++)
1191 PElmt[j] = PPlane[j]->GetPhys()
1193 gradElmt[j] = gradPlane[j]->GetPhys()
1195 CElmt[j] = CPlane[j]->GetPhys()
1198 JacElmt = JacPlane->GetPhys() + offset;
1213 bc = BndExp[n]->GetExp(i)->
1214 as<LocalRegions::Expansion1D> ();
1217 nbc = bc->GetTotPoints();
1220 normals = elmt->GetEdgeNormal(boundary);
1223 for(
int j = 0; j < nVel*nVel; ++j)
1227 elmt->GetEdgePhysVals(boundary,
1233 elmt->GetEdgePhysVals(boundary,
1238 elmt->GetEdgePhysVals(boundary,
1244 elmt->GetEdgePhysVals(boundary,
1254 bc = BndExp[n]->GetExp(i)->
1255 as<LocalRegions::Expansion2D> ();
1258 nbc = bc->GetTotPoints();
1261 normals = elmt->GetFaceNormal(boundary);
1264 for(
int j = 0; j < nVel*nVel; ++j)
1268 elmt->GetFacePhysVals(boundary,
1274 elmt->GetFacePhysVals(boundary,
1280 elmt->GetFacePhysVals(boundary,
1286 elmt->GetFacePhysVals(boundary,
1294 "Expansion not supported by FilterForces");
1301 for ( j = 0; j < nVel; j++)
1305 for ( k = 0; k < expdim; k++)
1316 for ( j = 0; j < nVel; j++ )
1319 for ( k = 0; k < expdim; k++ )
1334 for ( k = 0; k < nVel; k++ )
1343 for ( k = 0; k < nVel; k++ )
1346 for ( j = 0; j < nVel; j++ )
1354 for ( k = 0; k < nVel; k++ )
1359 for ( k = 0; k < nVel; k++ )
1362 for ( j = 0; j < nVel; j++ )
1370 for ( k = 0; k < nVel; k++ )
1376 for ( j = 0; j < expdim; j++)
1378 m_Fpplane[j][plane] += BndExp[n]->GetExp(i)->
1380 m_Fvplane[j][plane] += BndExp[n]->GetExp(i)->
1387 cnt += BndExp[n]->GetExpSize();
1396 for( i = 0; i < expdim; i++)
1406 for(plane = 0; plane <
m_nPlanes; plane++)
1410 for( i = 0; i < expdim; i++)
1412 for( j = 0; j < expdim; j++ )
1415 tmpV[i] += m_Fvplane[j][plane]*m_directions[i][j];
1419 for( i = 0; i < expdim; i++)
1422 m_Fvplane[i][plane] = tmpV[i];
1427 for(plane = 0; plane <
m_nPlanes; plane++)
1429 for( i = 0; i < expdim; i++)
1431 m_Ftplane[i][plane] =
m_Fpplane[i][plane] + m_Fvplane[i][plane];
1438 for (i = 0; i < pFields.num_elements(); ++i)
1440 pFields[i]->SetWaveSpace(
true);
1441 pFields[i]->HomogeneousFwdTrans(pFields[i]->GetPhys(),
1442 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
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)
Array< OneD, int > m_BCtoElmtID
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
std::vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
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::vector< unsigned int > m_boundaryRegionsIdList
ID's of boundary regions where we want the forces.
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
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