Nektar++
Functions
ExtractCriticalLayer.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList2D.h>
Include dependency graph for ExtractCriticalLayer.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)
 
int main (int argc, char *argv[])
 

Function Documentation

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

Definition at line 119 of file ExtractCriticalLayer.cpp.

Referenced by GetStreakLocation(), and main().

124 {
125  int i;
126  int npts = xc.num_elements();
127 
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);
133 
134  streak->BwdTrans(streak->GetCoeffs(),streak->UpdatePhys());
135  streak->PhysDeriv(MultiRegions::eY, streak->GetPhys(), derstreak);
136 
137  // set intiial xc to be equispaced over mesh and yc to be zero
138  NekDouble x_max = Vmath::Vmax(nq,x,1);
139  NekDouble x_min = Vmath::Vmin(nq,x,1);
140 
141  for(i = 0; i < npts; ++i)
142  {
143  xc[i] = x_min + (x_max - x_min)*i/((NekDouble)(npts-1));
144  yc[i] = 0.0;
145  }
146 
147 
148  int elmtid, offset,cnt;
149  NekDouble U,dU;
150  NekDouble F;
151  NekDouble ConvTol = 1e-9;
152  NekDouble CoordTol = 1e-5;
153  int maxiter = 100;
154  Array<OneD, NekDouble> coord(2);
155  NekDouble ytmp;
156  // Do Newton iteration on y direction
157  cerr << "[";
158  for(int e=0; e<npts; e++)
159  {
160  coord[0] = xc[e];
161  coord[1] = yc[e];
162 cout<<e<<endl;
163  if(!(e%10))
164  {
165  cerr << ".";
166  }
167 
168  F = 1000;
169  cnt = 0;
170  int its=0;
171  int attempt=0;
172 //cout<<"start guess: x="<<xc[e]<<" y="<<yc[e]<<endl;
173  while( abs(F)> 0.000000001)
174  {
175  ytmp = coord[1];
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;
181  F = U-cr;
182  ASSERTL0( coord[0]==xc[e], " x coordinate must remain the same");
183  //stvalues[e] = streak->GetExp(elmtid)->PhysEvaluate(coord, streak->GetPhys() +offset );
184 //cout<<"elmtid="<<elmtid<<" x="<<coord[0]<<" y="<<coord[1]<<" stvalue="<<U<<" dU="<<dU<<endl;
185  if(abs(coord[1])>1 )
186  {
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;
193 
194  if( (abs(Utmp-cr)>abs(F))||(abs(coord[1])>1) )
195  {
196  coord[1] = ytmp -0.01;
197  }
198 
199  attempt++;
200  }
201  else
202  {
203  ASSERTL0(abs(coord[1])<= 1, " y value out of bound +/-1");
204  }
205 
206  its++;
207  if(its>1000 && abs(F)< 0.0001)
208  {
209  cout<<"warning streak position obtained with precision:"<<F<<endl;
210  break;
211  }
212  else if(its>1000)
213  {
214  ASSERTL0(false, "no convergence after 1000 iterations");
215  }
216  }
217 
218  yc[e] = coord[1];
219  }
220  cerr << "]" << endl;
221 
222  // output to interface file
223  FILE *fp = fopen("interfacedat.geo","w");
224 
225  NekDouble y_max = Vmath::Vmax(nq,y,1);
226  NekDouble y_min = Vmath::Vmin(nq,y,1);
227 
228  cnt = 1;
229  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
230  cnt++,x_min,y_min);
231  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
232  cnt++,x_max,y_min);
233  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
234  cnt++,x_max,y_max);
235  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
236  cnt++,x_min,y_max);
237 
238  for(i = 0; i < npts; ++i)
239  {
240  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
241  cnt++,xc[i],yc[i]);
242  }
243 
244  fclose(fp);
245 
246  // output to interface_up file as bend of vertical shift and 45 degrees shift
247  fp = fopen("interfacedat_up.geo","w");
248 
249 
250  NekDouble nx,ny,norm;
251 
252  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
253 
254  for(i = 1; i < npts-1; ++i)
255  {
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;
259 
260  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
261  cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
262  }
263 
264  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
265 
266 
267  // output to interface_up file as bend of vertical shift and 45 degrees shift
268  fp = fopen("interfacedat_dn.geo","w");
269 
270  trans = -trans;
271 
272  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[0],yc[0]+trans);
273 
274  for(i = 1; i < npts-1; ++i)
275  {
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;
279 
280  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",
281  cnt++,xc[i]+nx*trans,yc[i]+ny*trans);
282  }
283 
284  fprintf(fp,"Point(%d)={%12.10lf,%12.10lf,0,1.0}; \n",cnt++,xc[npts-1],yc[npts-1]+trans);
285 
286 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Vmath::Vmax
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
Vmath::Vmin
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
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
int main ( int  argc,
char *  argv[] 
)

Definition at line 49 of file ExtractCriticalLayer.cpp.

References Computestreakpositions(), and npts.

50 {
51  int i,j;
52  NekDouble cr = 0;
53 
54  if(argc !=3)
55  {
56  fprintf(stderr,"Usage: ./ExtractCriticalLayer meshfile fieldfile \n");
57  exit(1);
58  }
59 
60  //------------------------------------------------------------
61  // Create Session file.
62  LibUtilities::SessionReaderSharedPtr vSession
63  = LibUtilities::SessionReader::CreateInstance(argc, argv);
64  //-----------------------------------------------------------
65 
66  //-------------------------------------------------------------
67  // Read in mesh from input file
68  SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(vSession);
69  //------------------------------------------------------------
70 
71  //-------------------------------------------------------------
72  // Define Streak Expansion
73  MultiRegions::ExpListSharedPtr streak;
74 
75  streak = MemoryManager<MultiRegions::ExpList2D>
76  ::AllocateSharedPtr(vSession,graphShPt);
77 
78  //---------------------------------------------------------------
79 
80  //----------------------------------------------
81  // Import field file.
82  string fieldfile(argv[argc-1]);
83  vector<SpatialDomains::FieldDefinitionsSharedPtr> fielddef;
84  vector<vector<NekDouble> > fielddata;
85  graphShPt->Import(fieldfile,fielddef,fielddata);
86  //----------------------------------------------
87 
88  //----------------------------------------------
89  // Copy data from field file
90  string streak_field("w");
91  for(unsigned int i = 0; i < fielddata.size(); ++i)
92  {
93  streak->ExtractDataToCoeffs(fielddef [i],
94  fielddata[i],
95  streak_field);
96  }
97  //----------------------------------------------
98 
99  int npts;
100  vSession->LoadParameter("NumCriticalLayerPts",npts,30);
101  Array<OneD, NekDouble> x_c(npts);
102  Array<OneD, NekDouble> y_c(npts);
103  NekDouble xc, yc;
104 
105  NekDouble trans;
106  vSession->LoadParameter("WidthOfLayers",trans,0.1);
107 
108  Computestreakpositions(streak,x_c, y_c,cr,trans);
109 
110  cout << "# x_c y_c" << endl;
111  for(i = 0; i < npts; ++i)
112  {
113  fprintf(stdout,"%12.10lf %12.10lf \n",x_c[i],y_c[i]);
114  //cout << x_c[i] << " " << y_c[i] << endl;
115  }
116 
117 }
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
Computestreakpositions
void Computestreakpositions(MultiRegions::ExpListSharedPtr &streak, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, NekDouble cr, NekDouble trans)
Definition: ExtractCriticalLayer.cpp:119
Nektar::MultiRegions::ExpListSharedPtr
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
Generated on Tue Dec 8 2015 23:45:31 for Nektar++ by   doxygen 1.8.9.1