41 #include <boost/format.hpp>
54 FilterHistoryPoints::FilterHistoryPoints(
59 ParamMap::const_iterator it;
62 it = pParams.find(
"OutputFile");
63 if (it == pParams.end())
69 ASSERTL0(it->second.length() > 0,
"Missing parameter 'OutputFile'.");
79 it = pParams.find(
"OutputFrequency");
80 if (it == pParams.end())
94 it = pParams.find(
"OutputPlane");
95 if (it == pParams.end())
105 it = pParams.find(
"WaveSpace");
106 if (it == pParams.end())
112 std::string sOption = it->second.c_str();
114 ( boost::iequals(sOption,
"yes"));
119 it = pParams.find(
"Points");
120 ASSERTL0(it != pParams.end(),
"Missing parameter 'Points'.");
143 "No history points in stream.");
148 vector<unsigned int> planeIDs;
151 int dim = pFields[0]->GetGraph()->GetSpaceDimension();
169 int nplanes = pFields[0]->GetHomogeneousBasis()
170 ->GetZ().num_elements();
171 NekDouble lhom = pFields[0]->GetHomoLen();
176 plane = floor((gloCoord[2]*nplanes)/lhom);
183 NekDouble Z = (pFields[0]->GetHomogeneousBasis()
188 cout <<
"Reseting History point from z = " << gloCoord[2]
189 <<
" to z = " << Z << endl;
192 planeIDs.push_back(plane);
198 gloCoord[1], gloCoord[2]);
209 int vRank = vComm->GetRowComm()->GetRank();
215 std::vector<Array<OneD, NekDouble> > LocCoords;
220 for (i = 0; i < vHP; ++i)
230 idList[i] = pFields[0]->GetPlane(0)->GetExpIndex(gloCoord,locCoords,
235 idList[i] = pFields[0]->GetExpIndex(gloCoord,locCoords,
240 LocCoords.push_back(locCoords);
248 pFields[0]->GetExp(idList[i])->GetGeom();
252 for (
int j = 0; j < g->GetCoordim(); ++j)
254 e->BwdTrans(g->GetCoeffs(j), coordVals);
255 NekDouble x = e->PhysEvaluate(locCoords, coordVals)
271 for (i = 0; i < vHP; ++i)
273 if (dist_loc[i] == dist[i])
293 = pFields[0]->GetZIDs();
294 for(j = 0; j < IDs.num_elements(); ++j)
296 if(IDs[j] == planeIDs[i])
302 if(j != IDs.num_elements())
317 for (i = 0; i < vHP; ++i)
321 if (procList[i] != vRank)
341 if (vComm->GetRank() == 0)
343 for (i = 0; i < vHP; ++i)
352 + boost::lexical_cast<std::string>(gloCoord[0])
354 + boost::lexical_cast<std::string>(gloCoord[1])
356 + boost::lexical_cast<std::string>(gloCoord[2])
357 +
" cannot be found in the mesh.");
363 cout <<
"Warning: History point " << i <<
" at ("
364 << gloCoord[0] <<
"," << gloCoord[1] <<
","
365 << gloCoord[2] <<
") lies a distance of "
366 << sqrt(dist[i]) <<
" from the manifold." << endl;
372 m_session->MatchSolverInfo(
"Driver",
"Adaptive",
384 for (i = 0; i < pFields.num_elements(); ++i)
398 for (i = 0; i < vHP; ++i)
437 int numFields = pFields.num_elements();
440 std::list<std::pair<SpatialDomains::PointGeomSharedPtr, Array<OneD, NekDouble> > >
::iterator x;
447 for (j = 0; j < numFields; ++j)
454 locCoord = (*x).second;
455 expId = (*x).first->GetVid();
461 ASSERTL0( pFields[j]->GetWaveSpace() ==
true,
462 "HistoryPoints in wavespace require that solution is in wavespace");
464 if ( pFields[j]->GetWaveSpace() ==
false ||
m_waveSpace)
468 physvals = pFields[j]->GetPlane(plane)->
469 UpdatePhys() + pFields[j]->GetPhys_Offset(expId);
472 if(pFields[j]->GetPhysState() ==
false)
474 pFields[j]->GetPlane(plane)->GetExp(expId)->
475 BwdTrans(pFields[j]->GetPlane(plane)->
476 GetCoeffs() + pFields[j]->
477 GetCoeff_Offset(expId),physvals);
480 value = pFields[j]->GetPlane(plane)->GetExp(expId)->
481 StdPhysEvaluate(locCoord,physvals);
487 std::vector<unsigned int> eIDs;
488 int nPlanes = pFields[j]->GetZIDs().num_elements();
489 int elmtsPerPlane = pFields[j]->GetExpSize()/nPlanes;
491 for (
int n = 0; n < nPlanes; n++)
493 eIDs.push_back(expId + n*elmtsPerPlane);
500 ASSERTL0(tmp,
"Failed to type cast expansion");
505 AllocateSharedPtr(*tmp, eIDs);
507 for (
int n = 0; n < nPlanes; n++)
510 exp->GetPlane(n)->UpdatePhys();
511 if(pFields[j]->GetPhysState())
513 int nq = exp->GetPlane(0)->GetTotPoints();
515 pFields[j]->GetPlane(n)->GetPhys() +
516 pFields[j]->GetPhys_Offset(expId);
522 pFields[j]->GetPlane(n)->GetCoeffs() +
523 pFields[j]->GetCoeff_Offset(expId);
524 exp->GetPlane(n)->GetExp(0)->
525 BwdTrans(fromCoeffs, toPhys);
528 exp->HomogeneousBwdTrans(exp->GetPhys(), exp->UpdatePhys());
532 physvals = exp->GetPlane(plane)->UpdatePhys();
534 value = exp->GetPlane(plane)->GetExp(0)->
535 StdPhysEvaluate(locCoord,physvals);
542 data[m_historyLocalPointMap[k]*numFields+j] = value;
550 locCoord = (*x).second;
551 expId = (*x).first->GetVid();
553 physvals = pFields[j]->UpdatePhys() + pFields[j]->GetPhys_Offset(expId);
556 if(pFields[j]->GetPhysState() ==
false)
558 pFields[j]->GetExp(expId)->BwdTrans(pFields[j]->GetCoeffs() + pFields[j]->GetCoeff_Offset(expId),physvals);
562 data[
m_historyLocalPointMap[k]*numFields+j] = pFields[j]->GetExp(expId)->StdPhysEvaluate(locCoord,physvals);
572 if (vComm->GetRank() == 0)
579 for (
int j = 0; j < numFields; ++j)
594 if (pFields[0]->GetComm()->GetRank() == 0)
std::ofstream m_outputStream
#define ASSERTL0(condition, msg)
std::map< int, int > m_historyLocalPointMap
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static const NekDouble kVertexTheSameDouble
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
std::stringstream m_historyPointStream
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
static const NekDouble kNekZeroTol
unsigned int m_outputPlane
plane to take history point from if using a homogeneous1D expansion
NekDouble Evaluate() const
SpatialDomains::PointGeomVector m_historyPoints
unsigned int m_outputFrequency
std::map< std::string, std::string > ParamMap
LibUtilities::SessionReaderSharedPtr m_session
SOLVER_UTILS_EXPORT ~FilterHistoryPoints()
virtual SOLVER_UTILS_EXPORT void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
virtual SOLVER_UTILS_EXPORT void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
virtual SOLVER_UTILS_EXPORT bool v_IsTimeDependent()
Array< OneD, int > m_planeIDs
static const NekDouble kGeomFactorsTol
FilterFactory & GetFilterFactory()
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
virtual SOLVER_UTILS_EXPORT void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
boost::shared_ptr< Geometry > GeometrySharedPtr
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
std::list< std::pair< SpatialDomains::PointGeomSharedPtr, Array< OneD, NekDouble > > > m_historyList
boost::shared_ptr< PointGeom > PointGeomSharedPtr
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.