126 int npts = xc.num_elements();
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);
134 streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
141 for(i = 0; i <
npts; ++i)
143 xc[i] = x_min + (x_max - x_min)*i/((
NekDouble)(npts-1));
148 int elmtid, offset,cnt;
154 Array<OneD, NekDouble> coord(2);
158 for(
int e=0; e<
npts; e++)
173 while( abs(F)> 0.000000001)
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;
182 ASSERTL0( coord[0]==xc[e],
" x coordinate must remain the same");
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;
194 if( (abs(Utmp-cr)>abs(F))||(abs(coord[1])>1) )
196 coord[1] = ytmp -0.01;
203 ASSERTL0(abs(coord[1])<= 1,
" y value out of bound +/-1");
207 if(its>1000 && abs(F)< 0.0001)
209 cout<<
"warning streak position obtained with precision:"<<F<<endl;
214 ASSERTL0(
false,
"no convergence after 1000 iterations");
223 FILE *fp = fopen(
"interfacedat.geo",
"w");
229 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
231 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
233 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
235 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
238 for(i = 0; i <
npts; ++i)
240 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
247 fp = fopen(
"interfacedat_up.geo",
"w");
252 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
254 for(i = 1; i < npts-1; ++i)
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;
260 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
261 cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
264 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
268 fp = fopen(
"interfacedat_dn.geo",
"w");
272 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
274 for(i = 1; i < npts-1; ++i)
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;
280 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
281 cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
284 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
#define ASSERTL0(condition, msg)
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.