Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExtractCriticalLayer.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExtractCriticalLayer.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Extract location of critical layer from streak file
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <cstdio>
37 #include <cstdlib>
38 
39 #include <MultiRegions/ExpList.h>
40 #include <MultiRegions/ExpList2D.h>
41 
42 
43 void Computestreakpositions(MultiRegions::ExpListSharedPtr &streak,
44  Array<OneD, NekDouble> &xc,
45  Array<OneD, NekDouble> &yc,
46  NekDouble cr,
47  NekDouble trans);
48 
49 int main(int argc, char *argv[])
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 }
118 
119 void Computestreakpositions(MultiRegions::ExpListSharedPtr &streak,
120  Array<OneD, NekDouble> &xc,
121  Array<OneD, NekDouble> &yc,
122  NekDouble cr,
123  NekDouble trans)
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 }
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2