Code-Eli  0.3.6
d1o3.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_d1o3_hpp
14 #define eli_mutil_fd_d1o3_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 d1o3
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  d+=a[i]*(*itphi);
50 
51  return d;
52  }
53 
54  public:
55  d1o3() : nnodes(4), n_order(3), st(LEFT_BIASED)
56  {
57  }
58 
59  d1o3(const stencil &s) : nnodes(4), n_order(3), st(s)
60  {
61  }
62 
63  d1o3(const d1o3<data__> &d) : nnodes(4), n_order(3), st(d.st)
64  {
65  }
66 
68  {
69  }
70 
71  void set_stencil(const stencil &s)
72  {
73  st=s;
74  }
75  const stencil & get_stencil() const
76  {
77  return st;
78  }
79 
80  int order(bool /*uniform*/) const
81  {
82  return n_order;
83  }
84 
85  size_t number_nodes() const
86  {
87  return nnodes;
88  }
89 
90  template<typename iti__> std::ptrdiff_t index(iti__ iti) const
91  {
92  size_t i0;
93 
94  switch(st)
95  {
96  case(LEFT):
97  {
98  i0=3;
99  (*iti)=-3;++iti;
100  (*iti)=-2;++iti;
101  (*iti)=-1;++iti;
102  (*iti)= 0;
103  break;
104  }
105  case(LEFT_BIASED):
106  {
107  i0=2;
108  (*iti)=-2;++iti;
109  (*iti)=-1;++iti;
110  (*iti)= 0;++iti;
111  (*iti)= 1;
112  break;
113  }
114  case(RIGHT_BIASED):
115  {
116  i0=1;
117  (*iti)=-1;++iti;
118  (*iti)= 0;++iti;
119  (*iti)= 1;++iti;
120  (*iti)= 2;
121  break;
122  }
123  case(RIGHT):
124  {
125  i0=0;
126  (*iti)=0;++iti;
127  (*iti)=1;++iti;
128  (*iti)=2;++iti;
129  (*iti)=3;
130  break;
131  }
132  default:
133  {
134  i0=4;
135  assert(false);
136  return i0;
137  }
138  }
139 
140  return i0;
141  }
142 
143  template<typename itphi__> int evaluate(data__ &d, itphi__ itphi, const data__ &dx) const
144  {
145  std::vector<data__> a(number_nodes());
146  int rtn;
147 
148  rtn=coefficients(a.begin(), dx);
149  if (rtn==0)
150  d=calculate_dot(a.begin(), itphi);
151 
152  return rtn;
153  }
154 
155  template<typename itphi__, typename itx__> int evaluate(data__ &d, itphi__ itphi, itx__ itx) const
156  {
157  std::vector<data__> a(number_nodes());
158  int rtn;
159 
160  rtn=coefficients(a.begin(), itx);
161  if (rtn==0)
162  d=calculate_dot(a.begin(), itphi);
163 
164  return rtn;
165  }
166 
167  template<typename itc__> int coefficients(itc__ itc, const data__ &dx) const
168  {
169  switch(st)
170  {
171  case(LEFT):
172  {
173  (*itc)=(static_cast<data__>(-1.0)/static_cast<data__>(3.0))/dx;++itc;
174  (*itc)=static_cast<data__>( 1.5)/dx;++itc;
175  (*itc)=static_cast<data__>(-3.0)/dx;++itc;
176  (*itc)=(static_cast<data__>( 11)/static_cast<data__>(6.0))/dx;
177 
178  break;
179  }
180  case(LEFT_BIASED):
181  {
182  (*itc)=(static_cast<data__>( 1.0)/static_cast<data__>(6.0))/dx;++itc;
183  (*itc)=static_cast<data__>(-1.0)/dx;++itc;
184  (*itc)=static_cast<data__>(0.5)/dx;++itc;
185  (*itc)=(static_cast<data__>( 1.0)/static_cast<data__>(3.0))/dx;
186 
187  break;
188  }
189  case(RIGHT_BIASED):
190  {
191  (*itc)=(static_cast<data__>(-1.0)/static_cast<data__>(3.0))/dx;++itc;
192  (*itc)=static_cast<data__>(-0.5)/dx;++itc;
193  (*itc)=static_cast<data__>(1.0)/dx;++itc;
194  (*itc)=(static_cast<data__>(-1.0)/static_cast<data__>(6.0))/dx;
195 
196  break;
197  }
198  case(RIGHT):
199  {
200  (*itc)=(static_cast<data__>(-11)/static_cast<data__>(6.0))/dx;++itc;
201  (*itc)=static_cast<data__>( 3.0)/dx;++itc;
202  (*itc)=static_cast<data__>(-1.5)/dx;++itc;
203  (*itc)=(static_cast<data__>( 1.0)/static_cast<data__>(3.0))/dx;
204 
205  break;
206  }
207  default:
208  {
209  assert(false);
210  return -1;
211  }
212  }
213 
214  return 0;
215  }
216 
217  template<typename itc__, typename itx__> int coefficients(itc__ itc, itx__ itx) const
218  {
219  std::vector<data__> x(number_nodes());
220 
221  // extract the x-locations
222  for (size_t i=0; i<number_nodes(); ++i, ++itx)
223  x[i]=(*itx);
224 
225  switch(st)
226  {
227  case(LEFT):
228  {
229  data__ alphaim1, betaim2, gammaim3;
230  alphaim1=x[3]-x[2];
231  betaim2=x[3]-x[1];
232  gammaim3=x[3]-x[0];
233 
234  (*itc)=-(alphaim1*betaim2)/(gammaim3*(gammaim3-betaim2)*(gammaim3-alphaim1));++itc;
235  (*itc)=+(alphaim1*gammaim3)/(betaim2*(gammaim3-betaim2)*(betaim2-alphaim1));++itc;
236  (*itc)=-(gammaim3*betaim2)/(alphaim1*(betaim2-alphaim1)*(gammaim3-alphaim1));++itc;
237  (*itc)=+(1/alphaim1+1/betaim2+1/gammaim3);
238 
239  break;
240  }
241  case(LEFT_BIASED):
242  {
243  data__ alphai, alphaim1, betaim2;
244  alphai=x[3]-x[2];
245  alphaim1=x[2]-x[1];
246  betaim2=x[2]-x[0];
247 
248  (*itc)=+(alphai*alphaim1)/(betaim2*(alphai+betaim2)*(betaim2-alphaim1));++itc;
249  (*itc)=-(betaim2*alphai)/(alphaim1*(alphai+alphaim1)*(betaim2-alphaim1));++itc;
250  (*itc)=+(1/alphaim1-1/alphai+1/betaim2);++itc;
251  (*itc)=+(alphaim1*betaim2)/(alphai*(alphai+alphaim1)*(alphai+betaim2));
252 
253  break;
254  }
255  case(RIGHT_BIASED):
256  {
257  data__ alphai, alphaim1, betai;
258  alphaim1=x[1]-x[0];
259  alphai=x[2]-x[1];
260  betai=x[3]-x[1];
261 
262  (*itc)=-(alphai*betai)/(alphaim1*(alphaim1+betai)*(alphai+alphaim1));++itc;
263  (*itc)=-(1/alphai-1/alphaim1+1/betai);++itc;
264  (*itc)=+(alphaim1*betai)/(alphai*(alphai+alphaim1)*(betai-alphai));++itc;
265  (*itc)=-(alphai*alphaim1)/(betai*(alphaim1+betai)*(betai-alphai));
266 
267  break;
268  }
269  case(RIGHT):
270  {
271  data__ alphai, betai, gammai;
272  alphai=x[1]-x[0];
273  betai=x[2]-x[0];
274  gammai=x[3]-x[0];
275 
276  (*itc)=-(1/alphai+1/betai+1/gammai);++itc;
277  (*itc)=+(betai*gammai)/(alphai*(gammai-alphai)*(betai-alphai));++itc;
278  (*itc)=-(alphai*gammai)/(betai*(gammai-betai)*(betai-alphai));++itc;
279  (*itc)=+(alphai*betai)/(gammai*(gammai-alphai)*(gammai-betai));
280 
281  break;
282  }
283  default:
284  {
285  assert(false);
286  return -1;
287  }
288  }
289 
290  return 0;
291  }
292 
293  int truncation_error(data__ &te, const data__ &phi4, const data__ &dx) const
294  {
295  switch(st)
296  {
297  case(LEFT):
298  {
299  te=phi4*dx*dx*dx/4;
300  break;
301  }
302  case(LEFT_BIASED):
303  {
304  te=-phi4*dx*dx*dx/12;
305  break;
306  }
307  case(RIGHT_BIASED):
308  {
309  te=phi4*dx*dx*dx/12;
310  break;
311  }
312  case(RIGHT):
313  {
314  te=-phi4*dx*dx*dx/4;
315  break;
316  }
317  default:
318  {
319  assert(false);
320  return -1;
321  }
322  }
323 
324  return 0;
325  }
326 
327  template<typename itx__> int truncation_error(data__ &te, const data__ &phi4, itx__ itx) const
328  {
329  std::vector<data__> x(number_nodes());
330 
331  // extract the x-locations
332  for (size_t i=0; i<number_nodes(); ++i, ++itx)
333  x[i]=(*itx);
334 
335  switch(st)
336  {
337  case(LEFT):
338  {
339  data__ alphaim1, betaim2, gammaim3;
340  alphaim1=x[3]-x[2];
341  betaim2=x[3]-x[1];
342  gammaim3=x[3]-x[0];
343 
344  te=phi4*alphaim1*betaim2*gammaim3/static_cast<data__>(24);
345  break;
346  }
347  case(LEFT_BIASED):
348  {
349  data__ alphai, alphaim1, betaim2;
350  alphai=x[3]-x[2];
351  alphaim1=x[2]-x[1];
352  betaim2=x[2]-x[0];
353 
354  te=-phi4*alphai*alphaim1*betaim2/static_cast<data__>(24);
355  break;
356  }
357  case(RIGHT_BIASED):
358  {
359  data__ alphai, alphaim1, betai;
360  alphaim1=x[1]-x[0];
361  alphai=x[2]-x[1];
362  betai=x[3]-x[1];
363 
364  te=phi4*alphai*alphaim1*betai/static_cast<data__>(24);
365  break;
366  }
367  case(RIGHT):
368  {
369  data__ alphai, betai, gammai;
370  alphai=x[1]-x[0];
371  betai=x[2]-x[0];
372  gammai=x[3]-x[0];
373 
374  te=-phi4*alphai*betai*gammai/static_cast<data__>(24);
375  break;
376  }
377  default:
378  {
379  assert(false);
380  return -1;
381  }
382  }
383 
384  return 0;
385  }
386  };
387  }
388  }
389 }
390 
391 #endif
int order(bool) const
Definition: d1o3.hpp:80
Definition: math.hpp:20
stencil
Definition: d1o3.hpp:30
const stencil & get_stencil() const
Definition: d1o3.hpp:75
~d1o3()
Definition: d1o3.hpp:67
Definition: d1o3.hpp:34
Definition: d1o3.hpp:27
size_t number_nodes() const
Definition: d1o3.hpp:85
Definition: d1o3.hpp:32
void set_stencil(const stencil &s)
Definition: d1o3.hpp:71
int coefficients(itc__ itc, const data__ &dx) const
Definition: d1o3.hpp:167
Definition: d1o3.hpp:33
int truncation_error(data__ &te, const data__ &phi4, const data__ &dx) const
Definition: d1o3.hpp:293
data__ calculate_dot(itc__ a, itphi__ itphi) const
Definition: d1o3.hpp:44
d1o3(const stencil &s)
Definition: d1o3.hpp:59
d1o3(const d1o3< data__ > &d)
Definition: d1o3.hpp:63
std::ptrdiff_t index(iti__ iti) const
Definition: d1o3.hpp:90
const size_t nnodes
Definition: d1o3.hpp:39
stencil st
Definition: d1o3.hpp:41
int truncation_error(data__ &te, const data__ &phi4, itx__ itx) const
Definition: d1o3.hpp:327
int coefficients(itc__ itc, itx__ itx) const
Definition: d1o3.hpp:217
int evaluate(data__ &d, itphi__ itphi, itx__ itx) const
Definition: d1o3.hpp:155
Definition: d1o3.hpp:35
int evaluate(data__ &d, itphi__ itphi, const data__ &dx) const
Definition: d1o3.hpp:143
d1o3()
Definition: d1o3.hpp:55
const int n_order
Definition: d1o3.hpp:40