MathUtils.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7
8 #include "MathUtils.h"
9 #include "libmesh/utility.h"
10
11 namespace MathUtils
12 {
13
14 Real
15 poly1Log(Real x, Real tol, int deriv)
16 {
17  Real c1 = 1.0 / tol;
18  Real c2 = std::log(tol) - 1.0;
19
20  Real value = 0.0;
21
22  if (deriv == 0)
23  {
24  if (x < tol)
25  value = c1 * x + c2;
26  else
27  value = std::log(x);
28  }
29  else if (deriv == 1)
30  {
31  if (x < tol)
32  value = c1;
33  else
34  value = 1.0 / x;
35  }
36  else if (deriv == 2)
37  {
38  if (x < tol)
39  value = 0.0;
40  else
41  value = -1.0 / (x * x);
42  }
43  else if (deriv == 3)
44  {
45  if (x < tol)
46  value = 0.0;
47  else
48  value = 2.0 / (x * x * x);
49  }
50
51  return value;
52 }
53
54 Real
55 poly2Log(Real x, Real tol, int deriv)
56 {
57  Real c1 = -0.5 / (tol * tol);
58  Real c2 = 2.0 / tol;
59  Real c3 = std::log(tol) - 3.0 / 2.0;
60
61  Real value = 0.0;
62
63  if (deriv == 0)
64  {
65  if (x < tol)
66  value = c1 * x * x + c2 * x + c3;
67  else
68  value = std::log(x);
69  }
70  else if (deriv == 1)
71  {
72  if (x < tol)
73  value = 2.0 * c1 * x + c2;
74  else
75  value = 1.0 / x;
76  }
77  else if (deriv == 2)
78  {
79  if (x < tol)
80  value = 2.0 * c1;
81  else
82  value = -1.0 / (x * x);
83  }
84  else if (deriv == 3)
85  {
86  if (x < tol)
87  value = 0.0;
88  else
89  value = 2.0 / (x * x * x);
90  }
91
92  return value;
93 }
94
95 Real
96 poly3Log(Real x, Real tol, int order)
97 {
98  Real c1 = 1.0 / (3.0 * tol * tol * tol);
99  Real c2 = -3.0 / (2.0 * tol * tol);
100  Real c3 = 3.0 / tol;
101  Real c4 = std::log(tol) - 11.0 / 6.0;
102
103  Real value = 0.0;
104
105  if (order == 0)
106  {
107  if (x < tol)
108  value = c1 * x * x * x + c2 * x * x + c3 * x + c4;
109  else
110  value = std::log(x);
111  }
112  else if (order == 1)
113  {
114  if (x < tol)
115  value = 3.0 * c1 * x * x + 2.0 * c2 * x + c3;
116  else
117  value = 1.0 / x;
118  }
119  else if (order == 2)
120  {
121  if (x < tol)
122  value = 6.0 * c1 * x + 2.0 * c2;
123  else
124  value = -1.0 / (x * x);
125  }
126  else if (order == 3)
127  {
128  if (x < tol)
129  value = 6.0 * c1;
130  else
131  value = 2.0 / (x * x * x);
132  }
133  return value;
134 }
135
136 Real
137 poly4Log(Real x, Real tol, int order)
138 {
139  Real value = 0.0;
140
141  if (order == 0)
142  {
143  if (x < tol)
144  value = std::log(tol) + (x - tol) / tol - (x - tol) * (x - tol) / (2.0 * tol * tol) +
145  (x - tol) * (x - tol) * (x - tol) / (3.0 * tol * tol * tol) -
146  (x - tol) * (x - tol) * (x - tol) * (x - tol) / (4.0 * tol * tol * tol * tol) +
147  (x - tol) * (x - tol) * (x - tol) * (x - tol) * (x - tol) /
148  (5.0 * tol * tol * tol * tol * tol) -
149  (x - tol) * (x - tol) * (x - tol) * (x - tol) * (x - tol) * (x - tol) /
150  (6.0 * tol * tol * tol * tol * tol * tol);
151  else
152  value = std::log(x);
153  }
154  else if (order == 1)
155  {
156  if (x < tol)
157  value = 1.0 / tol - 2.0 * (x - tol) / (2.0 * tol * tol) +
158  3.0 * (x - tol) * (x - tol) / (3.0 * tol * tol * tol) -
159  4.0 * (x - tol) * (x - tol) * (x - tol) / (4.0 * tol * tol * tol * tol) +
160  5.0 * (x - tol) * (x - tol) * (x - tol) * (x - tol) /
161  (5.0 * tol * tol * tol * tol * tol) -
162  6.0 * (x - tol) * (x - tol) * (x - tol) * (x - tol) * (x - tol) /
163  (6.0 * tol * tol * tol * tol * tol * tol);
164  else
165  value = 1.0 / x;
166  }
167  else if (order == 2)
168  {
169  if (x < tol)
170  value = -2.0 * 1.0 / (2.0 * tol * tol) + 3.0 * 2.0 * (x - tol) / (3.0 * tol * tol * tol) -
171  4.0 * 3.0 * (x - tol) * (x - tol) / (4.0 * tol * tol * tol * tol) +
172  5.0 * 4.0 * (x - tol) * (x - tol) * (x - tol) / (5.0 * tol * tol * tol * tol * tol) -
173  6.0 * 5.0 * (x - tol) * (x - tol) * (x - tol) * (x - tol) /
174  (6.0 * tol * tol * tol * tol * tol * tol);
175  else
176  value = -1.0 / (x * x);
177  }
178  else if (order == 3)
179  {
180  if (x < tol)
181  value = 3.0 * 2.0 * 1.0 / (3.0 * tol * tol * tol) -
182  4.0 * 3.0 * 2.0 * (x - tol) / (4.0 * tol * tol * tol * tol) +
183  5.0 * 4.0 * 3.0 * (x - tol) * (x - tol) / (5.0 * tol * tol * tol * tol * tol) -
184  6.0 * 5.0 * 4.0 * (x - tol) * (x - tol) * (x - tol) /
185  (6.0 * tol * tol * tol * tol * tol * tol);
186  else
187  value = 2.0 / (x * x * x);
188  }
189  return value;
190 }
191
193 Real
194 taylorLog(Real x)
195 {
196  Real y = (x - 1.0) / (x + 1.0);
197  Real val = 1.0;
198  for (unsigned int i = 0; i < 5; ++i)
199  {
200  Real exponent = i + 2.0;
201  val += 1.0 / (2.0 * (i + 1.0) + 1.0) * std::pow(y, exponent);
202  }
203
204  return val * 2.0 * y;
205 }
206
207 // fast positive integer powers
208 Real
209 pow(Real x, unsigned int e)
210 {
211  Real result = 1.0;
212  while (e)
213  {
214  // if bit 0 is set multiply the current power of two factor of the exponent
215  if (e & 1)
216  result *= x;
217
218  // x is incrementally set to consecutive powers of powers of two
219  x *= x;
220
221  // bit shift the exponent down
222  e >>= 1;
223  }
224
225  return result;
226 }
227
228 } // namespace MathUtils
Real poly1Log(Real x, Real tol, int deriv)
Definition: MathUtils.C:15
Real poly3Log(Real x, Real tol, int order)
Definition: MathUtils.C:96
Real pow(Real x, unsigned int e)
Definition: MathUtils.C:209
Real poly2Log(Real x, Real tol, int deriv)
Definition: MathUtils.C:55
Real poly4Log(Real x, Real tol, int order)
Definition: MathUtils.C:137
Real taylorLog(Real x)
Definition: MathUtils.C:194
static const double tol
Definition: XFEMFuncs.h:26
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)