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