Chombo + EB + MF  3.2
SmoothAbsoluteValue.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _SMOOTHABSOLUTEVALUE_H_
12 #define _SMOOTHABSOLUTEVALUE_H_
13 
14 #include "MayDay.H"
15 #include "RealVect.H"
16 #include "IntVect.H"
17 #include "BaseIF.H"
18 
19 #include "NamespaceHeader.H"
20 
21 ///
22 /**
23  Functions to take a smooth absolute value of the difference between two functions
24  A_e(a,b) = |a-b|
25 
26  max(f, g) = (a + b + |a - b|)/2
27  We need to make abs(a-b) smooth.
28 
29  We take the convolution of the absolute value with cos^4(pi w/(2 delta))
30  */
32 {
33 public:
34 
35  ///
36  /**
37  Set a known f and g (from short list I have programmed) to check answer against
38  1 = double ramp with
39  ramp_normal_0 = 2. 1. 0.
40  ramp_normal_1 = -2. 1. 0.
41  ramp_point_0 = 0. 1. 0.
42  ramp_point_1 = 1. 1. 0.
43 
44  2 = double sphere with
45  sphere_center_0 = 0 0 0
46  sphere_center_1 = 1 0 0
47  sphere_radius_0 = 0.75
48  sphere_radius_1 = 0.75
49  */
50  static void setKnownFunction(int a_whichfunc)
51  {
52  s_knownFunc = a_whichfunc;
53  }
54 
55  ///
57  const BaseIF* a_g,
58  const Real & a_delta)
59  {
60  m_f = a_f;
61  m_g = a_g;
62  m_d = a_delta;
63  m_pi = 4.*atan(1.0);
64  }
65 
66  ///
68  {; }
69 
70  ///
71  virtual Real smoothAbsFMinusG(const IntVect& a_deriv,
72  const RealVect& a_point) const;
73 
74  ///
75  /**
76  returns -1 if w < -delta, 1 if w > delta, 0 otherwise
77  reduces to regular |f-g| unless case == 0
78  */
79  void getWCase(int & a_case,
80  Real & a_wval,
81  const RealVect & a_point)const;
82 
83 protected:
84 
85  ///
86  /**
87  Here is the logic of this stuff. We have three cases that reduce
88  to two cases. w = f(x) - g(x)
89  case 1: (w > delta): ---- whole integral is above zero
90  answer = abs(w)
91  case -1: (w < - delta): ---- whole integral is below zero
92  answer = abs(w)
93  case 0: (-delta <= w <= delta) --- have to split integral into above and below
94  answer = functionAem();
95  */
96 
97 
98  virtual Real valueAem(const RealVect& a_point) const;
99 
100  ///
101  virtual Real firstDerivAem(const IntVect& a_deriv,
102  const RealVect& a_point) const;
103 
104  ///
105  virtual Real secondDerivAem(const IntVect& a_deriv,
106  const RealVect& a_point) const;
107 
108  ///
109  virtual Real thirdDerivAem(const IntVect& a_deriv,
110  const RealVect& a_point) const;
111 
112  ///
113  virtual Real fourthDerivAem(const IntVect& a_deriv,
114  const RealVect& a_point) const;
115 
116  ///just checks nan
117  bool isBogus(const Real& a_number) const;
118 
119  ///if s_knownFunc is set, check against known answer
120  void checkAgainstKnown(const Real & a_myAnswer,
121  const IntVect & a_deriv,
122  const RealVect & a_point) const;
123 
124 
125  //the two implicit functions are owned by others
126  const BaseIF* m_f;
127  const BaseIF* m_g;
128  //delta = the smoothing length scale
130 
131  //pi, you know, pi = 4atan(1)
133  ///for debugging against known functions
134  static int s_knownFunc;
135 private:
136 
138  {
139  MayDay::Abort("SmoothAbsoluteValue doesn't allow copy construction");
140  }
141 
143  {
144  MayDay::Abort("SmoothAbsoluteValue uses strong construction");
145  }
146 
147  void operator=(const SmoothAbsoluteValue& a_inputIF)
148  {
149  MayDay::Abort("SmoothAbsoluteValue doesn't allow assignment");
150  }
151 };
152 ///
153 /**
154  offset sphere test case
155  The mac daddy of test cases
156  Test case where
157  f: (x- 0 )^2 + (y- 0 )^2 + (z- 0 )^2 - (1/2)^2;
158  g: (x-(1/2))^2 + (y-(1/2))^2 + (z-(1/2))^2 - (1/4)^2;
159  sphere_center_0 = 0 0 0
160  sphere_center_1 = 0.5 0.5 0.5
161  sphere_radius_0 = 0.5
162  sphere_radius_1 = 0.25
163 **/
165 {
166 public:
167  OffsetSphereExact(const Real & a_delta,
168  const Real & a_pi)
169  {
170  m_d = a_delta;
171  m_pi = a_pi;
172  }
173 
175  {
176  }
177 
178  ///this one calls the other functions
179  Real value(const IntVect& a_deriv,
180  const RealVect& a_point) const;
181 
182 
183  ///
184  Real valueAem(const RealVect& a_point) const;
185 
186 
187  ///
188  Real firstDerivAem(const IntVect& a_deriv,
189  const RealVect& a_point) const;
190 
191 
192  ///
193  Real secondDerivAem(const IntVect& a_deriv,
194  const RealVect& a_point) const;
195 
196 
197  ///
198  Real thirdDerivAem(const IntVect& a_deriv,
199  const RealVect& a_point) const;
200 
201 
202  ///
203  Real fourthDerivAem(const IntVect& a_deriv,
204  const RealVect& a_point) const;
205 private:
206  void getXYIXIY(int & a_ix, int& a_iy,
207  Real& a_x, Real& a_y,
208  const int & a_xderivVal,
209  const IntVect& a_deriv,
210  const RealVect& a_point) const;
211 
213  {
214  MayDay::Error("this uses strong construction");
215  }
216 
219 
220 };
221 
222 
223 ///
224 /**
225  Double sphere test case
226  Test case where
227 f(x,y,z) = x^2 + y^2 + z^2 - 0.75^2
228 g(x,y,z) = (x-1)^2 + y^2 + z^2 - 0.75^2
229 g(x,y,z) = (-2*(x-1) + (y - 1))
230 f(x,y.z) - g(x,y,z) = x^2 - (x-1)^2
231 sphere_center_0 = 0 0 0
232 sphere_center_1 = 1 0 0
233 sphere_radius_0 = 0.75
234 sphere_radius_1 = 0.75
235 **/
237 {
238 public:
239  DoubleSphereExact(const Real & a_delta,
240  const Real & a_pi)
241  {
242  m_d = a_delta;
243  m_pi = a_pi;
244  }
245 
247  {
248  }
249 
250  ///this one calls the other functions
251  Real value(const IntVect& a_deriv,
252  const RealVect& a_point) const;
253 
254 
255  ///
256  Real valueAem(const RealVect& a_point) const;
257 
258 
259  ///
260  Real firstDerivAem(const IntVect& a_deriv,
261  const RealVect& a_point) const;
262 
263 
264  ///
265  Real secondDerivAem(const IntVect& a_deriv,
266  const RealVect& a_point) const;
267 
268 
269  ///
270  Real thirdDerivAem(const IntVect& a_deriv,
271  const RealVect& a_point) const;
272 
273 
274  ///
275  Real fourthDerivAem(const IntVect& a_deriv,
276  const RealVect& a_point) const;
277 private:
279  {
280  MayDay::Error("this uses strong construction");
281  }
282 
285 
286 };
287 
288 ///
289 /**
290  Double sphere test case
291  Test case where
292  f(x,y) = 2*(x-0) + (y - 1)
293  g(x,y) =-2*(x-1) + (y - 1)
294  in the input file:
295  ramp_normal_0 = 2. 1. 0.
296  ramp_normal_1 = -2. 1. 0.
297  ramp_point_0 = 0. 1. 0.
298  ramp_point_1 = 1. 1. 0.
299 **/
301 {
302 public:
303  DoubleRampExact(const Real & a_delta,
304  const Real & a_pi)
305  {
306  m_d = a_delta;
307  m_pi = a_pi;
308  }
309 
311  {
312  }
313 
314  ///this one calls the other functions
315  Real value(const IntVect& a_deriv,
316  const RealVect& a_point) const;
317 
318 
319  ///
320  Real valueAem(const RealVect& a_point) const;
321 
322 
323  ///
324  Real firstDerivAem(const IntVect& a_deriv,
325  const RealVect& a_point) const;
326 
327 
328  ///
329  Real secondDerivAem(const IntVect& a_deriv,
330  const RealVect& a_point) const;
331 
332 
333  ///
334  Real thirdDerivAem(const IntVect& a_deriv,
335  const RealVect& a_point) const;
336 
337 
338  ///
339  Real fourthDerivAem(const IntVect& a_deriv,
340  const RealVect& a_point) const;
341 private:
342 
344  {
345  MayDay::Error("this uses strong construction");
346  }
347 
350 
351 };
352 #include "NamespaceFooter.H"
353 #endif
Real m_d
Definition: SmoothAbsoluteValue.H:348
Real m_pi
Definition: SmoothAbsoluteValue.H:218
OffsetSphereExact()
Definition: SmoothAbsoluteValue.H:212
Real m_d
Definition: SmoothAbsoluteValue.H:283
~OffsetSphereExact()
Definition: SmoothAbsoluteValue.H:174
void operator=(const SmoothAbsoluteValue &a_inputIF)
Definition: SmoothAbsoluteValue.H:147
DoubleSphereExact()
Definition: SmoothAbsoluteValue.H:278
OffsetSphereExact(const Real &a_delta, const Real &a_pi)
Definition: SmoothAbsoluteValue.H:167
virtual Real secondDerivAem(const IntVect &a_deriv, const RealVect &a_point) const
DoubleSphereExact(const Real &a_delta, const Real &a_pi)
Definition: SmoothAbsoluteValue.H:239
const BaseIF * m_g
Definition: SmoothAbsoluteValue.H:127
Real m_pi
Definition: SmoothAbsoluteValue.H:132
Definition: SmoothAbsoluteValue.H:236
virtual Real firstDerivAem(const IntVect &a_deriv, const RealVect &a_point) const
~DoubleRampExact()
Definition: SmoothAbsoluteValue.H:310
Real m_d
Definition: SmoothAbsoluteValue.H:129
const BaseIF * m_f
Definition: SmoothAbsoluteValue.H:126
bool isBogus(const Real &a_number) const
just checks nan
void checkAgainstKnown(const Real &a_myAnswer, const IntVect &a_deriv, const RealVect &a_point) const
if s_knownFunc is set, check against known answer
SmoothAbsoluteValue()
Definition: SmoothAbsoluteValue.H:142
Real m_d
Definition: SmoothAbsoluteValue.H:217
virtual Real smoothAbsFMinusG(const IntVect &a_deriv, const RealVect &a_point) const
Definition: BaseIF.H:32
~DoubleSphereExact()
Definition: SmoothAbsoluteValue.H:246
Definition: SmoothAbsoluteValue.H:164
double Real
Definition: REAL.H:33
virtual Real fourthDerivAem(const IntVect &a_deriv, const RealVect &a_point) const
void getWCase(int &a_case, Real &a_wval, const RealVect &a_point) const
Definition: SmoothAbsoluteValue.H:31
virtual Real thirdDerivAem(const IntVect &a_deriv, const RealVect &a_point) const
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
virtual Real valueAem(const RealVect &a_point) const
A Real vector in SpaceDim-dimensional space.
Definition: RealVect.H:41
SmoothAbsoluteValue(const BaseIF *a_f, const BaseIF *a_g, const Real &a_delta)
Definition: SmoothAbsoluteValue.H:56
virtual ~SmoothAbsoluteValue()
Definition: SmoothAbsoluteValue.H:67
Real m_pi
Definition: SmoothAbsoluteValue.H:284
SmoothAbsoluteValue(const SmoothAbsoluteValue &a_inputIF)
Definition: SmoothAbsoluteValue.H:137
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
Real m_pi
Definition: SmoothAbsoluteValue.H:349
static int s_knownFunc
for debugging against known functions
Definition: SmoothAbsoluteValue.H:134
DoubleRampExact()
Definition: SmoothAbsoluteValue.H:343
DoubleRampExact(const Real &a_delta, const Real &a_pi)
Definition: SmoothAbsoluteValue.H:303
static void setKnownFunction(int a_whichfunc)
Definition: SmoothAbsoluteValue.H:50
Definition: SmoothAbsoluteValue.H:300
static void Abort(const char *const a_msg=m_nullString)
Print out message to cerr and exit via abort() (if serial) or MPI_Abort() (if parallel).