Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 using namespace std;
41 using namespace Nektar;
42 
46  NekDouble cr,
47  NekDouble trans)
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
STL namespace.
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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