41#include <boost/algorithm/string.hpp>
54 const std::weak_ptr<EquationSystem> &pEquation,
const ParamMap &pParams)
55 :
Filter(pSession, pEquation)
60 auto it = pParams.find(
"OutputFile");
61 if (it == pParams.end())
67 ASSERTL0(it->second.length() > 0,
"Empty parameter 'OutputFile'.");
73 it = pParams.find(
"Composites");
74 ASSERTL0(it != pParams.end(),
"Missing parameter 'Composites'.");
75 ASSERTL0(it->second.length() > 0,
"Empty parameter 'Composites'.");
77 std::vector<std::string> splitComposite;
79 boost::token_compress_on);
84 size_t first = comp.find_first_of(
'[') + 1;
85 size_t last = comp.find_last_of(
']') - 1;
86 auto tmpString = comp.substr(first, last - first + 1);
88 std::vector<unsigned int> tmpVec;
91 ASSERTL0(parseGood && !tmpVec.empty(),
92 "Unable to read composite regions index range for "
101 it = pParams.find(
"OutputPrecision");
102 if (it == pParams.end())
108 ASSERTL0(it->second.length() > 0,
"Empty parameter 'OutputPrecision'.");
109 precision = std::stoi(it->second);
113 auto equationSys =
m_equ.lock();
114 ASSERTL0(equationSys,
"Weak pointer expired");
118 m_comm = pSession->GetComm();
119 if (
m_comm->GetRank() == 0)
123 m_outFile.setf(std::ios::scientific, std::ios::floatfield);
131 std::string varName = equationSys->GetVariable(j);
132 m_outFile <<
" " + compName +
"_" + varName +
"_integral";
139 it = pParams.find(
"OutputFrequency");
140 if (it == pParams.end())
146 ASSERTL0(it->second.length() > 0,
"Empty parameter 'OutputFrequency'.");
161 auto expList = pFields[0]->GetExp();
162 auto meshGraph = pFields[0]->GetGraph();
164 for (
size_t i = 0; i < expList->size(); ++i)
166 auto exp = (*expList)[i];
171 std::map<size_t, size_t> geomIdToTraceId;
172 auto trace = pFields[0]->GetTrace()->GetExp();
173 for (
size_t i = 0; i < trace->size(); ++i)
175 auto exp = (*trace)[i];
176 geomIdToTraceId[exp->GetGeom()->GetGlobalID()] = i;
180 auto composites = pFields[0]->GetGraph()->GetComposites();
181 size_t meshDim = pFields[0]->GetGraph()->GetMeshDimension();
186 if (composites.find(
m_compVector[i][0]) == composites.end())
191 std::vector<std::shared_ptr<SpatialDomains::Geometry>> geomVec =
194 composites[
m_compVector[i][0]]->m_geomVec[0]->GetShapeDim();
197 std::vector<size_t> compGeomIds;
200 ASSERTL0(composites.find(compNum) != composites.end(),
201 "In FilterIntegral defined composite C[" +
202 std::to_string(compNum) +
203 "] does not exist in the mesh.")
205 auto compGeom = composites[compNum]->m_geomVec;
207 for (
auto &geom : compGeom)
209 compGeomIds.emplace_back(geom->GetGlobalID());
214 dim == compGeom[0]->GetShapeDim(),
215 "Differing geometry dimensions specified in FilterIntegral '" +
219 std::vector<std::pair<LocalRegions::ExpansionSharedPtr, int>>
220 tmpCompExp(compGeomIds.size());
224 if (dim == pFields[0]->GetShapeDimension())
226 for (
size_t j = 0; j < compGeomIds.size(); ++j)
228 tmpCompExp[j] = std::make_pair(
235 else if (meshDim > 1 && dim == meshDim - 1)
237 for (
size_t j = 0; j < compGeomIds.size(); ++j)
240 (*trace)[geomIdToTraceId[compGeomIds[j]]];
243 exp->GetLeftAdjacentElementExp();
244 int leftAdjElmtFace = exp->GetLeftAdjacentElementTrace();
246 tmpCompExp[j] = std::make_pair(leftAdjElmtExp, leftAdjElmtFace);
252 "FilterIntegral: Only composite dimensions equal to or "
253 "one lower than the mesh dimension are supported.")
274 if (
m_comm->GetRank() == 0)
284 phys = pFields[i]->GetPhys();
293 size_t dim = compExp[0].first->GetGeom()->GetShapeDim();
294 size_t meshDim = pFields[i]->GetGraph()->GetMeshDimension();
299 for (
auto &expPair : compExp)
302 auto exp = expPair.first;
306 size_t offset = pFields[i]->GetPhys_Offset(
308 input = exp->Integral(phys + offset);
310 else if (meshDim > 1 && dim == meshDim - 1)
313 exp->GetTracePhysVals(expPair.second, exp, phys,
318 ->GetExp(exp->GetGeom()->GetTid(expPair.second))
319 ->Integral(tracePhys);
324 "FilterIntegral: Only composite dimensions "
325 "equal to or one lower than the mesh "
326 "dimension are supported.")
330 c += fabs(sum) >= fabs(input) ? (sum - t) + input
339 if (
m_comm->GetRank() == 0)
346 if (
m_comm->GetRank() == 0)
360 if (
m_comm->GetRank() == 0)
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
NekDouble Evaluate() const
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static bool GenerateSeqVector(const std::string &str, std::vector< unsigned int > &out)
Takes a comma-separated compressed string and converts it to entries in a vector.
LibUtilities::SessionReaderSharedPtr m_session
const std::weak_ptr< EquationSystem > m_equ
std::map< std::string, std::string > ParamMap
std::map< int, std::vector< std::pair< LocalRegions::ExpansionSharedPtr, int > > > m_compExpMap
Map of composite ID to vector of expansions an face/edge local ID.
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) final
Calculate integral over requested composites.
std::map< size_t, size_t > m_geomElmtIdToExpId
Mapping from geometry ID to expansion ID.
size_t m_outputFrequency
Frequency to write to output data file in timesteps.
size_t m_numVariables
Number of fields to perform integral on.
std::ofstream m_outFile
Out file.
bool v_IsTimeDependent() final
Returns true as filter depends on time.
LibUtilities::CommSharedPtr m_comm
Global communicator.
static std::string className
Name of the class.
void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) final
Finalise the filter.
std::vector< std::vector< unsigned int > > m_compVector
Vector of vector of composites IDs as integers.
std::vector< std::string > m_splitCompString
Vector of composite IDs as a single string.
SOLVER_UTILS_EXPORT FilterIntegral(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< EquationSystem > &pEquation, const ParamMap &pParams)
Constructs the integral filter and parses filter options, opens file.
void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) final
Parse composite list and geometric entities to integrate over.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
FilterFactory & GetFilterFactory()