Overture  Version 25
OGPulseFunction.h
Go to the documentation of this file.
1 #ifndef OG_PULSE_FUNCTION_H
2 #define OG_PULSE_FUNCTION_H
3 
4 #include "OGFunction.h"
5 
6 //===========================================================================================
7 //
8 // Define "pulse" functions that are useful for testing AMR applications.
9 //
10 // Define a Function and it derivatives
11 // This function can be used to define the "exact solution" for
12 // an Overlapping Grid Appliciation (aka. TwilightZone Flow)
13 //
14 //
15 //===========================================================================================
16 class OGPulseFunction : public OGFunction
17 {
18 
19  public:
20 
21 
22  //
23  // Create a polynomial with the given degree, number Of Space Dimensions and for
24  // a maximum number of components
25  //
26  OGPulseFunction(int numberOfDimensions_ = 2,
27  int numberOfComponents_ =1,
28  real a0_ =1.,
29  real a1_ =5.,
30  real c0_ =0.,
31  real c1_ =0.,
32  real c2_ =0.,
33  real v0_ =1.,
34  real v1_ =1.,
35  real v2_ =1.,
36  real p_ =2.);
37 
38  OGPulseFunction(const OGPulseFunction & ogp );
39  OGPulseFunction & operator=(const OGPulseFunction & ogp );
40  ~OGPulseFunction(){}
41 
42  int setRadius( real radius );
43  int setShape( real p_ );
44  int setCentre( real c0_ =0.,
45  real c1_ =0.,
46  real c2_ =0. );
47  int setVelocity( real v0_ =1.,
48  real v1_ =1.,
49  real v2_ =1. );
50 
51  void setParameters(int numberOfDimensions_ = 2,
52  int numberOfComponents_ =1,
53  real a0_ =1.,
54  real a1_ =5.,
55  real c0_ =0.,
56  real c1_ =0.,
57  real c2_ =0.,
58  real v0_ =1.,
59  real v1_ =1.,
60  real v2_ =1.,
61  real p_ =2.);
62 
63  // Here are the member functions that you must define
64  // u(x,y,z,n) : function value at x,y,z for component n
65  // ux(x,y,z,n) : partial derivative of u with respect to x
66  // ...etc...
67 
68  inline real u(const real x, const real y, const real z, const int n, const real t=0.)
69  {
70  real b0=c0+v0*t;
71  real b1=c1+v1*t;
72  real b2=c2+v2*t;
73  real r;
74  if( numberOfDimensions==2 )
75  r= SQR(x-b0)+SQR(y-b1);
76  else if( numberOfDimensions==3 )
77  r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
78  else
79  r=SQR(x-b0);
80  real rp=pow(r,p);
81 
82  return a0*exp( - a1*rp );
83  }
84  inline real ut(const real x, const real y, const real z, const int n=0, const real t=0. )
85  {
86  real b0=c0+v0*t;
87  real b1=c1+v1*t;
88  real b2=c2+v2*t;
89  real r;
90  if( numberOfDimensions==2 )
91  r= SQR(x-b0)+SQR(y-b1);
92  else if( numberOfDimensions==3 )
93  r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
94  else
95  r=SQR(x-b0);
96  real rp=pow(r,p);
97  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
98  if( numberOfDimensions==2 )
99  return a0*exp( - a1*rp )*( 2.*a1*p*rp/r*( (x-b0)*v0 + (y-b1)*v1 ));
100  else if( numberOfDimensions==3 )
101  return a0*exp( - a1*rp )*( 2.*a1*p*rp/r*( (x-b0)*v0 + (y-b1)*v1 +(z-b2)*v2 ));
102  else
103  return a0*exp( - a1*rp )*( 2.*a1*p*rp/r*( (x-b0)*v0 ) );
104  }
105  inline real ux(const real x, const real y, const real z, const int n=0, const real t=0.)
106  {
107  real b0=c0+v0*t;
108  real b1=c1+v1*t;
109  real b2=c2+v2*t;
110 
111  real r;
112  if( numberOfDimensions==2 )
113  r= SQR(x-b0)+SQR(y-b1);
114  else if( numberOfDimensions==3 )
115  r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
116  else
117  r=SQR(x-b0);
118  real rp=pow(r,p);
119  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
120  return a0*exp(-a1*rp)*( -2.*a1*p*rp/r*( (x-b0)) );
121  }
122  inline real uy(const real x, const real y, const real z, const int n=0, const real t=0.)
123  {
124  if( numberOfDimensions==1 )
125  return 0.;
126  real b0=c0+v0*t;
127  real b1=c1+v1*t;
128  real b2=c2+v2*t;
129 
130  real r;
131  if( numberOfDimensions==2 )
132  r= SQR(x-b0)+SQR(y-b1);
133  else if( numberOfDimensions==3 )
134  r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
135  else
136  r=SQR(x-b0);
137  real rp=pow(r,p);
138  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
139  return a0*exp(-a1*rp)*( -2.*a1*p*rp/r*( (y-b1)) );
140  }
141  inline real uxx(const real x, const real y, const real z, const int n=0, const real t=0.)
142  {
143  real b0=c0+v0*t;
144  real b1=c1+v1*t;
145  real b2=c2+v2*t;
146 
147  real r;
148  if( numberOfDimensions==2 )
149  r= SQR(x-b0)+SQR(y-b1);
150  else if( numberOfDimensions==3 )
151  r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
152  else
153  r=SQR(x-b0);
154 
155  real rp=pow(r,p);
156  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
157  real g = -2.*a1*p*rp/r*(x-b0);
158  return a0*exp(-a1*rp)*(g*g-2.*a1*p*rp/r*(1.+2.*(p-1.)/r*SQR(x-b0)));
159  }
160 
161  inline real uxy(const real x, const real y, const real z, const int n=0, const real t=0.)
162  {
163  if( numberOfDimensions==1 )
164  return 0.;
165  real b0=c0+v0*t;
166  real b1=c1+v1*t;
167  real b2=c2+v2*t;
168 
169  real r;
170  if( numberOfDimensions==2 )
171  r= SQR(x-b0)+SQR(y-b1);
172  else if( numberOfDimensions==3 )
173  r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
174  else
175  r=SQR(x-b0);
176 
177  real rp=pow(r,p);
178  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
179  real gx = -2.*a1*p*rp/r*(x-b0);
180  real gy = -2.*a1*p*rp/r*(y-b1);
181  // N.B. write this as rp/r/r not rp/(r*r) to avoid nan's
182  return a0*exp(-a1*rp)*(gx*gy-4.*a1*p*(p-1.)*((rp/r)/r)*(x-b0)*(y-b1));
183  }
184  inline real uyy(const real x, const real y, const real z, const int n=0, const real t=0.)
185  {
186  if( numberOfDimensions==1 )
187  return 0.;
188  real b0=c0+v0*t;
189  real b1=c1+v1*t;
190  real b2=c2+v2*t;
191 
192  real r;
193  if( numberOfDimensions==2 )
194  r= SQR(x-b0)+SQR(y-b1);
195  else if( numberOfDimensions==3 )
196  r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
197  else
198  r=SQR(x-b0);
199 
200  real rp=pow(r,p);
201  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
202  real g = -2.*a1*p*rp/r*(y-b1);
203  return a0*exp(-a1*rp)*(g*g-2.*a1*p*rp/r*(1.+2.*(p-1.)/r*SQR(y-b1)));
204  }
205  inline real uz(const real x, const real y, const real z, const int n=0, const real t=0.)
206  {
207  if( numberOfDimensions<=2 )
208  return 0.;
209  real b0=c0+v0*t;
210  real b1=c1+v1*t;
211  real b2=c2+v2*t;
212 
213  real r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
214  real rp=pow(r,p);
215  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
216 
217  return a0*exp(-a1*rp)*( -2.*a1*p*rp/r*( (z-b2)) );
218  }
219  inline real uxz(const real x, const real y, const real z, const int n=0, const real t=0.)
220  {
221  if( numberOfDimensions<=2 )
222  return 0.;
223  real b0=c0+v0*t;
224  real b1=c1+v1*t;
225  real b2=c2+v2*t;
226 
227  real r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
228 
229  real rp=pow(r,p);
230  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
231  real gx = -2.*a1*p*rp/r*(x-b0);
232  real gz = -2.*a1*p*rp/r*(z-b2);
233  return a0*exp(-a1*rp)*(gx*gz-4.*a1*p*(p-1.)*rp/(r*r)*(x-b0)*(z-b2));
234  }
235  inline real uyz(const real x, const real y, const real z, const int n=0, const real t=0.)
236  {
237  if( numberOfDimensions<=2 )
238  return 0.;
239  real b0=c0+v0*t;
240  real b1=c1+v1*t;
241  real b2=c2+v2*t;
242 
243  real r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
244 
245  real rp=pow(r,p);
246  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
247  real gy = -2.*a1*p*rp/r*(y-b1);
248  real gz = -2.*a1*p*rp/r*(z-b2);
249  return a0*exp(-a1*rp)*(gy*gz-4.*a1*p*(p-1.)*rp/(r*r)*(y-b1)*(z-b2));
250  }
251  inline real uzz(const real x, const real y, const real z, const int n=0, const real t=0.)
252  {
253  if( numberOfDimensions<=2 )
254  return 0.;
255  real b0=c0+v0*t;
256  real b1=c1+v1*t;
257  real b2=c2+v2*t;
258 
259  real r= SQR(x-b0)+SQR(y-b1)+SQR(z-b2);
260 
261  real rp=pow(r,p);
262  r=max(r,REAL_MIN*100.); // avoid division by zero, note that rp/r*(x-b0) == 0 if r==0, p>.5
263  real g = -2.*a1*p*rp/r*(z-b2);
264  return a0*exp(-a1*rp)*(g*g-2.*a1*p*rp/r*(1.+2.*(p-1.)/r*SQR(z-b2)));
265  }
266 
267  // ========= Here are versions with a new naming convention ===========
268 
269  // (default arguments not allowed on operators)
270  virtual real operator()(const real , const real , const real);
271  virtual real operator()(const real , const real , const real , const int);
272  virtual real operator()(const real , const real , const real , const int, const real);
273 
274  virtual real t(const real , const real , const real , const int=0, const real=0.);
275  virtual real x(const real , const real , const real , const int=0, const real=0.);
276  virtual real y(const real , const real , const real , const int=0, const real=0.);
277  virtual real xx(const real , const real , const real , const int=0, const real=0.);
278  virtual real xy(const real , const real , const real , const int=0, const real=0.);
279  virtual real yy(const real , const real , const real , const int=0, const real=0.);
280  virtual real z(const real , const real , const real , const int=0, const real=0.);
281  virtual real xz(const real , const real , const real , const int=0, const real=0.);
282  virtual real yz(const real , const real , const real , const int=0, const real=0.);
283  virtual real zz(const real , const real , const real , const int=0, const real=0.);
284  virtual real xxx(const real , const real , const real , const int=0, const real=0.);
285  virtual real xxxx(const real , const real , const real , const int=0, const real=0.);
286  // obtain a general derivative
287  virtual real gd(const int & ntd, const int & nxd, const int & nyd, const int & nzd,
288  const real , const real , const real , const int=0, const real=0. );
289 
290  virtual RealDistributedArray operator()(const MappedGrid & c, const Index & I1, const Index & I2,
291  const Index & I3, const Index & N);
292  virtual RealDistributedArray operator()(const MappedGrid & c, const Index & I1, const Index & I2,
293  const Index & I3, const Index & N, const real t,
294  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
295 
296  virtual RealDistributedArray t(const MappedGrid & c, const Index & I1, const Index & I2,
297  const Index & I3, const Index & N, const real t=0.,
298  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
299  virtual RealDistributedArray x(const MappedGrid & c, const Index & I1, const Index & I2,
300  const Index & I3, const Index & N, const real t=0.,
301  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
302  virtual RealDistributedArray y(const MappedGrid & c, const Index & I1, const Index & I2,
303  const Index & I3, const Index & N, const real t=0.,
304  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
305  virtual RealDistributedArray z(const MappedGrid & c, const Index & I1, const Index & I2,
306  const Index & I3, const Index & N, const real t=0.,
307  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
308  virtual RealDistributedArray xx(const MappedGrid & c, const Index & I1, const Index & I2,
309  const Index & I3, const Index & N, const real t=0.,
310  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
311  virtual RealDistributedArray yy(const MappedGrid & c, const Index & I1, const Index & I2,
312  const Index & I3, const Index & N, const real t=0.,
313  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
314  virtual RealDistributedArray zz(const MappedGrid & c, const Index & I1, const Index & I2,
315  const Index & I3, const Index & N, const real t=0.,
316  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
317  virtual RealDistributedArray xy(const MappedGrid & c, const Index & I1, const Index & I2,
318  const Index & I3, const Index & N, const real t=0.,
319  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
320  virtual RealDistributedArray xz(const MappedGrid & c, const Index & I1, const Index & I2,
321  const Index & I3, const Index & N, const real t=0.,
322  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
323  virtual RealDistributedArray yz(const MappedGrid & c, const Index & I1, const Index & I2,
324  const Index & I3, const Index & N, const real t=0.,
325  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
326  virtual RealDistributedArray laplacian(const MappedGrid & c, const Index & I1, const Index & I2,
327  const Index & I3, const Index & N, const real t=0.,
328  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
329  virtual RealDistributedArray xxx(const MappedGrid & c, const Index & I1, const Index & I2,
330  const Index & I3, const Index & N, const real t=0.,
331  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
332  virtual RealDistributedArray xxxx(const MappedGrid & c, const Index & I1, const Index & I2,
333  const Index & I3, const Index & N, const real t=0.,
334  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering) ;
335  virtual RealDistributedArray gd(const int & ntd, const int & nxd, const int & nyd, const int & nzd,
336  const MappedGrid & c, const Index & I1, const Index & I2,
337  const Index & I3, const Index & N, const real t=0.,
338  const GridFunctionParameters::GridFunctionType & centering
339  =GridFunctionParameters::defaultCentering);
340 
341  // obtain a general derivative using serial arrays
342  virtual realSerialArray& gd( realSerialArray & result, // put result here
343  const realSerialArray & x, // coordinates to use if isRectangular==true
344  const int numberOfDimensions,
345  const bool isRectangular,
346  const int & ntd, const int & nxd, const int & nyd, const int & nzd,
347  const Index & I1, const Index & I2,
348  const Index & I3, const Index & N,
349  const real t=0., int option =0 );
350 
351  realCompositeGridFunction operator()(CompositeGrid & cg);
352  realCompositeGridFunction operator()(CompositeGrid & cg, const Index & N=nullIndex);
353  realCompositeGridFunction operator()(CompositeGrid & cg, const Index & N, const real t,
354  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
355  realCompositeGridFunction t(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
356  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
357  realCompositeGridFunction x(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
358  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
359  realCompositeGridFunction y(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
360  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
361  realCompositeGridFunction z(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
362  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
363  realCompositeGridFunction xx(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
364  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
365  realCompositeGridFunction xy(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
366  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
367  realCompositeGridFunction xz(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
368  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
369  realCompositeGridFunction yy(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
370  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
371  realCompositeGridFunction yz(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
372  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
373  realCompositeGridFunction zz(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
374  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
375  realCompositeGridFunction laplacian(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
376  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
377  realCompositeGridFunction xxx(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
378  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
379  realCompositeGridFunction xxxx(CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
380  const GridFunctionParameters::GridFunctionType & centering=GridFunctionParameters::defaultCentering);
381 
382  realCompositeGridFunction gd(const int & ntd, const int & nxd, const int & nyd, const int & nzd,
383  CompositeGrid & cg, const Index & N=nullIndex, const real t=0.,
384  const GridFunctionParameters::GridFunctionType & centering
385  =GridFunctionParameters::defaultCentering);
386 
387  private:
388  int numberOfComponents, numberOfDimensions;
389  real a0,a1,c0,c1,c2,c3,v0,v1,v2,p;
390 
391  void checkArguments(const Index & N);
392 
393 };
394 
395 
396 #undef TIME
397 
398 #endif
Generated on Fri Jan 4 2013 10:17:57 for Overture by   doxygen 1.8.3