43 #include <boost/math/special_functions/fpclassify.hpp>
52 ProcessInterpField::create,
53 "Interpolates one field to another, requires fromxml, "
54 "fromfld to be defined");
61 "Xml file form which to interpolate field");
63 "Fld file form which to interpolate field");
66 "Lower bound for interpolation value");
68 "Upper bound for interpolation value");
70 "Default value if point is outside domain");
82 cout <<
"Processing interpolation" << endl;
87 std::vector<std::string> files;
90 files.push_back(
m_config[
"fromxml"].as<string>());
98 int coordim =
m_f->m_exp[0]->GetCoordim(0);
99 int npts =
m_f->m_exp[0]->GetTotPoints();
103 for(
int i = 0; i < coordim; ++i)
108 for(
int i = coordim; i < 3; ++i)
113 m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
115 rng->m_checkShape =
false;
119 rng->m_doZrange =
true;
123 rng->m_doYrange =
true;
127 rng->m_doXrange =
true;
132 ASSERTL0(
false,
"too many values specfied in range");
145 if(!expansions.size())
151 SpatialDomains::ExpansionMap::const_iterator expIt;
154 for (expIt = expansions.begin(); expIt != expansions.end();
157 ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
160 string fromfld =
m_config[
"fromfld"].as<
string>();
166 int NumHomogeneousDir =
m_fromField->m_fielddef[0]->m_numHomogeneousDir;
172 int nfields =
m_fromField->m_fielddef[0]->m_fields.size();
177 m_f->m_exp.resize(nfields);
180 for(i = 1; i < nfields; ++i)
182 m_f->m_exp[i] =
m_f->AppendExpList(NumHomogeneousDir);
188 for(
int j = 0; j < nfields; ++j)
190 for (i = 0; i <
m_fromField->m_fielddef.size(); i++)
203 int nq1 =
m_f->m_exp[0]->GetTotPoints();
211 m_f->m_exp[0]->GetCoords(x1, y1);
213 else if (coordim == 3)
215 m_f->m_exp[0]->GetCoords(x1, y1, z1);
218 if(
m_f->m_session->GetComm()->TreatAsRankZero())
220 cout <<
"Interpolating [" << flush;
228 x1, y1, z1, clamp_low, clamp_up,def_value);
230 if(
m_f->m_session->GetComm()->TreatAsRankZero())
237 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
238 =
m_f->m_exp[0]->GetFieldDefinitions();
239 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
241 for (
int j = 0; j < nfields; ++j)
243 m_f->m_exp[j]->FwdTrans(
m_f->m_exp[j]->GetPhys(),
244 m_f->m_exp[j]->UpdateCoeffs());
245 for (i = 0; i < FieldDef.size(); ++i)
247 FieldDef[i]->m_fields.push_back(
m_fromField->m_fielddef[0]->m_fields[j]);
248 m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
252 m_f->m_fielddef = FieldDef;
253 m_f->m_data = FieldData;
257 vector<MultiRegions::ExpListSharedPtr> &field0,
258 vector<MultiRegions::ExpListSharedPtr> &field1,
266 int expdim = field0[0]->GetCoordim(0);
269 int nq1 = field1[0]->GetTotPoints();
272 static int intpts = 0;
274 ASSERTL0(field0.size() == field1.size(),
275 "Input field dimension must be same as output dimension");
277 for (r = 0; r < nq1; r++)
287 elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3,
true);
291 offset = field0[0]->GetPhys_Offset(field0[0]->
292 GetOffset_Elmt_Id(elmtid));
294 for (f = 0; f < field1.size(); ++f)
297 value = field0[f]->GetExp(elmtid)->
298 StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset);
300 if ((boost::math::isnan)(value))
302 ASSERTL0(
false,
"new value is not a number");
310 else if( value < clamp_low)
315 field1[f]->UpdatePhys()[r] = value;
321 for (f = 0; f < field1.size(); ++f)
323 field1[f]->UpdatePhys()[r] = def_value;
327 if (intpts%1000 == 0)
#define ASSERTL0(condition, msg)
virtual ~ProcessInterpField()
pair< ModuleType, string > ModuleKey
static Array< OneD, NekDouble > NullNekDouble1DArray
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
void InterpolateField(vector< MultiRegions::ExpListSharedPtr > &field0, vector< MultiRegions::ExpListSharedPtr > &field1, Array< OneD, NekDouble > x, Array< OneD, NekDouble > y, Array< OneD, NekDouble > z, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value)
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
map< string, ConfigOption > m_config
List of configuration values.
FieldSharedPtr m_f
Field object.
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
FieldSharedPtr m_fromField
boost::shared_ptr< DomainRange > DomainRangeShPtr
boost::shared_ptr< Field > FieldSharedPtr
Represents a command-line configuration option.
static FieldMetaDataMap NullFieldMetaDataMap
std::map< int, ExpansionShPtr > ExpansionMap
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.