Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
FldCalcBCs.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ContField1D.h>
#include <MultiRegions/ContField2D.h>
#include <MultiRegions/ContField3D.h>
#include <MultiRegions/ContField3DHomogeneous1D.h>
#include <MultiRegions/ContField3DHomogeneous2D.h>
#include <LocalRegions/SegExp.h>
#include <LocalRegions/Expansion2D.h>
#include <SpatialDomains/MeshGraph2D.h>
Include dependency graph for FldCalcBCs.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
void SetFields (SpatialDomains::MeshGraphSharedPtr &mesh, SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions, LibUtilities::SessionReaderSharedPtr &session, Array< OneD, MultiRegions::ExpListSharedPtr > &Exp, int nvariables)
Array< OneD, int > GetReflectionIndex (MultiRegions::ExpListSharedPtr Exp, int Ireg)
void Extractlayerdata (Array< OneD, int > Iregions, int coordim, SpatialDomains::MeshGraphSharedPtr &mesh, LibUtilities::SessionReaderSharedPtr &session, SpatialDomains::BoundaryConditions &bcs, Array< OneD, MultiRegions::ExpListSharedPtr > &infields, MultiRegions::ContField1DSharedPtr &outfieldx, MultiRegions::ContField1DSharedPtr &outfieldy, MultiRegions::ExpListSharedPtr &streak, bool symm, Array< OneD, int > Refindices, NekDouble alpha, NekDouble cr)
void Manipulate (Array< OneD, int > Iregions, int coordim, SpatialDomains::MeshGraphSharedPtr &mesh, SpatialDomains::BoundaryConditions &bcs, Array< OneD, MultiRegions::ExpListSharedPtr > &infields, Array< OneD, MultiRegions::ExpList1DSharedPtr > &outfieldx, Array< OneD, MultiRegions::ExpList1DSharedPtr > &outfieldy, MultiRegions::ExpListSharedPtr &streak)
void CalcNonLinearForcing (SpatialDomains::MeshGraphSharedPtr &mesh, LibUtilities::SessionReaderSharedPtr &session, string fieldfile, Array< OneD, MultiRegions::ExpListSharedPtr > &waveFields, MultiRegions::ExpListSharedPtr &streak, Array< OneD, int > Refindices, bool symm)
Array< OneD, int > GetReflectionIndex2D (MultiRegions::ExpListSharedPtr wavefield)
void WriteBcs (string variable, int region, string fieldfile, SpatialDomains::MeshGraphSharedPtr &mesh, MultiRegions::ContField1DSharedPtr &outregionfield)
void WriteFld (string outfile, SpatialDomains::MeshGraphSharedPtr &mesh, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)

Function Documentation

void CalcNonLinearForcing ( SpatialDomains::MeshGraphSharedPtr &  mesh,
LibUtilities::SessionReaderSharedPtr &  session,
string  fieldfile,
Array< OneD, MultiRegions::ExpListSharedPtr > &  waveFields,
MultiRegions::ExpListSharedPtr &  streak,
Array< OneD, int >  Refindices,
bool  symm 
)

Definition at line 1684 of file FldCalcBCs.cpp.

References Nektar::SpatialDomains::eNeumann, Vmath::Fill(), npts, Vmath::Smul(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vvtvp(), and Nektar::LibUtilities::Write().

Referenced by main().

{
int npts = waveFields[0]->GetPlane(0)->GetNpoints();
int ncoeffs = waveFields[0]->GetPlane(0)->GetNcoeffs();
Array<OneD, NekDouble> val(npts), der1(2*npts);
Array<OneD, NekDouble> der2 = der1 + npts;
Array<OneD, MultiRegions::ExpListSharedPtr> waveVelocities;
int nvel = waveFields.num_elements()-1;
waveVelocities = Array<OneD, MultiRegions::ExpListSharedPtr>(nvel);
//fill velocity fields (both coeffs, phys)
for( int i=0; i< nvel; i++)
{
waveVelocities[i] = waveFields[i];
Vmath::Vcopy(npts, waveFields[i]->GetPlane(0)->GetPhys(),1, waveVelocities[i]->GetPlane(0)->UpdatePhys(),1 );
Vmath::Vcopy(npts, waveFields[i]->GetPlane(1)->GetPhys(),1, waveVelocities[i]->GetPlane(1)->UpdatePhys(),1 );
Vmath::Vcopy(ncoeffs, waveFields[i]->GetPlane(0)->GetCoeffs(),1, waveVelocities[i]->GetPlane(0)->UpdateCoeffs(),1 );
Vmath::Vcopy(ncoeffs, waveFields[i]->GetPlane(1)->GetCoeffs(),1, waveVelocities[i]->GetPlane(1)->UpdateCoeffs(),1 );
}
MultiRegions::ExpListSharedPtr wavePressure;
wavePressure = waveFields[nvel];
//nvel=lastfield!!!
Vmath::Vcopy(npts, waveFields[nvel]->GetPlane(0)->GetPhys(),1, wavePressure->GetPlane(0)->UpdatePhys(),1 );
Vmath::Vcopy(npts, waveFields[nvel]->GetPlane(1)->GetPhys(),1, wavePressure->GetPlane(1)->UpdatePhys(),1 );
Vmath::Vcopy(ncoeffs, waveFields[nvel]->GetPlane(0)->GetCoeffs(),1, wavePressure->GetPlane(0)->UpdateCoeffs(),1 );
Vmath::Vcopy(ncoeffs, waveFields[nvel]->GetPlane(1)->GetCoeffs(),1, wavePressure->GetPlane(1)->UpdateCoeffs(),1 );
static int projectfield = -1;
// Set project field to be first field that has a Neumann
// boundary since this not impose any condition on the vertical boundaries
// Othersise set to zero.
if(projectfield == -1)
{
Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
for(int i = 0; i < waveVelocities.num_elements(); ++i)
{
BndConds = waveVelocities[i]->GetBndConditions();
for(int j = 0; j < BndConds.num_elements(); ++j)
{
if(BndConds[j]->GetBoundaryConditionType() == SpatialDomains::eNeumann)
{
projectfield = i;
break;
}
}
if(projectfield != -1)
{
break;
}
}
if(projectfield == -1)
{
cout << "using first field to project non-linear forcing which imposes a Dirichlet condition" << endl;
projectfield = 0;
}
}
// determine inverse of area normalised field.
wavePressure->GetPlane(0)->BwdTrans(wavePressure->GetPlane(0)->GetCoeffs(),
wavePressure->GetPlane(0)->UpdatePhys());
wavePressure->GetPlane(1)->BwdTrans(wavePressure->GetPlane(1)->GetCoeffs(),
wavePressure->GetPlane(1)->UpdatePhys());
// Determine normalisation of pressure so that |P|/A = 1
NekDouble norm = 0, l2;
l2 = wavePressure->GetPlane(0)->L2(wavePressure->GetPlane(0)->GetPhys());
norm = l2*l2;
l2 = wavePressure->GetPlane(1)->L2(wavePressure->GetPlane(0)->GetPhys());
norm += l2*l2;
Vmath::Fill(2*npts,1.0,der1,1);
NekDouble area = waveVelocities[0]->GetPlane(0)->PhysIntegral(der1);
norm = sqrt(area/norm);
// Get hold of arrays.
waveVelocities[0]->GetPlane(0)->BwdTrans(waveVelocities[0]->GetPlane(0)->GetCoeffs(),waveVelocities[0]->GetPlane(0)->UpdatePhys());
Array<OneD, NekDouble> u_real = waveVelocities[0]->GetPlane(0)->UpdatePhys();
Vmath::Smul(npts,norm,u_real,1,u_real,1);
waveVelocities[0]->GetPlane(1)->BwdTrans(waveVelocities[0]->GetPlane(1)->GetCoeffs(),waveVelocities[0]->GetPlane(1)->UpdatePhys());
Array<OneD, NekDouble> u_imag = waveVelocities[0]->GetPlane(1)->UpdatePhys();
Vmath::Smul(npts,norm,u_imag,1,u_imag,1);
waveVelocities[1]->GetPlane(0)->BwdTrans(waveVelocities[1]->GetPlane(0)->GetCoeffs(),waveVelocities[1]->GetPlane(0)->UpdatePhys());
Array<OneD, NekDouble> v_real = waveVelocities[1]->GetPlane(0)->UpdatePhys();
Vmath::Smul(npts,norm,v_real,1,v_real,1);
waveVelocities[1]->GetPlane(1)->BwdTrans(waveVelocities[1]->GetPlane(1)->GetCoeffs(),waveVelocities[1]->GetPlane(1)->UpdatePhys());
Array<OneD, NekDouble> v_imag = waveVelocities[1]->GetPlane(1)->UpdatePhys();
Vmath::Smul(npts,norm,v_imag,1,v_imag,1);
// Calculate non-linear terms for x and y directions
// d/dx(u u* + u* u)
Vmath::Vmul (npts,u_real,1,u_real,1,val,1);
Vmath::Vvtvp(npts,u_imag,1,u_imag,1,val,1,val,1);
Vmath::Smul (npts,2.0,val,1,val,1);
waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
// d/dy(v u* + v* u)
Vmath::Vmul (npts,u_real,1,v_real,1,val,1);
Vmath::Vvtvp(npts,u_imag,1,v_imag,1,val,1,val,1);
Vmath::Smul (npts,2.0,val,1,val,1);
waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
Vmath::Vadd(npts,der1,1,der2,1,der1,1);
NekDouble rho = session->GetParameter("RHO");
NekDouble Re = session->GetParameter("RE");
cout<<"Re="<<Re<<endl;
NekDouble waveForceMag = rho*rho*std::pow(Re,-1./3.);
Array<OneD, Array<OneD, NekDouble > > vwiForcing;
vwiForcing = Array<OneD, Array<OneD, NekDouble> > (2);
vwiForcing[0] = Array<OneD, NekDouble> (2*ncoeffs);
for(int i = 1; i < 2; ++i)
{
vwiForcing[i] = vwiForcing[i-1] + ncoeffs;
}
if(projectfield!=0)
{
waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,vwiForcing[0]);
}
else
{
waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,vwiForcing[0]);
}
cout<<"waveforcemag="<<waveForceMag<<endl;
Vmath::Smul(ncoeffs,-waveForceMag, vwiForcing[0],1, vwiForcing[0],1);
// d/dx(u v* + u* v)
waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
// d/dy(v v* + v* v)
Vmath::Vmul(npts,v_real,1,v_real,1,val,1);
Vmath::Vvtvp(npts,v_imag,1,v_imag,1,val,1,val,1);
Vmath::Smul (npts,2.0,val,1,val,1);
waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
Vmath::Vadd(npts,der1,1,der2,1,der1,1);
if(projectfield!=0)
{
waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,vwiForcing[1]);
}
else
{
waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,vwiForcing[1]);
}
Vmath::Smul(ncoeffs,- waveForceMag, vwiForcing[1],1, vwiForcing[1],1);
if(symm==true)
{
cout<<"symmetrization"<<endl;
waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(vwiForcing[0],der1);
for(int i = 0; i < npts; ++i)
{
val[i] = 0.5*(der1[i] - der1[Refindices[i]]);
}
waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, vwiForcing[0]);
waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(vwiForcing[1],der1);
for(int i = 0; i < npts; ++i)
{
val[i] = 0.5*(der1[i] - der1[Refindices[i]]);
}
waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, vwiForcing[1]);
}
// dump output
Array<OneD, std::string> variables(2);
Array<OneD, Array<OneD, NekDouble> > outfield(2);
variables[0] = "u";
variables[1] = "v";
outfield[0] = vwiForcing[0];
outfield[1] = vwiForcing[1];
string sessionName = fieldfile.substr(0,fieldfile.find_last_of("."));
std::string outname = sessionName + ".vwi";
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
= waveVelocities[0]->GetPlane(0)->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
//Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(outfield.num_elements());
string var;
for(int j=0; j< vwiForcing.num_elements(); ++j)
{
//fieldcoeffs =vwiForcing;
for(int i=0; i< FieldDef.size(); i++)
{
//var = vSession->GetVariable(j);
var = variables[j];
FieldDef[i]->m_fields.push_back(var);
waveVelocities[0]->GetPlane(0)->AppendFieldData(FieldDef[i], FieldData[i], vwiForcing[j]);
}
}
LibUtilities::Write(outname,FieldDef,FieldData);
}
void Extractlayerdata ( Array< OneD, int >  Iregions,
int  coordim,
SpatialDomains::MeshGraphSharedPtr &  mesh,
LibUtilities::SessionReaderSharedPtr &  session,
SpatialDomains::BoundaryConditions &  bcs,
Array< OneD, MultiRegions::ExpListSharedPtr > &  infields,
MultiRegions::ContField1DSharedPtr &  outfieldx,
MultiRegions::ContField1DSharedPtr &  outfieldy,
MultiRegions::ExpListSharedPtr &  streak,
bool  symm,
Array< OneD, int >  Refindices,
NekDouble  alpha,
NekDouble  cr 
)

