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().

1689  {
1690 
1691 
1692 
1693  int npts = waveFields[0]->GetPlane(0)->GetNpoints();
1694  int ncoeffs = waveFields[0]->GetPlane(0)->GetNcoeffs();
1695 
1696  Array<OneD, NekDouble> val(npts), der1(2*npts);
1697  Array<OneD, NekDouble> der2 = der1 + npts;
1698 
1699 
1700 
1701  Array<OneD, MultiRegions::ExpListSharedPtr> waveVelocities;
1702 
1703  int nvel = waveFields.num_elements()-1;
1704  waveVelocities = Array<OneD, MultiRegions::ExpListSharedPtr>(nvel);
1705  //fill velocity fields (both coeffs, phys)
1706  for( int i=0; i< nvel; i++)
1707  {
1708  waveVelocities[i] = waveFields[i];
1709  Vmath::Vcopy(npts, waveFields[i]->GetPlane(0)->GetPhys(),1, waveVelocities[i]->GetPlane(0)->UpdatePhys(),1 );
1710  Vmath::Vcopy(npts, waveFields[i]->GetPlane(1)->GetPhys(),1, waveVelocities[i]->GetPlane(1)->UpdatePhys(),1 );
1711 
1712  Vmath::Vcopy(ncoeffs, waveFields[i]->GetPlane(0)->GetCoeffs(),1, waveVelocities[i]->GetPlane(0)->UpdateCoeffs(),1 );
1713  Vmath::Vcopy(ncoeffs, waveFields[i]->GetPlane(1)->GetCoeffs(),1, waveVelocities[i]->GetPlane(1)->UpdateCoeffs(),1 );
1714 
1715  }
1716  MultiRegions::ExpListSharedPtr wavePressure;
1717  wavePressure = waveFields[nvel];
1718 
1719  //nvel=lastfield!!!
1720  Vmath::Vcopy(npts, waveFields[nvel]->GetPlane(0)->GetPhys(),1, wavePressure->GetPlane(0)->UpdatePhys(),1 );
1721  Vmath::Vcopy(npts, waveFields[nvel]->GetPlane(1)->GetPhys(),1, wavePressure->GetPlane(1)->UpdatePhys(),1 );
1722 
1723  Vmath::Vcopy(ncoeffs, waveFields[nvel]->GetPlane(0)->GetCoeffs(),1, wavePressure->GetPlane(0)->UpdateCoeffs(),1 );
1724  Vmath::Vcopy(ncoeffs, waveFields[nvel]->GetPlane(1)->GetCoeffs(),1, wavePressure->GetPlane(1)->UpdateCoeffs(),1 );
1725 
1726 
1727 
1728 
1729  static int projectfield = -1;
1730  // Set project field to be first field that has a Neumann
1731  // boundary since this not impose any condition on the vertical boundaries
1732  // Othersise set to zero.
1733  if(projectfield == -1)
1734  {
1735  Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
1736 
1737  for(int i = 0; i < waveVelocities.num_elements(); ++i)
1738  {
1739  BndConds = waveVelocities[i]->GetBndConditions();
1740  for(int j = 0; j < BndConds.num_elements(); ++j)
1741  {
1742  if(BndConds[j]->GetBoundaryConditionType() == SpatialDomains::eNeumann)
1743  {
1744  projectfield = i;
1745  break;
1746  }
1747  }
1748  if(projectfield != -1)
1749  {
1750  break;
1751  }
1752  }
1753  if(projectfield == -1)
1754  {
1755  cout << "using first field to project non-linear forcing which imposes a Dirichlet condition" << endl;
1756  projectfield = 0;
1757  }
1758  }
1759 
1760 
1761  // determine inverse of area normalised field.
1762  wavePressure->GetPlane(0)->BwdTrans(wavePressure->GetPlane(0)->GetCoeffs(),
1763  wavePressure->GetPlane(0)->UpdatePhys());
1764  wavePressure->GetPlane(1)->BwdTrans(wavePressure->GetPlane(1)->GetCoeffs(),
1765  wavePressure->GetPlane(1)->UpdatePhys());
1766 
1767  // Determine normalisation of pressure so that |P|/A = 1
1768  NekDouble norm = 0, l2;
1769  l2 = wavePressure->GetPlane(0)->L2(wavePressure->GetPlane(0)->GetPhys());
1770  norm = l2*l2;
1771  l2 = wavePressure->GetPlane(1)->L2(wavePressure->GetPlane(0)->GetPhys());
1772  norm += l2*l2;
1773  Vmath::Fill(2*npts,1.0,der1,1);
1774  NekDouble area = waveVelocities[0]->GetPlane(0)->PhysIntegral(der1);
1775  norm = sqrt(area/norm);
1776 
1777  // Get hold of arrays.
1778  waveVelocities[0]->GetPlane(0)->BwdTrans(waveVelocities[0]->GetPlane(0)->GetCoeffs(),waveVelocities[0]->GetPlane(0)->UpdatePhys());
1779  Array<OneD, NekDouble> u_real = waveVelocities[0]->GetPlane(0)->UpdatePhys();
1780  Vmath::Smul(npts,norm,u_real,1,u_real,1);
1781  waveVelocities[0]->GetPlane(1)->BwdTrans(waveVelocities[0]->GetPlane(1)->GetCoeffs(),waveVelocities[0]->GetPlane(1)->UpdatePhys());
1782  Array<OneD, NekDouble> u_imag = waveVelocities[0]->GetPlane(1)->UpdatePhys();
1783  Vmath::Smul(npts,norm,u_imag,1,u_imag,1);
1784  waveVelocities[1]->GetPlane(0)->BwdTrans(waveVelocities[1]->GetPlane(0)->GetCoeffs(),waveVelocities[1]->GetPlane(0)->UpdatePhys());
1785  Array<OneD, NekDouble> v_real = waveVelocities[1]->GetPlane(0)->UpdatePhys();
1786  Vmath::Smul(npts,norm,v_real,1,v_real,1);
1787  waveVelocities[1]->GetPlane(1)->BwdTrans(waveVelocities[1]->GetPlane(1)->GetCoeffs(),waveVelocities[1]->GetPlane(1)->UpdatePhys());
1788  Array<OneD, NekDouble> v_imag = waveVelocities[1]->GetPlane(1)->UpdatePhys();
1789  Vmath::Smul(npts,norm,v_imag,1,v_imag,1);
1790 
1791  // Calculate non-linear terms for x and y directions
1792  // d/dx(u u* + u* u)
1793  Vmath::Vmul (npts,u_real,1,u_real,1,val,1);
1794  Vmath::Vvtvp(npts,u_imag,1,u_imag,1,val,1,val,1);
1795  Vmath::Smul (npts,2.0,val,1,val,1);
1796  waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
1797 
1798 
1799  // d/dy(v u* + v* u)
1800  Vmath::Vmul (npts,u_real,1,v_real,1,val,1);
1801  Vmath::Vvtvp(npts,u_imag,1,v_imag,1,val,1,val,1);
1802  Vmath::Smul (npts,2.0,val,1,val,1);
1803  waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
1804 
1805  Vmath::Vadd(npts,der1,1,der2,1,der1,1);
1806 
1807  NekDouble rho = session->GetParameter("RHO");
1808  NekDouble Re = session->GetParameter("RE");
1809 cout<<"Re="<<Re<<endl;
1810  NekDouble waveForceMag = rho*rho*std::pow(Re,-1./3.);
1811 
1812  Array<OneD, Array<OneD, NekDouble > > vwiForcing;
1813  vwiForcing = Array<OneD, Array<OneD, NekDouble> > (2);
1814  vwiForcing[0] = Array<OneD, NekDouble> (2*ncoeffs);
1815  for(int i = 1; i < 2; ++i)
1816  {
1817  vwiForcing[i] = vwiForcing[i-1] + ncoeffs;
1818  }
1819 
1820  if(projectfield!=0)
1821  {
1822  waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,vwiForcing[0]);
1823  }
1824  else
1825  {
1826  waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,vwiForcing[0]);
1827  }
1828 
1829 cout<<"waveforcemag="<<waveForceMag<<endl;
1830  Vmath::Smul(ncoeffs,-waveForceMag, vwiForcing[0],1, vwiForcing[0],1);
1831 
1832  // d/dx(u v* + u* v)
1833  waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
1834 
1835  // d/dy(v v* + v* v)
1836  Vmath::Vmul(npts,v_real,1,v_real,1,val,1);
1837  Vmath::Vvtvp(npts,v_imag,1,v_imag,1,val,1,val,1);
1838  Vmath::Smul (npts,2.0,val,1,val,1);
1839  waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
1840 
1841  Vmath::Vadd(npts,der1,1,der2,1,der1,1);
1842 
1843  if(projectfield!=0)
1844  {
1845  waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,vwiForcing[1]);
1846  }
1847  else
1848  {
1849  waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,vwiForcing[1]);
1850  }
1851 
1852  Vmath::Smul(ncoeffs,- waveForceMag, vwiForcing[1],1, vwiForcing[1],1);
1853 
1854  if(symm==true)
1855  {
1856  cout<<"symmetrization"<<endl;
1857 
1858 
1859  waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(vwiForcing[0],der1);
1860  for(int i = 0; i < npts; ++i)
1861  {
1862  val[i] = 0.5*(der1[i] - der1[Refindices[i]]);
1863  }
1864 
1865  waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, vwiForcing[0]);
1866 
1867 
1868  waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(vwiForcing[1],der1);
1869  for(int i = 0; i < npts; ++i)
1870  {
1871  val[i] = 0.5*(der1[i] - der1[Refindices[i]]);
1872  }
1873 
1874  waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, vwiForcing[1]);
1875 
1876  }
1877 
1878  // dump output
1879  Array<OneD, std::string> variables(2);
1880  Array<OneD, Array<OneD, NekDouble> > outfield(2);
1881  variables[0] = "u";
1882  variables[1] = "v";
1883  outfield[0] = vwiForcing[0];
1884  outfield[1] = vwiForcing[1];
1885 
1886  string sessionName = fieldfile.substr(0,fieldfile.find_last_of("."));
1887  std::string outname = sessionName + ".vwi";
1888 
1889  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
1890  = waveVelocities[0]->GetPlane(0)->GetFieldDefinitions();
1891  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
1892  //Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(outfield.num_elements());
1893  string var;
1894  for(int j=0; j< vwiForcing.num_elements(); ++j)
1895  {
1896  //fieldcoeffs =vwiForcing;
1897  for(int i=0; i< FieldDef.size(); i++)
1898  {
1899  //var = vSession->GetVariable(j);
1900  var = variables[j];
1901  FieldDef[i]->m_fields.push_back(var);
1902  waveVelocities[0]->GetPlane(0)->AppendFieldData(FieldDef[i], FieldData[i], vwiForcing[j]);
1903  }
1904  }
1905  LibUtilities::Write(outname,FieldDef,FieldData);
1906 
1907 
1908  }
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::LibUtilities::Write
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap)
Write a field file in serial only.
Definition: FieldIO.cpp:81
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
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::iterator, Vmath::Smul(), Vmath::Vcopy(), Vmath::Vmul(), and Vmath::Zero().

Referenced by main().

701  {
702  //1 I regions is expected: (the layer is the last region)
703  ASSERTL0(Iregions.num_elements()==3, "something wrong with the number of I layers");
704  //int Iup =Iregions[0];
705  int Ireg =Iregions[2];
706 cout<<"layer region="<<Ireg<<endl;
707 
708 
709  //take the pressure from infields
710 
711 
712  int pfield;
713  if(infields.num_elements()==2)
714  {
715  pfield =1;
716  }
717  else
718  {
719  pfield =3;
720  }
721 
722  MultiRegions::ExpListSharedPtr pressure =infields[pfield];
723 
724  Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =infields[0]->GetBndCondExpansions();
725  MultiRegions::ExpListSharedPtr Ilayer = Iexp[Ireg];
726 
727  int nq1D = Ilayer->GetPlane(0)->GetTotPoints();
728 //decomment int nq1D = Ilayer->GetTotPoints();
729  Array<OneD,NekDouble> x0d(nq1D);
730  Array<OneD,NekDouble> x1d(nq1D);
731  Array<OneD,NekDouble> x2d(nq1D);
732 
733  Ilayer->GetPlane(0)->GetCoords(x0d,x1d,x2d);
734 //decomment Ilayer->GetCoords(x0d,x1d,x2d);
735 //start
736  //take the streak data along the layers (getbondcondexpansions does not work!!!)
737  //Array<OneD, MultiRegions::ExpListSharedPtr> Istreak = streak->GetBndCondExpansions();
738  //Array<OneD, MultiRegions::ExpListSharedPtr> Istreak =
739  // infields[1]->GetPlane(0)->GetBndCondExpansions();
740  //MultiRegions::ExpListSharedPtr Istreakup = Istreak[Iup];
741  //MultiRegions::ExpListSharedPtr Istreakdown = Istreak[Idown];
742 
743 
744 
745 
746  int np = streak->GetTotPoints();
747  Array<OneD, NekDouble> gradx(np);
748  Array<OneD, NekDouble> grady(np);
749 
750  int totcoeffs = streak->GetNcoeffs();
751  Array<OneD, NekDouble> gradxcoeffs(totcoeffs);
752  Array<OneD, NekDouble> gradycoeffs(totcoeffs);
753 
754  streak->PhysDeriv(streak->GetPhys(), gradx, grady);
755  streak->FwdTrans_IterPerExp(gradx, gradxcoeffs);
756  streak->FwdTrans_IterPerExp(grady, gradycoeffs);
757 
758 
759 
760 
761  Array<OneD, Array<OneD, NekDouble> > normals;
762  normals = Array<OneD, Array<OneD, NekDouble> >(2);
763  Array<OneD, Array<OneD, NekDouble> > tangents;
764  tangents = Array<OneD, Array<OneD, NekDouble> >(2);
765 
766 
767 
768 
769  Array<OneD, NekDouble> nx(nq1D);
770  Array<OneD, NekDouble> ny(nq1D);
771  Array<OneD, NekDouble> tx(nq1D);
772  Array<OneD, NekDouble> ty(nq1D);
773 
774 
775 
776 
777 
778 
779  const SpatialDomains::BoundaryRegionCollection &bregions
780  = bcs.GetBoundaryRegions();
781  SpatialDomains::Composite Icompreg;
782 
783  SpatialDomains::SegGeomSharedPtr segmentGeomreg;
784  SpatialDomains::ElementEdgeVectorSharedPtr elementreg;
785  StdRegions::Orientation orientreg;
786  SpatialDomains::BoundaryRegion::iterator regIt;
787  Array<OneD, unsigned int> bmapreg;
788  //fields to store the Re Im parts in physical space:
789  Array<OneD,NekDouble> Rephysreg (nq1D);
790  Array<OneD,NekDouble> Imphysreg (nq1D);
791 
792  //jacobian for the layer, std_deriv
793  Array<OneD, NekDouble> Jac (nq1D);
794  Array<OneD, NekDouble> dPrestd (nq1D);
795 
796  //Array<OneD, NekDouble> stphysup (nq1D);
797  Array<OneD, NekDouble> stphysreg (nq1D);
798  Array<OneD, NekDouble> stphysgradxreg (nq1D);
799  Array<OneD, NekDouble> stphysgradyreg (nq1D);
800 
801  //metrics definitions
802  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > deriv;
803  Array<OneD, NekDouble> dxde1 (nq1D);
804  Array<OneD, NekDouble> dyde1 (nq1D);
805  Array<OneD, NekDouble> dxde2 (nq1D);
806  Array<OneD, NekDouble> dyde2 (nq1D);
807  Array<OneD, NekDouble> gmat0 (nq1D);
808  Array<OneD, NekDouble> gmat3 (nq1D);
809  Array<OneD, NekDouble> gmat1 (nq1D);
810  Array<OneD, NekDouble> gmat2 (nq1D);
811 
812  Array<OneD, int> Elmtid(Ilayer->GetPlane(0)->GetExpSize());
813  Array<OneD, int> Edgeid(Ilayer->GetPlane(0)->GetExpSize());
814  //nummodes*nedges_layerdown
815  int Nregcoeffs;
816 
817  for(
818  regIt= bregions.find(Ireg)->second->begin();
819  regIt !=bregions.find(Ireg)->second->end();
820  ++regIt)
821  {
822  Icompreg=regIt->second;
823 
824 
825 
826 
827 
828  //assuming that ncoeffs is equal for all the edges (k=0)
829  int nummodes = Ilayer->GetExp(0)->GetNcoeffs();
830  Array<OneD, int> offsetup(Icompreg->size());
831  offsetup[0]=0;
832  int offsetregIExp;
833 
834  //array coeffs:
835  Nregcoeffs = Icompreg->size()*nummodes;
836  Array<OneD, NekDouble> Recoeffsreg (Nregcoeffs);
837  Array<OneD, NekDouble> Imcoeffsreg (Nregcoeffs);
838 
839  Array<OneD, NekDouble> stcoeffsreg (Nregcoeffs);
840  Array<OneD, NekDouble> stgradxcoeffsreg (Nregcoeffs);
841  Array<OneD, NekDouble> stgradycoeffsreg (Nregcoeffs);
842 
843 
844 
845  //define nq points for each edge
846  int nqedge =nq1D/Icompreg->size();
847 
848  Array<OneD,NekDouble> x0edge(nqedge);
849  Array<OneD,NekDouble> x1edge(nqedge);
850  Array<OneD, NekDouble> Rephysedgereg(nqedge);
851  Array<OneD, NekDouble> Imphysedgereg(nqedge);
852  Array<OneD, NekDouble> stgradxcoeffsedge(nqedge);
853  Array<OneD, NekDouble> stgradycoeffsedge(nqedge);
854 
855 
856  Array<TwoD, const NekDouble> gmat;
857  Array<OneD, NekDouble> gmat0_2D (nqedge*nqedge);
858  Array<OneD, NekDouble> gmat3_2D (nqedge*nqedge);
859  Array<OneD, NekDouble> gmat1_2D (nqedge*nqedge);
860  Array<OneD, NekDouble> gmat2_2D (nqedge*nqedge);
861 
862 
863 
864 
865  for(int f=0; f<2; ++f)
866  {
867  normals[f] = Array<OneD, NekDouble>(nqedge);
868  tangents[f] = Array<OneD, NekDouble>(nqedge);
869  }
870 
871  for(int k=0; k< Icompreg->size(); k++)//loop over segments of each layer
872  {
873 
874  if(
875  !(segmentGeomreg = boost::dynamic_pointer_cast<
876  SpatialDomains::SegGeom>((*Icompreg)[k]))
877  )
878  { ASSERTL0(false, "dynamic cast to a SegGeom failed"); }
879 
880  int EIDreg = segmentGeomreg->GetEid();
881  offsetregIExp = Ilayer->GetPlane(0)->GetCoeff_Offset(k);
882 //cout<<"edge="<<EIDreg<<" CHECK OFFSET="<<offsetregIExp<<endl;
883  elementreg = boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(mesh)
884  ->GetElementsFromEdge(segmentGeomreg);
885  int elmtidreg = ((*elementreg)[0]->m_Element)->GetGlobalID();
886  int dimelmtreg = ((*elementreg)[0]->m_Element)->GetNumEdges();
887  Elmtid[k] = elmtidreg;
888  Edgeid[k] = EIDreg;
889  orientreg =(boost::dynamic_pointer_cast<SpatialDomains::Geometry2D>
890  ((*elementreg)[0]->m_Element))->GetEorient((*elementreg)[0]
891  ->m_EdgeIndx);
892 
893  gmat = pressure->GetPlane(0)->GetExp(elmtidreg)->GetMetricInfo()->GetGmat(
894  pressure->GetPlane(0)->GetExp(elmtidreg)->GetPointsKeys());
895  int nq2D = nqedge*nqedge;
896 
897  // setup map between global and local ids
898  map<int, int>EdgeGIDreg;
899  int id,cnt,f;
900  for(cnt=f=0; f<dimelmtreg; f++)
901  {
902  id = (boost::dynamic_pointer_cast<SpatialDomains::Geometry2D>
903  ((*elementreg)[0]->m_Element))->GetEid(f);
904  EdgeGIDreg[id]=cnt++;
905  }
906 //cout<<" GIDreg="<<elmtidreg<<" num edges in this element="<<dimelmtreg<<endl;
907  int localidedge = EdgeGIDreg.find(EIDreg)->second;
908 //cout<<"locaidedge from m_EdgeInx ="<<localidedge<<endl;
909  localidedge = (*elementreg)[0]->m_EdgeIndx;
910 //cout<<"check map local id of Eid="<<EIDreg<<" is="<<localidedge<<endl;
911 
912 
913  //set bmaps -------------
914  Array<OneD, unsigned int > bmapedgereg;
915  Array<OneD, int > signreg;
916  pressure->GetExp(elmtidreg)->GetEdgeToElementMap(EdgeGIDreg.find(EIDreg)->second,
917  orientreg,bmapedgereg,signreg);
918 
919 //cout<<"SIGN down="<<signreg[0]<<endl;
920  //Extracting the coeffs...
921  Array<OneD, NekDouble> Recoeffs = pressure->GetPlane(0)->GetCoeffs();
922  Array<OneD, NekDouble> Imcoeffs = pressure->GetPlane(1)->GetCoeffs();
923  Array<OneD, NekDouble> stcoeffs = streak->GetCoeffs();
924  //OFFSET EDGE REFERING TO ELEMENT EXPLIST MISSING
925  int Reoffsetreg = pressure->GetPlane(0)->GetCoeff_Offset(elmtidreg);
926  int Imoffsetreg = pressure->GetPlane(1)->GetCoeff_Offset(elmtidreg);
927  int stoffsetreg = streak->GetCoeff_Offset(elmtidreg);
928 
929  //cout<<" Re offsetdown="<<Reoffsetreg<<" Im offset="<<Imoffsetreg<<endl;
930  //hypothesis: 2 planes (0 Re, 1 Im)
931  Array<OneD, NekDouble> Recoeffsedgereg(bmapedgereg.num_elements());
932  Array<OneD, NekDouble> Imcoeffsedgereg(bmapedgereg.num_elements());
933  Array<OneD, NekDouble> stcoeffsedgereg(bmapedgereg.num_elements());
934 
935 
936  //grad edge coeffs
937  //cout<<" Recoeffsedgeup num elements="<<Recoeffsedgereg.num_elements()
938  //<<" n bmapedge="<<bmapedgereg.num_elements()<<endl;
939 
940  for(int d=0; d<bmapedgereg.num_elements(); d++)
941  {
942  //pressure
943  Recoeffsedgereg[d]=Recoeffs[Reoffsetreg+bmapedgereg[d]];
944  Imcoeffsedgereg[d]=Imcoeffs[Imoffsetreg+bmapedgereg[d]];
945  Recoeffsreg[offsetregIExp +d] = Recoeffsedgereg[d];
946  Imcoeffsreg[offsetregIExp +d] = Imcoeffsedgereg[d];
947 
948  //streak coeffs
949  stcoeffsedgereg[d] = stcoeffs[stoffsetreg + bmapedgereg[d]];
950  stgradxcoeffsedge[d] = gradxcoeffs[stoffsetreg + bmapedgereg[d]];
951  stgradycoeffsedge[d] = gradycoeffs[stoffsetreg + bmapedgereg[d]];
952  stcoeffsreg[offsetregIExp +d] = stcoeffsedgereg[d];
953  stgradxcoeffsreg[offsetregIExp +d] = stgradxcoeffsedge[d];
954  stgradycoeffsreg[offsetregIExp +d] = stgradycoeffsedge[d];
955 //cout<<"results: Re="<<Recoeffsedgereg[d]<<" Im="<<Imcoeffsedgereg[d]<<endl;
956 
957  }
958 
959 
960 
961 
962 
963 
964 
965  //store the jacobian, std_deriv
966  Array<OneD, NekDouble> jacedge (nqedge);
967  Array<OneD, NekDouble> Pre_edge (nqedge);
968  Array<OneD, NekDouble> dPre_edge (nqedge);
969  jacedge = Ilayer->GetPlane(0)->GetExp(k)->GetMetricInfo()->GetJac(Ilayer->GetPlane(0)->GetExp(k)->GetPointsKeys());
970 
971 
972  //check if the metrix is the same for streak and Ilayer obj
973  StdRegions::StdExpansion1DSharedPtr edgestdexp;
974  edgestdexp = boost::dynamic_pointer_cast<StdRegions::StdExpansion1D> (
975  Ilayer->GetPlane(0)->GetExp(k) );
976  edgestdexp->BwdTrans(Recoeffsedgereg, Pre_edge);
977 //cout<<"call PhysTensorDERIV......"<<endl;
978  edgestdexp->StdPhysDeriv(Pre_edge,
979  dPre_edge);
980  for(int q=0; q<nqedge; q++)
981  {
982  Jac[k*nqedge +q]= jacedge[q];
983  dPrestd[k*nqedge +q] = dPre_edge[q];
984 
985  }
986 
987  //cout<<"Nlayercoeffs="<<Recoeffsreg.num_elements()<<" Nlayerphys="<<Rephysreg.num_elements()<<endl;
988 
989  //cout<<"extract normal for edge="<<k<<endl;
990  //cout<<"extract tangent for edge="<<k<<endl;
991  //k is the number of the edge according to the composite list
992  LocalRegions::Expansion1DSharedPtr edgeexp =
993  boost::dynamic_pointer_cast<LocalRegions::Expansion1D>
994  (Ilayer->GetPlane(0)->GetExp(k) );
995  int localEid = edgeexp->GetLeftAdjacentElementEdge();
996  normals = edgeexp->
997  GetLeftAdjacentElementExp()->GetEdgeNormal(localEid);
998  LocalRegions::SegExpSharedPtr bndSegExp =
999  boost::dynamic_pointer_cast<LocalRegions::SegExp>(Ilayer->GetPlane(0)->GetExp(k));
1000  // This function call has be deprecated -- cc
1001  //tangents = (bndSegExp)->GetMetricInfo()->GetEdgeTangent();
1002  int physoffsetregIExp = Ilayer->GetPlane(0)->GetPhys_Offset(k);
1003  for(int e=0; e< nqedge; e++)
1004  {
1005 
1006  //nx[k*nqedge +e] = normals[0][e];
1007  //ny[k*nqedge +e] = normals[1][e];
1008  //tx[k*nqedge +e] = tangents[0][e];
1009  //ty[k*nqedge +e] = tangents[1][e];
1010 //cout<<"offset="<<offsetregIExp<<" nx="<<normals[0][e]<<" ny="<<normals[1][e]<<" e="<<e<<endl;
1011  nx[physoffsetregIExp +e] = normals[0][e];
1012  ny[physoffsetregIExp +e] = normals[1][e];
1013  tx[physoffsetregIExp +e] = tangents[0][e];
1014  ty[physoffsetregIExp +e] = tangents[1][e];
1015 
1016 //cout<<"offset="<<offsetregIExp<<" nx="<<normals[0][e]<<" ny="<<normals[1][e]<<" e="<<e<<endl;
1017 //cout<<offsetregIExp +e<<" "<<x0d[offsetregIExp +e]<<" nx="<<nx[offsetregIExp +e]<<" ny="<<ny[offsetregIExp +e]<<" e="<<e<<endl;
1018 //cout<<"tannn tx="<<tangents[0][e]<<" ty="<<tangents[1][e]<<endl;
1019  }
1020 
1021 
1022 
1023  //extract metrics factors...
1024  //int nq2D= nqedge*nqedge;
1025  LibUtilities::PointsKeyVector ptsKeys = pressure->GetPlane(0)->GetExp(elmtidreg)->GetPointsKeys();
1026  deriv = pressure->GetPlane(0)->GetExp(elmtidreg)->GetMetricInfo()->GetDeriv(ptsKeys);
1027 
1028  Array<OneD, NekDouble> derivelmt(nq2D);
1029  int offsetregphys = Ilayer->GetPlane(0)->GetPhys_Offset(k);
1030  Array<OneD, NekDouble> deriv_i(nqedge);
1031  Vmath::Vcopy(nq2D, &(deriv[0][0][0]),1, &(derivelmt[0]),1);
1032  pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
1033  Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dxde1[offsetregphys]),1);
1034  Vmath::Zero(nqedge, deriv_i,1);
1035  Vmath::Zero(nq2D, derivelmt,1);
1036  Vmath::Vcopy(nq2D, &(deriv[1][1][0]),1, &(derivelmt[0]),1);
1037  pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
1038  Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dyde2[offsetregphys]),1);
1039  Vmath::Zero(nqedge, deriv_i,1);
1040  Vmath::Zero(nq2D, derivelmt,1);
1041  Vmath::Vcopy(nq2D, &(deriv[0][1][0]),1, &(derivelmt[0]),1);
1042  pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
1043  Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dyde1[offsetregphys]),1);
1044  Vmath::Zero(nqedge, deriv_i,1);
1045  Vmath::Zero(nq2D, derivelmt,1);
1046  Vmath::Vcopy(nq2D, &(deriv[1][0][0]),1, &(derivelmt[0]),1);
1047  pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge, derivelmt, deriv_i);
1048  Vmath::Vcopy(nqedge, &(deriv_i[0]),1, &(dxde2[offsetregphys]),1);
1049 
1050  Array<OneD, NekDouble> gmati (nqedge);
1051  Vmath::Vcopy(nq2D, &(gmat[0][0]),1, &(gmat0_2D[0]),1);
1052  Vmath::Vcopy(nq2D, &(gmat[3][0]),1, &(gmat3_2D[0]),1);
1053  Vmath::Vcopy(nq2D, &(gmat[1][0]),1, &(gmat1_2D[0]),1);
1054  Vmath::Vcopy(nq2D, &(gmat[2][0]),1, &(gmat2_2D[0]),1);
1055  pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat0_2D, gmati);
1056  Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat0[offsetregphys]),1);
1057  Vmath::Zero(nqedge, gmati,1);
1058  pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat1_2D, gmati);
1059  Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat1[offsetregphys]),1);
1060  Vmath::Zero(nqedge, gmati,1);
1061  pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat2_2D, gmati);
1062  Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat2[offsetregphys]),1);
1063  Vmath::Zero(nqedge, gmati,1);
1064  pressure->GetPlane(0)->GetExp(elmtidreg)->GetEdgePhysVals(localidedge,gmat3_2D, gmati);
1065  Vmath::Vcopy(nqedge, &(gmati[0]),1, &(gmat3[offsetregphys]),1);
1066 
1067  }
1068  //bwd transform:
1069  Ilayer->GetPlane(0)->BwdTrans(Recoeffsreg, Rephysreg);
1070  Ilayer->GetPlane(1)->BwdTrans(Imcoeffsreg, Imphysreg);
1071  //streak phys values:
1072  Ilayer->GetPlane(0)->BwdTrans(stcoeffsreg, stphysreg);
1073  Ilayer->GetPlane(0)->BwdTrans_IterPerExp(stgradxcoeffsreg, stphysgradxreg);
1074  Ilayer->GetPlane(0)->BwdTrans_IterPerExp(stgradycoeffsreg, stphysgradyreg);
1075  }
1076 
1077 
1078  //TANGENTS WRONG WITH 3DHOMOGENEOUS 1D...
1079  //calculate the curvature:
1080 
1081 
1082  //attempt to smooth the tangent components:
1083  Array<OneD, NekDouble> tcoeffs (Nregcoeffs,0.0);
1084  //outfieldx->FwdTrans(tx, tcoeffs);
1085  //outfieldx->BwdTrans(tcoeffs, tx);
1086  //Vmath::Zero(Nregcoeffs, tcoeffs,1);
1087  //outfieldx->FwdTrans(ty, tcoeffs);
1088  //outfieldx->BwdTrans(tcoeffs, ty);
1089  //Vmath::Zero(Nregcoeffs, tcoeffs,1);
1090  Array<OneD, NekDouble> dtx(nq1D,0.0);
1091  Array<OneD, NekDouble> dty(nq1D,0.0);
1092  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, tx, dtx);
1093  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, ty, dty);
1094  //attempt to smooth the tangent derivative:
1095  outfieldx->FwdTrans_IterPerExp(dtx, tcoeffs);
1096  outfieldx->BwdTrans_IterPerExp(tcoeffs, dtx);
1097  Vmath::Zero(Nregcoeffs, tcoeffs,1);
1098  outfieldx->FwdTrans_IterPerExp(dty, tcoeffs);
1099  outfieldx->BwdTrans_IterPerExp(tcoeffs, dty);
1100  Array<OneD, NekDouble> f_z(nq1D,0.0);
1101  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, x1d, f_z);
1102  Array<OneD, NekDouble> fcoeffs(Nregcoeffs,0.0);
1103  outfieldx->FwdTrans(f_z, fcoeffs);
1104  outfieldx->BwdTrans(fcoeffs, f_z);
1105  Vmath::Zero(Nregcoeffs, fcoeffs,1);
1106  Array<OneD, NekDouble> f_zz(nq1D,0.0);
1107 
1108  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, f_z, f_zz);
1109  outfieldx->FwdTrans(f_zz, fcoeffs);
1110  outfieldx->BwdTrans(fcoeffs, f_zz);
1111  Array<OneD, NekDouble> delta(nq1D);
1112  Array<OneD, NekDouble> curv_unsigned(nq1D,0.0);
1113  Array<OneD, NekDouble> curv(nq1D,0.0);
1114 
1115 
1116  //cout<<"x y tx ty dtx dty curv"<<endl;
1117  for(int t=0; t<nq1D; t++)
1118  {
1119 
1120  curv_unsigned[t]=sqrt(dtx[t]*dtx[t] +dty[t]*dty[t]);
1121  delta[t] = 1+f_z[t]*f_z[t];
1122  //Hall's definition...
1123  curv[t] = f_zz[t]/sqrt(delta[t]*delta[t]*delta[t]);
1124  //cout<<setw(13)<<x0d[t]<<" "<<setw(13)<<x1d[t]<<" "<<setw(13)<<tx[t]<<setw(13)<<" "<<ty[t]<<" "
1125  //<<nx[t]<<" "<<ny[t]<<" "<<setw(13)<<dtx[t]<<" "<<dty[t]<<" "<<curv[t]<<endl;
1126  }
1127 
1128 
1129  //attempt to smooth the curvature
1130  Array<OneD, NekDouble> curv_coeffs (Nregcoeffs);
1131  outfieldx->FwdTrans(curv, curv_coeffs);
1132  outfieldx->BwdTrans(curv_coeffs, curv);
1133  //P square before norm
1134  Array<OneD, NekDouble> P2reg (nq1D, 0.0);
1135  for(int e=0; e< nq1D; e++)
1136  {
1137  P2reg[e] = Rephysreg[e]*Rephysreg[e] + Imphysreg[e]*Imphysreg[e];
1138  }
1139 
1140 
1141  //calculate the s variable:
1142  Array<OneD, NekDouble> s0d(nq1D,0.0);
1143 /*
1144  Array<OneD, NekDouble> tmp(nq1D);
1145  for( int e=0; e<nq1D; e++)
1146  {
1147  tmp[e] = sqrt(1+f_z[e]*f_z[e]);
1148 
1149  }
1150 */
1151 
1152  //normalise the pressure norm*sqrt[ (\int Psquare)/Area]=1 =>
1153  NekDouble norm2D;
1154 //////////////////////////////////////////////////2D norm
1155 
1156  Array<OneD, NekDouble> P2_2D(np, 0.0);
1157  for(int h=0; h< np; h++)
1158  {
1159  P2_2D[h] = pressure->GetPlane(0)->GetPhys()[h]*pressure->GetPlane(0)->GetPhys()[h]
1160  + pressure->GetPlane(1)->GetPhys()[h]*pressure->GetPlane(1)->GetPhys()[h];
1161  }
1162  NekDouble intP2_2D = pressure->GetPlane(0)->PhysIntegral(P2_2D);
1163 
1164 
1165 
1166  NekDouble tmp = pressure->GetPlane(0)->L2(pressure->GetPlane(0)->GetPhys());
1167  norm2D = tmp*tmp;
1168  tmp = pressure->GetPlane(1)->L2(pressure->GetPlane(1)->GetPhys());
1169  norm2D += tmp*tmp;
1170 
1171  Array<OneD, NekDouble> I (2*np,1.0);
1172  //Vmath::Fill(2*np,1.0,I,1);
1173  NekDouble Area = pressure->GetPlane(0)->PhysIntegral(I);
1174  norm2D = sqrt(Area/norm2D);
1175 
1176  Array<OneD, NekDouble> I_int(np,1.0);
1177  NekDouble Area1 = pressure->GetPlane(0)->PhysIntegral(I_int);
1178  NekDouble normint2D = sqrt(Area1/intP2_2D);
1179 cout<<"norm2D="<<norm2D<<" area="<<Area1<<" intP2_2D="<<intP2_2D<<" normint2D="<<normint2D<<endl;
1180 
1181 ///////////////////////////////////////////////////////////1D norm
1182  NekDouble norm;
1183  Array<OneD, NekDouble> sqrtlen(nq1D);
1184  for(int u=0; u<nq1D; u++)
1185  {
1186  sqrtlen[u] = sqrt(1+f_z[u]*f_z[u]);
1187  s0d[u] = Ilayer->GetPlane(0)->PhysIntegral(sqrtlen);
1188  }
1189  NekDouble length = Ilayer->GetPlane(0)->PhysIntegral(sqrtlen);
1190  NekDouble int1D = Ilayer->GetPlane(0)->PhysIntegral(P2reg);
1191  norm = sqrt(length/int1D);
1192 
1193  NekDouble scal = Area/length;
1194 cout<<"norm1D="<<norm<<" norm2D/norm1D="<<norm2D/norm<<endl;
1195 cout<<"scal="<<scal<<endl;
1196  //norm*pressure
1197  Vmath::Smul(nq1D,norm,Rephysreg,1,Rephysreg,1);
1198  Vmath::Smul(nq1D,norm,Imphysreg,1,Imphysreg,1);
1199 
1200 
1201  //P square after norm
1202  Array<OneD, NekDouble> P2reg_aft (nq1D, 0.0);
1203  for(int e=0; e< nq1D; e++)
1204  {
1205  P2reg_aft[e] = Rephysreg[e]*Rephysreg[e] + Imphysreg[e]*Imphysreg[e];
1206  }
1207  //print out the fields
1208 
1209 
1210 //cout<<"derivative of the pressure"<<endl;
1211  Array<OneD, NekDouble> dP_re = Array<OneD, NekDouble>(nq1D,0.0);
1212  Array<OneD, NekDouble> dP_im = Array<OneD, NekDouble>(nq1D,0.0);
1213  Array<OneD, NekDouble> dP_square2 = Array<OneD, NekDouble>(nq1D,0.0);
1214 
1215  //attempt to smooth the normal
1216  //Array<OneD, NekDouble> ncoeffs(Nregcoeffs,0.0);
1217  //outfieldx->FwdTrans(nx, ncoeffs);
1218  //outfieldx->BwdTrans(ncoeffs, nx);
1219  //Vmath::Zero(Nregcoeffs, ncoeffs,1);
1220  //outfieldx->FwdTrans(ny, ncoeffs);
1221  //outfieldx->BwdTrans(ncoeffs, ny);
1222  //Vmath::Zero(Nregcoeffs, ncoeffs,1);
1223 
1224  Array<OneD, NekDouble> dUreg (nq1D);
1225  //Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eN, stphysreg,dUreg);
1226  for(int w=0; w<nq1D; w++)
1227  {
1228  dUreg[w] = nx[w]*stphysgradxreg[w] +ny[w]*stphysgradyreg[w];
1229  }
1230 
1231 
1232 
1233  Array<OneD, NekDouble> dUcoeffs(Nregcoeffs,0.0);
1234  //outfieldx->FwdTrans(dUreg, dUcoeffs);
1235  //outfieldx->BwdTrans(dUcoeffs, dUreg);
1236 
1237 
1238 
1239  Array<OneD, NekDouble> jaccoeffs(Nregcoeffs);
1240  Array<OneD, NekDouble> Jacreg (nq1D);
1241  Array<OneD, NekDouble> dPrecoeffs (Nregcoeffs);
1242  //attempt to smooth the jac and the dPre
1243  outfieldx->FwdTrans(Jac, jaccoeffs);
1244  outfieldx->BwdTrans(jaccoeffs, Jacreg);
1245  //outfieldx->FwdTrans(dPre, dPrecoeffs);
1246  //outfieldx->BwdTrans(dPrecoeffs, dPre);
1247 
1248 
1249  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS,Rephysreg,dP_re);
1250  Ilayer->GetPlane(1)->PhysDeriv(MultiRegions::eS,Imphysreg,dP_im);
1251  //attempt to smooth the derivative:
1252  Array<OneD, NekDouble> dP_re_coeffs (Nregcoeffs,0.0);
1253  Array<OneD, NekDouble> dP_im_coeffs (Nregcoeffs,0.0);
1254  Array<OneD, NekDouble> dP_imtest (nq1D);
1255  outfieldx->FwdTrans(dP_re, dP_re_coeffs);
1256  outfieldx->FwdTrans(dP_im, dP_im_coeffs);
1257  outfieldx->BwdTrans(dP_re_coeffs, dP_re);
1258  outfieldx->BwdTrans(dP_im_coeffs, dP_im);
1259 
1260 
1261 
1262  Array<OneD, NekDouble> dP_square = Array<OneD, NekDouble>(nq1D, 0.0);
1263 cout<<"x"<<" P_re"<<" dP_re"<<" streak"<<" dstreak"<<" pjump"<<endl;
1264  for(int s=0; s<nq1D; s++)
1265  {
1266  dP_square[s] = dP_re[s]*dP_re[s] +dP_im[s]*dP_im[s];
1267  }
1268 
1269 
1270 
1271 
1272  //attempt to smooth the pressure
1273  Array<OneD, NekDouble> dPsquare_coeffs (Nregcoeffs,0.0);
1274  outfieldx->FwdTrans(dP_square, dPsquare_coeffs);
1275  outfieldx->BwdTrans(dPsquare_coeffs, dP_square);
1276 
1277 
1278  Array<OneD, NekDouble> mu53 (nq1D,0.0);
1279  Array<OneD, NekDouble> d2v (nq1D,0.0);
1280  Array<OneD, NekDouble> prod (nq1D,0.0);
1281 
1282  double pow = 1.0/3.0;
1283  double base;
1284  for(int y=0; y<nq1D; y++)
1285  {
1286  base =dUreg[y];
1287  mu53[y] = std::pow ((base*base),pow);
1288  mu53[y] = 1/(base*mu53[y]);
1289  //prod[y] = mu53[y]*dP_square[y];
1290  }
1291  //attempting to smooth mu53
1292  Array<OneD, NekDouble> mucoeffs(Nregcoeffs,0.0);
1293  outfieldx->FwdTrans(mu53, mucoeffs);
1294  outfieldx->BwdTrans(mucoeffs, mu53);
1295  Vmath::Vmul(nq1D, mu53,1, dP_square,1, prod,1);
1296 
1297  //attempting to smooth field...
1298  Array<OneD, NekDouble> prod_coeffs (Nregcoeffs,0.0);
1299  outfieldx->FwdTrans(prod, prod_coeffs);
1300  outfieldx->BwdTrans(prod_coeffs, prod);
1301  Vmath::Zero(Nregcoeffs,prod_coeffs,1);
1302 
1303  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, prod ,d2v);
1304  //attempting to smooth the vjump
1305  Vmath::Zero( Nregcoeffs, prod_coeffs,1);
1306  outfieldx->FwdTrans(d2v, prod_coeffs);
1307  outfieldx->BwdTrans(prod_coeffs, d2v);
1308  //test prod deriv... add background..
1309  //Array<OneD, NekDouble> prod_b(nq1D);
1310  //Array<OneD, NekDouble> dersprod_b(nq1D);
1311  //Vmath::Sadd(nq1D,0.001, prod,1,prod_b,1);
1312  //Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, prod_b ,dersprod_b);
1313 
1314 
1315  Array<OneD, NekDouble> pjump(nq1D);
1316  Array<OneD, NekDouble> vjump(nq1D);
1317 
1318 
1319  //test a sin curve....
1320 /*
1321  Array<OneD, NekDouble> fun1D(nq1D);
1322  Array<OneD, NekDouble> derfun1D(nq1D);
1323  NekDouble pi =3.14159265;
1324  for(int e=0; e<nq1D; e++)
1325  {
1326  fun1D[e]= std::sin((x0d[e]));
1327  }
1328  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, fun1D, derfun1D);
1329 
1330  //add noise to a function
1331  Array<OneD, NekDouble> fun1D_b(nq1D);
1332  Array<OneD, NekDouble> dersfun1D_b(nq1D);
1333  int cnt=0;
1334  int cnt1=1;
1335 
1336  for(int v=0; v<nq1D/2; v++)
1337  {
1338 
1339  fun1D_b[cnt] = fun1D[cnt] +0.000000;
1340 //cout<<cnt<<"x="<<x0d[cnt]<<" fun1D_b="<<fun1D_b[cnt]<<" fun1D="<<fun1D[cnt]<<" diff="<<
1341 //fun1D[cnt]-fun1D_b[cnt]<<endl;
1342  //ASSERTL0((fun1D_b[cnt]-fun1D[cnt])==0.0001, "problem");
1343  cnt = cnt +2;
1344  fun1D_b[cnt1] = fun1D[cnt1] -0.000000;
1345 //cout<<cnt1<<"x="<<x0d[cnt1]<<" fun1D_b="<<fun1D_b[cnt1]<<" fun1D="<<fun1D[cnt1]<<" diff="<<
1346 //fun1D[cnt1]-fun1D_b[cnt1]<<endl;
1347  //ASSERTL0((fun1D_b[cnt1]-fun1D[cnt1])==-0.0001, "problem");
1348  cnt1 = cnt1 +2;
1349  if( (v==((nq1D/2) -1)) && (nq1D%2==1)
1350  )
1351  {
1352  fun1D_b[cnt] = fun1D[cnt] +0.000000;
1353 //cout<<cnt<<"x="<<x0d[cnt]<<" fun1D_b="<<fun1D_b[cnt]<<" fun1D="<<fun1D[cnt]<<" diff="<<
1354 //fun1D[cnt]-fun1D_b[cnt]<<endl;
1355  }
1356  }
1357 
1358  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, fun1D_b ,dersfun1D_b);
1359 */
1360 /*
1361 
1362  int nqedge =nq1D/Icompreg->size();
1363  Array<OneD, NekDouble> phys_edge (nqedge);
1364  Array<OneD, NekDouble> der_edge (nqedge);
1365  Array<OneD, NekDouble> der_fun1Db_std (nq1D);
1366  //extract stdderiv....
1367 
1368  for(int w=0; w<Icompreg->size() ; w++)
1369  {
1370  Vmath::Vcopy(nqedge,&(fun1D_b[w*nqedge]),1,&(phys_edge[0]),1);
1371  //check if the metrix is the same for streak and Ilayer obj
1372  StdRegions::StdExpansion1DSharedPtr edgestdexp;
1373  edgestdexp = boost::dynamic_pointer_cast<StdRegions::StdExpansion1D> (
1374  Ilayer->GetPlane(0)->GetExp(w) );
1375  edgestdexp->StdPhysDeriv(phys_edge,
1376  der_edge);
1377  Vmath::Vcopy(nqedge, &(der_edge[0]),1, &(der_fun1Db_std[w*nqedge]),1);
1378  Vmath::Zero(nqedge, phys_edge,1);
1379  Vmath::Zero(nqedge, der_edge,1);
1380 
1381  Vmath::Vcopy(nqedge, &(Rephysreg[w*nqedge]),1, &(phys_edge[0]),1);
1382  edgestdexp->StdPhysDeriv(phys_edge, der_edge);
1383  Vmath::Vcopy(nqedge, &(der_edge[0]),1, &(der_Pre_std[w*nqedge]),1);
1384  Vmath::Zero(nqedge, phys_edge,1);
1385  Vmath::Zero(nqedge, der_edge,1);
1386  Vmath::Vcopy(nqedge, &(fun1D[w*nqedge]),1, &(phys_edge[0]),1);
1387  edgestdexp->StdPhysDeriv(phys_edge,
1388  der_edge);
1389  Vmath::Vcopy(nqedge, &(der_edge[0]),1, &(derfun1D[w*nqedge]),1);
1390  Vmath::Zero(nqedge, phys_edge,1);
1391  Vmath::Zero(nqedge, der_edge,1);
1392 
1393  }
1394 
1395 */
1396 
1397  //final jump conditions
1398  //n0=2*pi*(2/3)^(2/3)*(-2/3)! = 2*pi*(2/3)^(2/3)*gamma(1/3)
1399  NekDouble n0 = 12.845424015;
1400  //use the session to get the values of rho,alpha...
1401  NekDouble rho ;
1402  //rho = 0.08;
1403  rho = session->GetParameter("RHO");
1404  //load alpha only if not manually specified
1405  if(alpha==0)
1406  {
1407  alpha = session->GetParameter("ALPHA");
1408  }
1409 cout<<"alpha="<<alpha<<endl;
1410 
1411  ASSERTL0(alpha!=0, "alpha cannot be 0");
1412  NekDouble alpha53;
1413  //alpha53=1;
1414  alpha53 = std::pow ((alpha*alpha),pow);
1415  alpha53 = 1/(alpha*alpha53);
1416  for(int c=0; c<nq1D; c++)
1417  {
1418  pjump[c] = -n0*alpha53*curv[c]*prod[c] ;
1419  }
1420  Vmath::Smul(nq1D, n0,d2v,1, vjump,1);
1421  Vmath::Smul(nq1D, alpha53,vjump,1, vjump,1);
1422 
1423 cout<<"alpha^-5/3="<<alpha53<<endl;
1424 
1425  Array<OneD, NekDouble> txcoeffs (Nregcoeffs);
1426  Array<OneD, NekDouble> tycoeffs (Nregcoeffs);
1427  //attempt to smooth the tangent components:
1428  outfieldx->FwdTrans(tx, txcoeffs);
1429  outfieldx->FwdTrans(ty, tycoeffs);
1430  outfieldy->BwdTrans(txcoeffs, tx);
1431  outfieldy->BwdTrans(tycoeffs, ty);
1432 cout<<"RHO=="<<rho<<endl;
1433 
1434 
1435 
1436 //end
1437 //PAY ATTENTION to the sign (vjump -/+ pjump)*t!!!!!!!!!!
1438  for(int j=0; j<nq1D; j++)
1439  {
1440 //start
1441 
1442  (outfieldx->UpdatePhys())[j] =
1443  rho*rho*(vjump[j]*tx[j]-pjump[j]);
1444  (outfieldy->UpdatePhys())[j] =
1445  rho*rho*(vjump[j]*ty[j]-pjump[j]);
1446 //cout<<x0d[j]<<" "<<curv[j]*prod[j]<<" "<<(outfieldx->GetPhys())[j]<<" "<<(outfieldy->GetPhys())[j]<<endl;
1447 //end
1448 
1449 //decomment
1450 /*
1451  (outfieldx->UpdatePhys())[j] =
1452  20*sin(2*x0d[j]);
1453  (outfieldy->UpdatePhys())[j] =
1454  20*2*cos(2*x0d[j])/3.14;
1455 */
1456  }
1457 
1458 
1459  //FINAL REFINEMENT:::
1460 //start
1461 
1462  Array<OneD, NekDouble> finalcoeffs (Nregcoeffs);
1463  outfieldx->FwdTrans(outfieldx->GetPhys(), finalcoeffs);
1464  outfieldx->BwdTrans(finalcoeffs, outfieldx->UpdatePhys());
1465  Vmath::Zero(Nregcoeffs,finalcoeffs,1);
1466  outfieldy->FwdTrans(outfieldy->GetPhys(), finalcoeffs);
1467  outfieldy->BwdTrans(finalcoeffs, outfieldy->UpdatePhys());
1468 
1469  //symmetrize the jumps conditions
1470  if(symm ==true)
1471  {
1472  Array<OneD, NekDouble> tmpx(nq1D);
1473  Array<OneD, NekDouble> tmpy(nq1D);
1474 
1475 cout<<"symmetrise the jump conditions"<<endl;
1476 
1477  for(int i = 0; i < nq1D; ++i)
1478  {
1479  tmpx[i] =0.5*(outfieldx->GetPhys()[i] - outfieldx->GetPhys()[Refindices[i]]);
1480  tmpy[i] =0.5*(outfieldy->GetPhys()[i] - outfieldy->GetPhys()[Refindices[i]]);
1481  }
1482  Vmath::Vcopy(nq1D, tmpx,1, outfieldx->UpdatePhys(),1);
1483  Vmath::Vcopy(nq1D, tmpy,1, outfieldy->UpdatePhys(),1);
1484  }
1485 
1486 /*
1487 
1488  //calculate J,K
1489  Array<OneD, NekDouble> lambda(nq1D);
1490  Array<OneD, NekDouble> a(nq1D);
1491  Array<OneD, NekDouble> dPre_dz(nq1D);
1492  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, Rephysreg, dPre_dz);
1493  Array<OneD, NekDouble> dPim_dz(nq1D);
1494  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, Imphysreg, dPim_dz);
1495  Array<OneD, NekDouble> dPre_dz2(nq1D);
1496  Vmath::Vmul(nq1D, dPre_dz,1, dPre_dz,1, dPre_dz2,1);
1497  Array<OneD, NekDouble> dPim_dz2(nq1D);
1498  Vmath::Vmul(nq1D, dPim_dz,1, dPim_dz,1, dPim_dz2,1);
1499  Array<OneD, NekDouble> dP_dz2(nq1D);
1500  Vmath::Vadd(nq1D, dPre_dz2,1,dPim_dz2,1,dP_dz2,1);
1501  Array<OneD, NekDouble> ddP_dz2z(nq1D);
1502  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, dP_dz2, ddP_dz2z);
1503 
1504 
1505  Array<OneD, NekDouble> delta5(nq1D);
1506  Array<OneD, NekDouble> a53(nq1D);
1507  pow = 1./3.;
1508  for(int e=0; e<nq1D; e++)
1509  {
1510  //delta[e] = 1+ f_z[e]*f_z[e];
1511  delta5[e] = delta[e]*delta[e]*delta[e]*delta[e]*delta[e];
1512  lambda[e] = dUreg[e]/sqrt(delta[e]);
1513  a[e] = lambda[e]*alpha/delta[e];
1514  a53[e] = std::pow( (a[e]*a[e]), pow);
1515 //cout<<"a^2/3="<<a53[e]<<endl;
1516  a53[e] = a[e]*a53[e];
1517  }
1518 
1519  Array<OneD, NekDouble> delta_z(nq1D);
1520  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, delta, delta_z);
1521  Array<OneD, NekDouble> a_z(nq1D);
1522  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eX, a, a_z);
1523 
1524  Array<OneD, NekDouble> J(nq1D);
1525  Array<OneD, NekDouble> K(nq1D);
1526  Array<OneD, NekDouble> f1(nq1D);
1527  Array<OneD, NekDouble> f2(nq1D);
1528  //NekDouble a53,delta5;
1529  for(int f=0; f<nq1D; f++)
1530  {
1531  //a53 = std::pow(a[f],5./3.);
1532  //delta5 = std::pow( delta[f], 5);
1533  J[f] =( (n0*rho*rho)/(a53[f]*delta5[f]) )*
1534  (
1535  ( (-7*delta_z[f])/(2*delta[f]) - (5*a_z[f] )/(3*a[f]) )*dP_dz2[f] +
1536  ddP_dz2z[f]
1537  );
1538  K[f] = -(n0*rho*rho)/(a53[f]*delta5[f])*f_zz[f]*dP_dz2[f];
1539  f1[f] = lambda[f]*(K[f]-delta[f]*f_z[f]*J[f]);
1540  f2[f] = -lambda[f]*(f_z[f]*K[f] +delta[f]*J[f]);
1541  }
1542 
1543 
1544 
1545 
1546  //test dermu
1547  Array<OneD, NekDouble> dermu53(nq1D);
1548  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, mu53, dermu53);
1549 
1550  //test d(dP_square)ds
1551  Array<OneD, NekDouble> ddP_square_ds(nq1D);
1552  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, dP_square, ddP_square_ds);
1553 */
1554 
1555 
1556  cout<<"length layer="<<length<<endl;
1557  //calc dUds
1558  Array<OneD, NekDouble> dUds(nq1D);
1559  Ilayer->GetPlane(0)->PhysDeriv(MultiRegions::eS, stphysreg, dUds);
1560 
1561  Array<OneD, int> indexjacwarn(nq1D,-1);
1562  NekDouble invjac2D;
1563  NekDouble jactest;
1564  NekDouble laytest=0;
1565  NekDouble basetest;
1566  for(int g=0; g<nq1D; g++)
1567  {
1568 
1569 
1570  //NBBB dPre is the stdDERIV!!!
1571  invjac2D = (gmat0[g]*tx[g] +gmat2[g]*ty[g]);
1572  jactest = (invjac2D + (1/Jac[g]))/2.;
1573  jactest = (invjac2D- (1/Jac[g]))/jactest;
1574  basetest = abs(curv[g]);
1575 /*
1576  cout<<x0d[g]<<" "<<
1577 //stphysreg[g]<<" "<<lambda[g]<<" "<<
1578 //delta[g]<<" "<<
1579 //f_z[g]<<" "<<a53[g]<<" "<<
1580 //delta5[g]<<" "<<delta_z[g]<<" "<<a_z[g]<<" "<<
1581 //" "<<ddP_dz2z[g]<<" "<<dP_dz2[g]<<" "<<
1582 //dP_square[g]<<" "<<dP_dz2[g]/delta[g]<<" "<<
1583 //ddP_square_ds[g]<<" "<<ddP_dz2z[g]/(delta[g]*sqrt(delta[g]))<<" "<<
1584 //ddP_dz2z[g]/delta[g]<<" "<<
1585 //a[g]<<" "<<a53[g]<<" "<<
1586 vjump[g]*tx[g]*sqrt(delta[g])<<" "<<
1587 J[g]<<" "<<K[g]<<" "<<pjump[g]<<" "<<
1588 f1[g]<<" "<<f2[g]<<" "<<
1589 f_zz[g]/(sqrt(delta5[g]))<<" "<<
1590 //f_zz[g]/sqrt(delta[g]*delta[g]*delta[g])<<" "
1591 curv[g]<<" "<<
1592 outfieldx->GetPhys()[g]<<" "
1593 <<outfieldy->GetPhys()[g]
1594 <<endl;
1595 */
1596 
1597  laytest = (abs(stphysreg[g])-cr);
1598 
1599 
1600 cout<<setw(14)<<x0d[g]<<" "<<setw(13)<<x1d[g]<<
1601 " "<<Rephysreg[g]<<" "<<dP_re[g]<<" "
1602 <<stphysreg[g]<<" "<<dUreg[g]<<" "<<mu53[g]<<" "<<curv[g]
1603 
1604 //<<" "<<-1.5*std::pow(abs(curv[g]),-1./5.)<<" "
1605 //<<" "<<std::pow(delta[g],-0.5)-0.45*ty[g]+(f_z[g]*f_z[g])/(std::pow(delta[g],0.5))<<" "
1606 
1607 <<" "<<f_z[g]<<" "
1608 <<f_zz[g]<<" "
1609 //<<-1.5*std::pow(abs(curv[g]/delta[g]),-1./5.)<<" "
1610 <<s0d[g]
1611 <<setw(13)<<" "<<dP_square[g]<<" "<<setw(13)<<prod[g]
1612 <<" "<<
1613 
1614 //(gmat0[g]*tx[g] +gmat2[g]*ty[g])<<" "<<1/Jac[g]
1615 //<<" "<<dermu53[g]
1616 //<<" "<<derfun1D[g]
1617 //nx[g]<<" "<<ny[g]<<" "<<
1618 tx[g]<<" "<<ty[g]<<" "<<
1619 //f1[g]<<" "<<f2[g]
1620 //prod_b[g]<<" "<<dersprod_b[g]
1621 //fun1D[g]<<" "<<derfun1D[g]<<" "<<
1622 //der_fun1Db_std[g]<<" "<<
1623 //fun1D_b[g]<<" "<<dersfun1D_b[g]
1624 pjump[g]<<setw(13)<<" "<<d2v[g]<<" "
1625 <<outfieldx->GetPhys()[g]<<" "
1626 <<outfieldy->GetPhys()[g]<<" "<<Imphysreg[g]<<" "<<dP_im[g]
1627 <<" "<<"AA"<<endl;
1628 
1629 //cout<<"laytest="<<laytest<<endl;
1630 //cout<<"(gmat0[g]*tx[g] +gmat2[g]*ty[g])="<<(gmat0[g]*tx[g] +gmat2[g]*ty[g])<<
1631 //" 1/jac1D="<<1./Jac[g]<<endl;
1632  ASSERTL0(invjac2D*Jac[g]>0, " sign jac problem..");
1633  //30% error is allowed!!!
1634  ASSERTL0(abs(jactest)<0.3, "jac 1D problem..");
1635  //save info of point where the error >20%
1636  if(jactest>=0.2)
1637  {
1638  indexjacwarn[g] =g;
1639  }
1640 
1641  }
1642  if(cr==0)
1643  {
1644  ASSERTL0(abs(laytest)<0.002, "critical layer wrong");
1645  }
1646  else
1647  {
1648  cout<<"WARNING: crit may be wrong"<<endl;
1649  }
1650 
1651 // gamma(1/3)= 2.6789385347077476337
1652  // need the tangents related to the expList1D outfieldx[region]
1653 
1654 
1655  for(int m=0; m<nq1D; m++)
1656  {
1657  if(indexjacwarn[m]!=-1)
1658  {
1659  cout<<"warning: point with jacerr>20% index="<<m<<" x="<<x0d[m]<<" y="<<x1d[m]<<endl;
1660  }
1661  }
1662 
1663  for(int a=0; a<Elmtid.num_elements(); a++)
1664  {
1665 cout<<"elmt id="<<Elmtid[a]<<" edge id="<<Edgeid[a]<<endl;
1666  }
1667 //end
1668  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:220
Nektar::LocalRegions::SegExpSharedPtr
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::SpatialDomains::SegGeomSharedPtr
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
Nektar::SpatialDomains::BoundaryRegionCollection
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:206
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::StdRegions::StdExpansion1DSharedPtr
boost::shared_ptr< StdExpansion1D > StdExpansion1DSharedPtr
Definition: StdExpansion1D.h:112
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SpatialDomains::Composite
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
Nektar::SpatialDomains::ElementEdgeVectorSharedPtr
boost::shared_ptr< ElementEdgeVector > ElementEdgeVectorSharedPtr
Definition: MeshGraph.h:134
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
Nektar::SpatialDomains::Geometry2D
2D geometry information
Definition: Geometry2D.h:65
Nektar::StdRegions::StdExpansion::BwdTrans
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space...
Definition: StdExpansion.h:525
Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:227
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::LocalRegions::Expansion1DSharedPtr
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
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().

638  {
639  int i,j;
640  //remember Exp is the streak so GetPlane(0) is not needed
641  Array<OneD, MultiRegions::ExpListSharedPtr> Iexp =Exp->GetBndCondExpansions();
642  MultiRegions::ExpListSharedPtr Ilayer = Iexp[Ireg];
643  int npts = Ilayer->GetNpoints();
644  Array<OneD, int> index(npts);
645 
646  Array<OneD, NekDouble> coord(2);
647  Array<OneD, NekDouble> coord_x(npts);
648  Array<OneD, NekDouble> coord_y(npts);
649 
650  //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
651  Ilayer->GetCoords(coord_x,coord_y);
652  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
653  NekDouble tol = NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
654  NekDouble xnew,ynew;
655 
656  int start = npts-1;
657 //cout<<"xmax="<<xmax<<endl;
658  for(i = 0; i < npts; ++i)
659  {
660 //cout<<"Ref index for x="<<coord_x[i]<<endl;
661  xnew = - coord_x[i] + xmax;
662  ynew = - coord_y[i];
663 
664  for(j = start; j >=0 ; --j)
665  {
666  if((coord_x[j]-xnew)*(coord_x[j]-xnew) < tol)
667  {
668  index[i] = j;
669  start = j;
670  break;
671  }
672  }
673 
674  if(j == -1)
675  {
676 
677  for(j = npts-1; j > start; --j)
678  {
679 
680  if((coord_x[j]-xnew)*(coord_x[j]-xnew) < tol)
681  {
682  index[i] = j;
683  break;
684  }
685  }
686  ASSERTL0(j != start,"Failsed to find matching point");
687  }
688  }
689  return index;
690  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::NekConstants::kGeomFactorsTol
static const NekDouble kGeomFactorsTol
Definition: NektarUnivConsts.hpp:50
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().

1913  {
1914  int i,j;
1915  int npts = wavefield->GetPlane(0)->GetNpoints();
1916  Array<OneD, int> index(npts);
1917 
1918  Array<OneD, NekDouble> coord(2);
1919  Array<OneD, NekDouble> coord_x(npts);
1920  Array<OneD, NekDouble> coord_y(npts);
1921 
1922  //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
1923  wavefield->GetPlane(0)->GetCoords(coord_x,coord_y);
1924  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
1925  NekDouble tol = NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
1926  NekDouble xnew,ynew;
1927 
1928  int start = npts-1;
1929  for(i = 0; i < npts; ++i)
1930  {
1931  xnew = - coord_x[i] + xmax;
1932  ynew = - coord_y[i];
1933 
1934  for(j = start; j >=0 ; --j)
1935  {
1936  if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1937  {
1938  index[i] = j;
1939  start = j;
1940  break;
1941  }
1942  }
1943 
1944  if(j == -1)
1945  {
1946 
1947  for(j = npts-1; j > start; --j)
1948  {
1949 
1950  if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1951  {
1952  index[i] = j;
1953  break;
1954  }
1955  }
1956  ASSERTL0(j != start,"Failsed to find matching point");
1957  }
1958  }
1959  return index;
1960  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::NekConstants::kGeomFactorsTol
static const NekDouble kGeomFactorsTol
Definition: NektarUnivConsts.hpp:50
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::LibUtilities::ePolyEvenlySpaced, Extractlayerdata(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), GetReflectionIndex(), GetReflectionIndex2D(), Nektar::LibUtilities::Import(), Nektar::NekConstants::kNekUnsetDouble, Manipulate(), Nektar::SpatialDomains::MeshGraph::Read(), SetFields(), WriteBcs(), and WriteFld().

110 {
111  void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh,
112  SpatialDomains::BoundaryConditionsSharedPtr& boundaryConditions,
113  LibUtilities::SessionReaderSharedPtr &session,
114  Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,int nvariables);
115  Array<OneD, int> GetReflectionIndex(MultiRegions::ExpListSharedPtr Exp,
116  int Ireg);
117  Array<OneD, int> GetReflectionIndex2D(MultiRegions::ExpListSharedPtr wavefield);
118  void Extractlayerdata(Array<OneD, int> Iregions, int coordim,
119  SpatialDomains::MeshGraphSharedPtr &mesh,
120  LibUtilities::SessionReaderSharedPtr &session,
121  SpatialDomains::BoundaryConditions &bcs,
122  Array<OneD,MultiRegions::ExpListSharedPtr> &infields,
123  MultiRegions::ContField1DSharedPtr &outfieldx,
124  MultiRegions::ContField1DSharedPtr &outfieldy,
125  MultiRegions::ExpListSharedPtr &streak, bool symm,
126  Array<OneD, int> Refindices, NekDouble alpha, NekDouble cr);
127  void Manipulate(Array<OneD, int> Iregions, int coordim, SpatialDomains::MeshGraphSharedPtr &mesh,
128  SpatialDomains::BoundaryConditions &bcs,
129  Array<OneD,MultiRegions::ExpList1DSharedPtr> &infields,
130  Array<OneD,MultiRegions::ExpList1DSharedPtr> &outfieldx,
131  Array<OneD,MultiRegions::ExpList1DSharedPtr> &outfieldy,
132  MultiRegions::ExpListSharedPtr &streak);
133 
134  void CalcNonLinearForcing(SpatialDomains::MeshGraphSharedPtr &mesh,
135  LibUtilities::SessionReaderSharedPtr &session, string fieldfile,
136  Array<OneD,MultiRegions::ExpListSharedPtr> &waveFields,
137  MultiRegions::ExpListSharedPtr &streak,
138  Array<OneD, int> Refindices, bool symm );
139  void WriteBcs(string variable, int region, string fieldfile, SpatialDomains::MeshGraphSharedPtr &mesh,
140  MultiRegions::ContField1DSharedPtr &outregionfield);
141  void WriteFld(string outfile, SpatialDomains::MeshGraphSharedPtr &mesh, Array<OneD,
142  MultiRegions::ExpListSharedPtr> &fields );
143 
144 //REMARK:
145 /**
146 To create another bcs file comment everything from '//start' to '//end'
147 decomment the lines which follow '//decomment'
148  and run with the command: ./FldCalcBCs meshfile fieldfile fieldfile
149 **/
150 
151 
152 
153 
154  if(argc >= 6 || argc <4)
155  {
156  fprintf(stderr,"Usage: ./FldCalcBCs meshfile fieldfile streakfile (optional)alpha\n");
157  exit(1);
158  }
159 cout<<"argc="<<argc<<endl;
160  //change argc from 4 to 5 allow the loading of alpha to be optional
161  if(argc==4){ argc=5;}
162  //----------------------------------------------
163  string meshfile(argv[argc-4]);
164 
165  //setting argc=2 will take consider only the first two words argv[0],argv[1]
166  LibUtilities::SessionReaderSharedPtr vSession
167  = LibUtilities::SessionReader::CreateInstance(2, argv);
168 
169  // Read in mesh from input file
170  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);//meshfile);
171  //----------------------------------------------
172 
173  // Also read and store the boundary conditions
174  SpatialDomains::BoundaryConditionsSharedPtr boundaryConditions;
175  boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
176  ::AllocateSharedPtr(vSession, graphShPt);
177  SpatialDomains::BoundaryConditions bcs(vSession, graphShPt);
178  //----------------------------------------------
179 
180  //load alpha is provided
181  NekDouble alpha=0;
182  string alp_s;
183 
184  if(argv[argc-1])
185  {
186  alp_s = argv[argc-1];
187  alpha = boost::lexical_cast<double>(alp_s);
188  }
189  //----------------------------------------------
190 
191  //load the phase:
192  NekDouble cr=0.0;
193  if( vSession->DefinesSolverInfo("INTERFACE")
194  &&vSession->GetSolverInfo("INTERFACE")=="phase" )
195  {
196  vSession->LoadParameter("phase",cr,NekConstants::kNekUnsetDouble);
197  }
198 cout<<"cr="<<cr<<endl;
199  //-----------------------------------------------
200 
201 
202  //determine if the symmetrization is on:
203  bool symm =true;
204  vSession->MatchSolverInfo("Symmetrization","True",symm,true);
205  if( symm == true )
206  {
207  cout<<"symmetrization is active"<<endl;
208  }
209 
210 
211  // Import field file.
212  string fieldfile(argv[argc-3]);
213  vector<LibUtilities::FieldDefinitionsSharedPtr> fielddef;
214  vector<vector<NekDouble> > fielddata;
215  LibUtilities::Import(fieldfile,fielddef,fielddata);
216  //----------------------------------------------
217 
218  // Define Expansion
219  std::string solvtype = vSession->GetSolverInfo("SOLVERTYPE");
220  int nfields;
221  nfields = fielddef[0]->m_fields.size();
222  Array<OneD, MultiRegions::ExpListSharedPtr> fields;
223  if(solvtype == "CoupledLinearisedNS" && nfields==1)
224  {
225  nfields++;
226  }
227  fields= Array<OneD, MultiRegions::ExpListSharedPtr>(nfields);
228 
229  int lastfield = 0;
230  if(solvtype == "CoupledLinearisedNS" && nfields!=2)
231  {
232  SetFields(graphShPt,boundaryConditions,vSession,fields,nfields-1);
233  //decomment
234  //nfields = nfields-1;
235  //start
236  lastfield = nfields-1;
237  cout<<"Set pressure: "<<lastfield<<endl;
238  int nplanes = fielddef[0]->m_numModes[2];
239  const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::ePolyEvenlySpaced);
240  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
241  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
242  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
243  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,false,false,graphShPt,fielddef[0]->m_fields[0]);
244  fields[lastfield] = Exp3DH1;
245  //end
246  }
247  else if(solvtype == "CoupledLinearisedNS" && nfields==2)
248  {
249  cout<<"Set pressure split"<<endl;
250 
251  int nplanes = fielddef[0]->m_numModes[2];
252  const LibUtilities::PointsKey Pkey(nplanes,LibUtilities::ePolyEvenlySpaced);
253  const LibUtilities::BasisKey Bkey(fielddef[0]->m_basis[2],nplanes,Pkey);
254  NekDouble lz = fielddef[0]->m_homogeneousLengths[0];
255  cout<<"set ppp"<<endl;
256  /*
257  //to use GetBndCondExpansions() a contfield is needed and you j=have to
258  //call it "u" to get the bndconds
259  MultiRegions::ContField3DHomogeneous1DSharedPtr Cont3DH1;
260  Cont3DH1 = MemoryManager<MultiRegions::ContField3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,false,false,graphShPt, "u");
261  fields[0] = Cont3DH1;
262 */
263  fields[0] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
264  ::AllocateSharedPtr (vSession,Bkey,lz,false,false,graphShPt,vSession->GetVariable(0));
265 
266 
267 
268  //pressure is field 1
269  lastfield = nfields-1;
270  MultiRegions::ExpList3DHomogeneous1DSharedPtr Exp3DH1;
271  Exp3DH1 = MemoryManager<MultiRegions::ExpList3DHomogeneous1D>::AllocateSharedPtr(vSession,Bkey,lz,false,false,graphShPt,fielddef[0]->m_fields[0]);
272  fields[lastfield] = Exp3DH1;
273  }
274  else
275  {
276 
277  SetFields(graphShPt,boundaryConditions,vSession,fields,nfields);
278  }
279  //----------------------------------------------
280 
281  // Copy data from file:fill fields with the fielddata
282  if(lastfield==1)
283  {
284 
285  for(int i = 0; i < fielddata.size(); ++i)
286  {
287  fields[0]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[0]->m_fields[0], fields[0]->UpdateCoeffs());
288  }
289 
290  fields[0]->BwdTrans_IterPerExp(fields[0]->GetCoeffs(),fields[0]->UpdatePhys());
291 cout<<"field:"<<fielddef[0]->m_fields[0]<<endl;
292 /*
293  //bwd plane 0
294  fields[0]->GetPlane(0)->BwdTrans_IterPerExp(fields[0]->GetPlane(0)->GetCoeffs(),
295  fields[0]->GetPlane(0)->UpdatePhys() );
296 cout<<"hjhj"<<endl;
297  //bwd plane 1
298  fields[0]->GetPlane(1)->BwdTrans_IterPerExp(fields[0]->GetPlane(1)->GetCoeffs(),
299  fields[0]->GetPlane(1)->UpdatePhys() );
300 */
301 
302  for(int i = 0; i < fielddata.size(); ++i)
303  {
304  fields[lastfield]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[i]->m_fields[0], fields[lastfield]->UpdateCoeffs());
305  }
306  fields[lastfield]->BwdTrans_IterPerExp(fields[lastfield]->GetCoeffs(),fields[lastfield]->UpdatePhys());
307 /*
308  //bwd plane 0
309  fields[lastfield]->GetPlane(0)->BwdTrans_IterPerExp(fields[lastfield]->GetPlane(0)->GetCoeffs(),
310  fields[lastfield]->GetPlane(0)->UpdatePhys() );
311 
312  //bwd plane 1
313  fields[lastfield]->GetPlane(1)->BwdTrans_IterPerExp(fields[lastfield]->GetPlane(1)->GetCoeffs(),
314  fields[lastfield]->GetPlane(1)->UpdatePhys() );
315 */
316  }
317  else
318  {
319  for(int j = 0; j < nfields; ++j)
320  {
321  for(int i = 0; i < fielddata.size(); ++i)
322  {
323  fields[j]->ExtractDataToCoeffs(fielddef[i],fielddata[i],fielddef[i]->m_fields[j], fields[j]->UpdateCoeffs());
324  }
325  fields[j]->BwdTrans_IterPerExp(fields[j]->GetCoeffs(),fields[j]->UpdatePhys());
326  }
327  }
328 
329  //----------------------------------------------
330  //the phys values are the same (WRONG) but the coeffs are CORRECT!!!
331 /*
332  for(int g=0; g<fields[0]->GetPlane(1)->GetNcoeffs(); g++)
333  {
334 cout<<"g="<<g<<" coeff f0="<<fields[lastfield]->GetPlane(0)->GetCoeff(g)<<" f1="<<fields[lastfield]->GetPlane(1)->GetCoeff(g)<<endl;
335  }
336 */
337 
338  // import the streak
339  //import the streak field
340  //add y to obtain the streak w=w'+y and write fld file
341  MultiRegions::ExpListSharedPtr streak;
342 //start
343  streak = MemoryManager<MultiRegions::ContField2D>
344  ::AllocateSharedPtr(vSession, graphShPt, "w",true);
345  //streak = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr
346  // (vSession, graphShPt, true, "w");
347 
348  string streakfile(argv[argc-2]);
349  vector<LibUtilities::FieldDefinitionsSharedPtr> streakdef;
350  vector<vector<NekDouble> > streakdata;
351  LibUtilities::Import(streakfile,streakdef,streakdata);
352  //attention: streakdata.size==2 because the session file is 3DHomo1D but in
353  //reality the streak is real quantity
354 
355 
356  // Copy data from file:fill streak with the streakdata
357 
358  for(int i = 0; i < streakdata.size(); ++i)
359  {
360  streak->ExtractDataToCoeffs(streakdef[i],streakdata[i],streakdef[i]->m_fields[0], streak->UpdateCoeffs());
361  }
362  streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
363 
364  //end
365  //---------------------------------------------------------------
366  // determine the I regions
367  //hypothesis: the number of I regions is the same for all the variables
368  //hypothesis: all the I regions have the same nq points
369  int nIregions, lastIregion;
370  const Array<OneD, SpatialDomains::BoundaryConditionShPtr> bndConditions = streak->GetBndConditions();
371 
372  Array<OneD, int> Iregions =Array<OneD, int>(bndConditions.num_elements(),-1);
373  //3 internal layers:
374  Array<OneD, int> Ilayers =Array<OneD, int>(3,-1);
375  nIregions=0;
376  int nbnd= bndConditions.num_elements();
377  for(int r=0; r<nbnd; r++)
378  {
379  if(bndConditions[r]->GetUserDefined()==SpatialDomains::eCalcBC)
380  {
381  lastIregion=r;
382  Iregions[r]=r;
383  Ilayers[nIregions]=r;
384  nIregions++;
385 //cout<<"Iregion="<<r<<endl;
386  }
387  }
388  ASSERTL0(nIregions>0,"there is any boundary region with the tag USERDEFINEDTYPE=""CalcBC"" specified");
389  //-----------------------------------------------------------------
390 
391  //set 1D output fields:
392  //initialise fields
393  const SpatialDomains::BoundaryRegionCollection &bregions = bcs.GetBoundaryRegions();
394 
395  MultiRegions::ExpList1DSharedPtr outfieldxtmp;
396  MultiRegions::ExpList1DSharedPtr outfieldytmp;
397  MultiRegions::ContField1DSharedPtr outfieldx;
398  MultiRegions::ContField1DSharedPtr outfieldy;
399  //initialie fields
400  //const SpatialDomains::BoundaryRegionCollection &bregions = bcs.GetBoundaryRegions();
401  //declarecoeffsphysarray?????????!!!! true?!
402 
403  outfieldxtmp = MemoryManager<MultiRegions::ExpList1D>
404  ::AllocateSharedPtr(*(bregions.find(lastIregion)->second), graphShPt, true);
405  outfieldytmp = MemoryManager<MultiRegions::ExpList1D>
406  ::AllocateSharedPtr(*(bregions.find(lastIregion)->second), graphShPt, true);
407 
408  outfieldx = MemoryManager<MultiRegions::ContField1D>
409  ::AllocateSharedPtr(vSession, *outfieldxtmp);
410 
411  outfieldy = MemoryManager<MultiRegions::ContField1D>
412  ::AllocateSharedPtr(vSession, *outfieldytmp);
413 
414 //YOU NEED TO DEFINE THE OUTFIELDS AS ContField2D to use them in the 2D NS!!!!!!!!!!!
415 //(not 3DHomogeneous1D!!!!!!!
416 
417 /*
418  //define the outfields:
419  int nplanes = fielddef[0]->m_numModes[1];
420  int nzlines = fielddef[0]->m_numModes[2];
421  //ATTENTO Polyevenlyspaced NON Fourier!!!
422  const LibUtilities::PointsKey PkeyZ(nzlines,LibUtilities::ePolyEvenlySpaced);
423  const LibUtilities::BasisKey BkeyZ(fielddef[0]->m_basis[2],nzlines,PkeyZ);
424  NekDouble lz = fielddef[0]->m_homogeneousLengths[1];
425  for(int r=0; r<nlayers; r++)
426  {
427 
428  bndfieldx[r]= MemoryManager<MultiRegions::ExpList2DHomogeneous1D>
429  ::AllocateSharedPtr(vSession,BkeyZ,lz,false,graphShPt);
430 cout<<"OOOK"<<endl;
431  bndfieldy[r]= MemoryManager<MultiRegions::ExpList2DHomogeneous1D>
432  ::AllocateSharedPtr(vSession,BkeyZ,lz,false,graphShPt);
433  }
434  //--------------------------------------------------------
435 */
436 
437  bool coarseVWI=false;
438 
439  if(coarseVWI==false)
440  {
441  //manipulate data
442  //for 2 variables(u,v) only:
443  int coordim = graphShPt->GetMeshDimension();
444  //remark Ilayers[2] is the critical layer
445  static Array<OneD, int> Refindices = GetReflectionIndex(streak, Ilayers[2]);
446  Extractlayerdata(Ilayers,coordim, graphShPt,vSession, bcs,
447  fields, outfieldx,outfieldy,streak,symm,Refindices, alpha,cr);
448 
449  //--------------------------------------------------------------------------------------
450 
451 
452  //write bcs files: one for each I region and each variable
453  string var="u";
454  WriteBcs(var,lastIregion, fieldfile,graphShPt,outfieldx);
455  var="v";
456  WriteBcs(var,lastIregion, fieldfile,graphShPt,outfieldy);
457 
458  //--------------------------------------------------------------------------------
459  }
460  else
461  {
462  cout<<"CalcNonLinearForcing"<<endl;
463  static Array<OneD, int> Refindices;
464  if(symm==true)
465  {
466  Refindices = GetReflectionIndex2D(fields[0]);
467  }
468  CalcNonLinearForcing(graphShPt,vSession, fieldfile, fields, streak, Refindices,symm ) ;
469  }
470 
471 }
WriteFld
void WriteFld(string outfile, SpatialDomains::MeshGraphSharedPtr &mesh, Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Definition: FldCalcBCs.cpp:1990
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::MultiRegions::ContField1DSharedPtr
boost::shared_ptr< ContField1D > ContField1DSharedPtr
Definition: ContField1D.h:237
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
GetReflectionIndex2D
Array< OneD, int > GetReflectionIndex2D(MultiRegions::ExpListSharedPtr wavefield)
Definition: FldCalcBCs.cpp:1911
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::LibUtilities::Import
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > ElementiDs)
Imports an FLD file.
Definition: FieldIO.cpp:115
Nektar::LibUtilities::ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:63
Nektar::MultiRegions::ExpList1DSharedPtr
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
SetFields
void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh, SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions, LibUtilities::SessionReaderSharedPtr &session, Array< OneD, MultiRegions::ExpListSharedPtr > &Exp, int nvariables)
Definition: FldCalcBCs.cpp:474
Manipulate
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: FldCalcBCs.cpp:1671
Nektar::SpatialDomains::BoundaryRegionCollection
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:206
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
CalcNonLinearForcing
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: FldCalcBCs.cpp:1684
WriteBcs
void WriteBcs(string variable, int region, string fieldfile, SpatialDomains::MeshGraphSharedPtr &mesh, MultiRegions::ContField1DSharedPtr &outregionfield)
Definition: FldCalcBCs.cpp:1962
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
GetReflectionIndex
Array< OneD, int > GetReflectionIndex(MultiRegions::ExpListSharedPtr Exp, int Ireg)
Definition: FldCalcBCs.cpp:637
Nektar::NekConstants::kNekUnsetDouble
static const NekDouble kNekUnsetDouble
Definition: NektarUnivConsts.hpp:47
Extractlayerdata
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: FldCalcBCs.cpp:692
Nektar::SpatialDomains::BoundaryConditionsSharedPtr
boost::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:271
Nektar::MultiRegions::ExpList3DHomogeneous1DSharedPtr
boost::shared_ptr< ExpList3DHomogeneous1D > ExpList3DHomogeneous1DSharedPtr
Shared pointer to an ExpList3DHomogeneous1D object.
Definition: ExpList3DHomogeneous1D.h:50
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
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().

1678  {
1679 
1680  }
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().

478  {
479  // Setting parameteres for homogenous problems
480  NekDouble LhomX; ///< physical length in X direction (if homogeneous)
481  NekDouble LhomY; ///< physical length in Y direction (if homogeneous)
482  NekDouble LhomZ; ///< physical length in Z direction (if homogeneous)
483 
484  bool DeclareCoeffPhysArrays = true;
485  int npointsX; ///< number of points in X direction (if homogeneous)
486  int npointsY; ///< number of points in Y direction (if homogeneous)
487  int npointsZ; ///< number of points in Z direction (if homogeneous)
488  int HomoDirec = 0;
489  bool useFFT = false;
490  bool deal = false;
491  ///Parameter for homogeneous expansions
492  enum HomogeneousType
493  {
494  eHomogeneous1D,
495  eHomogeneous2D,
496  eHomogeneous3D,
497  eNotHomogeneous
498  };
499 
500  enum HomogeneousType HomogeneousType = eNotHomogeneous;
501 
502  if(session->DefinesSolverInfo("HOMOGENEOUS"))
503  {
504  std::string HomoStr = session->GetSolverInfo("HOMOGENEOUS");
505  //m_spacedim = 3;
506 
507  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
508  (HomoStr == "1D")||(HomoStr == "Homo1D"))
509  {
510  HomogeneousType = eHomogeneous1D;
511  npointsZ = session->GetParameter("HomModesZ");
512  LhomZ = session->GetParameter("LZ");
513  HomoDirec = 1;
514  }
515 
516  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
517  (HomoStr == "2D")||(HomoStr == "Homo2D"))
518  {
519  HomogeneousType = eHomogeneous2D;
520  npointsY = session->GetParameter("HomModesY");
521  LhomY = session->GetParameter("LY");
522  npointsZ = session->GetParameter("HomModesZ");
523  LhomZ = session->GetParameter("LZ");
524  HomoDirec = 2;
525  }
526 
527  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
528  (HomoStr == "3D")||(HomoStr == "Homo3D"))
529  {
530  HomogeneousType = eHomogeneous3D;
531  npointsX = session->GetParameter("HomModesX");
532  LhomX = session->GetParameter("LX");
533  npointsY = session->GetParameter("HomModesY");
534  LhomY = session->GetParameter("LY");
535  npointsZ = session->GetParameter("HomModesZ");
536  LhomZ = session->GetParameter("LZ");
537  HomoDirec = 3;
538  }
539 
540  if(session->DefinesSolverInfo("USEFFT"))
541  {
542  useFFT = true;
543  }
544  }
545 
546  int i;
547  int expdim = mesh->GetMeshDimension();
548  //Exp= Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
549  // I can always have 3 variables in a 2D mesh (oech vel component i a function which can depend on 1-3 var)
550  // Continuous Galerkin projection
551 
552  switch(expdim)
553  {
554  case 1:
555  {
556  if(HomogeneousType == eHomogeneous2D)
557  {
558  const LibUtilities::PointsKey PkeyY(npointsY,LibUtilities::eFourierEvenlySpaced);
559  const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,npointsY,PkeyY);
560  const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
561  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
562 
563  for(i = 0 ; i < nvariables; i++)
564  {
565  Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous2D>
566  ::AllocateSharedPtr(session,BkeyY,BkeyZ,LhomY,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
567  }
568  }
569  else
570  {
571  for(i = 0 ; i < nvariables; i++)
572  {
573  Exp[i] = MemoryManager<MultiRegions::ContField1D>
574  ::AllocateSharedPtr(session,mesh,session->GetVariable(i));
575  }
576  }
577 
578  break;
579  }
580  case 2:
581  {
582  if(HomogeneousType == eHomogeneous1D)
583  {
584  const LibUtilities::PointsKey PkeyZ(npointsZ,LibUtilities::eFourierEvenlySpaced);
585  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,npointsZ,PkeyZ);
586  for(i = 0 ; i < nvariables; i++)
587  {
588  Exp[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
589  ::AllocateSharedPtr(session,BkeyZ,LhomZ,useFFT,deal,mesh,session->GetVariable(i));
590  }
591  }
592  else
593  {
594  i = 0;
595  MultiRegions::ContField2DSharedPtr firstfield;
596  firstfield = MemoryManager<MultiRegions::ContField2D>
597  ::AllocateSharedPtr(session,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
598 
599  Exp[0] = firstfield;
600  for(i = 1 ; i < nvariables; i++)
601  {
602  Exp[i] = MemoryManager<MultiRegions::ContField2D>
603  ::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i),DeclareCoeffPhysArrays);
604  }
605  }
606  break;
607  }
608  case 3:
609  {
610  if(HomogeneousType == eHomogeneous3D)
611  {
612  ASSERTL0(false,"3D fully periodic problems not implemented yet");
613  }
614  else
615  {
616  i = 0;
617  MultiRegions::ContField3DSharedPtr firstfield =
618  MemoryManager<MultiRegions::ContField3D>
619  ::AllocateSharedPtr(session,mesh,session->GetVariable(i));
620 
621  Exp[0] = firstfield;
622  for(i = 1 ; i < nvariables; i++)
623  {
624  Exp[i] = MemoryManager<MultiRegions::ContField3D>
625  ::AllocateSharedPtr(*firstfield,mesh,session->GetVariable(i));
626  }
627  }
628  break;
629  }
630  default:
631  ASSERTL0(false,"Expansion dimension not recognised");
632  break;
633  }
634  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::MultiRegions::ContField2DSharedPtr
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:293
Nektar::LibUtilities::eFourier
Fourier Expansion .
Definition: BasisType.h:52
Nektar::LibUtilities::eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:64
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::MultiRegions::ContField3DSharedPtr
boost::shared_ptr< ContField3D > ContField3DSharedPtr
Definition: ContField3D.h:191
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
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().

1964  {
1965  string outfile = fieldfile.substr(0,fieldfile.find_last_of("."));
1966  outfile +="_"+variable+"_";
1967  char ibnd[16]="";
1968  sprintf(ibnd,"%d",region);
1969  outfile +=ibnd;
1970  string endfile(".bc");
1971  outfile += endfile;
1972  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);
1973  outregionfield->FwdTrans_IterPerExp(outregionfield->GetPhys(),outregionfield->UpdateCoeffs());
1974  fieldcoeffs[0] = outregionfield->UpdateCoeffs();
1975  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
1976  = outregionfield->GetFieldDefinitions();
1977  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
1978  // copy Data into FieldData and set variable
1979  //for(int i = 0; i < FieldDef.size(); ++i)
1980  //{
1981  //FieldDef[i]->m_fields.push_back(variable);
1982  FieldDef[0]->m_fields.push_back(variable);
1983  //outregionfield->AppendFieldData(FieldDef[i], FieldData[i]);
1984  outregionfield->AppendFieldData(FieldDef[0], FieldData[0]);
1985  //outregionfield->AppendFieldData(FieldDef[i], FieldData[i],fieldcoeffs[0]);
1986  //}
1987  LibUtilities::Write(outfile,FieldDef,FieldData);
1988  }
Nektar::LibUtilities::Write
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap)
Write a field file in serial only.
Definition: FieldIO.cpp:81
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().

1991  {
1992  //variables array dimension is defined by fields!!!
1993  Array<OneD, std::string> variables(fields.num_elements());
1994  if( variables.num_elements()==2)
1995  {
1996  variables[0]="u";
1997  variables[1]="v";
1998  }
1999  else if( variables.num_elements()==3)
2000  {
2001  variables[0]="u";
2002  variables[1]="v";
2003  variables[2]="w";
2004  }
2005  else
2006  {
2007  ASSERTL0(false, " something goes wrong... ");
2008  }
2009  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(fields.num_elements());
2010 
2011  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
2012  = fields[0]->GetFieldDefinitions();
2013 
2014  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
2015  //this cycle is needed in order to use the the appendfielddata with coeffs function!!!
2016  for(int i = 0; i < fields.num_elements(); ++i)
2017  {
2018  if (fields[i]->GetPhysState()==true)
2019  {
2020  //fields[i]->FwdTrans_IterPerExp(fields[i]->GetPhys(),fields[i]->UpdateCoeffs());
2021  fields[i]->FwdTrans(fields[i]->GetPhys(),fields[i]->UpdateCoeffs());
2022 
2023  }
2024  fieldcoeffs[i] = fields[i]->UpdateCoeffs();
2025  }
2026 
2027  // copy Data into FieldData and set variable
2028  for(int j = 0; j < fieldcoeffs.num_elements(); ++j)
2029  {
2030  //fieldcoeffs[j] = fields[j]->UpdateCoeffs();
2031  for(int i = 0; i < FieldDef.size(); ++i)
2032  {
2033  // Could do a search here to find correct variable
2034  FieldDef[i]->m_fields.push_back(variables[j]);
2035  //remark: appendfielddata without fieldcoeffs input gives always the first field!!!
2036  //fields[0]->AppendFieldData(FieldDef[i], FieldData[i]);
2037  fields[0]->AppendFieldData(FieldDef[i], FieldData[i],fieldcoeffs[j]);
2038  }
2039  }
2040  LibUtilities::Write(outfile,FieldDef,FieldData);
2041 
2042  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::LibUtilities::Write
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap)
Write a field file in serial only.
Definition: FieldIO.cpp:81
Generated on Sat Jun 4 2016 00:02:28 for Nektar++ by   doxygen 1.8.8