Nektar++
ExtractCriticalLayerFunctions.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ComputeCriticalLayer.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: Compute location of critical layer from streak file
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <cstdio>
37 
38 #include <MultiRegions/ExpList.h>
39 
40 
44  NekDouble cr,
45  NekDouble trans)
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
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
void Computestreakpositions(MultiRegions::ExpListSharedPtr &streak, Array< OneD, NekDouble > &xc, Array< OneD, NekDouble > &yc, NekDouble cr, NekDouble trans)
static const NekDouble kNekUnsetDouble