Definition at line 692 of file FldCalcBCs.cpp.

References ASSERTL0, Nektar::StdRegions::StdExpansion::BwdTrans(), Nektar::MultiRegions::eS, Nektar::MultiRegions::eX, Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), Nektar::LocalRegions::Expansion1D::GetLeftAdjacentElementEdge(), Nektar::StdRegions::StdExpansion::GetPointsKeys(), Nektar::iterator, Vmath::Smul(), Vmath::Vcopy(), Vmath::Vmul(), and Vmath::Zero().

Referenced by main().

{
//1 I regions is expected: (the layer is the last region)
ASSERTL0(Iregions.num_elements()==3, "something wrong with the number of I layers");
//int Iup =Iregions[0];
int Ireg =Iregions[2];
cout<<"layer region="<<Ireg<<endl;
//take the pressure from infields
int pfield;
if(infields.num_elements()==2)
{
pfield =1;
}
else
{
pfield =3;
}
MultiRegions::ExpListSharedPtr pressure =infields[pfield];
Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =infields[0]->GetBndCondExpansions();
MultiRegions::ExpListSharedPtr Ilayer = Iexp[Ireg];
int nq1D = Ilayer->GetPlane(0)->GetTotPoints();
//decomment int nq1D = Ilayer->GetTotPoints();
Array<OneD,NekDouble> x0d(nq1D);
Array<OneD,NekDouble> x1d(nq1D);
Array<OneD,NekDouble> x2d(nq1D);
Ilayer->GetPlane(0)->GetCoords(x0d,x1d,x2d);
//decomment Ilayer->GetCoords(x0d,x1d,x2d);
//start
//take the streak data along the layers (getbondcondexpansions does not work!!!)
//Array<OneD, MultiRegions::ExpListSharedPtr> Istreak = streak->GetBndCondExpansions();
//Array<OneD, MultiRegions::ExpListSharedPtr> Istreak =
// infields[1]->GetPlane(0)->GetBndCondExpansions();
//MultiRegions::ExpListSharedPtr Istreakup = Istreak[Iup];
//MultiRegions::ExpListSharedPtr Istreakdown = Istreak[Idown];
int np = streak->GetTotPoints();
Array<OneD, NekDouble> gradx(np);
Array<OneD, NekDouble> grady(np);
int totcoeffs = streak->GetNcoeffs();
Array<OneD, NekDouble> gradxcoeffs(totcoeffs);
Array<OneD, NekDouble> gradycoeffs(totcoeffs);
streak->PhysDeriv(streak->GetPhys(), gradx, grady);
streak->FwdTrans_IterPerExp(gradx, gradxcoeffs);
streak->FwdTrans_IterPerExp(grady, gradycoeffs);
Array<OneD, Array<OneD, NekDouble> > normals;
normals = Array<OneD, Array<OneD, NekDouble> >(2);
Array<OneD, Array<OneD, NekDouble> > tangents;
tangents = Array<OneD, Array<OneD, NekDouble> >(2);
Array<OneD, NekDouble> nx(nq1D);
Array<OneD, NekDouble> ny(nq1D);
Array<OneD, NekDouble> tx(nq1D);
Array<OneD, NekDouble> ty(nq1D);
const SpatialDomains::BoundaryRegionCollection &bregions
= bcs.GetBoundaryRegions();
SpatialDomains::Composite Icompreg;
SpatialDomains::SegGeomSharedPtr segmentGeomreg;
SpatialDomains::ElementEdgeVectorSharedPtr elementreg;
StdRegions::Orientation orientreg;
SpatialDomains::BoundaryRegion::iterator regIt;
Array<OneD, unsigned int> bmapreg;
//fields to store the Re Im parts in physical space:
Array<OneD,NekDouble> Rephysreg (nq1D);
Array<OneD,NekDouble> Imphysreg (nq1D);
//jacobian for the layer, std_deriv
Array<OneD, NekDouble> Jac (nq1D);
Array<OneD, NekDouble> dPrestd (nq1D);
//Array<OneD, NekDouble> stphysup (nq1D);
Array<OneD, NekDouble> stphysreg (nq1D);
Array<OneD, NekDouble> stphysgradxreg (nq1D);
Array<OneD, NekDouble> stphysgradyreg (nq1D);
//metrics definitions
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > deriv;
Array<OneD, NekDouble> dxde1 (nq1D);
Array<OneD, NekDouble> dyde1 (nq1D);
Array<OneD, NekDouble> dxde2 (nq1D);
Array<OneD, NekDouble> dyde2 (nq1D);
Array<OneD, NekDouble> gmat0 (nq1D);
Array<OneD, NekDouble> gmat3 (nq1D);
Array<OneD, NekDouble> gmat1 (nq1D);
Array<OneD, NekDouble> gmat2 (nq1D);
Array<OneD, int> Elmtid(Ilayer->GetPlane(0)->GetExpSize());
Array<OneD, int> Edgeid(Ilayer->GetPlane(0)->GetExpSize());
//nummodes*nedges_layerdown
int Nregcoeffs;
for(
regIt= bregions.find(Ireg)->second->begin();
regIt !=bregions.find(Ireg)->second->end();
++regIt)
{
Icompreg=regIt->second;
//assuming that ncoeffs is equal for all the edges (k=0)
int nummodes = Ilayer->GetExp(0)->GetNcoeffs();
Array<OneD, int> offsetup(Icompreg->size());
offsetup[0]=0;
int offsetregIExp;
//array coeffs:
Nregcoeffs = Icompreg->size()*nummodes;
Array<OneD, NekDouble> Recoeffsreg (Nregcoeffs);
Array<OneD, NekDouble> Imcoeffsreg (Nregcoeffs);
Array<OneD, NekDouble> stcoeffsreg (Nregcoeffs);
Array<OneD, NekDouble> stgradxcoeffsreg (Nregcoeffs);
Array<OneD, NekDouble> stgradycoeffsreg (Nregcoeffs);
//define nq points for each edge
int nqedge =nq1D/Icompreg->size();
Array<OneD,NekDouble> x0edge(nqedge);
Array<OneD,NekDouble> x1edge(nqedge);
Array<OneD, NekDouble> Rephysedgereg(nqedge);
Array<OneD, NekDouble> Imphysedgereg(nqedge);
Array<OneD, NekDouble> stgradxcoeffsedge(nqedge);
Array<OneD, NekDouble> stgradycoeffsedge(nqedge);
Array<TwoD, const NekDouble> gmat;
Array<OneD, NekDouble> gmat0_2D (nqedge*nqedge);
Array<OneD, NekDouble> gmat3_2D (nqedge*nqedge);
Array<OneD, NekDouble> gmat1_2D (nqedge*nqedge);
Array<OneD, NekDouble> gmat2_2D (nqedge*nqedge);
for(int f=0; f<2; ++f)
{
normals[f] = Array<OneD, NekDouble>(nqedge);
tangents[f] = Array<OneD, NekDouble>(nqedge);
}
for(int k=0; k< Icompreg->size(); k++)//loop over segments of each layer
{
if(
!(segmentGeomreg = boost::dynamic_pointer_cast<
SpatialDomains::SegGeom>((*Icompreg)[k]))
)
{ ASSERTL0(false, "dynamic cast to a SegGeom failed"); }
int EIDreg = segmentGeomreg->GetEid();
offsetregIExp = Ilayer->GetPlane(0)->GetCoeff_Offset(k);
//cout<<"edge="<<EIDreg<<" CHECK OFFSET="<<offsetregIExp<<endl;
elementreg = boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(mesh)
->GetElementsFromEdge(segmentGeomreg);
int elmtidreg = ((*elementreg)[0]->m_Element)->GetGlobalID();
int dimelmtreg = ((*elementreg)[0]->m_Element)->GetNumEdges();
Elmtid[k] = elmtidreg;
Edgeid[k] = EIDreg;
orientreg =(boost::dynamic_pointer_cast<SpatialDomains::Geometry2D>
((*elementreg)[0]->m_Element))->GetEorient((*elementreg)[0]
->m_EdgeIndx);
gmat = pressure->GetPlane(0)->GetExp(elmtidreg)->GetMetricInfo()->GetGmat(
pressure->GetPlane(0)->GetExp(elmtidreg)->GetPointsKeys());
int nq2D = nqedge*nqedge;
// setup map between global and local ids
map<int, int>EdgeGIDreg;
int id,cnt,f;
for(cnt=f=0; f<dimelmtreg; f++)
{
id = (boost::dynamic_pointer_cast<SpatialDomains::Geometry2D>
((*elementreg)[0]->m_Element))->GetEid(f);
EdgeGIDreg[id]=cnt++;
}
//cout<<" GIDreg="<<elmtidreg<<" num edges in this element="<<dimelmtreg<<endl;
int localidedge = EdgeGIDreg.find(EIDreg)->second;
//cout<<"locaidedge from m_EdgeInx ="<<localidedge<<endl;
localidedge = (*elementreg)[0]->m_EdgeIndx;
//cout<<"check map local id of Eid="<<EIDreg<<" is="<<localidedge<<endl;
//set bmaps -------------
Array<OneD, unsigned int > bmapedgereg;
Array<OneD, int > signreg;
pressure->GetExp(elmtidreg)->GetEdgeToElementMap(EdgeGIDreg.find(EIDreg)->second,
orientreg,bmapedgereg,signreg);
//cout<<"SIGN down="<<signreg[0]<<endl;
//Extracting the coeffs...
Array<OneD, NekDouble> Recoeffs = pressure->GetPlane(0)->GetCoeffs();
Array<OneD, NekDouble> Imcoeffs = pressure->GetPlane(1)->GetCoeffs();
Array<OneD, NekDouble> stcoeffs = streak->GetCoeffs();
//OFFSET EDGE REFERING TO ELEMENT EXPLIST MISSING
int Reoffsetreg = pressure->GetPlane(0)->GetCoeff_Offset(elmtidreg);
int Imoffsetreg = pressure->GetPlane(1)->GetCoeff_Offset(elmtidreg);
int stoffsetreg = streak->GetCoeff_Offset(elmtidreg);
//cout<<" Re offsetdown="<<Reoffsetreg<<" Im offset="<<Imoffsetreg<<endl;
//hypothesis: 2 planes (0 Re, 1 Im)
Array<OneD, NekDouble> Recoeffsedgereg(bmapedgereg.num_elements());
Array<OneD, NekDouble> Imcoeffsedgereg(bmapedgereg.num_elements());
Array<OneD, NekDouble> stcoeffsedgereg(bmapedgereg.num_elements());
//grad edge coeffs
//cout<<" Recoeffsedgeup num elements="<<Recoeffsedgereg.num_elements()
//<<" n bmapedge="<<bmapedgereg.num_elements()<<endl;
for(int d=0; d<bmapedgereg.num_elements(); d++)
{
//pressure
Recoeffsedgereg[d]=Recoeffs[Reoffsetreg+bmapedgereg[d]];
Imcoeffsedgereg[d]=Imcoeffs[Imoffsetreg+bmapedgereg[d]];
Recoeffsreg[offsetregIExp +d] = Recoeffsedgereg[d];
Imcoeffsreg[offsetregIExp +d] = Imcoeffsedgereg[d];
//streak coeffs
stcoeffsedgereg[d] = stcoeffs[stoffsetreg + bmapedgereg[d]];
stgradxcoeffsedge[d] = gradxcoeffs[stoffsetreg + bmapedgereg[d]];
stgradycoeffsedge[d] = gradycoeffs[stoffsetreg + bmapedgereg[d]];
stcoeffsreg[offsetregIExp +d] = stcoeffsedgereg[d];
stgradxcoeffsreg[offsetregIExp +d] = stgradxcoeffsedge[d];
stgradycoeffsreg[offsetregIExp +d] = stgradycoeffsedge[d];
//cout<<"results: Re="<<Recoeffsedgereg[d]<<" Im="<<Imcoeffsedgereg[d]<<endl;
}
//store the jacobian, std_deriv
Array<OneD, NekDouble> jacedge (nqedge);
Array<OneD, NekDouble> Pre_edge (nqedge);
Array<OneD, NekDouble> dPre_edge (nqedge);
jacedge = Ilayer->GetPlane(0)->GetExp(k)->GetMetricInfo()->GetJac(Ilayer->GetPlane(0)->GetExp(k)->GetPointsKeys());
//check if the metrix is the same for streak and Ilayer obj
StdRegions::StdExpansion1DSharedPtr edgestdexp;
edgestdexp = boost::dynamic_pointer_cast<StdRegions::StdExpansion1D> (
Ilayer->GetPlane(0)->GetExp(k) );
edgestdexp->BwdTrans(Recoeffsedgereg, Pre_edge);
//cout<<"call PhysTensorDERIV......"<<endl;
edgestdexp->StdPhysDeriv(Pre_edge,
dPre_edge);
for(int q=0; q<nqedge; q++)
{
Jac[k*nqedge +q]= jacedge[q];
dPrestd[k*nqedge +q] = dPre_edge[q];
}
//cout<<"Nlayercoeffs="<<Recoeffsreg.num_elements()<<" Nlayerphys="<<Rephysreg.num_elements()<<endl;
//cout<<"extract normal for edge="<<k<<endl;
//cout<<"extract tangent for edge="<<k<<endl;
//k is the number of the edge according to the composite list
LocalRegions::Expansion1DSharedPtr edgeexp =
boost::dynamic_pointer_cast<LocalRegions::Expansion1D>
(Ilayer->GetPlane(0)->GetExp(k) );
int localEid = edgeexp->GetLeftAdjacentElementEdge();
normals = edgeexp->
GetLeftAdjacentElementExp()->GetEdgeNormal(localEid);
LocalRegions::SegExpSharedPtr bndSegExp =
boost::dynamic_pointer_cast<LocalRegions::SegExp>(Ilayer->GetPlane(0)->GetExp(k));
// This function call has be deprecated -- cc
//tangents = (bndSegExp)->GetMetricInfo()->GetEdgeTangent();
int physoffsetregIExp = Ilayer->GetPlane(0)->GetPhys_Offset(k);
for(int e=0; e< nqedge; e++)
{
//nx[k*nqedge +e] = normals[0][e];
//ny[k*nqedge +e] = normals[1][e];
//tx[k*nqedge +e] = tangents[0][e];
//ty[k*nqedge +e] = tangents[1][e];
//cout<<"offset="<<offsetregIExp<<" nx="<<normals[0][e]<<" ny="<<normals[1][e]<<" e="<<e<<endl;
nx[physoffsetregIExp +e] = normals[0][e];
ny[physoffsetregIExp +e] = normals[1][e];
tx[physoffsetregIExp +e] = tangents[0][e];
ty[physoffsetregIExp +e] = tangents[1][e];
//cout<<"offset="<<offsetregIExp<<" nx="<<normals[0][e]<<" ny="<<normals[1][e]<<" e="<<e<<endl;
//cout<<offsetregIExp +e<<" "<<x0d[offsetregIExp +e]<<" nx="<<nx[offsetregIExp +e]<<" ny="<<ny[offsetregIExp +e]<<" e="<<e<<endl;
//cout<<"tannn tx="<<tangents[0][e]<<" ty="<<tangents[1][e]<<endl;
}
//extract metrics factors...
//int nq2D= nqedge*nqedge;
LibUtilities::PointsKeyVector ptsKeys = pressure->GetPlane(0)->GetExp(elmtidreg)->GetPointsKeys();
deriv = pressure->GetPlane(0)->GetExp(elmtidreg)->GetMetricInfo()->GetDeriv(ptsKeys);
Array<OneD, NekDouble> derivelmt(nq2D);
int offsetregphys = Ilayer->GetPlane(0)->GetPhys_Offset(k);
Array<OneD, NekDouble> deriv_i(nqedge);
Vmath::Vcopy(nq2D, &(deriv[0][0][0]),1, &(derivelmt[0]),1);
pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dxde1[offsetregphys]),1);
Vmath::Zero(nqedge, deriv_i,1);
Vmath::Zero(nq2D, derivelmt,1);
Vmath::Vcopy(nq2D, &(deriv[1][1][0]),1, &(derivelmt[0]),1);
pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dyde2[offsetregphys]),1);
Vmath::Zero(nqedge, deriv_i,1);
Vmath::Zero(nq2D, derivelmt,1);
Vmath::Vcopy(nq2D, &(deriv[0][1][0]),1, &(derivelmt[0]),1);
pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dyde1[offsetregphys]),1);
Vmath::Zero(nqedge, deriv_i,1);
Vmath::Zero(nq2D, derivelmt,1);
Vmath::Vcopy(nq2D, &(deriv[1][0][0]),1, &(derivelmt[0]),1);
pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dxde2[offsetregphys]),1);
Array<OneD, NekDouble> gmati (nqedge);
Vmath::Vcopy(nq2D, &(gmat[0][0]),1, &(gmat0_2D[0]),1);
Vmath::Vcopy(nq2D, &(gmat[3][0]),1, &(gmat3_2D[0]),1);
Vmath::Vcopy(nq2D, &(gmat[1][0]),1, &(gmat1_2D[0]),1);
Vmath::Vcopy(nq2D, &(gmat[2][0]),1, &(gmat2_2D[0]),1);
pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat0_2D, gmati);
Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat0[offsetregphys]),1);
Vmath::Zero(nqedge, gmati,1);
pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat1_2D, gmati);
Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat1[offsetregphys]),1);
Vmath::Zero(nqedge, gmati,1);
pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat2_2D, gmati);
Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat2[offsetregphys]),1);
Vmath::Zero(nqedge, gmati,1);
pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat3_2D, gmati);
Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat3[offsetregphys]),1);
}
//bwd transform:
Ilayer->GetPlane(0)->BwdTrans(Recoeffsreg, Rephysreg);
Ilayer->GetPlane(1)->BwdTrans(Imcoeffsreg, Imphysreg);
//streak phys values:
Ilayer->GetPlane(0)->BwdTrans(stcoeffsreg, stphysreg);
Ilayer->GetPlane(0)->BwdTrans_IterPerExp(stgradxcoeffsreg, stphysgradxreg);
Ilayer->GetPlane(0)->BwdTrans_IterPerExp(stgradycoeffsreg, stphysgradyreg);
}
//TANGENTS WRONG WITH 3DHOMOGENEOUS 1D...
//calculate the curvature:
//attempt to smooth the tangent components:
Array<OneD, NekDouble> tcoeffs (Nregcoeffs,0.0);
//outfieldx->FwdTrans(tx, tcoeffs);
//outfieldx->BwdTrans(tcoeffs, tx);
//Vmath::Zero(Nregcoeffs, tcoeffs,1);
//outfieldx->FwdTrans(ty, tcoeffs);
//outfieldx->BwdTrans(tcoeffs, ty);
//Vmath::Zero(Nregcoeffs, tcoeffs,1);
Array<OneD, NekDouble> dtx(nq1D,0.0);
Array<OneD, NekDouble> dty(nq1D,0.0);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, tx, dtx);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, ty, dty);
//attempt to smooth the tangent derivative:
outfieldx->FwdTrans_IterPerExp(dtx, tcoeffs);
outfieldx->BwdTrans_IterPerExp(tcoeffs, dtx);
Vmath::Zero(Nregcoeffs, tcoeffs,1);
outfieldx->FwdTrans_IterPerExp(dty, tcoeffs);
outfieldx->BwdTrans_IterPerExp(tcoeffs, dty);
Array<OneD, NekDouble> f_z(nq1D,0.0);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, x1d, f_z);
Array<OneD, NekDouble> fcoeffs(Nregcoeffs,0.0);
outfieldx->FwdTrans(f_z, fcoeffs);
outfieldx->BwdTrans(fcoeffs, f_z);
Vmath::Zero(Nregcoeffs, fcoeffs,1);
Array<OneD, NekDouble> f_zz(nq1D,0.0);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, f_z, f_zz);
outfieldx->FwdTrans(f_zz, fcoeffs);
outfieldx->BwdTrans(fcoeffs, f_zz);
Array<OneD, NekDouble> delta(nq1D);
Array<OneD, NekDouble> curv_unsigned(nq1D,0.0);
Array<OneD, NekDouble> curv(nq1D,0.0);
//cout<<"x y tx ty dtx dty curv"<<endl;
for(int t=0; t<nq1D; t++)
{
curv_unsigned[t]=sqrt(dtx[t]*dtx[t] +dty[t]*dty[t]);
delta[t] = 1+f_z[t]*f_z[t];
//Hall's definition...
curv[t] = f_zz[t]/sqrt(delta[t]*delta[t]*delta[t]);
//cout<<setw(13)<<x0d[t]<<" "<<setw(13)<<x1d[t]<<" "<<setw(13)<<tx[t]<<setw(13)<<" "<<ty[t]<<" "
//<<nx[t]<<" "<<ny[t]<<" "<<setw(13)<<dtx[t]<<" "<<dty[t]<<" "<<curv[t]<<endl;
}
//attempt to smooth the curvature
Array<OneD, NekDouble> curv_coeffs (Nregcoeffs);
outfieldx->FwdTrans(curv, curv_coeffs);
outfieldx->BwdTrans(curv_coeffs, curv);
//P square before norm
Array<OneD, NekDouble> P2reg (nq1D, 0.0);
for(int e=0; e< nq1D; e++)
{
P2reg[e] = Rephysreg[e]*Rephysreg[e] + Imphysreg[e]*Imphysreg[e];
}
//calculate the s variable:
Array<OneD, NekDouble> s0d(nq1D,0.0);
/*
Array<OneD, NekDouble> tmp(nq1D);
for( int e=0; e<nq1D; e++)
{
tmp[e] = sqrt(1+f_z[e]*f_z[e]);
}
*/
//normalise the pressure norm*sqrt[ (\int Psquare)/Area]=1 =>
NekDouble norm2D;
//////////////////////////////////////////////////2D norm
Array<OneD, NekDouble> P2_2D(np, 0.0);
for(int h=0; h< np; h++)
{
P2_2D[h] = pressure->GetPlane(0)->GetPhys()[h]*pressure->GetPlane(0)->GetPhys()[h]
+ pressure->GetPlane(1)->GetPhys()[h]*pressure->GetPlane(1)->GetPhys()[h];
}
NekDouble intP2_2D = pressure->GetPlane(0)->PhysIntegral(P2_2D);
NekDouble tmp = pressure->GetPlane(0)->L2(pressure->GetPlane(0)->GetPhys());
norm2D = tmp*tmp;
tmp = pressure->GetPlane(1)->L2(pressure->GetPlane(1)->GetPhys());
norm2D += tmp*tmp;
Array<OneD, NekDouble> I (2*np,1.0);
//Vmath::Fill(2*np,1.0,I,1);
NekDouble Area = pressure->GetPlane(0)->PhysIntegral(I);
norm2D = sqrt(Area/norm2D);
Array<OneD, NekDouble> I_int(np,1.0);
NekDouble Area1 = pressure->GetPlane(0)->PhysIntegral(I_int);
NekDouble normint2D = sqrt(Area1/intP2_2D);
cout<<"norm2D="<<norm2D<<" area="<<Area1<<" intP2_2D="<<intP2_2D<<" normint2D="<<normint2D<<endl;
///////////////////////////////////////////////////////////1D norm
NekDouble norm;
Array<OneD, NekDouble> sqrtlen(nq1D);
for(int u=0; u<nq1D; u++)
{
sqrtlen[u] = sqrt(1+f_z[u]*f_z[u]);
s0d[u] = Ilayer->GetPlane(0)->PhysIntegral(sqrtlen);
}
NekDouble length = Ilayer->GetPlane(0)->PhysIntegral(sqrtlen);
NekDouble int1D = Ilayer->GetPlane(0)->PhysIntegral(P2reg);
norm = sqrt(length/int1D);
NekDouble scal = Area/length;
cout<<"norm1D="<<norm<<" norm2D/norm1D="<<norm2D/norm<<endl;
cout<<"scal="<<scal<<endl;
//norm*pressure
Vmath::Smul(nq1D,norm,Rephysreg,1,Rephysreg,1);
Vmath::Smul(nq1D,norm,Imphysreg,1,Imphysreg,1);
//P square after norm
Array<OneD, NekDouble> P2reg_aft (nq1D, 0.0);
for(int e=0; e< nq1D; e++)
{
P2reg_aft[e] = Rephysreg[e]*Rephysreg[e] + Imphysreg[e]*Imphysreg[e];
}
//print out the fields
//cout<<"derivative of the pressure"<<endl;
Array<OneD, NekDouble> dP_re = Array<OneD, NekDouble>(nq1D,0.0);
Array<OneD, NekDouble> dP_im = Array<OneD, NekDouble>(nq1D,0.0);
Array<OneD, NekDouble> dP_square2 = Array<OneD, NekDouble>(nq1D,0.0);
//attempt to smooth the normal
//Array<OneD, NekDouble> ncoeffs(Nregcoeffs,0.0);
//outfieldx->FwdTrans(nx, ncoeffs);
//outfieldx->BwdTrans(ncoeffs, nx);
//Vmath::Zero(Nregcoeffs, ncoeffs,1);
//outfieldx->FwdTrans(ny, ncoeffs);
//outfieldx->BwdTrans(ncoeffs, ny);
//Vmath::Zero(Nregcoeffs, ncoeffs,1);
Array<OneD, NekDouble> dUreg (nq1D);
//Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eN, stphysreg,dUreg);
for(int w=0; w<nq1D; w++)
{
dUreg[w] = nx[w]*stphysgradxreg[w] +ny[w]*stphysgradyreg[w];
}
Array<OneD, NekDouble> dUcoeffs(Nregcoeffs,0.0);
//outfieldx->FwdTrans(dUreg, dUcoeffs);
//outfieldx->BwdTrans(dUcoeffs, dUreg);
Array<OneD, NekDouble> jaccoeffs(Nregcoeffs);
Array<OneD, NekDouble> Jacreg (nq1D);
Array<OneD, NekDouble> dPrecoeffs (Nregcoeffs);
//attempt to smooth the jac and the dPre
outfieldx->FwdTrans(Jac, jaccoeffs);
outfieldx->BwdTrans(jaccoeffs, Jacreg);
//outfieldx->FwdTrans(dPre, dPrecoeffs);
//outfieldx->BwdTrans(dPrecoeffs, dPre);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS,Rephysreg,dP_re);
Ilayer->GetPlane(1)->PhysDeriv(MultiRegions::eS,Imphysreg,dP_im);
//attempt to smooth the derivative:
Array<OneD, NekDouble> dP_re_coeffs (Nregcoeffs,0.0);
Array<OneD, NekDouble> dP_im_coeffs (Nregcoeffs,0.0);
Array<OneD, NekDouble> dP_imtest (nq1D);
outfieldx->FwdTrans(dP_re, dP_re_coeffs);
outfieldx->FwdTrans(dP_im, dP_im_coeffs);
outfieldx->BwdTrans(dP_re_coeffs, dP_re);
outfieldx->BwdTrans(dP_im_coeffs, dP_im);
Array<OneD, NekDouble> dP_square = Array<OneD, NekDouble>(nq1D, 0.0);
cout<<"x"<<" P_re"<<" dP_re"<<" streak"<<" dstreak"<<" pjump"<<endl;
for(int s=0; s<nq1D; s++)
{
dP_square[s] = dP_re[s]*dP_re[s] +dP_im[s]*dP_im[s];
}
//attempt to smooth the pressure
Array<OneD, NekDouble> dPsquare_coeffs (Nregcoeffs,0.0);
outfieldx->FwdTrans(dP_square, dPsquare_coeffs);
outfieldx->BwdTrans(dPsquare_coeffs, dP_square);
Array<OneD, NekDouble> mu53 (nq1D,0.0);
Array<OneD, NekDouble> d2v (nq1D,0.0);
Array<OneD, NekDouble> prod (nq1D,0.0);
double pow = 1.0/3.0;
double base;
for(int y=0; y<nq1D; y++)
{
base =dUreg[y];
mu53[y] = std::pow ((base*base),pow);
mu53[y] = 1/(base*mu53[y]);
//prod[y] = mu53[y]*dP_square[y];
}
//attempting to smooth mu53
Array<OneD, NekDouble> mucoeffs(Nregcoeffs,0.0);
outfieldx->FwdTrans(mu53, mucoeffs);
outfieldx->BwdTrans(mucoeffs, mu53);
Vmath::Vmul(nq1D, mu53,1, dP_square,1, prod,1);
//attempting to smooth field...
Array<OneD, NekDouble> prod_coeffs (Nregcoeffs,0.0);
outfieldx->FwdTrans(prod, prod_coeffs);
outfieldx->BwdTrans(prod_coeffs, prod);
Vmath::Zero(Nregcoeffs,prod_coeffs,1);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, prod ,d2v);
//attempting to smooth the vjump
Vmath::Zero( Nregcoeffs, prod_coeffs,1);
outfieldx->FwdTrans(d2v, prod_coeffs);
outfieldx->BwdTrans(prod_coeffs, d2v);
//test prod deriv... add background..
//Array<OneD, NekDouble> prod_b(nq1D);
//Array<OneD, NekDouble> dersprod_b(nq1D);
//Vmath::Sadd(nq1D,0.001, prod,1,prod_b,1);
//Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, prod_b ,dersprod_b);
Array<OneD, NekDouble> pjump(nq1D);
Array<OneD, NekDouble> vjump(nq1D);
//test a sin curve....
/*
Array<OneD, NekDouble> fun1D(nq1D);
Array<OneD, NekDouble> derfun1D(nq1D);
NekDouble pi =3.14159265;
for(int e=0; e<nq1D; e++)
{
fun1D[e]= std::sin((x0d[e]));
}
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, fun1D, derfun1D);
//add noise to a function
Array<OneD, NekDouble> fun1D_b(nq1D);
Array<OneD, NekDouble> dersfun1D_b(nq1D);
int cnt=0;
int cnt1=1;
for(int v=0; v<nq1D/2; v++)
{
fun1D_b[cnt] = fun1D[cnt] +0.000000;
//cout<<cnt<<"x="<<x0d[cnt]<<" fun1D_b="<<fun1D_b[cnt]<<" fun1D="<<fun1D[cnt]<<" diff="<<
//fun1D[cnt]-fun1D_b[cnt]<<endl;
//ASSERTL0((fun1D_b[cnt]-fun1D[cnt])==0.0001, "problem");
cnt = cnt +2;
fun1D_b[cnt1] = fun1D[cnt1] -0.000000;
//cout<<cnt1<<"x="<<x0d[cnt1]<<" fun1D_b="<<fun1D_b[cnt1]<<" fun1D="<<fun1D[cnt1]<<" diff="<<
//fun1D[cnt1]-fun1D_b[cnt1]<<endl;
//ASSERTL0((fun1D_b[cnt1]-fun1D[cnt1])==-0.0001, "problem");
cnt1 = cnt1 +2;
if( (v==((nq1D/2) -1)) && (nq1D%2==1)
)
{
fun1D_b[cnt] = fun1D[cnt] +0.000000;
//cout<<cnt<<"x="<<x0d[cnt]<<" fun1D_b="<<fun1D_b[cnt]<<" fun1D="<<fun1D[cnt]<<" diff="<<
//fun1D[cnt]-fun1D_b[cnt]<<endl;
}
}
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, fun1D_b ,dersfun1D_b);
*/
/*
int nqedge =nq1D/Icompreg->size();
Array<OneD, NekDouble> phys_edge (nqedge);
Array<OneD, NekDouble> der_edge (nqedge);
Array<OneD, NekDouble> der_fun1Db_std (nq1D);
//extract stdderiv....
for(int w=0; w<Icompreg->size() ; w++)
{
Vmath::Vcopy(nqedge,&(fun1D_b[w*nqedge]),1,&(phys_edge[0]),1);
//check if the metrix is the same for streak and Ilayer obj
StdRegions::StdExpansion1DSharedPtr edgestdexp;
edgestdexp = boost::dynamic_pointer_cast<StdRegions::StdExpansion1D> (
Ilayer->GetPlane(0)->GetExp(w) );
edgestdexp->StdPhysDeriv(phys_edge,
der_edge);
Vmath::Vcopy(nqedge, &(der_edge[0]),1, &(der_fun1Db_std[w*nqedge]),1);
Vmath::Zero(nqedge, phys_edge,1);
Vmath::Zero(nqedge, der_edge,1);
Vmath::Vcopy(nqedge, &(Rephysreg[w*nqedge]),1, &(phys_edge[0]),1);
edgestdexp->StdPhysDeriv(phys_edge, der_edge);
Vmath::Vcopy(nqedge, &(der_edge[0]),1, &(der_Pre_std[w*nqedge]),1);
Vmath::Zero(nqedge, phys_edge,1);
Vmath::Zero(nqedge, der_edge,1);
Vmath::Vcopy(nqedge, &(fun1D[w*nqedge]),1, &(phys_edge[0]),1);
edgestdexp->StdPhysDeriv(phys_edge,
der_edge);
Vmath::Vcopy(nqedge, &(der_edge[0]),1, &(derfun1D[w*nqedge]),1);
Vmath::Zero(nqedge, phys_edge,1);
Vmath::Zero(nqedge, der_edge,1);
}
*/
//final jump conditions
//n0=2*pi*(2/3)^(2/3)*(-2/3)! = 2*pi*(2/3)^(2/3)*gamma(1/3)
NekDouble n0 = 12.845424015;
//use the session to get the values of rho,alpha...
NekDouble rho ;
//rho = 0.08;
rho = session->GetParameter("RHO");
//load alpha only if not manually specified
if(alpha==0)
{
alpha = session->GetParameter("ALPHA");
}
cout<<"alpha="<<alpha<<endl;
ASSERTL0(alpha!=0, "alpha cannot be 0");
NekDouble alpha53;
//alpha53=1;
alpha53 = std::pow ((alpha*alpha),pow);
alpha53 = 1/(alpha*alpha53);
for(int c=0; c<nq1D; c++)
{
pjump[c] = -n0*alpha53*curv[c]*prod[c] ;
}
Vmath::Smul(nq1D, n0,d2v,1, vjump,1);
Vmath::Smul(nq1D, alpha53,vjump,1, vjump,1);
cout<<"alpha^-5/3="<<alpha53<<endl;
Array<OneD, NekDouble> txcoeffs (Nregcoeffs);
Array<OneD, NekDouble> tycoeffs (Nregcoeffs);
//attempt to smooth the tangent components:
outfieldx->FwdTrans(tx, txcoeffs);
outfieldx->FwdTrans(ty, tycoeffs);
outfieldy->BwdTrans(txcoeffs, tx);
outfieldy->BwdTrans(tycoeffs, ty);
cout<<"RHO=="<<rho<<endl;
//end
//PAY ATTENTION to the sign (vjump -/+ pjump)*t!!!!!!!!!!
for(int j=0; j<nq1D; j++)
{
//start
(outfieldx->UpdatePhys())[j] =
rho*rho*(vjump[j]*tx[j]-pjump[j]);
(outfieldy->UpdatePhys())[j] =
rho*rho*(vjump[j]*ty[j]-pjump[j]);
//cout<<x0d[j]<<" "<<curv[j]*prod[j]<<" "<<(outfieldx->GetPhys())[j]<<" "<<(outfieldy->GetPhys())[j]<<endl;
//end
//decomment
/*
(outfieldx->UpdatePhys())[j] =
20*sin(2*x0d[j]);
(outfieldy->UpdatePhys())[j] =
20*2*cos(2*x0d[j])/3.14;
*/
}
//FINAL REFINEMENT:::
//start
Array<OneD, NekDouble> finalcoeffs (Nregcoeffs);
outfieldx->FwdTrans(outfieldx->GetPhys(), finalcoeffs);
outfieldx->BwdTrans(finalcoeffs, outfieldx->UpdatePhys());
Vmath::Zero(Nregcoeffs,finalcoeffs,1);
outfieldy->FwdTrans(outfieldy->GetPhys(), finalcoeffs);
outfieldy->BwdTrans(finalcoeffs, outfieldy->UpdatePhys());
//symmetrize the jumps conditions
if(symm ==true)
{
Array<OneD, NekDouble> tmpx(nq1D);
Array<OneD, NekDouble> tmpy(nq1D);
cout<<"symmetrise the jump conditions"<<endl;
for(int i = 0; i < nq1D; ++i)
{
tmpx[i] =0.5*(outfieldx->GetPhys()[i] - outfieldx->GetPhys()[Refindices[i]]);
tmpy[i] =0.5*(outfieldy->GetPhys()[i] - outfieldy->GetPhys()[Refindices[i]]);
}
Vmath::Vcopy(nq1D, tmpx,1, outfieldx->UpdatePhys(),1);
Vmath::Vcopy(nq1D, tmpy,1, outfieldy->UpdatePhys(),1);
}
/*
//calculate J,K
Array<OneD, NekDouble> lambda(nq1D);
Array<OneD, NekDouble> a(nq1D);
Array<OneD, NekDouble> dPre_dz(nq1D);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, Rephysreg, dPre_dz);
Array<OneD, NekDouble> dPim_dz(nq1D);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, Imphysreg, dPim_dz);
Array<OneD, NekDouble> dPre_dz2(nq1D);
Vmath::Vmul(nq1D, dPre_dz,1, dPre_dz,1, dPre_dz2,1);
Array<OneD, NekDouble> dPim_dz2(nq1D);
Vmath::Vmul(nq1D, dPim_dz,1, dPim_dz,1, dPim_dz2,1);
Array<OneD, NekDouble> dP_dz2(nq1D);
Vmath::Vadd(nq1D, dPre_dz2,1,dPim_dz2,1,dP_dz2,1);
Array<OneD, NekDouble> ddP_dz2z(nq1D);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, dP_dz2, ddP_dz2z);
Array<OneD, NekDouble> delta5(nq1D);
Array<OneD, NekDouble> a53(nq1D);
pow = 1./3.;
for(int e=0; e<nq1D; e++)
{
//delta[e] = 1+ f_z[e]*f_z[e];
delta5[e] = delta[e]*delta[e]*delta[e]*delta[e]*delta[e];
lambda[e] = dUreg[e]/sqrt(delta[e]);
a[e] = lambda[e]*alpha/delta[e];
a53[e] = std::pow( (a[e]*a[e]), pow);
//cout<<"a^2/3="<<a53[e]<<endl;
a53[e] = a[e]*a53[e];
}
Array<OneD, NekDouble> delta_z(nq1D);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, delta, delta_z);
Array<OneD, NekDouble> a_z(nq1D);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, a, a_z);
Array<OneD, NekDouble> J(nq1D);
Array<OneD, NekDouble> K(nq1D);
Array<OneD, NekDouble> f1(nq1D);
Array<OneD, NekDouble> f2(nq1D);
//NekDouble a53,delta5;
for(int f=0; f<nq1D; f++)
{
//a53 = std::pow(a[f],5./3.);
//delta5 = std::pow( delta[f], 5);
J[f] =( (n0*rho*rho)/(a53[f]*delta5[f]) )*
(
( (-7*delta_z[f])/(2*delta[f]) - (5*a_z[f] )/(3*a[f]) )*dP_dz2[f] +
ddP_dz2z[f]
);
K[f] = -(n0*rho*rho)/(a53[f]*delta5[f])*f_zz[f]*dP_dz2[f];
f1[f] = lambda[f]*(K[f]-delta[f]*f_z[f]*J[f]);
f2[f] = -lambda[f]*(f_z[f]*K[f] +delta[f]*J[f]);
}
//test dermu
Array<OneD, NekDouble> dermu53(nq1D);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, mu53, dermu53);
//test d(dP_square)ds
Array<OneD, NekDouble> ddP_square_ds(nq1D);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, dP_square, ddP_square_ds);
*/
cout<<"length layer="<<length<<endl;
//calc dUds
Array<OneD, NekDouble> dUds(nq1D);
Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, stphysreg, dUds);
Array<OneD, int> indexjacwarn(nq1D,-1);
NekDouble invjac2D;
NekDouble jactest;
NekDouble laytest=0;
NekDouble basetest;
for(int g=0; g<nq1D; g++)
{
//NBBB dPre is the stdDERIV!!!
invjac2D = (gmat0[g]*tx[g] +gmat2[g]*ty[g]);
jactest = (invjac2D + (1/Jac[g]))/2.;
jactest = (invjac2D- (1/Jac[g]))/jactest;
basetest = abs(curv[g]);
/*
cout<<x0d[g]<<" "<<
//stphysreg[g]<<" "<<lambda[g]<<" "<<
//delta[g]<<" "<<
//f_z[g]<<" "<<a53[g]<<" "<<
//delta5[g]<<" "<<delta_z[g]<<" "<<a_z[g]<<" "<<
//" "<<ddP_dz2z[g]<<" "<<dP_dz2[g]<<" "<<
//dP_square[g]<<" "<<dP_dz2[g]/delta[g]<<" "<<
//ddP_square_ds[g]<<" "<<ddP_dz2z[g]/(delta[g]*sqrt(delta[g]))<<" "<<
//ddP_dz2z[g]/delta[g]<<" "<<
//a[g]<<" "<<a53[g]<<" "<<
vjump[g]*tx[g]*sqrt(delta[g])<<" "<<
J[g]<<" "<<K[g]<<" "<<pjump[g]<<" "<<
f1[g]<<" "<<f2[g]<<" "<<
f_zz[g]/(sqrt(delta5[g]))<<" "<<
//f_zz[g]/sqrt(delta[g]*delta[g]*delta[g])<<" "
curv[g]<<" "<<
outfieldx->GetPhys()[g]<<" "
<<outfieldy->GetPhys()[g]
<<endl;
*/
laytest = (abs(stphysreg[g])-cr);
cout<<setw(14)<<x0d[g]<<" "<<setw(13)<<x1d[g]<<
" "<<Rephysreg[g]<<" "<<dP_re[g]<<" "
<<stphysreg[g]<<" "<<dUreg[g]<<" "<<mu53[g]<<" "<<curv[g]
//<<" "<<-1.5*std::pow(abs(curv[g]),-1./5.)<<" "
//<<" "<<std::pow(delta[g],-0.5)-0.45*ty[g]+(f_z[g]*f_z[g])/(std::pow(delta[g],0.5))<<" "
<<" "<<f_z[g]<<" "
<<f_zz[g]<<" "
//<<-1.5*std::pow(abs(curv[g]/delta[g]),-1./5.)<<" "
<<s0d[g]
<<setw(13)<<" "<<dP_square[g]<<" "<<setw(13)<<prod[g]
<<" "<<
//(gmat0[g]*tx[g] +gmat2[g]*ty[g])<<" "<<1/Jac[g]
//<<" "<<dermu53[g]
//<<" "<<derfun1D[g]
//nx[g]<<" "<<ny[g]<<" "<<
tx[g]<<" "<<ty[g]<<" "<<
//f1[g]<<" "<<f2[g]
//prod_b[g]<<" "<<dersprod_b[g]
//fun1D[g]<<" "<<derfun1D[g]<<" "<<
//der_fun1Db_std[g]<<" "<<
//fun1D_b[g]<<" "<<dersfun1D_b[g]
pjump[g]<<setw(13)<<" "<<d2v[g]<<" "
<<outfieldx->GetPhys()[g]<<" "
<<outfieldy->GetPhys()[g]<<" "<<Imphysreg[g]<<" "<<dP_im[g]
<<" "<<"AA"<<endl;
//cout<<"laytest="<<laytest<<endl;
//cout<<"(gmat0[g]*tx[g] +gmat2[g]*ty[g])="<<(gmat0[g]*tx[g] +gmat2[g]*ty[g])<<
//" 1/jac1D="<<1./Jac[g]<<endl;
ASSERTL0(invjac2D*Jac[g]>0, " sign jac problem..");
//30% error is allowed!!!
ASSERTL0(abs(jactest)<0.3, "jac 1D problem..");
//save info of point where the error >20%
if(jactest>=0.2)
{
indexjacwarn[g] =g;
}
}
if(cr==0)
{
ASSERTL0(abs(laytest)<0.002, "critical layer wrong");
}
else
{
cout<<"WARNING: crit may be wrong"<<endl;
}
// gamma(1/3)= 2.6789385347077476337
// need the tangents related to the expList1D outfieldx[region]
for(int m=0; m<nq1D; m++)
{
if(indexjacwarn[m]!=-1)
{
cout<<"warning: point with jacerr>20% index="<<m<<" x="<<x0d[m]<<" y="<<x1d[m]<<endl;
}
}
for(int a=0; a<Elmtid.num_elements(); a++)
{
cout<<"elmt id="<<Elmtid[a]<<" edge id="<<Edgeid[a]<<endl;
}
//end
}
Array<OneD, int> GetReflectionIndex ( MultiRegions::ExpListSharedPtr  Exp,
int  Ireg 
)

