Nektar++
Functions
ExtractCriticalLayerFunctions.cpp File Reference
#include <cstdio>
#include <MultiRegions/ExpList.h>
Include dependency graph for ExtractCriticalLayerFunctions.cpp:

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)
 

Function Documentation

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

Definition at line 41 of file ExtractCriticalLayerFunctions.cpp.

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

46 {
47  int i;
48  int npts = xc.num_elements();
49 
50  int nq = streak->GetTotPoints();
51  Array<OneD, NekDouble> derstreak(nq);
54  streak->GetCoords(x,y);
55 
56  streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
57  streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);
58 
59  // set intiial xc to be equispaced over mesh and yc to be zero
60  NekDouble x_max = Vmath::Vmax(nq,x,1);
61  NekDouble x_min = Vmath::Vmin(nq,x,1);
62 
63  for(i = 0; i < npts; ++i)
64  {
65  xc[i] = x_min + (x_max - x_min)*i/((NekDouble)(npts-1));
66  yc[i] = 0.0;
67  }
68 
69 
70  int elmtid, offset,cnt;
71  NekDouble U,dU;
72  NekDouble F;
73  NekDouble ConvTol = 1e-9;
74  NekDouble CoordTol = 1e-5;
75  int maxiter = 100;
76  Array<OneD, NekDouble> coord(2);
77 
78  // Do Newton iteration on y direction
79  cerr << "[";
80  for(int e=0; e<npts; e++)
81  {
82  coord[0] = xc[e];
83  coord[1] = yc[e];
84 
85  if(!(e%10))
86  {
87  cerr << ".";
88  }
89 
90  F = 1000;
91  cnt = 0;
92  while((abs(F)> ConvTol)&&(cnt < maxiter))
93  {
94  elmtid = streak->GetExpIndex(coord,CoordTol);
95  offset = streak->GetPhys_Offset(elmtid);
96 
97  U = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() + offset);
98  dU = streak->GetExp(elmtid)->PhysEvaluate(coord, derstreak + offset);
99 
100  coord[1] = coord[1] - (U-cr)/dU;
101 
102  F = U-cr;
103  cnt++;
104  }
105  ASSERTL0(cnt < maxiter, "Failed to converge Newton iteration");
106 
107  yc[e] = coord[1];
108  }
109  cerr << "]" << endl;
110 
111  if(trans != NekConstants::kNekUnsetDouble)
112  {
113  // output to interface file
114  FILE *fp = fopen("interfacedat.geo","w");
115 
116  NekDouble y_max = Vmath::Vmax(nq,y,1);
117  NekDouble y_min = Vmath::Vmin(nq,y,1);
118 
119  cnt = 1;
120  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
121  cnt++,x_min,y_min);
122  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
123  cnt++,x_max,y_min);
124  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
125  cnt++,x_max,y_max);
126  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
127  cnt++,x_min,y_max);
128 
129  for(i = 0; i < npts; ++i)
130  {
131  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
132  cnt++,xc[i],yc[i]);
133  }
134 
135  fclose(fp);
136 
137 
138  // output to interface_up file as bend of vertical shift and 45 degrees shift
139  fp = fopen("interfacedat_up.geo","w");
140 
141 
142  NekDouble nx,ny,norm;
143 
144  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
145 
146  for(i = 1; i < npts-1; ++i)
147  {
148  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]));
149  nx = (yc[i-1]-yc[i+1])/norm;
150  ny = (xc[i+1]-xc[i-1])/norm;
151 
152  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
153  cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
154  }
155 
156  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
157 
158 
159  // output to interface_up file as bend of vertical shift and 45 degrees shift
160  fp = fopen("interfacedat_dn.geo","w");
161 
162  trans = -trans;
163 
164  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
165 
166  for(i = 1; i < npts-1; ++i)
167  {
168  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]));
169  nx = (yc[i-1]-yc[i+1])/norm;
170  ny = (xc[i+1]-xc[i-1])/norm;
171 
172  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
173  cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
174  }
175 
176  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
177  }
178 
179 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
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:756
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:848
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
static const NekDouble kNekUnsetDouble