Write mesh to output file.
113 int rank =
m_f->m_comm->GetRank();
114 int nprocs =
m_f->m_comm->GetSize();
116 if((
m_f->m_verbose)&&(rank == 0))
118 cout <<
"Processing point interpolation" << endl;
126 if(
m_config[
"line"].as<string>().compare(
"NotSet") != 0)
128 string help =
m_config[
"line"].as<
string>();
129 vector<NekDouble> values;
131 m_config[
"line"].as<string>().c_str(),values),
132 "Failed to interpret line string");
135 "line string should contain 2 Dim+1 values "
136 "N,x0,y0,z0,x1,y1,z1");
139 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N is not an integer");
140 ASSERTL0(values[0] > 1,
"N is not a valid number");
142 int dim = (values.size()-1)/2;
143 int npts = values[0];
144 Array<OneD, Array<OneD, NekDouble> > pts(dim);
146 for(
int i = 0; i < dim; ++i)
148 pts[i] = Array<OneD,NekDouble>(
npts);
151 for(
int i = 0; i <
npts; ++i)
153 pts[0][i] = values[1]
154 + i/((
NekDouble)(npts-1))*(values[dim+1] - values[1]);
157 pts[1][i] = values[2]
158 + i/((
NekDouble)(npts-1))*(values[dim+2] - values[2]);
162 pts[2][i] = values[3]
163 + i/((
NekDouble)(npts-1))*(values[dim+3] - values[3]);
172 m_f->m_fieldPts->SetPointsPerEdge(ppe);
175 else if(
m_config[
"plane"].as<string>().compare(
"NotSet") != 0)
177 string help =
m_config[
"plane"].as<
string>();
178 vector<NekDouble> values;
180 m_config[
"plane"].as<string>().c_str(),values),
181 "Failed to interpret plane string");
184 "plane string should contain 4 Dim+2 values "
185 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
188 ASSERTL0(std::modf(values[0], &tmp) == 0.0,
"N1 is not an integer");
189 ASSERTL0(std::modf(values[1], &tmp) == 0.0,
"N2 is not an integer");
191 ASSERTL0(values[0] > 1,
"N1 is not a valid number");
192 ASSERTL0(values[1] > 1,
"N2 is not a valid number");
194 int dim = (values.size()-2)/4;
196 int npts1 = values[0];
197 int npts2 = values[1];
199 Array<OneD, Array<OneD, NekDouble> > pts(dim);
201 int totpts = npts1*npts2;
202 int nlocpts = totpts/nprocs;
206 for(
int i = 0; i < dim; ++i)
208 pts[i] = Array<OneD,NekDouble>(nlocpts);
214 for(
int j = 0; j < npts2; ++j)
216 for(
int i = 0; i < npts1; ++i)
219 if((cnt >= rank*nlocpts)&&(cnt < (rank+1)*nlocpts))
222 (values[2] + i/((
NekDouble)(npts1-1))*(values[dim+2] - values[2]))*(1.0-j/((
NekDouble)(npts2-1))) +
223 (values[3*dim+2] + i/((
NekDouble)(npts1-1))*(values[2*dim+2] - values[3*dim+2]))*(j/((
NekDouble)(npts2-1)));
226 (values[3] + i/((
NekDouble)(npts1-1))*(values[dim+3] - values[3]))*(1.0-j/((
NekDouble)(npts2-1))) +
227 (values[3*dim+3] + i/((
NekDouble)(npts1-1))*(values[2*dim+3] - values[3*dim+3]))*(j/((
NekDouble)(npts2-1)));
232 (values[4] + i/((
NekDouble)(npts1-1))*(values[dim+4] - values[4]))*(1.0-j/((
NekDouble)(npts2-1))) +
233 (values[3*dim+4] + i/((
NekDouble)(npts1-1))*(values[2*dim+4] - values[3*dim+4]))*(j/((
NekDouble)(npts2-1)));
243 totpts = totpts - rank*nlocpts;
245 for(
int i = 0; i < dim; ++i)
247 pts[i] = Array<OneD,NekDouble>(totpts);
253 for(
int j = 0; j < npts2; ++j)
255 for(
int i = 0; i < npts1; ++i)
258 if(cnt >= rank*nlocpts)
261 (values[2] + i/((
NekDouble)(npts1-1))*(values[dim+2] - values[2]))*(1.0-j/((
NekDouble)(npts2-1))) +
262 (values[3*dim+2] + i/((
NekDouble)(npts1-1))*(values[2*dim+2] - values[3*dim+2]))*(j/((
NekDouble)(npts2-1)));
265 (values[3] + i/((
NekDouble)(npts1-1))*(values[dim+3] - values[3]))*(1.0-j/((
NekDouble)(npts2-1))) +
266 (values[3*dim+3] + i/((
NekDouble)(npts1-1))*(values[2*dim+3] - values[3*dim+3]))*(j/((
NekDouble)(npts2-1)));
271 (values[4] + i/((
NekDouble)(npts1-1))*(values[dim+4] - values[4]))*(1.0-j/((
NekDouble)(npts2-1))) +
272 (values[3*dim+4] + i/((
NekDouble)(npts1-1))*(values[2*dim+4] - values[3*dim+4]))*(j/((
NekDouble)(npts2-1)));
282 ppe.push_back(npts1);
283 ppe.push_back(npts2);
286 m_f->m_fieldPts->SetPointsPerEdge(ppe);
289 else if(
m_config[
"box"].as<string>().compare(
"NotSet") != 0)
291 string help =
m_config[
"box"].as<
string>();
292 vector<NekDouble> values;
294 m_config[
"box"].as<string>().c_str(),values),
295 "Failed to interpret box string");
298 "box string should contain 9 values "
299 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
303 int npts1 = values[0];
304 int npts2 = values[1];
305 int npts3 = values[2];
307 Array<OneD, Array<OneD, NekDouble> > pts(dim);
309 int totpts = npts1*npts2*npts3;
310 int nlocpts = totpts/nprocs;
316 for(
int i = 0; i < dim; ++i)
318 pts[i] = Array<OneD,NekDouble>(totpts);
324 for(
int k = 0; k < npts3; ++k)
326 for(
int j = 0; j < npts2; ++j)
328 for(
int i = 0; i < npts1; ++i)
330 if((cnt >= rank*nlocpts)&&(cnt < (rank+1)*nlocpts))
332 pts[0][cntloc] = values[3] +
333 i/((
NekDouble)(npts1-1))*(values[4]-values[3]);
334 pts[1][cntloc] = values[5] +
335 j/((
NekDouble)(npts2-1))*(values[6]-values[5]);
336 pts[2][cntloc] = values[7] +
337 k/((
NekDouble)(npts3-1))*(values[8]-values[7]);
347 totpts = totpts - rank*nlocpts;
349 for(
int i = 0; i < dim; ++i)
351 pts[i] = Array<OneD,NekDouble>(totpts);
357 for(
int k = 0; k < npts3; ++k)
359 for(
int j = 0; j < npts2; ++j)
361 for(
int i = 0; i < npts1; ++i)
363 if(cnt >= rank*nlocpts)
365 pts[0][cntloc] = values[3] +
366 i/((
NekDouble)(npts1-1))*(values[4]-values[3]);
367 pts[1][cntloc] = values[5] +
368 j/((
NekDouble)(npts2-1))*(values[6]-values[5]);
369 pts[2][cntloc] = values[7] +
370 k/((
NekDouble)(npts3-1))*(values[8]-values[7]);
380 ppe.push_back(npts1);
381 ppe.push_back(npts2);
382 ppe.push_back(npts3);
385 m_f->m_fieldPts->SetPointsPerEdge(ppe);
386 vector<NekDouble> boxdim;
387 boxdim.assign(&values[3],&values[3]+6);
388 m_f->m_fieldPts->SetBoxSize(boxdim);
395 std::vector<std::string> files;
397 files.push_back(
m_config[
"fromxml"].as<string>());
404 int coordim =
m_f->m_fieldPts->GetDim();
405 int npts =
m_f->m_fieldPts->GetNpoints();
406 Array<OneD, Array<OneD, NekDouble> > pts;
407 m_f->m_fieldPts->GetPts(pts);
409 rng->m_checkShape =
false;
417 rng->m_doZrange =
true;
420 if(rng->m_zmax == rng->m_zmin)
426 rng->m_doYrange =
true;
430 rng->m_doXrange =
true;
435 ASSERTL0(
false,
"too many values specfied in range");
448 Array<OneD,int> ElementGIDs(expansions.size());
449 SpatialDomains::ExpansionMap::const_iterator expIt;
452 for (expIt = expansions.begin(); expIt != expansions.end();
455 ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
460 ASSERTL0(i > 0,
"No elements are set. Are the interpolated points "
461 "wihtin the domain given by the xml files?");
463 string fromfld =
m_config[
"fromfld"].as<
string>();
464 fromField->m_fld->Import(fromfld,fromField->m_fielddef,
469 int NumHomogeneousDir = fromField->m_fielddef[0]->m_numHomogeneousDir;
473 fromField->m_graph->SetExpansions(fromField->m_fielddef);
475 int nfields = fromField->m_fielddef[0]->m_fields.size();
477 fromField->m_exp.resize(nfields);
478 fromField->m_exp[0] = fromField->SetUpFirstExpList(NumHomogeneousDir,
true);
480 m_f->m_exp.resize(nfields);
483 for(i = 1; i < nfields; ++i)
485 fromField->m_exp[i] = fromField->AppendExpList(NumHomogeneousDir);
489 for(
int j = 0; j < nfields; ++j)
491 for (i = 0; i < fromField->m_fielddef.size(); i++)
493 fromField->m_exp[j]->ExtractDataToCoeffs(
494 fromField->m_fielddef[i],
495 fromField->m_data[i],
496 fromField->m_fielddef[0]->m_fields[j],
497 fromField->m_exp[j]->UpdateCoeffs());
499 fromField->m_exp[j]->BwdTrans(fromField->m_exp[j]->GetCoeffs(),
500 fromField->m_exp[j]->UpdatePhys());
503 m_f->m_fieldPts->AddField(newPts, fromField->m_fielddef[0]->m_fields[j]);
508 cout <<
"Interpolating on proc 0 [" << flush;
516 clamp_low, clamp_up, def_value, !rank);
#define ASSERTL0(condition, msg)
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.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
void InterpolateFieldToPts(vector< MultiRegions::ExpListSharedPtr > &field0, Array< OneD, Array< OneD, NekDouble > > &pts, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value, bool isRoot)
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.
boost::shared_ptr< DomainRange > DomainRangeShPtr
boost::shared_ptr< Field > FieldSharedPtr
static PtsFieldSharedPtr NullPtsField
static bool GenerateUnOrderedVector(const char *const str, std::vector< NekDouble > &vec)
static FieldMetaDataMap NullFieldMetaDataMap
std::map< int, ExpansionShPtr > ExpansionMap