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 expansion"
56 "basis and output fld file. Requires .pts .xml and .fld files.");
65 "reset any nan value to prescribed value");
67 if((f->m_inputfiles.count(
"pts") == 0))
69 cout << endl <<
"A pts input file must be specified for the boundary extraction module" << endl;
71 cout <<
"Usage: Fieldconvert -m pointdatatofld file.pts file.xml out.fld" << endl;
75 if((f->m_inputfiles.count(
"xml") == 0)&&(f->m_inputfiles.count(
"xml.gz") == 0))
77 cout << endl <<
"An xml or xml.gz input file must be specified for the boundary extraction module" << endl;
78 cout <<
"Usage: Fieldconvert -m pointdatatofld file.pts file.xml out.fld" << endl;
91 bool setnantovalue =
false;
94 if(!boost::iequals(
m_config[
"setnantovalue"].as<string>(),
"NotSet"))
102 "No input points found");
104 int nFields =
m_f->m_fieldPts->GetNFields();
106 "No field values provided in input");
108 int dim =
m_f->m_fieldPts->GetDim();
111 m_f->m_exp.resize(nFields);
112 for(i = 1; i < nFields; ++i)
114 m_f->m_exp[i] =
m_f->AppendExpList(0);
118 m_f->m_fieldPts->GetPts(pts);
123 for(
int i = 0; i < pts[0].num_elements(); ++i)
125 for(
int j = 0; j < nFields; ++j)
127 if ((boost::math::isnan)(pts[j+dim][i]))
129 pts[j+dim][i] = defvalue;
137 int totcoeffs =
m_f->m_exp[0]->GetNcoeffs();
139 ASSERTL0(pts[0].num_elements() != totcoeffs,
"length of points in .pts file is different "
140 "to the number of coefficients in expansion ");
142 for(
int i = 0; i < nFields; ++i)
146 for(
int e = 0; e <
m_f->m_exp[0]->GetNumElmts(); ++e)
148 int ncoeffs =
m_f->m_exp[i]->GetExp(e)->GetNcoeffs();
149 int offset =
m_f->m_exp[i]->GetCoeff_Offset(e);
150 Vmath::Vcopy(ncoeffs,&pts[i+dim][cnt],1,&coeffs[offset],1);
152 m_f->m_exp[i]->GetExp(e)->EquiSpacedToCoeffs(coeffs+offset,tmp = coeffs+offset);
159 int totpoints =
m_f->m_exp[0]->GetTotPoints();
166 m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
169 if(pts[0].num_elements() != totpoints)
171 WARNINGL0(
false,
"length of points in .pts file is different to the number of quadrature points in xml file");
172 totpoints = min(totpoints,(
int)pts[0].num_elements());
176 for(i = 0; i < totpoints; ++i)
178 for(j = 0; j < dim; ++j)
180 if(fabs(coords[j][i]-pts[j][i]) > 1e-4)
182 string outstring =
"Coordinates do not match within 1e-4: "
183 + boost::lexical_cast<
string>(coords[j][i])
185 + boost::lexical_cast<string>(pts[j][i])
186 +
" diff " + boost::lexical_cast<
string>(fabs(coords[j][i]-pts[j][i]));;
192 for(j = 0;j < nFields; ++j)
194 m_f->m_exp[j]->SetPhys(i, pts[j+dim][i]);
198 if(
m_f->m_session->GetComm()->GetRank() == 0)
204 for(i = 0; i < nFields; ++i)
206 m_f->m_exp[i]->FwdTrans_IterPerExp(
m_f->m_exp[i]->GetPhys(),
207 m_f->m_exp[i]->UpdateCoeffs());
212 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
213 =
m_f->m_exp[0]->GetFieldDefinitions();
214 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
216 for (j = 0; j < nFields; ++j)
218 for (i = 0; i < FieldDef.size(); ++i)
220 FieldDef[i]->m_fields.push_back(
m_f->m_fieldPts->GetFieldName(j));
222 m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
226 m_f->m_fielddef = FieldDef;
227 m_f->m_data = FieldData;
#define ASSERTL0(condition, msg)
pair< ModuleType, string > ModuleKey
map< string, ConfigOption > m_config
List of configuration values.
FieldSharedPtr m_f
Field object.
virtual ~ProcessPointDataToFld()
#define WARNINGL0(condition, msg)
boost::shared_ptr< Field > FieldSharedPtr
Represents a command-line configuration option.
static PtsFieldSharedPtr NullPtsField
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
ModuleFactory & GetModuleFactory()
Abstract base class for processing modules.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.