Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FldCalcBCs.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 #include <cmath>
4 #include <iomanip>
5 #include <LibUtilities/BasicUtils/SessionReader.h>
6 #include <MultiRegions/ExpList.h>
7 #include <MultiRegions/ExpList1D.h>
8 #include <MultiRegions/ContField1D.h>
9 #include <MultiRegions/ContField2D.h>
10 #include <MultiRegions/ContField3D.h>
11 #include <MultiRegions/ContField3DHomogeneous1D.h>
12 #include <MultiRegions/ContField3DHomogeneous2D.h>
13 #include <LocalRegions/SegExp.h>
14 #include <LocalRegions/Expansion2D.h>
15 #include <SpatialDomains/MeshGraph2D.h>
16 
17 
18 using namespace Nektar;
19 
20 /**
21  * Below is the old function from GeomFactors1D::v_ComputeEdgeTangents. This
22  * However, the vectors were computed in the function below, but is never called
23  * and therefore was never working.
24  *
25  * This utility therefore does not presently work, but the (now removed)
26  * function is kept here for reference.
27  *
28  * -- Chris Cantwell
29  *
30  void GeomFactors1D::v_ComputeEdgeTangents(
31  const GeometrySharedPtr &geom,
32  const int edge,
33  const LibUtilities::PointsKey &to_key)
34  {
35  int k;
36  int nquad= to_key.GetNumPoints();
37  Geometry2DSharedPtr g;
38  ASSERTL0(g= boost::dynamic_pointer_cast<Geometry2D>(geom),
39  "FAIL");
40 
41  GeomFactorsSharedPtr gf = geom->GetMetricInfo();
42  //cannot use m_type here
43  //GeomType gtype = gf->GetGtype();
44  GeomType gtype= m_type;
45  m_tangent =Array<OneD, Array<OneD, NekDouble> >(m_coordDim);
46  for( k=0; k< m_coordDim; ++k)
47  {
48  m_tangent[k] = Array<OneD, NekDouble>(nquad);
49  }
50 
51  int i;
52 
53  DerivStorage deriv = GetDeriv(m_pointsKey);
54  //FillDeriv(deriv, m_pointsKey);
55 
56  NekDouble fac;
57  // Regular geometry case
58  if((gtype == eRegular)||(gtype == eMovingRegular))
59  {
60 
61  for(i = 0; i < m_coordDim; ++i)
62  {
63  Vmath::Fill(nquad, deriv[0][i][0],m_tangent[i],1);
64  }
65 
66  // normalise
67  fac = 0.0;
68  for(i =0 ; i < m_coordDim; ++i)
69  {
70  fac += m_tangent[i][0]*m_tangent[i][0];
71  }
72  fac = 1.0/sqrt(fac);
73  for (i = 0; i < m_coordDim; ++i)
74  {
75  Vmath::Smul(nquad,fac,m_tangent[i],1,m_tangent[i],1);
76  }
77  }
78 
79  else // Set up deformed tangents
80  {
81 
82  Array<OneD, NekDouble> jac(m_coordDim*nquad);
83 
84  for(i = 0; i < m_coordDim; ++i)
85  {
86  for(int j=0; j<nquad; j++)
87  {
88  m_tangent[i][j] = deriv[0][i][j];
89  }
90  }
91  //normalise normal vectors
92  Array<OneD,NekDouble> work(nquad,0.0);
93  for(i = 0; i < m_coordDim; ++i)
94  {
95  Vmath::Vvtvp(nquad, m_tangent[i],1, m_tangent[i],1,work,1,work,1);
96  }
97 
98  Vmath::Vsqrt(nquad,work,1,work,1);
99  Vmath::Sdiv(nquad,1.0,work,1,work,1);
100 
101  for(i = 0; i < m_coordDim; ++i)
102  {
103  Vmath::Vmul(nquad, m_tangent[i],1,work,1,m_tangent[i],1);
104  }
105  }
106  }
107  *
108  */
109 int main(int argc, char *argv[])
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 }
472 
473  // Define Expansion
474  void SetFields(SpatialDomains::MeshGraphSharedPtr &mesh,
475  SpatialDomains::BoundaryConditionsSharedPtr &boundaryConditions,
476  LibUtilities::SessionReaderSharedPtr &session,
477  Array<OneD,MultiRegions::ExpListSharedPtr> &Exp,int nvariables)
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  }
635 
636 
637  Array<OneD, int> GetReflectionIndex( MultiRegions::ExpListSharedPtr Exp, int Ireg)
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  }
691 
692  void Extractlayerdata(Array<OneD, int> Iregions, int coordim,
693  SpatialDomains::MeshGraphSharedPtr &mesh,
694  LibUtilities::SessionReaderSharedPtr &session,
695  SpatialDomains::BoundaryConditions &bcs,
696  Array<OneD,MultiRegions::ExpListSharedPtr> &infields,
697  MultiRegions::ContField1DSharedPtr &outfieldx,
698  MultiRegions::ContField1DSharedPtr &outfieldy,
699  MultiRegions::ExpListSharedPtr &streak, bool symm,
700  Array<OneD, int> Refindices, NekDouble alpha,NekDouble cr)
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  }
1669 
1670 
1671  void Manipulate(Array<OneD, int> Iregions,int coordim, SpatialDomains::MeshGraphSharedPtr &mesh,
1672  SpatialDomains::BoundaryConditions &bcs,
1673  Array<OneD,MultiRegions::ExpListSharedPtr> &infields,
1674  Array<OneD,MultiRegions::ExpList1DSharedPtr> &outfieldx,
1675  Array<OneD,MultiRegions::ExpList1DSharedPtr> &outfieldy,
1676  MultiRegions::ExpListSharedPtr &streak)
1677 
1678  {
1679 
1680  }
1681 
1682 
1683 
1684  void CalcNonLinearForcing(SpatialDomains::MeshGraphSharedPtr &mesh,
1685  LibUtilities::SessionReaderSharedPtr &session, string fieldfile,
1686  Array<OneD,MultiRegions::ExpListSharedPtr> &waveFields,
1687  MultiRegions::ExpListSharedPtr &streak,
1688  Array<OneD, int> Refindices, bool symm )
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  }
1909 
1910 
1911  Array<OneD, int> GetReflectionIndex2D(
1912  MultiRegions::ExpListSharedPtr wavefield)
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  }
1961 
1962  void WriteBcs(string variable, int region, string fieldfile,SpatialDomains::MeshGraphSharedPtr &mesh,
1963  MultiRegions::ContField1DSharedPtr &outregionfield)
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  }
1989 
1990  void WriteFld(string outfile,SpatialDomains::MeshGraphSharedPtr &mesh, Array<OneD, MultiRegions::ExpListSharedPtr> &fields)
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  }
2043 
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2