Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Static Public Member Functions | Static Public Attributes | Private Member Functions | Private Attributes | List of all members
Nektar::Utilities::ProcessInterpPoints Class Reference

This processing module interpolates one field to another. More...

#include <ProcessInterpPoints.h>

Inheritance diagram for Nektar::Utilities::ProcessInterpPoints:
Inheritance graph
[legend]
Collaboration diagram for Nektar::Utilities::ProcessInterpPoints:
Collaboration graph
[legend]

Public Member Functions

 ProcessInterpPoints (FieldSharedPtr f)
virtual ~ProcessInterpPoints ()
virtual void Process (po::variables_map &vm)
 Write mesh to output file.
- Public Member Functions inherited from Nektar::Utilities::ProcessModule
 ProcessModule (FieldSharedPtr p_f)
 ProcessModule (MeshSharedPtr p_m)
- Public Member Functions inherited from Nektar::Utilities::Module
 Module (FieldSharedPtr p_f)
void RegisterConfig (string key, string value)
 Register a configuration option with a module.
void PrintConfig ()
 Print out all configuration options for a module.
void SetDefaults ()
 Sets default configuration options for those which have not been set.
bool GetRequireEquiSpaced (void)
void SetRequireEquiSpaced (bool pVal)
 Module (MeshSharedPtr p_m)
virtual void Process ()=0
void RegisterConfig (string key, string value)
void PrintConfig ()
void SetDefaults ()
MeshSharedPtr GetMesh ()
virtual void ProcessVertices ()
 Extract element vertices.

Static Public Member Functions

static boost::shared_ptr< Module > create (FieldSharedPtr f)
 Creates an instance of this class.

Static Public Attributes

static ModuleKey className

Private Member Functions

void InterpolateFieldToPts (vector< MultiRegions::ExpListSharedPtr > &field0, Array< OneD, Array< OneD, NekDouble > > &pts, NekDouble clamp_low, NekDouble clamp_up, NekDouble def_value)

Private Attributes

FieldSharedPtr m_fromField

Additional Inherited Members

- Protected Member Functions inherited from Nektar::Utilities::Module
virtual void ProcessEdges (bool ReprocessEdges=true)
 Extract element edges.
virtual void ProcessFaces (bool ReprocessFaces=true)
 Extract element faces.
virtual void ProcessElements ()
 Generate element IDs.
virtual void ProcessComposites ()
 Generate composites.
void ReorderPrisms (PerMap &perFaces)
 Reorder node IDs so that prisms and tetrahedra are aligned correctly.
void PrismLines (int prism, PerMap &perFaces, set< int > &prismsDone, vector< ElementSharedPtr > &line)
- Protected Attributes inherited from Nektar::Utilities::Module
FieldSharedPtr m_f
 Field object.
map< string, ConfigOption > m_config
 List of configuration values.
bool m_requireEquiSpaced
MeshSharedPtr m_mesh
 Mesh object.

Detailed Description

This processing module interpolates one field to another.

Definition at line 48 of file ProcessInterpPoints.h.

Constructor & Destructor Documentation

Nektar::Utilities::ProcessInterpPoints::ProcessInterpPoints ( FieldSharedPtr  f)

Definition at line 57 of file ProcessInterpPoints.cpp.

References ASSERTL0, and Nektar::Utilities::Module::m_config.

: ProcessModule(f)
{
m_config["fromxml"] =
ConfigOption(false, "NotSet",
"Xml file from which to interpolate field");
ASSERTL0(m_config["fromxml"].as<string>().compare("NotSet") != 0,
"Need to specify fromxml=file.xml");
m_config["fromfld"] =
ConfigOption(false,"NotSet",
"Fld file from which to interpolate field");
ASSERTL0(m_config["fromfld"].as<string>().compare("NotSet") != 0,
"Need to specify fromfld=file.fld ");
m_config["clamptolowervalue"] =
ConfigOption(false, "-10000000",
"Lower bound for interpolation value");
m_config["clamptouppervalue"] =
ConfigOption(false, "10000000",
"Upper bound for interpolation value");
m_config["defaultvalue"] =
ConfigOption(false, "0",
"Default value if point is outside domain");
m_config["line"] =
ConfigOption(false, "NotSet",
"Specify a line of N points using "
"line=N,x0,y0,z0,z1,y1,z1");
m_config["plane"] =
ConfigOption(false, "NotSet",
"Specify a plane of N1 x N1 points using "
"plane=N1,N2,x0,y0,z0,z1,y1,z1,x2,y2,z2,x3,"
"y3,z3");
}
Nektar::Utilities::ProcessInterpPoints::~ProcessInterpPoints ( )
virtual

