Code-Eli  0.3.6
d1o2.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_d1o2_hpp
14 #define eli_mutil_fd_d1o2_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 d1o2
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  d1o2() : nnodes(3), n_order(2), st(CENTER)
55  {
56  }
57 
58  d1o2(const stencil &s) : nnodes(3), n_order(2), st(s)
59  {
60  }
61 
62  d1o2(const d1o2<data__> &d) : nnodes(3), n_order(2), 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  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=2;
99  (*iti)=-2;++iti;
100  (*iti)=-1;++iti;
101  (*iti)=0;
102  break;
103  }
104  case(CENTER):
105  {
106  i0=1;
107  (*iti)=-1;++iti;
108  (*iti)=0;++iti;
109  (*iti)=1;
110  break;
111  }
112  case(RIGHT):
113  {
114  i0=0;
115  (*iti)=0;++iti;
116  (*iti)=1;++iti;
117  (*iti)=2;
118  break;
119  }
120  default:
121  {
122  i0=3;
123  assert(false);
124  return i0;
125  }
126  }
127 
128  return i0;
129  }
130 
131  template<typename itphi__> int evaluate(data__ &d, itphi__ itphi, const data__ &dx) const
132  {
133  std::vector<data__> a(number_nodes());
134  int rtn;
135 
136  rtn=coefficients(a.begin(), dx);
137  if (rtn==0)
138  d=calculate_dot(a.begin(), itphi);
139 
140  return rtn;
141  }
142 
143  template<typename itphi__, typename itx__> int evaluate(data__ &d, itphi__ itphi, itx__ itx) const
144  {
145  std::vector<data__> a(number_nodes());
146  int rtn;
147 
148  rtn=coefficients(a.begin(), itx);
149  if (rtn==0)
150  d=calculate_dot(a.begin(), itphi);
151 
152  return rtn;
153  }
154 
155  template<typename itc__> int coefficients(itc__ itc, const data__ &dx) const
156  {
157  switch(st)
158  {
159  case(LEFT):
160  {
161  (*itc)=static_cast<data__>( 0.5)/dx;++itc;
162  (*itc)=static_cast<data__>(-2.0)/dx;++itc;
163  (*itc)=static_cast<data__>( 1.5)/dx;
164 
165  break;
166  }
167  case(CENTER):
168  {
169  (*itc)=static_cast<data__>(-0.5)/dx;++itc;
170  (*itc)=static_cast<data__>( 0.0)/dx;++itc;
171  (*itc)=static_cast<data__>( 0.5)/dx;
172 
173  break;
174  }
175  case(RIGHT):
176  {
177  (*itc)=static_cast<data__>(-1.5)/dx;++itc;
178  (*itc)=static_cast<data__>( 2.0)/dx;++itc;
179  (*itc)=static_cast<data__>(-0.5)/dx;
180 
181  break;
182  }
183  default:
184  {
185  assert(false);
186  return -1;
187  }
188  }
189 
190  return 0;
191  }
192 
193  template<typename itc__, typename itx__> int coefficients(itc__ itc, itx__ itx) const
194  {
195  std::vector<data__> x(number_nodes());
196 
197  // extract the x-locations
198  x[0]=(*itx); ++itx;
199  x[1]=(*itx); ++itx;
200  x[2]=(*itx);
201 
202  switch(st)
203  {
204  case(LEFT):
205  {
206  data__ alphaim1, betaim2, alphaim2;
207  alphaim1=x[2]-x[1];
208  betaim2=x[2]-x[0];
209  alphaim2=x[1]-x[0];
210 
211  (*itc)=alphaim1/(betaim2*alphaim2);++itc;
212  (*itc)=-betaim2/(alphaim1*alphaim2);++itc;
213  (*itc)=(alphaim1+betaim2)/(alphaim1*betaim2);
214 
215  break;
216  }
217  case(CENTER):
218  {
219  data__ deli, delim1, gamip1;
220  gamip1=x[2]-x[0];
221  delim1=x[1]-x[0];
222  deli=x[2]-x[1];
223 
224  (*itc)=-deli/(delim1*gamip1);++itc;
225  (*itc)=(deli-delim1)/(deli*delim1);++itc;
226  (*itc)=delim1/(deli*gamip1);
227 
228  break;
229  }
230  case(RIGHT):
231  {
232  data__ alphai, betai, alphaip1;
233  alphai=x[1]-x[0];
234  betai=x[2]-x[0];
235  alphaip1=x[2]-x[1];
236 
237  (*itc)=-(alphai+betai)/(alphai*betai);++itc;
238  (*itc)=betai/(alphai*alphaip1);++itc;
239  (*itc)=-alphai/(betai*alphaip1);
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*dx/3;
260  break;
261  }
262  case(CENTER):
263  {
264  te=-phi3*dx*dx/6;
265  break;
266  }
267  case(RIGHT):
268  {
269  te=phi3*dx*dx/3;
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__ delim1, betaim2;
296 
297  delim1=x[2]-x[1];
298  betaim2=x[2]-x[0];
299 
300  te=phi3*delim1*betaim2/6;
301  break;
302  }
303  case(CENTER):
304  {
305  data__ deli, delim1;
306 
307  deli=x[2]-x[1];
308  delim1=x[1]-x[0];
309 
310  te=-phi3*deli*delim1/6;
311  break;
312  }
313  case(RIGHT):
314  {
315  data__ alphai, betai;
316 
317  alphai=x[1]-x[0];
318  betai=x[2]-x[0];
319 
320  te=phi3*alphai*betai/6;
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
int order(bool) const
Definition: d1o2.hpp:80
Definition: math.hpp:20
int truncation_error(data__ &te, const data__ &phi3, itx__ itx) const
Definition: d1o2.hpp:282
stencil st
Definition: d1o2.hpp:40
void set_stencil(const stencil &s)
Definition: d1o2.hpp:70
const stencil & get_stencil() const
Definition: d1o2.hpp:75
stencil
Definition: d1o2.hpp:30
Definition: d1o2.hpp:34
int coefficients(itc__ itc, const data__ &dx) const
Definition: d1o2.hpp:155
~d1o2()
Definition: d1o2.hpp:66
Definition: d1o2.hpp:32
data__ calculate_dot(itc__ a, itphi__ itphi) const
Definition: d1o2.hpp:43
int evaluate(data__ &d, itphi__ itphi, itx__ itx) const
Definition: d1o2.hpp:143
const int n_order
Definition: d1o2.hpp:39
Definition: d1o2.hpp:27
int coefficients(itc__ itc, itx__ itx) const
Definition: d1o2.hpp:193
d1o2()
Definition: d1o2.hpp:54
int evaluate(data__ &d, itphi__ itphi, const data__ &dx) const
Definition: d1o2.hpp:131
size_t number_nodes() const
Definition: d1o2.hpp:85
const size_t nnodes
Definition: d1o2.hpp:38
d1o2(const d1o2< data__ > &d)
Definition: d1o2.hpp:62
Definition: d1o2.hpp:33
d1o2(const stencil &s)
Definition: d1o2.hpp:58
int truncation_error(data__ &te, const data__ &phi3, const data__ &dx) const
Definition: d1o2.hpp:253
std::ptrdiff_t index(iti__ iti) const
Definition: d1o2.hpp:90