44 Array<OneD, NekDouble> &xc,
45 Array<OneD, NekDouble> &yc,
49 int main(
int argc,
char *argv[])
56 fprintf(stderr,
"Usage: ./ExtractCriticalLayer meshfile fieldfile \n");
63 = LibUtilities::SessionReader::CreateInstance(argc, argv);
75 streak = MemoryManager<MultiRegions::ExpList2D>
76 ::AllocateSharedPtr(vSession,graphShPt);
82 string fieldfile(argv[argc-1]);
83 vector<SpatialDomains::FieldDefinitionsSharedPtr> fielddef;
84 vector<vector<NekDouble> > fielddata;
85 graphShPt->Import(fieldfile,fielddef,fielddata);
90 string streak_field(
"w");
91 for(
unsigned int i = 0; i < fielddata.size(); ++i)
93 streak->ExtractDataToCoeffs(fielddef [i],
100 vSession->LoadParameter(
"NumCriticalLayerPts",npts,30);
101 Array<OneD, NekDouble> x_c(npts);
102 Array<OneD, NekDouble> y_c(npts);
106 vSession->LoadParameter(
"WidthOfLayers",trans,0.1);
110 cout <<
"# x_c y_c" << endl;
111 for(i = 0; i <
npts; ++i)
113 fprintf(stdout,
"%12.10lf %12.10lf \n",x_c[i],y_c[i]);
120 Array<OneD, NekDouble> &xc,
121 Array<OneD, NekDouble> &yc,
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.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr