Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 43 of file ExtractCriticalLayerFunctions.cpp.

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

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