Code-Eli  0.3.6
d2o2.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * CopyRIGHT (c) 2013 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All RIGHTs reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef eli_mutil_fd_d2o2_hpp
14 #define eli_mutil_fd_d2o2_hpp
15 
16 #include <vector>
17 
18 #include "eli/code_eli.hpp"
19 
20 namespace eli
21 {
22  namespace mutil
23  {
24  namespace fd
25  {
26  template<typename data__>
27  class d2o2
28  {
29  public:
30  enum stencil
31  {
32  LEFT=0,
36  };
37 
38  private:
39  const size_t nnodes;
40  const int n_order;
42 
43  protected:
44  template<typename __itc, typename __itphi> data__ calculate_dot(__itc a, __itphi itphi) const
45  {
46  data__ d(static_cast<data__>(0));
47 
48  for (size_t i=0; i<number_nodes(); ++i, ++itphi)
49  {
50  if (a[i]!=0)
51  d+=a[i]*(*itphi);
52  }
53 
54  return d;
55  }
56 
57  public:
58  d2o2() : nnodes(4), n_order(2), st(d2o2<data__>::RIGHT_BIASED)
59  {
60  }
61 
62  d2o2(const stencil &s) : nnodes(4), n_order(2), st(s)
63  {
64  }
65 
66  d2o2(const d2o2<data__> &d) : nnodes(4), n_order(2), st(d.st)
67  {
68  }
69 
71  {
72  }
73 
74  void set_stencil(const stencil &s)
75  {
76  st=s;
77  }
78 
79  const stencil & get_stencil() const
80  {
81  return st;
82  }
83 
84  int order(bool /*uniform*/) const
85  {
86  return n_order;
87  }
88 
89  size_t number_nodes() const
90  {
91  return nnodes;
92  }
93 
94  template<typename __iti> std::ptrdiff_t index(__iti iti) const
95  {
96  size_t i0;
97 
98  switch(st)
99  {
100  case(LEFT):
101  {
102  i0=3;
103  (*iti)=-3;++iti;
104  (*iti)=-2;++iti;
105  (*iti)=-1;++iti;
106  (*iti)= 0;
107  break;
108  }
109  case(LEFT_BIASED):
110  {
111  i0=2;
112  (*iti)=-2;++iti;
113  (*iti)=-1;++iti;
114  (*iti)= 0;++iti;
115  (*iti)= 1;
116  break;
117  }
118  case(RIGHT_BIASED):
119  {
120  i0=1;
121  (*iti)=-1;++iti;
122  (*iti)= 0;++iti;
123  (*iti)= 1;++iti;
124  (*iti)= 2;
125  break;
126  }
127  case(RIGHT):
128  {
129  i0=0;
130  (*iti)=0;++iti;
131  (*iti)=1;++iti;
132  (*iti)=2;++iti;
133  (*iti)=3;
134  break;
135  }
136  default:
137  {
138  i0=4;
139  assert(false);
140  return i0;
141  }
142  }
143 
144  return i0;
145  }
146 
147  template<typename __itphi> int evaluate(data__ &d, __itphi itphi, const data__ &dx) const
148  {
149  std::vector<data__> a(number_nodes());
150  int rtn;
151 
152  rtn=coefficients(a.begin(), dx);
153  if (rtn==0)
154  d=calculate_dot(a.begin(), itphi);
155 
156  return rtn;
157  }
158 
159  template<typename __itphi, typename __itx> int evaluate(data__ &d, __itphi itphi, __itx itx) const
160  {
161  std::vector<data__> a(number_nodes());
162  int rtn;
163 
164  rtn=coefficients(a.begin(), itx);
165  if (rtn==0)
166  d=calculate_dot(a.begin(), itphi);
167 
168  return rtn;
169  }
170 
171  template<typename __itc> int coefficients(__itc itc, const data__ &dx) const
172  {
173  switch(st)
174  {
175  case(LEFT):
176  {
177  (*itc)=static_cast<data__>(-1)/dx/dx;++itc;
178  (*itc)=static_cast<data__>( 4)/dx/dx;++itc;
179  (*itc)=static_cast<data__>(-5)/dx/dx;++itc;
180  (*itc)=static_cast<data__>( 2)/dx/dx;
181 
182  break;
183  }
184  case(LEFT_BIASED):
185  {
186  (*itc)=static_cast<data__>( 0)/dx/dx;++itc;
187  (*itc)=static_cast<data__>( 1)/dx/dx;++itc;
188  (*itc)=static_cast<data__>(-2)/dx/dx;++itc;
189  (*itc)=static_cast<data__>( 1)/dx/dx;
190 
191  break;
192  }
193  case(RIGHT_BIASED):
194  {
195  (*itc)=static_cast<data__>( 1)/dx/dx;++itc;
196  (*itc)=static_cast<data__>(-2)/dx/dx;++itc;
197  (*itc)=static_cast<data__>( 1)/dx/dx;++itc;
198  (*itc)=static_cast<data__>( 0)/dx/dx;
199 
200  break;
201  }
202  case(RIGHT):
203  {
204  (*itc)=static_cast<data__>( 2)/dx/dx;++itc;
205  (*itc)=static_cast<data__>(-5)/dx/dx;++itc;
206  (*itc)=static_cast<data__>( 4)/dx/dx;++itc;
207  (*itc)=static_cast<data__>(-1)/dx/dx;
208 
209  break;
210  }
211  default:
212  {
213  assert(false);
214  return -1;
215  }
216  }
217 
218  return 0;
219  }
220 
221  template<typename __itc, typename __itx> int coefficients(__itc itc, __itx itx) const
222  {
223  std::vector<data__> x(number_nodes());
224 
225  // extract the x-locations
226  for (size_t i=0; i<number_nodes(); ++i, ++itx)
227  x[i]=(*itx);
228 
229  switch(st)
230  {
231  case(LEFT):
232  {
233  data__ alphaim1, betaim2, gammaim3;
234  alphaim1=x[3]-x[2];
235  betaim2=x[3]-x[1];
236  gammaim3=x[3]-x[0];
237 
238  (*itc)=-2*(alphaim1+betaim2)/(gammaim3*(gammaim3-betaim2)*(gammaim3-alphaim1));++itc;
239  (*itc)=+2*(alphaim1+gammaim3)/(betaim2*(gammaim3-betaim2)*(betaim2-alphaim1));++itc;
240  (*itc)=-2*(gammaim3+betaim2)/(alphaim1*(betaim2-alphaim1)*(gammaim3-alphaim1));++itc;
241  (*itc)=+2*(alphaim1+betaim2+gammaim3)/(alphaim1*betaim2*gammaim3);
242 
243  break;
244  }
245  case(LEFT_BIASED):
246  {
247  data__ alphai, alphaim1, betaim2;
248  alphai=x[3]-x[2];
249  alphaim1=x[2]-x[1];
250  betaim2=x[2]-x[0];
251 
252  (*itc)=+2*(alphai-alphaim1)/(betaim2*(alphai+betaim2)*(betaim2-alphaim1));++itc;
253  (*itc)=+2*(betaim2-alphai)/(alphaim1*(alphai+alphaim1)*(betaim2-alphaim1));++itc;
254  (*itc)=-2*(alphaim1-alphai+betaim2)/(alphaim1*alphai*betaim2);++itc;
255  (*itc)=+2*(alphaim1+betaim2)/(alphai*(alphai+alphaim1)*(alphai+betaim2));
256 
257  break;
258  }
259  case(RIGHT_BIASED):
260  {
261  data__ alphai, alphaim1, betai;
262  alphaim1=x[1]-x[0];
263  alphai=x[2]-x[1];
264  betai=x[3]-x[1];
265 
266  (*itc)=+2*(alphai+betai)/(alphaim1*(alphaim1+betai)*(alphai+alphaim1));++itc;
267  (*itc)=-2*(alphai-alphaim1+betai)/(alphai*alphaim1*betai);++itc;
268  (*itc)=+2*(betai-alphaim1)/(alphai*(alphai+alphaim1)*(betai-alphai));++itc;
269  (*itc)=-2*(alphai-alphaim1)/(betai*(alphaim1+betai)*(betai-alphai));
270 
271  break;
272  }
273  // FIX: Change for 3rd order
274  case(RIGHT):
275  {
276  data__ alphai, betai, gammai;
277  alphai=x[1]-x[0];
278  betai=x[2]-x[0];
279  gammai=x[3]-x[0];
280 
281  (*itc)=+2*(alphai+betai+gammai)/(alphai*betai*gammai);++itc;
282  (*itc)=-2*(betai+gammai)/(alphai*(gammai-alphai)*(betai-alphai));++itc;
283  (*itc)=+2*(alphai+gammai)/(betai*(gammai-betai)*(betai-alphai));++itc;
284  (*itc)=-2*(alphai+betai)/(gammai*(gammai-alphai)*(gammai-betai));
285 
286  break;
287  }
288  default:
289  {
290  assert(false);
291  return -1;
292  }
293  }
294 
295  return 0;
296  }
297 
298  int truncation_error(data__ &te, const data__ &phi4, const data__ &dx) const
299  {
300  switch(st)
301  {
302  case(LEFT):
303  {
304  te=phi4*11*dx*dx/12;
305  break;
306  }
307  case(LEFT_BIASED):
308  {
309  te=-phi4*dx*dx/12;
310  break;
311  }
312  case(RIGHT_BIASED):
313  {
314  te=-phi4*dx*dx/12;
315  break;
316  }
317  case(RIGHT):
318  {
319  te=phi4*11*dx*dx/12;
320  break;
321  }
322  default:
323  {
324  assert(false);
325  return -1;
326  }
327  }
328 
329  return 0;
330  }
331 
332  template<typename __itx> int truncation_error(data__ &te, const data__ &phi4, __itx itx) const
333  {
334  std::vector<data__> x(number_nodes());
335 
336  // extract the x-locations
337  for (size_t i=0; i<number_nodes(); ++i, ++itx)
338  x[i]=(*itx);
339 
340  switch(st)
341  {
342  case(LEFT):
343  {
344  data__ alphaim1, betaim2, gammaim3;
345  alphaim1=x[3]-x[2];
346  betaim2=x[3]-x[1];
347  gammaim3=x[3]-x[0];
348 
349  te=phi4*(alphaim1*gammaim3+betaim2*(alphaim1+gammaim3))/static_cast<data__>(12);
350  break;
351  }
352  // FIX: These need to be fixed for d2o2
353  case(LEFT_BIASED):
354  {
355  data__ alphai, alphaim1, betaim2;
356  alphai=x[3]-x[2];
357  alphaim1=x[2]-x[1];
358  betaim2=x[2]-x[0];
359 
360  te=-phi4*(alphai*(alphaim1+betaim2)-alphaim1*betaim2)/static_cast<data__>(12);
361  break;
362  }
363  case(RIGHT_BIASED):
364  {
365  data__ alphai, alphaim1, betai;
366  alphaim1=x[1]-x[0];
367  alphai=x[2]-x[1];
368  betai=x[3]-x[1];
369 
370  te=-phi4*(alphaim1*(alphai+betai)-alphai*betai)/static_cast<data__>(12);
371  break;
372  }
373  case(RIGHT):
374  {
375  data__ alphai, betai, gammai;
376  alphai=x[1]-x[0];
377  betai=x[2]-x[0];
378  gammai=x[3]-x[0];
379 
380  te=phi4*(alphai*gammai+betai*(alphai+gammai))/static_cast<data__>(12);
381  break;
382  }
383  default:
384  {
385  assert(false);
386  return -1;
387  }
388  }
389 
390  return 0;
391  }
392  };
393  }
394  }
395 }
396 
397 #endif
stencil
Definition: d2o2.hpp:30
Definition: math.hpp:20
std::ptrdiff_t index(__iti iti) const
Definition: d2o2.hpp:94
const int n_order
Definition: d2o2.hpp:40
data__ calculate_dot(__itc a, __itphi itphi) const
Definition: d2o2.hpp:44
void set_stencil(const stencil &s)
Definition: d2o2.hpp:74
int order(bool) const
Definition: d2o2.hpp:84
Definition: d2o2.hpp:35
int evaluate(data__ &d, __itphi itphi, const data__ &dx) const
Definition: d2o2.hpp:147
~d2o2()
Definition: d2o2.hpp:70
const size_t nnodes
Definition: d2o2.hpp:39
int coefficients(__itc itc, const data__ &dx) const
Definition: d2o2.hpp:171
Definition: d2o2.hpp:32
size_t number_nodes() const
Definition: d2o2.hpp:89
Definition: d2o2.hpp:27
int truncation_error(data__ &te, const data__ &phi4, __itx itx) const
Definition: d2o2.hpp:332
int truncation_error(data__ &te, const data__ &phi4, const data__ &dx) const
Definition: d2o2.hpp:298
Definition: d2o2.hpp:34
d2o2()
Definition: d2o2.hpp:58
stencil st
Definition: d2o2.hpp:41
d2o2(const stencil &s)
Definition: d2o2.hpp:62
Definition: d2o2.hpp:33
d2o2(const d2o2< data__ > &d)
Definition: d2o2.hpp:66
int evaluate(data__ &d, __itphi itphi, __itx itx) const
Definition: d2o2.hpp:159
const stencil & get_stencil() const
Definition: d2o2.hpp:79
int coefficients(__itc itc, __itx itx) const
Definition: d2o2.hpp:221