Nektar++
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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description:
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <LibUtilities/Kernel/kernel.h>
36 
37 namespace Nektar
38 {
39  namespace LibUtilities
40  {
41  Kernel::Kernel(int order)
42  {
43  UpdateKernelOrder(order);
44  UpdateKernelNumOfCoeffs();
45  UpdateKernelWidth();
46  UpdateKernelBspline();
47  UpdateKernelCoeffs();
48  k_center = 0.0;
49  }
50 
51  void Kernel::UpdateKernelBspline()
52  {
53  Array<TwoD,NekDouble> out_final(k_order,k_order);
54 
55  int i,j;
56 
57  if(k_order == 1.0)
58  {
59  b_spline[0][0] = 1.0;
60 
61  }else if(k_order == 2)
62  {
63  NekDouble out[2][2] = {{1.0,0.0},
64  {-1.0,1.0}};
65  for(i = 0; i < k_order; i++)
66  {
67  for(j = 0; j < k_order; j++)
68  {
69  out_final[i][j] = out[i][j];
70  }
71  }
72 
73  }else if(k_order == 3)
74  {
75  NekDouble out[3][3] = {{0.5, 0.0, 0.0},
76  {-1.0, 1.0, 0.5},
77  {0.5, -1.0, 0.5}};
78  for(i = 0; i < k_order; i++)
79  {
80  for(j = 0; j < k_order; j++)
81  {
82  out_final[i][j] = out[i][j];
83  }
84  }
85 
86  }else if(k_order == 4)
87  {
88  NekDouble out[4][4] = {{1.0/6.0, 0, 0, 0},
89  {-0.5, 0.5, 0.5, 1.0/6.0},
90  {0.5, -1.0, 0, 2.0/3},
91  {-1.0/6.0, 0.5, -0.5, 1.0/6.0}};
92  for(i = 0; i < k_order; i++)
93  {
94  for(j = 0; j < k_order; j++)
95  {
96  out_final[i][j] = out[i][j];
97  }
98  }
99 
100  }else if(k_order == 5)
101  {
102  NekDouble out[5][5] = {{1.0 / 24.0, 0, 0, 0, 0},
103  {-1.0/6, 1.0/6, 0.25, 1.0/6, 1.0/24},
104  {0.25, -0.5, -0.25, 0.5, 11.0/24},
105  {-1.0/6, 0.5, -0.25, -0.5, 11.0/24},
106  {1.0/24, -1.0/6.0, 0.25, -1.0/6, 1.0/24}};
107  for(i = 0; i < k_order; i++)
108  {
109  for(j = 0; j < k_order; j++)
110  {
111  out_final[i][j] = out[i][j];
112  }
113  }
114 
115  }else if(k_order == 6)
116  {
117  NekDouble out[6][6] = {{1.0/1.020, 0, 0, 0, 0, 0},
118  {-1.0/24, 1.0/24, 1.0/12, 1.0/12, 1.0/24, 1.0/120},
119  {1.0/12, -1.0/6, -1.0/6, 1.0/6, 5.0/12, 13.0/60},
120  {-1.0/12, 0.25, 0, -1.0/2, 0, 11.0/20},
121  {1.0/24, -1.0/6, 1.0/6, 1.0/6, -5.0/12, 13.0/60},
122  {-1.0/120, 1.0/24, -1.0/12, 1.0/12, -1.0/24, 1.0/120}};
123  for(i = 0; i < k_order; i++)
124  {
125  for(j = 0; j < k_order; j++)
126  {
127  out_final[i][j] = out[i][j];
128  }
129  }
130 
131  }else if(k_order == 7)
132  {
133  NekDouble out[7][7] = {{1.0/720, 0, 0, 0, 0, 0, 0},
134  {-1.0/120, 1.0/120, 1.0/48, 1.0/36, 1.0/48, 1.0/120, 1.0/720},
135  {1.0/48, -1.0/24, -1.0/16, 1.0/36, 3.0/16, 5.0/24, 19.0/240},
136  {-1.0/36, 1.0/1.02, 1.0/24, -2.0/9, -5.0/24, 1.0/3, 151.0/360},
137  {1.0/48, -1.0/12, 1.0/24, 2.0/9, -5.0/24, -1.0/3, 151.0/360},
138  {-1.0/120, 1.0/24, -1.0/16, -1.0/36, 3.0/16, -5.0/24, 19.0/240},
139  {1.0/720, -1.0/120, 1.0/48, -1.0/36, 1.0/48, -1.0/120, 1.0/720}};
140  for(i = 0; i < k_order; i++)
141  {
142  for(j = 0; j < k_order; j++)
143  {
144  out_final[i][j] = out[i][j];
145  }
146  }
147 
148  }else if(k_order == 8)
149  {
150  NekDouble out[8][8] = {{1.0/5040, 0, 0, 0, 0, 0, 0, 0},
151  {-1.0/720, 1.0/720, 1.0/240, 1.0/144, 1.0/144, 1.0/240, 1.0/720, 1.0/5040 },
152  {1.0/240, -1.0/120, -1.0/60, 0, 1.0/18, 1.0/10, 7.0/90, 1.0/42},
153  {-1.0/144, 1.0/48, 1.0/48, -1.0/16, -19.0/144, 1.0/16, 49.0/144, 397.0/1680},
154  {1.0/144, -1.0/36, 0, 1.0/9, 0, -1.0/3, 0, 151.0/315},
155  {-1.0/240, 1.0/48, -1.0/48, -1.0/16, 19.0/144, 1.0/16, -49.0/144, 397.0/1680},
156  { 1.0/720, -1.0/120, 1.0/60, 0, -1.0/18, 1.0/10, -7.0/90, 1.0/42},
157  { -1.0/5040, 1.0/720, -1.0/240, 1.0/144, -1.0/144, 1.0/240, -1.0/720, 1.0/5040}};
158  for(i = 0; i < k_order; i++)
159  {
160  for(j = 0; j < k_order; j++)
161  {
162  out_final[i][j] = out[i][j];
163  }
164  }
165 
166  }
167 
168  b_spline = out_final;
169 
170  ASSERTL0(k_order <= 8, "Order is not supported");
171 
172 
173  }
174 
175  void Kernel::UpdateKernelCoeffs()
176  {
177  int i;
178  Array<OneD,NekDouble> out_final(k_ncoeffs);
179 
180  if (k_order == 2)
181  {
182  NekDouble out[3] = {-1.0/12, 7.0/6, -1.0/12};
183  for(i = 0; i < k_ncoeffs; i++)
184  {
185  out_final[i] = out[i];
186  }
187  }else if (k_order == 3)
188  {
189  NekDouble out[5] = {37.0/1920, -97.0/480, 437.0/320, -97.0/480, 37.0/1920};
190  for(i = 0; i < k_ncoeffs; i++)
191  {
192  out_final[i] = out[i];
193  }
194  }else if (k_order == 4)
195  {
196  NekDouble out[7] = {-41.0/7560, 311.0/5040, -919.0/2520, 12223.0/7560, -919.0/2520, 311.0/5040, -41.0/7560};
197  for(i = 0; i < k_ncoeffs; i++)
198  {
199  out_final[i] = out[i];
200  }
201  }else if (k_order == 5)
202  {
203  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};
204  for(i = 0; i < k_ncoeffs; i++)
205  {
206  out_final[i] = out[i];
207  }
208  }else if (k_order == 6)
209  {
210  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};
211  for(i = 0; i < k_ncoeffs; i++)
212  {
213  out_final[i] = out[i];
214  }
215  }else if (k_order == 7)
216  {
217  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};
218  for(i = 0; i < k_ncoeffs; i++)
219  {
220  out_final[i] = out[i];
221  }
222  }else if (k_order == 8)
223  {
224  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};
225  for(i = 0; i < k_ncoeffs; i++)
226  {
227  out_final[i] = out[i];
228  }
229  }
230 
231  k_coeffs = out_final;
232  ASSERTL0(k_order <= 8, "Order is not supported");
233 
234  }
235 
236  void Kernel::UpdateKernelBreaks(NekDouble h)
237  {
238  int i;
239  Array<OneD,NekDouble> temp(k_width+1);
240  temp[0] = -(k_width/2.0)*h; // it needs to get scaled by h
241  for(i = 1; i < k_width+1; i++)
242  {
243  temp[i] = temp[i-1]+h;
244  }
245  k_breaks = temp;
246  }
247 
248  void Kernel::MoveKernelCenter(NekDouble x_value, Array<OneD,NekDouble> &outarray)
249  {
250  int i;
251  for(i = 0; i < k_width+1; i++)
252  {
253  outarray[i] = k_breaks[i] + x_value;
254  }
255 
256  // Update the center of the kernel
257  k_center = x_value;
258  }
259 
260  void Kernel::FindMeshUnderKernel(Array<OneD,NekDouble> &inarray, NekDouble h,
261  Array<OneD,NekDouble> &outarray)
262  {
263  int j;
264  NekDouble first = ceil(inarray[0]/h)*h;
265  int index = k_width;
266  NekDouble last = floor(inarray[index]/h)*h;
267  int count = (int)((last-first)/h)+1; // number of mesh breaks under the kernel support
268  Array<OneD,NekDouble> mesh_breaks(count);
269  mesh_breaks[0] = first;
270  for(j = 1; j < count; j++)
271  {
272  mesh_breaks[j] = mesh_breaks[j-1]+h;
273  }
274  outarray = mesh_breaks;
275 
276  }
277 
278  void Kernel::EvaluateKernel(Array<OneD,NekDouble> inarray,NekDouble h,
279  Array<OneD,NekDouble> &outarray)
280  {
281  int gamma,i;
282  int degree = k_order - 1;
283  int nvalues = inarray.num_elements();
284  Array<OneD,NekDouble> bs_values(nvalues);
285 
286  for(i = 0; i < nvalues; i++ )
287  {
288  outarray[i] = 0.0;
289  }
290 
291  for(gamma = -degree; gamma <= degree; gamma++)
292  {
293  int cIndex = gamma+degree;
294 
295  // Evaluate the bSpline values
296  EvaluateBspline(inarray,h,k_center+(gamma*h),bs_values);
297 
298  for(i = 0; i < nvalues; i++ )
299  {
300  outarray[i] += k_coeffs[cIndex]*bs_values[i];
301  }
302  }
303 
304  }
305 
306  void Kernel::EvaluateBspline(Array<OneD,NekDouble> inarray, NekDouble h,
307  NekDouble offset, Array<OneD,NekDouble> &outarray)
308  {
309  int i;
310  NekDouble min_value = -k_order/2.0;
311  NekDouble max_value = k_order/2.0;
312 
313  int nvalues = inarray.num_elements();
314 
315  // Make a copy for further modifications
316  Array<OneD,NekDouble> inarray_cp(nvalues);
317 
318  for(i = 0; i < nvalues; i++)
319  {
320  inarray_cp[i] = inarray[i] - offset;
321  inarray_cp[i] = inarray_cp[i]/h;
322  int interval = (int)floor(inarray_cp[i] - min_value); // determines to which interval of the bspline the value belongs
323 
324  if(inarray_cp[i] >= min_value && inarray_cp[i] <= max_value)
325  {
326  if(interval >= k_order)
327  {
328  interval -= 1;
329  }
330  NekDouble shift = min_value + interval;
331  inarray_cp[i] -= shift;
332  outarray[i] = EvaluateBsplinePoly(inarray_cp[i],interval);
333  }else
334  {
335  outarray[i] = 0.0;
336  }
337 
338  }
339  }
340 
341  NekDouble Kernel::EvaluateBsplinePoly(NekDouble x_value,int interval)
342  {
343  int i;
344  int deg = k_order - 1;
345  NekDouble poly_value = b_spline[interval][0];
346 
347  for(i = 0; i < deg; i++)
348  {
349  poly_value = poly_value*x_value + b_spline[interval][i+1];
350 
351  }
352 
353  return poly_value;
354  }
355 
356  void Kernel::Sort(Array<OneD,NekDouble> &inarray1, Array<OneD,NekDouble> &inarray2,
357  Array<OneD,NekDouble> &outarray)
358  {
359  int j;
360  int kIndex = 0; // Keeps track of the kernel breaks
361  int mIndex = 0; // Keeps track of the mesh breaks
362  for(j = 0; j < outarray.num_elements(); j++)
363  {
364  if(mIndex >= inarray2.num_elements())
365  {
366  outarray[j] = inarray1[kIndex];
367  kIndex++;
368  }else if(kIndex >= inarray1.num_elements())
369  {
370  outarray[j] = inarray2[mIndex];
371  mIndex++;
372 
373  }else if(inarray1[kIndex] < inarray2[mIndex])
374  {
375  outarray[j] = inarray1[kIndex];
376  kIndex++;
377  }else
378  {
379  outarray[j] = inarray2[mIndex];
380  mIndex++;
381  }
382 
383  }
384  }
385 
386  }// end of namespace
387 }// end of namespace
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::LibUtilities::Kernel::EvaluateKernel
void EvaluateKernel(Array< OneD, NekDouble > inarray, NekDouble h, Array< OneD, NekDouble > &outarray)
This funciton evaluates the kernel at input values.
Definition: kernel.cpp:278
Nektar::LibUtilities::Kernel::Kernel
Kernel()
The default constructor.
Nektar::LibUtilities::Kernel::b_spline
Array< TwoD, NekDouble > b_spline
Definition: kernel.h:200
Nektar::LibUtilities::Kernel::UpdateKernelNumOfCoeffs
void UpdateKernelNumOfCoeffs()
This funciton sets the k_ncoeffs variable.
Definition: kernel.h:113
Nektar::LibUtilities::Kernel::EvaluateBsplinePoly
NekDouble EvaluateBsplinePoly(NekDouble x_value, int interval)
This funciton evaluates the piecewise bspline polynomial.
Definition: kernel.cpp:341
Nektar::LibUtilities::Kernel::EvaluateBspline
void EvaluateBspline(Array< OneD, NekDouble > inarray, NekDouble h, NekDouble offset, Array< OneD, NekDouble > &outarray)
This function evaluates the bspline at input values.
Definition: kernel.cpp:306
Nektar::LibUtilities::Kernel::FindMeshUnderKernel
void FindMeshUnderKernel(Array< OneD, NekDouble > &inarray, NekDouble h, Array< OneD, NekDouble > &outarray)
This funciton calculates the mesh breaks under the kernel support.
Definition: kernel.cpp:260
Nektar::LibUtilities::Kernel::UpdateKernelOrder
void UpdateKernelOrder(int order)
This funciton sets the k_order variable.
Definition: kernel.h:105
Nektar::LibUtilities::Kernel::k_coeffs
Array< OneD, NekDouble > k_coeffs
Definition: kernel.h:201
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::LibUtilities::Kernel::UpdateKernelCoeffs
void UpdateKernelCoeffs()
This funciton updates the kernel coefficients.
Definition: kernel.cpp:175
Nektar::LibUtilities::Kernel::MoveKernelCenter
void MoveKernelCenter(NekDouble x_value, Array< OneD, NekDouble > &outarray)
This funciton moves the center of the kernel to the.
Definition: kernel.cpp:248
Nektar::LibUtilities::Kernel::UpdateKernelBreaks
void UpdateKernelBreaks(NekDouble h)
This funciton updates the kernel breaks.
Definition: kernel.cpp:236
Nektar::LibUtilities::Kernel::k_breaks
Array< OneD, NekDouble > k_breaks
Definition: kernel.h:202
Nektar::LibUtilities::Kernel::Sort
void Sort(Array< OneD, NekDouble > &inarray1, Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
This funciton performs the ordered merge of.
Definition: kernel.cpp:356
Nektar::LibUtilities::Kernel::UpdateKernelWidth
void UpdateKernelWidth()
This funciton sets the kernel width size.
Definition: kernel.h:121
Nektar::LibUtilities::Kernel::UpdateKernelBspline
void UpdateKernelBspline()
The default destructor.
Definition: kernel.cpp:51
Generated on Wed Jan 20 2021 03:33:36 for Nektar++ by   doxygen 1.8.13