Definition at line 637 of file FldCalcBCs.cpp.

References ASSERTL0, Nektar::NekConstants::kGeomFactorsTol, npts, and Vmath::Vmax().

Referenced by main().

{
int i,j;
//remember Exp is the streak so GetPlane(0) is not needed
Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =Exp->GetBndCondExpansions();
MultiRegions::ExpListSharedPtr Ilayer = Iexp[Ireg];
int npts = Ilayer->GetNpoints();
Array<OneD, int> index(npts);
Array<OneD, NekDouble> coord(2);
Array<OneD, NekDouble> coord_x(npts);
Array<OneD, NekDouble> coord_y(npts);
//-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
Ilayer->GetCoords(coord_x,coord_y);
NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
NekDouble tol = NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
NekDouble xnew,ynew;
int start = npts-1;
//cout<<"xmax="<<xmax<<endl;
for(i = 0; i < npts; ++i)
{
//cout<<"Ref index for x="<<coord_x[i]<<endl;
xnew = - coord_x[i] + xmax;
ynew = - coord_y[i];
for(j = start; j >=0 ; --j)
{
if((coord_x[j]-xnew)*(coord_x[j]-xnew) < tol)
{
index[i] = j;
start = j;
break;
}
}
if(j == -1)
{
for(j = npts-1; j > start; --j)
{
if((coord_x[j]-xnew)*(coord_x[j]-xnew) < tol)
{
index[i] = j;
break;
}
}
ASSERTL0(j != start,"Failsed to find matching point");
}
}
return index;
}
Array<OneD, int> GetReflectionIndex2D ( MultiRegions::ExpListSharedPtr  wavefield)

