50 int npts = xc.num_elements();
52 int nq = streak->GetTotPoints();
56 streak->GetCoords(x,y);
58 streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
65 for(i = 0; i <
npts; ++i)
67 xc[i] = x_min + (x_max - x_min)*i/((
NekDouble)(npts-1));
72 int elmtid, offset,cnt;
82 for(
int e=0; e<
npts; e++)
94 while((abs(F)> ConvTol)&&(cnt < maxiter))
96 elmtid = streak->GetExpIndex(coord,CoordTol);
97 offset = streak->GetPhys_Offset(elmtid);
99 U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
100 dU = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
102 coord[1] = coord[1] - (U-cr)/dU;
107 ASSERTL0(cnt < maxiter,
"Failed to converge Newton iteration");
116 FILE *fp = fopen(
"interfacedat.geo",
"w");
122 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
124 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
126 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
128 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
131 for(i = 0; i <
npts; ++i)
133 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
141 fp = fopen(
"interfacedat_up.geo",
"w");
146 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
148 for(i = 1; i < npts-1; ++i)
150 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]));
151 nx = (yc[i-1]-yc[i+1])/norm;
152 ny = (xc[i+1]-xc[i-1])/norm;
154 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
155 cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
158 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
162 fp = fopen(
"interfacedat_dn.geo",
"w");
166 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
168 for(i = 1; i < npts-1; ++i)
170 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]));
171 nx = (yc[i-1]-yc[i+1])/norm;
172 ny = (xc[i+1]-xc[i-1])/norm;
174 fprintf(fp,
"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
175 cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
178 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< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
static const NekDouble kNekUnsetDouble