Definition at line 95 of file ProcessInterpPoints.cpp.

{
}

Member Function Documentation

static boost::shared_ptr<Module> Nektar::Utilities::ProcessInterpPoints::create ( FieldSharedPtr  f)
inlinestatic

Creates an instance of this class.

Definition at line 52 of file ProcessInterpPoints.h.

{
return MemoryManager<ProcessInterpPoints>::AllocateSharedPtr(f);
}
void Nektar::Utilities::ProcessInterpPoints::InterpolateFieldToPts ( vector< MultiRegions::ExpListSharedPtr > &  field0,
Array< OneD, Array< OneD, NekDouble > > &  pts,
NekDouble  clamp_low,
NekDouble  clamp_up,
NekDouble  def_value 
)
private

Definition at line 325 of file ProcessInterpPoints.cpp.

References ASSERTL0, and Nektar::Utilities::Module::m_f.

Referenced by Process().

{
int expdim = field0[0]->GetCoordim(0);
Array<OneD, NekDouble> coords(expdim), Lcoords(expdim);
int nq1 = pts[0].num_elements();
int elmtid, offset;
int r, f;
static int intpts = 0;
// resize data field
m_f->m_data.resize(field0.size());
for (f = 0; f < field0.size(); ++f)
{
m_f->m_data[f].resize(nq1);
}
for (r = 0; r < nq1; r++)
{
coords[0] = pts[0][r];
coords[1] = pts[1][r];
if (expdim == 3)
{
coords[2] = pts[2][r];
}
// Obtain Element and LocalCoordinate to interpolate
elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3);
if(elmtid >= 0)
{
offset = field0[0]->GetPhys_Offset(field0[0]->
GetOffset_Elmt_Id(elmtid));
for (f = 0; f < field0.size(); ++f)
{
NekDouble value;
value = field0[f]->GetExp(elmtid)->
StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset);
if ((boost::math::isnan)(value))
{
ASSERTL0(false, "new value is not a number");
}
else
{
value = (value > clamp_up)? clamp_up :
((value < clamp_low)? clamp_low :
value);
m_f->m_data[f][r] = value;
}
}
}
else
{
for (f = 0; f < field0.size(); ++f)
{
m_f->m_data[f][r] = def_value;
}
}
if (intpts%1000 == 0)
{
cout <<"." << flush;
}
intpts ++;
}
}
void Nektar::Utilities::ProcessInterpPoints::Process ( po::variables_map &  vm)
virtual

Write mesh to output file.

Implements Nektar::Utilities::Module.

Definition at line 99 of file ProcessInterpPoints.cpp.

References ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::Utilities::ePtsLine, Nektar::Utilities::ePtsPlane, Nektar::ParseUtils::GenerateUnOrderedVector(), InterpolateFieldToPts(), Nektar::Utilities::Module::m_config, Nektar::Utilities::Module::m_f, m_fromField, npts, Nektar::LibUtilities::NullFieldMetaDataMap, Nektar::Utilities::NullFieldPts, Nektar::SpatialDomains::MeshGraph::Read(), Vmath::Vmax(), and Vmath::Vmin().

{
if(m_f->m_verbose)
{
cout << "Processing point interpolation" << endl;
}
// Check for command line point specification if no .pts file specified
if(m_f->m_fieldPts == NullFieldPts)
{
int dim = 0;
if(m_config["line"].as<string>().compare("NotSet") != 0)
{
string help = m_config["line"].as<string>();
vector<NekDouble> values;
ASSERTL0(ParseUtils::GenerateUnOrderedVector(
m_config["line"].as<string>().c_str(),values),
"Failed to interpret line string");
ASSERTL0(values.size() > 3,
"line string should contain 2Dim+1 values "
"N,x0,y0,z0,x1,y1,z1");
m_f->m_fieldPts = MemoryManager<FieldPts>::AllocateSharedPtr();
dim = m_f->m_fieldPts->m_ptsDim = (values.size()-1)/2;
int npts = values[0];
m_f->m_fieldPts->m_npts.push_back(values[0]);
m_f->m_fieldPts->m_ptype = Utilities::ePtsLine;
m_f->m_fieldPts->m_pts = Array<OneD, Array<OneD, NekDouble> > (m_f->m_fieldPts->m_ptsDim);
for(int i = 0; i < dim; ++i)
{
m_f->m_fieldPts->m_pts[i] = Array<OneD,NekDouble>(npts);
}
for(int i = 0; i < npts; ++i)
{
m_f->m_fieldPts->m_pts[0][i] = values[1]
+ i/((NekDouble)(npts-1))*(values[dim+1] - values[1]);
if(dim > 1)
{
m_f->m_fieldPts->m_pts[1][i] = values[2]
+ i/((NekDouble)(npts-1))*(values[dim+2] - values[2]);
if(dim > 2)
{
m_f->m_fieldPts->m_pts[2][i] = values[3]
+ i/((NekDouble)(npts-1))*(values[dim+3] - values[3]);
}
}
}
}
else if(m_config["plane"].as<string>().compare("NotSet") != 0)
{
string help = m_config["plane"].as<string>();
vector<NekDouble> values;
ASSERTL0(ParseUtils::GenerateUnOrderedVector(
m_config["plane"].as<string>().c_str(),values),
"Failed to interpret plane string");
ASSERTL0(values.size() > 3,
"line string should contain 2Dim+1 values "
"N,x0,y0,z0,x1,y1,z1");
m_f->m_fieldPts = MemoryManager<FieldPts>::AllocateSharedPtr();
dim = m_f->m_fieldPts->m_ptsDim = (values.size()-1)/4;
int npts1 = values[0];
int npts2 = values[1];
m_f->m_fieldPts->m_npts.push_back(npts1);
m_f->m_fieldPts->m_npts.push_back(npts2);
m_f->m_fieldPts->m_ptype = Utilities::ePtsPlane;
m_f->m_fieldPts->m_pts = Array<OneD, Array<OneD, NekDouble> > (m_f->m_fieldPts->m_ptsDim);
for(int i = 0; i < dim; ++i)
{
m_f->m_fieldPts->m_pts[i] = Array<OneD,NekDouble>(npts1*npts2);
}
for(int j = 0; j < npts2; ++j)
{
for(int i = 0; i < npts1; ++i)
{
m_f->m_fieldPts->m_pts[0][i+j*npts1] =
(values[2] + i/((NekDouble)(npts1-1))*(values[dim+2] - values[2]))*(1.0-j/((NekDouble)(npts2-1))) +
(values[3*dim+2] + i/((NekDouble)(npts1-1))*(values[2*dim+2] - values[3*dim+2]))*(j/((NekDouble)(npts2-1)));
if(dim > 1)
{
m_f->m_fieldPts->m_pts[1][i+j*npts1] =
(values[3] + i/((NekDouble)(npts1-1))*(values[dim+3] - values[3]))*(1.0-j/((NekDouble)(npts2-1))) +
(values[3*dim+3] + i/((NekDouble)(npts1-1))*(values[2*dim+3] - values[3*dim+3]))*(j/((NekDouble)(npts2-1)));
if(dim > 2)
{
m_f->m_fieldPts->m_pts[2][i+j*npts1] =
(values[4] + i/((NekDouble)(npts1-1))*(values[dim+4] - values[4]))*(1.0-j/((NekDouble)(npts2-1))) +
(values[3*dim+4] + i/((NekDouble)(npts1-1))*(values[2*dim+4] - values[3*dim+4]))*(j/((NekDouble)(npts2-1)));
}
}
}
}
}
}
m_fromField = boost::shared_ptr<Field>(new Field());
std::vector<std::string> files;
// set up session file for from field
files.push_back(m_config["fromxml"].as<string>());
m_fromField->m_session = LibUtilities::SessionReader::
CreateInstance(0, 0, files);
// Set up range based on min and max of local parallel partition
SpatialDomains::DomainRangeShPtr rng = MemoryManager<SpatialDomains::DomainRange>::AllocateSharedPtr();
int coordim = m_f->m_fieldPts->m_ptsDim;
int npts = m_f->m_fieldPts->m_pts[0].num_elements();
Array<OneD, Array<OneD, NekDouble> > coords = m_f->m_fieldPts->m_pts;
switch(coordim)
{
case 3:
rng->doZrange = true;
rng->zmin = Vmath::Vmin(npts,coords[2],1);
rng->zmax = Vmath::Vmax(npts,coords[2],1);
case 2:
rng->doYrange = true;
rng->ymin = Vmath::Vmin(npts,coords[1],1);
rng->ymax = Vmath::Vmax(npts,coords[1],1);
case 1:
rng->doXrange = true;
rng->xmin = Vmath::Vmin(npts,coords[0],1);
rng->xmax = Vmath::Vmax(npts,coords[0],1);
break;
default:
ASSERTL0(false,"too many values specfied in range");
}
// setup rng parameters.
m_fromField->m_graph = SpatialDomains::MeshGraph::Read(m_fromField->m_session,rng);
// Read in local from field partitions
const SpatialDomains::ExpansionMap &expansions = m_fromField->m_graph->GetExpansions();
m_fromField->m_fld = MemoryManager<LibUtilities::FieldIO>
::AllocateSharedPtr(m_fromField->m_session->GetComm());
Array<OneD,int> ElementGIDs(expansions.size());
SpatialDomains::ExpansionMap::const_iterator expIt;
int i = 0;
for (expIt = expansions.begin(); expIt != expansions.end();
++expIt)
{
ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
}
string fromfld = m_config["fromfld"].as<string>();
m_fromField->m_fld->Import(fromfld,m_fromField->m_fielddef,
m_fromField->m_data,
LibUtilities::NullFieldMetaDataMap,
ElementGIDs);
int NumHomogeneousDir = m_fromField->m_fielddef[0]->m_numHomogeneousDir;
//----------------------------------------------
// Set up Expansion information to use mode order from field
m_fromField->m_graph->SetExpansions(m_fromField->m_fielddef);
int nfields = m_fromField->m_fielddef[0]->m_fields.size();
m_fromField->m_exp.resize(nfields);
m_fromField->m_exp[0] = m_fromField->SetUpFirstExpList(NumHomogeneousDir,true);
m_f->m_exp.resize(nfields);
// declare auxiliary fields.
for(i = 1; i < nfields; ++i)
{
m_fromField->m_exp[i] = m_fromField->AppendExpList(NumHomogeneousDir);
}
// load field into expansion in fromfield.
for(int j = 0; j < nfields; ++j)
{
for (i = 0; i < m_fromField->m_fielddef.size(); i++)
{
m_fromField->m_exp[j]->ExtractDataToCoeffs(
m_fromField->m_fielddef[i],
m_fromField->m_data[i],
m_fromField->m_fielddef[0]->m_fields[j],
m_fromField->m_exp[j]->UpdateCoeffs());
}
m_fromField->m_exp[j]->BwdTrans(m_fromField->m_exp[j]->GetCoeffs(),
m_fromField->m_exp[j]->UpdatePhys());
m_f->m_fieldPts->m_fields.push_back(m_fromField->m_fielddef[0]->m_fields[j]);
}
if(m_fromField->m_session->GetComm()->GetRank() == 0)
{
cout << "Interpolating [" << flush;
}
NekDouble clamp_low = m_config["clamptolowervalue"].as<NekDouble>();
NekDouble clamp_up = m_config["clamptouppervalue"].as<NekDouble>();
NekDouble def_value = m_config["defaultvalue"].as<NekDouble>();
InterpolateFieldToPts(m_fromField->m_exp, m_f->m_fieldPts->m_pts,
clamp_low, clamp_up, def_value);
if(m_fromField->m_session->GetComm()->GetRank() == 0)
{
cout << "]" << endl;
}
}

Member Data Documentation

ModuleKey Nektar::Utilities::ProcessInterpPoints::className
static
Initial value:
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eProcessModule, "interppoints"),
ProcessInterpPoints::create,
"Interpolates a set of points to another, requires fromfld and "
"fromxml to be defined, a line or plane of points can be "
"defined")

Definition at line 55 of file ProcessInterpPoints.h.

FieldSharedPtr Nektar::Utilities::ProcessInterpPoints::m_fromField
private

Definition at line 64 of file ProcessInterpPoints.h.

Referenced by Process().

Generated on Sun Mar 15 2015 00:12:10 for Nektar++ by   doxygen 1.8.1.2