45 #include <boost/math/special_functions/fpclassify.hpp>
51 ModuleKey ProcessPointDataToFld::className =
54 ProcessPointDataToFld::create,
55 "Given discrete data at quadrature points project them onto an "
57 "basis and output fld file. Requires .pts .xml and .fld files.");
65 false,
"NotSet",
"reset any nan value to prescribed value");
67 if ((f->m_inputfiles.count(
"pts") == 0))
70 <<
"A pts input file must be specified for the boundary "
75 <<
"Usage: Fieldconvert -m pointdatatofld file.pts file.xml out.fld"
80 if ((f->m_inputfiles.count(
"xml") == 0) &&
81 (f->m_inputfiles.count(
"xml.gz") == 0))
84 <<
"An xml or xml.gz input file must be specified for the "
85 "boundary extraction module"
88 <<
"Usage: Fieldconvert -m pointdatatofld file.pts file.xml out.fld"
102 if (
m_f->m_comm->TreatAsRankZero())
104 cout <<
"ProcessPointDataToFld: projecting data to expansion..."
110 bool setnantovalue =
false;
113 if (!boost::iequals(
m_config[
"setnantovalue"].as<string>(),
"NotSet"))
115 setnantovalue =
true;
121 "No input points found");
123 int nFields =
m_f->m_fieldPts->GetNFields();
124 ASSERTL0(nFields > 0,
"No field values provided in input");
126 int dim =
m_f->m_fieldPts->GetDim();
129 m_f->m_exp.resize(nFields);
130 for (i = 1; i < nFields; ++i)
132 m_f->m_exp[i] =
m_f->AppendExpList(0);
136 m_f->m_fieldPts->GetPts(pts);
141 for (
int i = 0; i < pts[0].num_elements(); ++i)
143 for (
int j = 0; j < nFields; ++j)
145 if ((boost::math::isnan)(pts[j + dim][i]))
147 pts[j + dim][i] = defvalue;
155 int totcoeffs =
m_f->m_exp[0]->GetNcoeffs();
157 ASSERTL0(pts[0].num_elements() != totcoeffs,
158 "length of points in .pts file is different "
159 "to the number of coefficients in expansion ");
161 for (
int i = 0; i < nFields; ++i)
165 for (
int e = 0; e <
m_f->m_exp[0]->GetNumElmts(); ++e)
167 int ncoeffs =
m_f->m_exp[i]->GetExp(e)->GetNcoeffs();
168 int offset =
m_f->m_exp[i]->GetCoeff_Offset(e);
169 Vmath::Vcopy(ncoeffs, &pts[i + dim][cnt], 1, &coeffs[offset],
172 m_f->m_exp[i]->GetExp(e)->EquiSpacedToCoeffs(
173 coeffs + offset, tmp = coeffs + offset);
180 int totpoints =
m_f->m_exp[0]->GetTotPoints();
187 m_f->m_exp[0]->GetCoords(coords[0], coords[1], coords[2]);
189 if (pts[0].num_elements() != totpoints)
191 WARNINGL0(
false,
"length of points in .pts file is different to "
192 "the number of quadrature points in xml file");
193 totpoints = min(totpoints, (
int)pts[0].num_elements());
197 for (i = 0; i < totpoints; ++i)
199 for (j = 0; j < dim; ++j)
201 if (fabs(coords[j][i] - pts[j][i]) > 1e-4)
204 "Coordinates do not match within 1e-4: " +
205 boost::lexical_cast<
string>(coords[j][i]) +
" versus " +
206 boost::lexical_cast<string>(pts[j][i]) +
" diff " +
207 boost::lexical_cast<
string>(
208 fabs(coords[j][i] - pts[j][i]));
215 for (j = 0; j < nFields; ++j)
217 m_f->m_exp[j]->SetPhys(i, pts[j + dim][i]);
221 if (
m_f->m_session->GetComm()->GetRank() == 0)
227 for (i = 0; i < nFields; ++i)
229 m_f->m_exp[i]->FwdTrans_IterPerExp(
m_f->m_exp[i]->GetPhys(),
230 m_f->m_exp[i]->UpdateCoeffs());
235 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
236 m_f->m_exp[0]->GetFieldDefinitions();
237 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
239 for (j = 0; j < nFields; ++j)
241 for (i = 0; i < FieldDef.size(); ++i)
243 FieldDef[i]->m_fields.push_back(
m_f->m_fieldPts->GetFieldName(j));
245 m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
249 m_f->m_fielddef = FieldDef;
250 m_f->m_data = FieldData;
map< string, ConfigOption > m_config
List of configuration values.
#define ASSERTL0(condition, msg)
Represents a command-line configuration option.
virtual void Process(po::variables_map &vm)
Write mesh to output file.
virtual ~ProcessPointDataToFld()
pair< ModuleType, string > ModuleKey
boost::shared_ptr< Field > FieldSharedPtr
#define WARNINGL0(condition, msg)
static PtsFieldSharedPtr NullPtsField
Abstract base class for processing modules.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
ModuleFactory & GetModuleFactory()
FieldSharedPtr m_f
Field object.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.