Definition at line 1911 of file FldCalcBCs.cpp.

References ASSERTL0, Nektar::NekConstants::kGeomFactorsTol, npts, and Vmath::Vmax().

Referenced by main().

{
int i,j;
int npts = wavefield->GetPlane(0)->GetNpoints();
Array<OneD, int> index(npts);
Array<OneD, NekDouble> coord(2);
Array<OneD, NekDouble> coord_x(npts);
Array<OneD, NekDouble> coord_y(npts);
//-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
wavefield->GetPlane(0)->GetCoords(coord_x,coord_y);
NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
NekDouble tol = NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
NekDouble xnew,ynew;
int start = npts-1;
for(i = 0; i < npts; ++i)
{
xnew = - coord_x[i] + xmax;
ynew = - coord_y[i];
for(j = start; j >=0 ; --j)
{
if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
{
index[i] = j;
start = j;
break;
}
}
if(j == -1)
{
for(j = npts-1; j > start; --j)
{
if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
{
index[i] = j;
break;
}
}
ASSERTL0(j != start,"Failsed to find matching point");
}
}
return index;
}
int main ( int  argc,
char *  argv[] 
)

Below is the old function from GeomFactors1D::v_ComputeEdgeTangents. This However, the vectors were computed in the function below, but is never called and therefore was never working.

This utility therefore does not presently work, but the (now removed) function is kept here for reference.

– Chris Cantwell 

 void GeomFactors1D::v_ComputeEdgeTangents(
             const GeometrySharedPtr &geom,
             const int edge,
             const LibUtilities::PointsKey &to_key)
 {
     int k;
     int nquad= to_key.GetNumPoints();
     Geometry2DSharedPtr g;
     ASSERTL0(g= boost::dynamic_pointer_cast<Geometry2D>(geom),
              "FAIL");

     GeomFactorsSharedPtr gf = geom->GetMetricInfo();

cannot use m_type here GeomType gtype = gf->GetGtype(); GeomType gtype= m_type; m_tangent =Array<OneD, Array<OneD, NekDouble> >(m_coordDim); for( k=0; k< m_coordDim; ++k) { m_tangent[k] = Array<OneD, NekDouble>(nquad); }

int i;

DerivStorage deriv = GetDeriv(m_pointsKey); FillDeriv(deriv, m_pointsKey); 

        NekDouble fac;

Regular geometry case if((gtype == eRegular)||(gtype == eMovingRegular)) { 

for(i = 0; i < m_coordDim; ++i)
{
        Vmath::Fill(nquad, deriv[0][i][0],m_tangent[i],1);
}

normalise fac = 0.0; for(i =0 ; i < m_coordDim; ++i) { fac += m_tangent[i][0]*m_tangent[i][0]; } fac = 1.0/sqrt(fac); for (i = 0; i < m_coordDim; ++i) { Vmath::Smul(nquad,fac,m_tangent[i],1,m_tangent[i],1); } }

else // Set up deformed tangents {

Array<OneD, NekDouble> jac(m_coordDim*nquad); 

for(i = 0; i < m_coordDim; ++i)
{
    for(int j=0; j<nquad; j++)
        {
    m_tangent[i][j] = deriv[0][i][j];
    }
}

normalise normal vectors Array<OneD,NekDouble> work(nquad,0.0); for(i = 0; i < m_coordDim; ++i) { Vmath::Vvtvp(nquad, m_tangent[i],1, m_tangent[i],1,work,1,work,1); }

Vmath::Vsqrt(nquad,work,1,work,1); Vmath::Sdiv(nquad,1.0,work,1,work,1);

for(i = 0; i < m_coordDim; ++i) { Vmath::Vmul(nquad, m_tangent[i],1,work,1,m_tangent[i],1); } } }

To create another bcs file comment everything from '//start' to '//end' decomment the lines which follow '//decomment' and run with the command: ./FldCalcBCs meshfile fieldfile fieldfile

Definition at line 109 of file FldCalcBCs.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, CalcNonLinearForcing(), Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::SpatialDomains::eCalcBC, Nektar::LibUtilities::ePolyEvenlySpaced, Extractlayerdata(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), GetReflectionIndex(), GetReflectionIndex2D(), Nektar::LibUtilities::Import(), Nektar::NekConstants::kNekUnsetDouble, Manipulate(), Nektar::SpatialDomains::MeshGraph::Read(), SetFields(), WriteBcs(), and WriteFld().

{
void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh,
SpatialDomains::BoundaryConditionsSharedPtr& boundaryConditions,
LibUtilities::SessionReaderSharedPtr &session,
Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,int nvariables);
Array<OneD, int> GetReflectionIndex(MultiRegions::ExpListSharedPtr Exp,
int Ireg);
Array<OneD, int> GetReflectionIndex2D(MultiRegions::ExpListSharedPtr wavefield);
void Extractlayerdata(Array<OneD, int> Iregions, int coordim,
SpatialDomains::MeshGraphSharedPtr &mesh,
LibUtilities::SessionReaderSharedPtr &session,
SpatialDomains::BoundaryConditions &bcs,
Array<OneD,MultiRegions::ExpListSharedPtr> &infields,
MultiRegions::ContField1DSharedPtr &outfieldx,
MultiRegions::ContField1DSharedPtr &outfieldy,
MultiRegions::ExpListSharedPtr &streak, bool symm,
Array<OneD, int> Refindices, NekDouble alpha, NekDouble cr);
void Manipulate(Array<OneD, int> Iregions, int coordim, SpatialDomains::MeshGraphSharedPtr &mesh,
SpatialDomains::BoundaryConditions &bcs,
Array<OneD,MultiRegions::ExpList1DSharedPtr> &infields,
Array<OneD,MultiRegions::ExpList1DSharedPtr> &outfieldx,
Array<OneD,MultiRegions::ExpList1DSharedPtr> &outfieldy,
MultiRegions::ExpListSharedPtr &streak);
void CalcNonLinearForcing(SpatialDomains::MeshGraphSharedPtr &mesh,
LibUtilities::SessionReaderSharedPtr &session, string fieldfile,
Array<OneD,MultiRegions::ExpListSharedPtr> &waveFields,
MultiRegions::ExpListSharedPtr &streak,
Array<OneD, int> Refindices, bool symm );
void WriteBcs(string variable, int region, string fieldfile, SpatialDomains::MeshGraphSharedPtr &mesh,
MultiRegions::ContField1DSharedPtr &outregionfield);
void WriteFld(string outfile, SpatialDomains::MeshGraphSharedPtr &mesh, Array<OneD,
MultiRegions::ExpListSharedPtr> &fields );
//REMARK:
/**
To create another bcs file comment everything from '//start' to '//end'
decomment the lines which follow '//decomment'
and run with the command: ./FldCalcBCs meshfile fieldfile fieldfile
**/
if(argc >= 6 || argc <4)
{
fprintf(stderr,"Usage: ./FldCalcBCs meshfile fieldfile streakfile (optional)alpha\n");
exit(1);
}
cout<<"argc="<<argc<<endl;
//change argc from 4 to 5 allow the loading of alpha to be optional
if(argc==4){ argc=5;}
//----------------------------------------------
string meshfile(argv[argc-4]);
//setting argc=2 will take consider only the first two words argv[0],argv[1]
LibUtilities::SessionReaderSharedPtr vSession
= LibUtilities::SessionReader::CreateInstance(2, argv);
// Read in mesh from input file
SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
//----------------------------------------------
// Also read and store the boundary conditions
SpatialDomains::BoundaryConditionsSharedPtr boundaryConditions;
boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
::AllocateSharedPtr(vSession, graphShPt);
SpatialDomains::BoundaryConditions bcs(vSession, graphShPt);
//----------------------------------------------
//load alpha is provided
NekDouble alpha=0;
string alp_s;
if(argv[argc-1])
{
alp_s = argv[argc-1];
alpha = boost::lexical_cast<double>(alp_s);
}
//----------------------------------------------
//load the phase:
NekDouble cr=0.0;
if( vSession->DefinesSolverInfo("INTERFACE")
&&vSession->GetSolverInfo("INTERFACE")=="phase" )
{
vSession->LoadParameter("phase",cr,NekConstants::kNekUnsetDouble);
}
cout<<"cr="<<cr<<endl;
//-----------------------------------------------
//determine if the symmetrization is on:
bool symm =true;
vSession->MatchSolverInfo("Symmetrization","True",symm,true);
if( symm == true )
{
cout<<"symmetrization is active"<<endl;
}
// Import field file.
string fieldfile(argv[argc-3]);
vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
vector<vector<NekDouble> > fielddata;
LibUtilities::Import(fieldfile,fielddef,fielddata);
//----------------------------------------------
// Define Expansion
std::string solvtype = vSession->GetSolverInfo("SOLVERTYPE");
int nfields;
nfields = fielddef[0]->m_fields.size();
Array<OneD, MultiRegions::ExpListSharedPtr> fields;
if(solvtype == "CoupledLinearisedNS" && nfields==1)
{
nfields++;
}
fields= Array<OneD, MultiRegions::ExpListSharedPtr>(nfields);
int lastfield = 0;
if(solvtype == "CoupledLinearisedNS" && nfields!=2)
{
SetFields(graphShPt,boundaryConditions,vSession,fields,nfields-1);
//decomment
//nfields = nfields-1;
//start
lastfield = nfields-1;
cout<<"Set pressure: "<<lastfield<<endl;
int nplanes = fielddef[0]->m_numModes[2];
const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::ePolyEvenlySpaced);
const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,false,false,graphShPt,fielddef[0]->m_fields[0]);
fields[lastfield] = Exp3DH1;
//end
}
else if(solvtype == "CoupledLinearisedNS" && nfields==2)
{
cout<<"Set pressure split"<<endl;
int nplanes = fielddef[0]->m_numModes[2];
const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::ePolyEvenlySpaced);
const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
cout<<"set ppp"<<endl;
/*
//to use GetBndCondExpansions() a contfield is needed and you j=have to
//call it "u" to get the bndconds
MultiRegions::ContField3DHomogeneous1DSharedPtr Cont3DH1;
Cont3DH1 = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,false,false,graphShPt, "u");
fields[0] = Cont3DH1;
*/
fields[0] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
::AllocateSharedPtr (vSession,Bkey,lz,false,false,graphShPt,vSession->GetVariable(0));
//pressure is field 1
lastfield = nfields-1;
MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,false,false,graphShPt,fielddef[0]->m_fields[0]);
fields[lastfield] = Exp3DH1;
}
else
{
SetFields(graphShPt,boundaryConditions,vSession,fields,nfields);
}
//----------------------------------------------
// Copy data from file:fill fields with the fielddata
if(lastfield==1)
{
for(int i = 0; i < fielddata.size(); ++i)
{
fields[0]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[0]->m_fields[0], fields[0]->UpdateCoeffs());
}
fields[0]->BwdTrans_IterPerExp(fields[0]->GetCoeffs(),fields[0]->UpdatePhys());
cout<<"field:"<<fielddef[0]->m_fields[0]<<endl;
/*
//bwd plane 0
fields[0]->GetPlane(0)->BwdTrans_IterPerExp(fields[0]->GetPlane(0)->GetCoeffs(),
fields[0]->GetPlane(0)->UpdatePhys() );
cout<<"hjhj"<<endl;
//bwd plane 1
fields[0]->GetPlane(1)->BwdTrans_IterPerExp(fields[0]->GetPlane(1)->GetCoeffs(),
fields[0]->GetPlane(1)->UpdatePhys() );
*/
for(int i = 0; i < fielddata.size(); ++i)
{
fields[lastfield]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[i]->m_fields[0], fields[lastfield]->UpdateCoeffs());
}
fields[lastfield]->BwdTrans_IterPerExp(fields[lastfield]->GetCoeffs(),fields[lastfield]->UpdatePhys());
/*
//bwd plane 0
fields[lastfield]->GetPlane(0)->BwdTrans_IterPerExp(fields[lastfield]->GetPlane(0)->GetCoeffs(),
fields[lastfield]->GetPlane(0)->UpdatePhys() );
//bwd plane 1
fields[lastfield]->GetPlane(1)->BwdTrans_IterPerExp(fields[lastfield]->GetPlane(1)->GetCoeffs(),
fields[lastfield]->GetPlane(1)->UpdatePhys() );
*/
}
else
{
for(int j = 0; j < nfields; ++j)
{
for(int i = 0; i < fielddata.size(); ++i)
{
fields[j]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[i]->m_fields[j], fields[j]->UpdateCoeffs());
}
fields[j]->BwdTrans_IterPerExp(fields[j]->GetCoeffs(),fields[j]->UpdatePhys());
}
}
//----------------------------------------------
//the phys values are the same (WRONG) but the coeffs are CORRECT!!!
/*
for(int g=0; g<fields[0]->GetPlane(1)->GetNcoeffs(); g++)
{
cout<<"g="<<g<<" coeff f0="<<fields[lastfield]->GetPlane(0)->GetCoeff(g)<<" f1="<<fields[lastfield]->GetPlane(1)->GetCoeff(g)<<endl;
}
*/
// import the streak
//import the streak field
//add y to obtain the streak w=w'+y and write fld file
MultiRegions::ExpListSharedPtr streak;
//start
streak = MemoryManager<MultiRegions::ContField2D>
::AllocateSharedPtr(vSession, graphShPt, "w",true);
//streak = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr
// (vSession, graphShPt, true, "w");
string streakfile(argv[argc-2]);
vector<LibUtilities::FieldDefinitionsSharedPtr> streakdef;
vector<vector<NekDouble> > streakdata;
LibUtilities::Import(streakfile,streakdef,streakdata);
//attention: streakdata.size==2 because the session file is 3DHomo1D but in
//reality the streak is real quantity
// Copy data from file:fill streak with the streakdata
for(int i = 0; i < streakdata.size(); ++i)
{
streak->ExtractDataToCoeffs(streakdef[i],streakdata[i],streakdef[i]->m_fields[0], streak->UpdateCoeffs());
}
streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
//end
//---------------------------------------------------------------
// determine the I regions
//hypothesis: the number of I regions is the same for all the variables
//hypothesis: all the I regions have the same nq points
int nIregions, lastIregion;
const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConditions = streak->GetBndConditions();
Array<OneD, int> Iregions =Array<OneD, int>(bndConditions.num_elements(),-1);
//3 internal layers:
Array<OneD, int> Ilayers =Array<OneD, int>(3,-1);
nIregions=0;
int nbnd= bndConditions.num_elements();
for(int r=0; r<nbnd; r++)
{
if(bndConditions[r]->GetUserDefined()==SpatialDomains::eCalcBC)
{
lastIregion=r;
Iregions[r]=r;
Ilayers[nIregions]=r;
nIregions++;
//cout<<"Iregion="<<r<<endl;
}
}
ASSERTL0(nIregions>0,"there is any boundary region with the tag USERDEFINEDTYPE=""CalcBC"" specified");
//-----------------------------------------------------------------
//set 1D output fields:
//initialise fields
const SpatialDomains::BoundaryRegionCollection &bregions = bcs.GetBoundaryRegions();
MultiRegions::ExpList1DSharedPtr outfieldxtmp;
MultiRegions::ExpList1DSharedPtr outfieldytmp;
MultiRegions::ContField1DSharedPtr outfieldx;
MultiRegions::ContField1DSharedPtr outfieldy;
//initialie fields
//const SpatialDomains::BoundaryRegionCollection &bregions = bcs.GetBoundaryRegions();
//declarecoeffsphysarray?????????!!!! true?!
outfieldxtmp = MemoryManager<MultiRegions::ExpList1D>
::AllocateSharedPtr(*(bregions.find(lastIregion)->second), graphShPt, true);
outfieldytmp = MemoryManager<MultiRegions::ExpList1D>
::AllocateSharedPtr(*(bregions.find(lastIregion)->second), graphShPt, true);
outfieldx = MemoryManager<MultiRegions::ContField1D>
::AllocateSharedPtr(vSession, *outfieldxtmp);
outfieldy = MemoryManager<MultiRegions::ContField1D>
::AllocateSharedPtr(vSession, *outfieldytmp);
//YOU NEED TO DEFINE THE OUTFIELDS AS ContField2D to use them in the 2D NS!!!!!!!!!!!
//(not 3DHomogeneous1D!!!!!!!
/*
//define the outfields:
int nplanes = fielddef[0]->m_numModes[1];
int nzlines = fielddef[0]->m_numModes[2];
//ATTENTO Polyevenlyspaced NON Fourier!!!
const LibUtilities::PointsKey PkeyZ(nzlines,LibUtilities::ePolyEvenlySpaced);
const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
for(int r=0; r<nlayers; r++)
{
bndfieldx[r]= MemoryManager<MultiRegions::ExpList2DHomogeneous1D>
::AllocateSharedPtr(vSession,BkeyZ,lz,false,graphShPt);
cout<<"OOOK"<<endl;
bndfieldy[r]= MemoryManager<MultiRegions::ExpList2DHomogeneous1D>
::AllocateSharedPtr(vSession,BkeyZ,lz,false,graphShPt);
}
//--------------------------------------------------------
*/
bool coarseVWI=false;
if(coarseVWI==false)
{
//manipulate data
//for 2 variables(u,v) only:
int coordim = graphShPt->GetMeshDimension();
//remark Ilayers[2] is the critical layer
static Array<OneD, int> Refindices = GetReflectionIndex(streak, Ilayers[2]);
Extractlayerdata(Ilayers,coordim, graphShPt,vSession, bcs,
fields, outfieldx,outfieldy,streak,symm,Refindices, alpha,cr);
//--------------------------------------------------------------------------------------
//write bcs files: one for each I region and each variable
string var="u";
WriteBcs(var,lastIregion, fieldfile,graphShPt,outfieldx);
var="v";
WriteBcs(var,lastIregion, fieldfile,graphShPt,outfieldy);
//--------------------------------------------------------------------------------
}
else
{
cout<<"CalcNonLinearForcing"<<endl;
static Array<OneD, int> Refindices;
if(symm==true)
{
Refindices = GetReflectionIndex2D(fields[0]);
}
CalcNonLinearForcing(graphShPt,vSession, fieldfile, fields, streak, Refindices,symm ) ;
}
}
void Manipulate ( Array< OneD, int >  Iregions,
int  coordim,
SpatialDomains::MeshGraphSharedPtr &  mesh,
SpatialDomains::BoundaryConditions &  bcs,
Array< OneD, MultiRegions::ExpListSharedPtr > &  infields,
Array< OneD, MultiRegions::ExpList1DSharedPtr > &  outfieldx,
Array< OneD, MultiRegions::ExpList1DSharedPtr > &  outfieldy,
MultiRegions::ExpListSharedPtr &  streak 
)

Definition at line 1671 of file FldCalcBCs.cpp.

Referenced by main().

{
}
void SetFields ( SpatialDomains::MeshGraphSharedPtr &  mesh,
SpatialDomains::BoundaryConditionsSharedPtr &  boundaryConditions,
LibUtilities::SessionReaderSharedPtr &  session,
Array< OneD, MultiRegions::ExpListSharedPtr > &  Exp,
int  nvariables 
)

< physical length in X direction (if homogeneous)

< physical length in Y direction (if homogeneous)

< physical length in Z direction (if homogeneous)

< number of points in X direction (if homogeneous)

< number of points in Y direction (if homogeneous)

< number of points in Z direction (if homogeneous)

Parameter for homogeneous expansions

Definition at line 474 of file FldCalcBCs.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::eFourier, and Nektar::LibUtilities::eFourierEvenlySpaced.

Referenced by main().

{
// Setting parameteres for homogenous problems
NekDouble LhomX; ///< physical length in X direction (if homogeneous)
NekDouble LhomY; ///< physical length in Y direction (if homogeneous)
NekDouble LhomZ; ///< physical length in Z direction (if homogeneous)
bool DeclareCoeffPhysArrays = true;
int npointsX; ///< number of points in X direction (if homogeneous)
int npointsY; ///< number of points in Y direction (if homogeneous)
int npointsZ; ///< number of points in Z direction (if homogeneous)
int HomoDirec = 0;
bool useFFT = false;
bool deal = false;
///Parameter for homogeneous expansions
enum HomogeneousType
{
eHomogeneous1D,
eHomogeneous2D,
eHomogeneous3D,
eNotHomogeneous
};
enum HomogeneousType HomogeneousType = eNotHomogeneous;
if(session->DefinesSolverInfo("HOMOGENEOUS"))
{
std::string HomoStr = session->GetSolverInfo("HOMOGENEOUS");
//m_spacedim = 3;
if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
(HomoStr == "1D")||(HomoStr == "Homo1D"))
{
HomogeneousType = eHomogeneous1D;
npointsZ = session->GetParameter("HomModesZ");
LhomZ = session->GetParameter("LZ");
HomoDirec = 1;
}
if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
(HomoStr == "2D")||(HomoStr == "Homo2D"))
{
HomogeneousType = eHomogeneous2D;
npointsY = session->GetParameter("HomModesY");
LhomY = session->GetParameter("LY");
npointsZ = session->GetParameter("HomModesZ");
LhomZ = session->GetParameter("LZ");
HomoDirec = 2;
}
if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
(HomoStr == "3D")||(HomoStr == "Homo3D"))
{
HomogeneousType = eHomogeneous3D;
npointsX = session->GetParameter("HomModesX");
LhomX = session->GetParameter("LX");
npointsY = session->GetParameter("HomModesY");
LhomY = session->GetParameter("LY");
npointsZ = session->GetParameter("HomModesZ");
LhomZ = session->GetParameter("LZ");
HomoDirec = 3;
}
if(session->DefinesSolverInfo("USEFFT"))
{
useFFT = true;
}
}
int i;
int expdim = mesh->GetMeshDimension();
//Exp= Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
// I can always have 3 variables in a 2D mesh (oech vel component i a function which can depend on 1-3 var)
// Continuous Galerkin projection
switch(expdim)
{
case 1:
{
if(HomogeneousType == eHomogeneous2D)
{
const LibUtilities::PointsKey PkeyY(npointsY,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,npointsY,PkeyY);
const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
for(i = 0 ; i < nvariables; i++)
{
Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous2D>
::AllocateSharedPtr(session,BkeyY,BkeyZ,LhomY,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
}
}
else
{
for(i = 0 ; i < nvariables; i++)
{
Exp[i] = MemoryManager<MultiRegions::ContField1D>
::AllocateSharedPtr(session,mesh,session->GetVariable(i));
}
}
break;
}
case 2:
{
if(HomogeneousType == eHomogeneous1D)
{
const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
for(i = 0 ; i < nvariables; i++)
{
Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
::AllocateSharedPtr(session,BkeyZ,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
}
}
else
{
i = 0;
MultiRegions::ContField2DSharedPtr firstfield;
firstfield = MemoryManager<MultiRegions::ContField2D>
::AllocateSharedPtr(session,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
Exp[0] = firstfield;
for(i = 1 ; i < nvariables; i++)
{
Exp[i] = MemoryManager<MultiRegions::ContField2D>
::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
}
}
break;
}
case 3:
{
if(HomogeneousType == eHomogeneous3D)
{
ASSERTL0(false,"3D fully periodic problems not implemented yet");
}
else
{
i = 0;
MultiRegions::ContField3DSharedPtr firstfield =
MemoryManager<MultiRegions::ContField3D>
::AllocateSharedPtr(session,mesh,session->GetVariable(i));
Exp[0] = firstfield;
for(i = 1 ; i < nvariables; i++)
{
Exp[i] = MemoryManager<MultiRegions::ContField3D>
::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i));
}
}
break;
}
default:
ASSERTL0(false,"Expansion dimension not recognised");
break;
}
}
void WriteBcs ( string  variable,
int  region,
string  fieldfile,
SpatialDomains::MeshGraphSharedPtr &  mesh,
MultiRegions::ContField1DSharedPtr &  outregionfield 
)

Definition at line 1962 of file FldCalcBCs.cpp.

References Nektar::LibUtilities::Write().

Referenced by main().

{
string outfile = fieldfile.substr(0,fieldfile.find_last_of("."));
outfile +="_"+variable+"_";
char ibnd[16]="";
sprintf(ibnd,"%d",region);
outfile +=ibnd;
string endfile(".bc");
outfile += endfile;
Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);
outregionfield->FwdTrans_IterPerExp(outregionfield->GetPhys(),outregionfield->UpdateCoeffs());
fieldcoeffs[0] = outregionfield->UpdateCoeffs();
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
= outregionfield->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
// copy Data into FieldData and set variable
//for(int i = 0; i < FieldDef.size(); ++i)
//{
//FieldDef[i]->m_fields.push_back(variable);
FieldDef[0]->m_fields.push_back(variable);
//outregionfield->AppendFieldData(FieldDef[i], FieldData[i]);
outregionfield->AppendFieldData(FieldDef[0], FieldData[0]);
//outregionfield->AppendFieldData(FieldDef[i], FieldData[i],fieldcoeffs[0]);
//}
LibUtilities::Write(outfile,FieldDef,FieldData);
}
void WriteFld ( string  outfile,
SpatialDomains::MeshGraphSharedPtr &  mesh,
Array< OneD, MultiRegions::ExpListSharedPtr > &  fields 
)

Definition at line 1990 of file FldCalcBCs.cpp.

References ASSERTL0, and Nektar::LibUtilities::Write().

Referenced by main().

{
//variables array dimension is defined by fields!!!
Array<OneD, std::string> variables(fields.num_elements());
if( variables.num_elements()==2)
{
variables[0]="u";
variables[1]="v";
}
else if( variables.num_elements()==3)
{
variables[0]="u";
variables[1]="v";
variables[2]="w";
}
else
{
ASSERTL0(false, " something goes wrong... ");
}
Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(fields.num_elements());
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
= fields[0]->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
//this cycle is needed in order to use the the appendfielddata with coeffs function!!!
for(int i = 0; i < fields.num_elements(); ++i)
{
if (fields[i]->GetPhysState()==true)
{
//fields[i]->FwdTrans_IterPerExp(fields[i]->GetPhys(),fields[i]->UpdateCoeffs());
fields[i]->FwdTrans(fields[i]->GetPhys(),fields[i]->UpdateCoeffs());
}
fieldcoeffs[i] = fields[i]->UpdateCoeffs();
}
// copy Data into FieldData and set variable
for(int j = 0; j < fieldcoeffs.num_elements(); ++j)
{
//fieldcoeffs[j] = fields[j]->UpdateCoeffs();
for(int i = 0; i < FieldDef.size(); ++i)
{
// Could do a search here to find correct variable
FieldDef[i]->m_fields.push_back(variables[j]);
//remark: appendfielddata without fieldcoeffs input gives always the first field!!!
//fields[0]->AppendFieldData(FieldDef[i], FieldData[i]);
fields[0]->AppendFieldData(FieldDef[i], FieldData[i],fieldcoeffs[j]);
}
}
LibUtilities::Write(outfile,FieldDef,FieldData);
}
Generated on Sun Mar 15 2015 00:11:39 for Nektar++ by   doxygen 1.8.1.2