Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
kernel.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File kernel.h
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Scientific Computing and Imaging Institute,
10 // University of Utah (USA) and Department of Aeronautics, Imperial
11 // College London (UK).
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:
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/Kernel/kernel.h>
37 
38 namespace Nektar
39 {
40  namespace LibUtilities
41  {
42  Kernel::Kernel(int order)
43  {
44  UpdateKernelOrder(order);
45  UpdateKernelNumOfCoeffs();
46  UpdateKernelWidth();
47  UpdateKernelBspline();
48  UpdateKernelCoeffs();
49  k_center = 0.0;
50  }
51 
52  void Kernel::UpdateKernelBspline()
53  {
54  Array<TwoD,NekDouble> out_final(k_order,k_order);
55 
56  int i,j;
57 
58  if(k_order == 1.0)
59  {
60  b_spline[0][0] = 1.0;
61 
62  }else if(k_order == 2)
63  {
64  NekDouble out[2][2] = {{1.0,0.0},
65  {-1.0,1.0}};
66  for(i = 0; i < k_order; i++)
67  {
68  for(j = 0; j < k_order; j++)
69  {
70  out_final[i][j] = out[i][j];
71  }
72  }
73 
74  }else if(k_order == 3)
75  {
76  NekDouble out[3][3] = {{0.5, 0.0, 0.0},
77  {-1.0, 1.0, 0.5},
78  {0.5, -1.0, 0.5}};
79  for(i = 0; i < k_order; i++)
80  {
81  for(j = 0; j < k_order; j++)
82  {
83  out_final[i][j] = out[i][j];
84  }
85  }
86 
87  }else if(k_order == 4)
88  {
89  NekDouble out[4][4] = {{1.0/6.0, 0, 0, 0},
90  {-0.5, 0.5, 0.5, 1.0/6.0},
91  {0.5, -1.0, 0, 2.0/3},
92  {-1.0/6.0, 0.5, -0.5, 1.0/6.0}};
93  for(i = 0; i < k_order; i++)
94  {
95  for(j = 0; j < k_order; j++)
96  {
97  out_final[i][j] = out[i][j];
98  }
99  }
100 
101  }else if(k_order == 5)
102  {
103  NekDouble out[5][5] = {{1.0 / 24.0, 0, 0, 0, 0},
104  {-1.0/6, 1.0/6, 0.25, 1.0/6, 1.0/24},
105  {0.25, -0.5, -0.25, 0.5, 11.0/24},
106  {-1.0/6, 0.5, -0.25, -0.5, 11.0/24},
107  {1.0/24, -1.0/6.0, 0.25, -1.0/6, 1.0/24}};
108  for(i = 0; i < k_order; i++)
109  {
110  for(j = 0; j < k_order; j++)
111  {
112  out_final[i][j] = out[i][j];
113  }
114  }
115 
116  }else if(k_order == 6)
117  {
118  NekDouble out[6][6] = {{1.0/1.020, 0, 0, 0, 0, 0},
119  {-1.0/24, 1.0/24, 1.0/12, 1.0/12, 1.0/24, 1.0/120},
120  {1.0/12, -1.0/6, -1.0/6, 1.0/6, 5.0/12, 13.0/60},
121  {-1.0/12, 0.25, 0, -1.0/2, 0, 11.0/20},
122  {1.0/24, -1.0/6, 1.0/6, 1.0/6, -5.0/12, 13.0/60},
123  {-1.0/120, 1.0/24, -1.0/12, 1.0/12, -1.0/24, 1.0/120}};
124  for(i = 0; i < k_order; i++)
125  {
126  for(j = 0; j < k_order; j++)
127  {
128  out_final[i][j] = out[i][j];
129  }
130  }
131 
132  }else if(k_order == 7)
133  {
134  NekDouble out[7][7] = {{1.0/720, 0, 0, 0, 0, 0, 0},
135  {-1.0/120, 1.0/120, 1.0/48, 1.0/36, 1.0/48, 1.0/120, 1.0/720},
136  {1.0/48, -1.0/24, -1.0/16, 1.0/36, 3.0/16, 5.0/24, 19.0/240},
137  {-1.0/36, 1.0/1.02, 1.0/24, -2.0/9, -5.0/24, 1.0/3, 151.0/360},
138  {1.0/48, -1.0/12, 1.0/24, 2.0/9, -5.0/24, -1.0/3, 151.0/360},
139  {-1.0/120, 1.0/24, -1.0/16, -1.0/36, 3.0/16, -5.0/24, 19.0/240},
140  {1.0/720, -1.0/120, 1.0/48, -1.0/36, 1.0/48, -1.0/120, 1.0/720}};
141  for(i = 0; i < k_order; i++)
142  {
143  for(j = 0; j < k_order; j++)
144  {
145  out_final[i][j] = out[i][j];
146  }
147  }
148 
149  }else if(k_order == 8)
150  {
151  NekDouble out[8][8] = {{1.0/5040, 0, 0, 0, 0, 0, 0, 0},
152  {-1.0/720, 1.0/720, 1.0/240, 1.0/144, 1.0/144, 1.0/240, 1.0/720, 1.0/5040 },
153  {1.0/240, -1.0/120, -1.0/60, 0, 1.0/18, 1.0/10, 7.0/90, 1.0/42},
154  {-1.0/144, 1.0/48, 1.0/48, -1.0/16, -19.0/144, 1.0/16, 49.0/144, 397.0/1680},
155  {1.0/144, -1.0/36, 0, 1.0/9, 0, -1.0/3, 0, 151.0/315},
156  {-1.0/240, 1.0/48, -1.0/48, -1.0/16, 19.0/144, 1.0/16, -49.0/144, 397.0/1680},
157  { 1.0/720, -1.0/120, 1.0/60, 0, -1.0/18, 1.0/10, -7.0/90, 1.0/42},
158  { -1.0/5040, 1.0/720, -1.0/240, 1.0/144, -1.0/144, 1.0/240, -1.0/720, 1.0/5040}};
159  for(i = 0; i < k_order; i++)
160  {
161  for(j = 0; j < k_order; j++)
162  {
163  out_final[i][j] = out[i][j];
164  }
165  }
166 
167  }
168 
169  b_spline = out_final;
170 
171  ASSERTL0(k_order <= 8, "Order is not supported");
172 
173 
174  }
175 
176  void Kernel::UpdateKernelCoeffs()
177  {
178  int i;
179  Array<OneD,NekDouble> out_final(k_ncoeffs);
180 
181  if (k_order == 2)
182  {
183  NekDouble out[3] = {-1.0/12, 7.0/6, -1.0/12};
184  for(i = 0; i < k_ncoeffs; i++)
185  {
186  out_final[i] = out[i];
187  }
188  }else if (k_order == 3)
189  {
190  NekDouble out[5] = {37.0/1920, -97.0/480, 437.0/320, -97.0/480, 37.0/1920};
191  for(i = 0; i < k_ncoeffs; i++)
192  {
193  out_final[i] = out[i];
194  }
195  }else if (k_order == 4)
196  {
197  NekDouble out[7] = {-41.0/7560, 311.0/5040, -919.0/2520, 12223.0/7560, -919.0/2520, 311.0/5040, -41.0/7560};
198  for(i = 0; i < k_ncoeffs; i++)
199  {
200  out_final[i] = out[i];
201  }
202  }else if (k_order == 5)
203  {
204  NekDouble out[9] = {153617.0/92897280, -35411.0/1658880, 3153959.0/23224320, -6803459.0/11612160, 18017975.0/9289728, -6803459.0/11612160, 3153959.0/23224320, -35411.0/1658880, 153617.0/92897280};
205  for(i = 0; i < k_ncoeffs; i++)
206  {
207  out_final[i] = out[i];
208  }
209  }else if (k_order == 6)
210  {
211  NekDouble out[11] = {-4201.0/7983360, 30773.0/3991680, -20813.0/380160, 2825.0/11088, -1179649.0/1330560, 1569217.0/665280, -1179649.0/1330560, 2825.0/11088, -20813.0/380160, 30773.0/3991680, -4201.0/7983360};
212  for(i = 0; i < k_ncoeffs; i++)
213  {
214  out_final[i] = out[i];
215  }
216  }else if (k_order == 7)
217  {
218  NekDouble out[13] = {13154671847.0/76517631590400.0, -18073154507.0/6376469299200.0, 287360344573.0/12752938598400.0, -2217732343517.0/19129407897600.0, 1240941746699.0/2833986355200.0, -275386671493.0/212548976640.0, 2648644782397.0/910924185600.0, -275386671493.0/212548976640.0, 1240941746699.0/2833986355200.0, -2217732343517.0/19129407897600.0, 287360344573.0/12752938598400.0, -18073154507.0/6376469299200.0, 13154671847.0/76517631590400.0};
219  for(i = 0; i < k_ncoeffs; i++)
220  {
221  out_final[i] = out[i];
222  }
223  }else if (k_order == 8)
224  {
225  NekDouble out[15] = {-800993.0/14010796800.0, 73587167.0/70053984000.0, -651305719.0/70053984000.0, 3714581677.0/70053984000.0, -3085236289.0/14010796800.0, 1426328231.0/2001542400.0, -43268401973.0/23351328000.0, 42401344373.0/11675664000.0, -43268401973.0/23351328000.0, 1426328231.0/2001542400.0, -3085236289.0/14010796800.0, 3714581677.0/70053984000.0, -651305719.0/70053984000.0, 73587167.0/70053984000.0, -800993.0/14010796800.0};
226  for(i = 0; i < k_ncoeffs; i++)
227  {
228  out_final[i] = out[i];
229  }
230  }
231 
232  k_coeffs = out_final;
233  ASSERTL0(k_order <= 8, "Order is not supported");
234 
235  }
236 
237  void Kernel::UpdateKernelBreaks(NekDouble h)
238  {
239  int i;
240  Array<OneD,NekDouble> temp(k_width+1);
241  temp[0] = -(k_width/2.0)*h; // it needs to get scaled by h
242  for(i = 1; i < k_width+1; i++)
243  {
244  temp[i] = temp[i-1]+h;
245  }
246  k_breaks = temp;
247  }
248 
249  void Kernel::MoveKernelCenter(NekDouble x_value, Array<OneD,NekDouble> &outarray)
250  {
251  int i;
252  for(i = 0; i < k_width+1; i++)
253  {
254  outarray[i] = k_breaks[i] + x_value;
255  }
256 
257  // Update the center of the kernel
258  k_center = x_value;
259  }
260 
261  void Kernel::FindMeshUnderKernel(Array<OneD,NekDouble> &inarray, NekDouble h,
262  Array<OneD,NekDouble> &outarray)
263  {
264  int j;
265  NekDouble first = ceil(inarray[0]/h)*h;
266  int index = k_width;
267  NekDouble last = floor(inarray[index]/h)*h;
268  int count = (int)((last-first)/h)+1; // number of mesh breaks under the kernel support
269  Array<OneD,NekDouble> mesh_breaks(count);
270  mesh_breaks[0] = first;
271  for(j = 1; j < count; j++)
272  {
273  mesh_breaks[j] = mesh_breaks[j-1]+h;
274  }
275  outarray = mesh_breaks;
276 
277  }
278 
279  void Kernel::EvaluateKernel(Array<OneD,NekDouble> inarray,NekDouble h,
280  Array<OneD,NekDouble> &outarray)
281  {
282  int gamma,i;
283  int degree = k_order - 1;
284  int nvalues = inarray.num_elements();
285  Array<OneD,NekDouble> bs_values(nvalues);
286 
287  for(i = 0; i < nvalues; i++ )
288  {
289  outarray[i] = 0.0;
290  }
291 
292  for(gamma = -degree; gamma <= degree; gamma++)
293  {
294  int cIndex = gamma+degree;
295 
296  // Evaluate the bSpline values
297  EvaluateBspline(inarray,h,k_center+(gamma*h),bs_values);
298 
299  for(i = 0; i < nvalues; i++ )
300  {
301  outarray[i] += k_coeffs[cIndex]*bs_values[i];
302  }
303  }
304 
305  }
306 
307  void Kernel::EvaluateBspline(Array<OneD,NekDouble> inarray, NekDouble h,
308  NekDouble offset, Array<OneD,NekDouble> &outarray)
309  {
310  int i;
311  NekDouble min_value = -k_order/2.0;
312  NekDouble max_value = k_order/2.0;
313 
314  int nvalues = inarray.num_elements();
315 
316  // Make a copy for further modifications
317  Array<OneD,NekDouble> inarray_cp(nvalues);
318 
319  for(i = 0; i < nvalues; i++)
320  {
321  inarray_cp[i] = inarray[i] - offset;
322  inarray_cp[i] = inarray_cp[i]/h;
323  int interval = (int)floor(inarray_cp[i] - min_value); // determines to which interval of the bspline the value belongs
324 
325  if(inarray_cp[i] >= min_value && inarray_cp[i] <= max_value)
326  {
327  if(interval >= k_order)
328  {
329  interval -= 1;
330  }
331  NekDouble shift = min_value + interval;
332  inarray_cp[i] -= shift;
333  outarray[i] = EvaluateBsplinePoly(inarray_cp[i],interval);
334  }else
335  {
336  outarray[i] = 0.0;
337  }
338 
339  }
340  }
341 
342  NekDouble Kernel::EvaluateBsplinePoly(NekDouble x_value,int interval)
343  {
344  int i;
345  int deg = k_order - 1;
346  NekDouble poly_value = b_spline[interval][0];
347 
348  for(i = 0; i < deg; i++)
349  {
350  poly_value = poly_value*x_value + b_spline[interval][i+1];
351 
352  }
353 
354  return poly_value;
355  }
356 
357  void Kernel::Sort(Array<OneD,NekDouble> &inarray1, Array<OneD,NekDouble> &inarray2,
358  Array<OneD,NekDouble> &outarray)
359  {
360  int j;
361  int kIndex = 0; // Keeps track of the kernel breaks
362  int mIndex = 0; // Keeps track of the mesh breaks
363  for(j = 0; j < outarray.num_elements(); j++)
364  {
365  if(mIndex >= inarray2.num_elements())
366  {
367  outarray[j] = inarray1[kIndex];
368  kIndex++;
369  }else if(kIndex >= inarray1.num_elements())
370  {
371  outarray[j] = inarray2[mIndex];
372  mIndex++;
373 
374  }else if(inarray1[kIndex] < inarray2[mIndex])
375  {
376  outarray[j] = inarray1[kIndex];
377  kIndex++;
378  }else
379  {
380  outarray[j] = inarray2[mIndex];
381  mIndex++;
382  }
383 
384  }
385  }
386 
387  }// end of namespace
388 }// end of namespace
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2