Commit a33189fd authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Added Filters base class and factory

Added ThresholdMax filter.
Moved activation time calculation in UnsteadySystem to be done by filter.
Closes ticket #65


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@3658 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent e8320133
......@@ -834,6 +834,15 @@ namespace Nektar
}
/**
*
*/
const FilterMap &SessionReader::GetFilters() const
{
return m_filters;
}
/**
*
*/
......@@ -925,6 +934,10 @@ namespace Nektar
e = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
ReadGeometricInfo(e);
e = docHandle.FirstChildElement("NEKTAR").FirstChildElement("FILTERS").Element();
ReadFilters(e);
}
......@@ -1491,6 +1504,53 @@ namespace Nektar
}
/**
*
*/
void SessionReader::ReadFilters(TiXmlElement *filters)
{
if (!filters)
{
return;
}
m_filters.clear();
TiXmlElement *filter = filters->FirstChildElement("FILTER");
while (filter)
{
std::map<std::string, std::string> vParams;
TiXmlElement *type = filter->FirstChildElement("TYPE");
ASSERTL0(type, "No type specified for filter.");
ASSERTL0(type->GetText(), "Empty type string specified.");
std::string typeStr = type->GetText();
TiXmlElement *param = filter->FirstChildElement("PARAM");
while (param)
{
TiXmlElement *name = param->FirstChildElement("NAME");
ASSERTL0(name, "No name specified for parameter.");
ASSERTL0(name->GetText(), "Empty name string for param.");
std::string nameStr = name->GetText();
TiXmlElement *value = param->FirstChildElement("VALUE");
ASSERTL0(value, "No value specified for parameter.");
ASSERTL0(value->GetText(), "Empty value string for param.");
std::string valueStr = value->GetText();
vParams[nameStr] = valueStr;
param = param->NextSiblingElement("PARAM");
}
m_filters.insert(std::pair<std::string, std::map<std::string, std::string> >(typeStr, vParams));
filter = filter->NextSiblingElement("FILTER");
}
}
/**
*
*/
......
......@@ -62,6 +62,7 @@ namespace Nektar
typedef std::vector<std::string> VariableList;
typedef std::map<std::string, EquationSharedPtr> EquationMap;
typedef std::map<std::string, std::string> TagMap;
typedef std::multimap<std::string, std::map<std::string, std::string> > FilterMap;
typedef std::map<std::string, int> EnumMap;
typedef std::map<std::string, EnumMap> EnumMapList;
......@@ -234,6 +235,9 @@ namespace Nektar
/// Returns the value of a specified tag.
LIB_UTILITIES_EXPORT const std::string &GetTag(const std::string& pName) const;
/* ------ FILTERS ------ */
LIB_UTILITIES_EXPORT const FilterMap& GetFilters() const;
/// Substitutes expressions defined in the XML document.
LIB_UTILITIES_EXPORT void SubstituteExpressions(std::string &expr);
......@@ -263,6 +267,8 @@ namespace Nektar
VariableList m_variables;
/// Custom tags.
TagMap m_tags;
/// Filters map.
FilterMap m_filters;
/// Be verbose
bool m_verbose;
......@@ -303,6 +309,8 @@ namespace Nektar
LIB_UTILITIES_EXPORT void ReadVariables(TiXmlElement *conditions);
/// Reads the FUNCTIONS section of the XML document.
LIB_UTILITIES_EXPORT void ReadFunctions(TiXmlElement *conditions);
/// Reads the FILTERS section of the XML document.
LIB_UTILITIES_EXPORT void ReadFilters(TiXmlElement *filters);
/// Perform a case-insensitive string comparison.
LIB_UTILITIES_EXPORT int NoCaseStringCompare(const std::string &s1, const std::string &s2) const;
......
......@@ -734,6 +734,11 @@ namespace Nektar
return m_comm;
}
SpatialDomains::MeshGraphSharedPtr GetGraph()
{
return m_graph;
}
// Wrapper functions for Homogeneous Expansions
inline LibUtilities::BasisSharedPtr GetHomogeneousBasis(void)
{
......
......@@ -15,6 +15,8 @@ SET(ADRSolverSource
../Auxiliary/DriverStandard.cpp
../Auxiliary/DriverArnoldi.cpp
../Auxiliary/DriverModifiedArnoldi.cpp
../Auxiliary/Filters/Filter.cpp
../Auxiliary/Filters/FilterThresholdMax.cpp
../Auxiliary/UnsteadySystem.cpp
../Auxiliary/EquationSystem.cpp)
......
/*
* Filter.cpp
*
* Created on: 24 Feb 2012
* Author: cc
*/
#include <Auxiliary/Filters/Filter.h>
namespace Nektar
{
FilterFactory& GetFilterFactory()
{
typedef Loki::SingletonHolder<FilterFactory,
Loki::CreateUsingNew,
Loki::NoDestroy > Type;
return Type::Instance();
}
Filter::Filter()
{
}
Filter::~Filter()
{
}
}
/*
* Filter.h
*
* Created on: 24 Feb 2012
* Author: cc
*/
#ifndef FILTER_H_
#define FILTER_H_
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <MultiRegions/ExpList.h>
namespace Nektar
{
class Filter;
/// A shared pointer to a Driver object
typedef boost::shared_ptr<Filter> FilterSharedPtr;
/// Datatype of the NekFactory used to instantiate classes derived from
/// the Driver class.
typedef LibUtilities::NekFactory<
std::string, Filter,
const std::map<std::string, std::string>&
> FilterFactory;
FilterFactory& GetFilterFactory();
class Filter
{
public:
Filter();
~Filter();
inline void Initialise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time);
inline void Update(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time);
inline void Finalise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time);
inline bool IsTimeDependent();
protected:
virtual void v_Initialise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time) = 0;
virtual void v_Update(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time) = 0;
virtual void v_Finalise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time) = 0;
virtual bool v_IsTimeDependent() = 0;
};
inline void Filter::Initialise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time)
{
v_Initialise(pFields, time);
}
inline void Filter::Update(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time)
{
v_Update(pFields, time);
}
inline void Filter::Finalise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time)
{
v_Finalise(pFields, time);
}
inline bool Filter::IsTimeDependent()
{
return v_IsTimeDependent();
}
}
#endif /* FILTER_H_ */
/*
* FilterThresholdMax.cpp
*
* Created on: 24 Feb 2012
* Author: cc
*/
#include <Auxiliary/Filters/FilterThresholdMax.h>
namespace Nektar
{
std::string FilterThresholdMax::className = GetFilterFactory().RegisterCreatorFunction("ThresholdMax", FilterThresholdMax::create);
FilterThresholdMax::FilterThresholdMax(const std::map<std::string, std::string> &pParams) :
Filter()
{
ASSERTL0(pParams.find("ThresholdValue") != pParams.end(),
"Missing parameter 'ThresholdValue'.");
m_thresholdValue = atof(pParams.find("ThresholdValue")->second.c_str());
ASSERTL0(pParams.find("InitialValue") != pParams.end(),
"Missing parameter 'InitialValue'.");
m_initialValue = atof(pParams.find("InitialValue")->second.c_str());
ASSERTL0(!(pParams.find("OutputFile")->second.empty()),
"Missing parameter 'OutputFile'.");
m_outputFile = pParams.find("OutputFile")->second;
}
FilterThresholdMax::~FilterThresholdMax()
{
}
void FilterThresholdMax::v_Initialise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time)
{
m_threshold = Array<OneD, NekDouble> (pFields[0]->GetNpoints(), m_initialValue);
}
void FilterThresholdMax::v_Update(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time)
{
int i;
NekDouble timestep = pFields[0]->GetSession()->GetParameter("TimeStep");
for (i = 0; i < pFields[0]->GetNpoints(); ++i)
{
if (m_threshold[i] < timestep && pFields[0]->GetPhys()[i] > m_thresholdValue)
{
m_threshold[i] = time;
}
}
}
void FilterThresholdMax::v_Finalise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time)
{
SpatialDomains::MeshGraphSharedPtr vGraph = pFields[0]->GetGraph();
std::vector<SpatialDomains::FieldDefinitionsSharedPtr> FieldDef
= pFields[0]->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
Array<OneD, NekDouble> vCoeffs(pFields[0]->GetNcoeffs());
pFields[0]->FwdTrans_IterPerExp(m_threshold, vCoeffs);
// copy Data into FieldData and set variable
for(int i = 0; i < FieldDef.size(); ++i)
{
// Could do a search here to find correct variable
FieldDef[i]->m_fields.push_back("m");
pFields[0]->AppendFieldData(FieldDef[i], FieldData[i], vCoeffs);
}
vGraph->Write(m_outputFile,FieldDef,FieldData);
}
bool FilterThresholdMax::v_IsTimeDependent()
{
return true;
}
}
/*
* FilterThresholdMax.h
*
* Created on: 24 Feb 2012
* Author: cc
*/
#ifndef FILTERTHRESHOLDMAX_H_
#define FILTERTHRESHOLDMAX_H_
#include <Auxiliary/Filters/Filter.h>
namespace Nektar
{
class FilterThresholdMax : public Filter
{
public:
friend class MemoryManager<FilterThresholdMax>;
/// Creates an instance of this class
static FilterSharedPtr create(const std::map<std::string, std::string> &pParams) {
FilterSharedPtr p = MemoryManager<FilterThresholdMax>::AllocateSharedPtr(pParams);
//p->InitObject();
return p;
}
///Name of the class
static std::string className;
FilterThresholdMax(const std::map<std::string, std::string> &pParams);
~FilterThresholdMax();
protected:
virtual void v_Initialise(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time);
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);
virtual bool v_IsTimeDependent();
private:
Array<OneD, NekDouble> m_threshold;
NekDouble m_thresholdValue;
NekDouble m_initialValue;
std::string m_outputFile;
};
}
#endif /* FILTERTHRESHOLDMAX_H_ */
......@@ -101,6 +101,14 @@ namespace Nektar
m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
m_session->LoadParameter("IO_HistorySteps", m_historysteps, 0);
m_session->LoadParameter("CFL", m_cfl, 0.0);
// Set up filters
LibUtilities::FilterMap::const_iterator x;
LibUtilities::FilterMap f = m_session->GetFilters();
for (x = f.begin(); x != f.end(); ++x)
{
m_filters.push_back(GetFilterFactory().CreateInstance(x->first, x->second));
}
}
......@@ -257,6 +265,12 @@ namespace Nektar
hisFile.open(outname.c_str());
}
std::vector<FilterSharedPtr>::iterator x;
for (x = m_filters.begin(); x != m_filters.end(); ++x)
{
(*x)->Initialise(m_fields, m_time);
}
const Array<OneD,int> ExpOrder = GetNumExpModesPerExp();
NekDouble TimeStability;
......@@ -411,8 +425,6 @@ namespace Nektar
else
{
m_actTime = Array<OneD, NekDouble> (m_fields[0]->GetNpoints(), 0.0);
m_actThreshold = 0.0;
for(n = 0; n < m_steps; ++n)
{
Timer timer;
......@@ -455,15 +467,10 @@ namespace Nektar
}
}
if (m_session->DefinesSolverInfo("ComputeActivationTime"))
std::vector<FilterSharedPtr>::iterator x;
for (x = m_filters.begin(); x != m_filters.end(); ++x)
{
for (i = 0; i < m_fields[0]->GetNpoints(); ++i)
{
if (m_actTime[i] < m_timestep && m_fields[0]->GetPhys()[i] > m_actThreshold)
{
m_actTime[i] = m_time;
}
}
(*x)->Update(m_fields, m_time);
}
// Write out history data to file
......@@ -493,32 +500,11 @@ namespace Nektar
m_fields[m_intVariables[i]]->UpdatePhys() = fields[i];
}
if (m_session->DefinesSolverInfo("ComputeActivationTime"))
{
NekDouble tmp;
// swap in
for (i = 0; i < m_fields[0]->GetNpoints(); ++i)
{
tmp = m_fields[0]->GetPhys()[i];
m_fields[0]->UpdatePhys()[i] = m_actTime[i];
m_actTime[i] = tmp;
m_fields[0]->FwdTrans_IterPerExp(m_fields[0]->GetPhys(), m_fields[0]->UpdateCoeffs());
m_fields[0]->SetPhysState(false);
}
// write out
Checkpoint_Output(nchk);
// swap out
for (i = 0; i < m_fields[0]->GetNpoints(); ++i)
{
tmp = m_fields[0]->GetPhys()[i];
m_fields[0]->UpdatePhys()[i] = m_actTime[i];
m_actTime[i] = tmp;
m_fields[0]->FwdTrans_IterPerExp(m_fields[0]->GetPhys(), m_fields[0]->UpdateCoeffs());
m_fields[0]->SetPhysState(false);
}
}
std::vector<FilterSharedPtr>::iterator x;
for (x = m_filters.begin(); x != m_filters.end(); ++x)
{
(*x)->Finalise(m_fields, m_time);
}
}
}
......
......@@ -37,6 +37,7 @@
#define NEKTAR_SOLVERS_ADRSOLVER_EQUATIONSYSTEMS_UNSTEADYSYSTEM_H
#include <Auxiliary/EquationSystem.h>
#include <Auxiliary/Filters/Filter.h>
#include <time.h>
namespace Nektar
......@@ -89,11 +90,10 @@ namespace Nektar
/// Indicates if explicit or implicit treatment of reaction is used.
bool m_explicitReaction;
Array<OneD, NekDouble> m_actTime;
NekDouble m_actThreshold;
std::vector<int> m_intVariables;
std::vector<FilterSharedPtr> m_filters;
/// Initialises UnsteadySystem class members.
UnsteadySystem(const LibUtilities::SessionReaderSharedPtr& pSession);
......
......@@ -15,7 +15,10 @@ SET(CardiacEPSolverSource
../Auxiliary/UnsteadySystem.cpp
../Auxiliary/EquationSystem.cpp
../Auxiliary/Driver.cpp
../Auxiliary/DriverStandard.cpp)
../Auxiliary/DriverStandard.cpp
../Auxiliary/Filters/Filter.cpp
../Auxiliary/Filters/FilterThresholdMax.cpp
)
ADD_SOLVER_EXECUTABLE(CardiacEPSolver solvers-extra
${CardiacEPSolverSource})
......@@ -6,6 +6,8 @@ SET(PulseWaveSolverSource
../Auxiliary/DriverStandard.cpp
../Auxiliary/DriverArnoldi.cpp
../Auxiliary/DriverModifiedArnoldi.cpp
../Auxiliary/Filters/Filter.cpp
../Auxiliary/Filters/FilterThresholdMax.cpp
../Auxiliary/UnsteadySystem.cpp
../Auxiliary/EquationSystem.cpp)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment