Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
ExtractCriticalLayerFunctions.h File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList2D.h>
Include dependency graph for ExtractCriticalLayerFunctions.h:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

void Computestreakpositions (MultiRegions::ExpListSharedPtr &streak, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, NekDouble cr, NekDouble trans=NekConstants::kNekUnsetDouble)
 

Function Documentation

void Computestreakpositions ( MultiRegions::ExpListSharedPtr &  streak,
Array< OneD, NekDouble > &  xc,
Array< OneD, NekDouble > &  yc,
NekDouble  cr,
NekDouble  trans = NekConstants::kNekUnsetDouble 
)

Definition at line 119 of file ExtractCriticalLayer.cpp.

References ASSERTL0, Nektar::MultiRegions::eY, npts, Vmath::Vmax(), and Vmath::Vmin().

Referenced by GetStreakLocation(), and main().

124 {
125  int i;
126  int npts = xc.num_elements();
127 
128  int nq = streak->GetTotPoints();
129  Array<OneD, NekDouble> derstreak(nq);
130  Array<OneD, NekDouble> x(nq);
131  Array<OneD, NekDouble> y(nq);
132  streak->GetCoords(x,y);
133 
134  streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
135  streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);
136 
137  // set intiial xc to be equispaced over mesh and yc to be zero
138  NekDouble x_max = Vmath::Vmax(nq,x,1);
139  NekDouble x_min = Vmath::Vmin(nq,x,1);
140 
141  for(i = 0; i < npts; ++i)
142  {
143  xc[i] = x_min + (x_max - x_min)*i/((NekDouble)(npts-1));
144  yc[i] = 0.0;
145  }
146 
147 
148  int elmtid, offset,cnt;
149  NekDouble U,dU;
150  NekDouble F;
151  NekDouble ConvTol = 1e-9;
152  NekDouble CoordTol = 1e-5;
153  int maxiter = 100;
154  Array<OneD, NekDouble> coord(2);
155  NekDouble ytmp;
156  // Do Newton iteration on y direction
157  cerr << "[";
158  for(int e=0; e<npts; e++)
159  {
160  coord[0] = xc[e];
161  coord[1] = yc[e];
162 cout<<e<<endl;
163  if(!(e%10))
164  {
165  cerr << ".";
166  }
167 
168  F = 1000;
169  cnt = 0;
170  int its=0;
171  int attempt=0;
172 //cout<<"start guess: x="<<xc[e]<<" y="<<yc[e]<<endl;
173  while( abs(F)> 0.000000001)
174  {
175  ytmp = coord[1];
176  elmtid = streak->GetExpIndex(coord,0.00001);
177  offset = streak->GetPhys_Offset(elmtid);
178  U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
179  dU = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
180  coord[1] = coord[1] - (U-cr)/dU;
181  F = U-cr;
182  ASSERTL0( coord[0]==xc[e], " x coordinate must remain the same");
183  //stvalues[e] = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() +offset );
184 //cout<<"elmtid="<<elmtid<<" x="<<coord[0]<<" y="<<coord[1]<<" stvalue="<<U<<" dU="<<dU<<endl;
185  if(abs(coord[1])>1 )
186  {
187  coord[1] = ytmp +0.01;
188  elmtid = streak->GetExpIndex(coord,0.00001);
189  offset = streak->GetPhys_Offset(elmtid);
190  NekDouble Utmp = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
191  NekDouble dUtmp = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
192  coord[1] = coord[1] - (Utmp-cr)/dUtmp;
193 
194  if( (abs(Utmp-cr)>abs(F))||(abs(coord[1])>1) )
195  {
196  coord[1] = ytmp -0.01;
197  }
198 
199  attempt++;
200  }
201  else
202  {
203  ASSERTL0(abs(coord[1])<= 1, " y value out of bound +/-1");
204  }
205 
206  its++;
207  if(its>1000 && abs(F)< 0.0001)
208  {
209  cout<<"warning streak position obtained with precision:"<<F<<endl;
210  break;
211  }
212  else if(its>1000)
213  {
214  ASSERTL0(false, "no convergence after 1000 iterations");
215  }
216  }
217 
218  yc[e] = coord[1];
219  }
220  cerr << "]" << endl;
221 
222  // output to interface file
223  FILE *fp = fopen("interfacedat.geo","w");
224 
225  NekDouble y_max = Vmath::Vmax(nq,y,1);
226  NekDouble y_min = Vmath::Vmin(nq,y,1);
227 
228  cnt = 1;
229  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
230  cnt++,x_min,y_min);
231  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
232  cnt++,x_max,y_min);
233  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
234  cnt++,x_max,y_max);
235  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
236  cnt++,x_min,y_max);
237 
238  for(i = 0; i < npts; ++i)
239  {
240  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
241  cnt++,xc[i],yc[i]);
242  }
243 
244  fclose(fp);
245 
246  // output to interface_up file as bend of vertical shift and 45 degrees shift
247  fp = fopen("interfacedat_up.geo","w");
248 
249 
250  NekDouble nx,ny,norm;
251 
252  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
253 
254  for(i = 1; i < npts-1; ++i)
255  {
256  norm = sqrt((xc[i+1]-xc[i-1])*(xc[i+1]-xc[i-1])+(yc[i+1]-yc[i-1])*(yc[i+1]-yc[i-1]));
257  nx = (yc[i-1]-yc[i+1])/norm;
258  ny = (xc[i+1]-xc[i-1])/norm;
259 
260  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
261  cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
262  }
263 
264  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
265 
266 
267  // output to interface_up file as bend of vertical shift and 45 degrees shift
268  fp = fopen("interfacedat_dn.geo","w");
269 
270  trans = -trans;
271 
272  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
273 
274  for(i = 1; i < npts-1; ++i)
275  {
276  norm = sqrt((xc[i+1]-xc[i-1])*(xc[i+1]-xc[i-1])+(yc[i+1]-yc[i-1])*(yc[i+1]-yc[i-1]));
277  nx = (yc[i-1]-yc[i+1])/norm;
278  ny = (xc[i+1]-xc[i-1])/norm;
279 
280  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
281  cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
282  }
283 
284  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
285 
286 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:857
static std::string npts
Definition: InputFld.cpp:43
double NekDouble