OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
mathfem.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2013 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "mathfem.h"
36 #include "floatarray.h"
37 
38 namespace oofem {
39 // measure dependent constant
40 #define CUBIC_ZERO 1.0e-100
41 
42 
43 void cubic(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num)
44 //
45 // solves cubic equation for real roots
46 //
47 // input:
48 // a,b,c,d - coefficients of equation in form:
49 // ax^3 + bx^2 + cx + d = 0
50 //
51 // output:
52 // r1,r2,r3 - roots (only first num roots is valid)
53 // num - number of roots resolved
54 //
55 {
56  double aa, p, q, D, u, v, phi;
57  double help;
58 
59  double norm = 1e-6 * ( fabs(a) + fabs(b) + fabs(c) ) + CUBIC_ZERO;
60 
61  if ( fabs(a) <= norm ) {
62  if ( fabs(b) <= norm ) {
63  if ( fabs(c) <= norm ) {
64  if ( fabs(d) <= norm ) {
65  * r1 = 0.0;
66  * num = 1;
67  return;
68  } else {
69  * num = 0;
70  return;
71  }
72  } else {
73  * r1 = -d / c;
74  * num = 1;
75  return;
76  }
77  } else {
78  if ( ( D = c * c - 4.0 * b * d ) < 0.0 ) {
79  * num = 0;
80  return;
81  } else {
82  //*r1 = (-c + sqrt(D)) / 2.0 / b;
83  //*r2 = (-c - sqrt(D)) / 2.0 / b;
84  if ( fabs(c) < norm ) {
85  help = -d / b;
86  if ( help > 0. ) {
87  * r1 = sqrt(help);
88  * r2 = -sqrt(help);
89  * num = 2;
90  return;
91  } else {
92  * num = 0;
93  return;
94  }
95  } else {
96  help = -( c + sgn(c) * sqrt(D) ) / 2.0;
97  * r1 = help / b;
98  * r2 = d / help;
99  * num = 2;
100  }
101  }
102  }
103  } else {
104  aa = a;
105  a = b / ( aa * 3.0 );
106  b = c / aa;
107  c = d / aa;
108  p = b - a * a * 3.0;
109  q = 2.0 * a * a * a - a * b + c;
110  D = q * q / 4.0 + p * p * p / 27.0;
111  if ( fabs(D) < CUBIC_ZERO ) {
112  if ( fabs(p * q) < CUBIC_ZERO ) {
113  * r1 = 0.0 - a;
114  * r2 = * r1;
115  * r3 = * r1;
116  * num = 3;
117  } else {
118  * r2 = cbrt(q / 2.0) - a;
119  * r1 = -2.0 * * r2 - a;
120  * num = 2;
121  }
122  } else {
123  if ( D > 0.0 ) {
124  u = -q / 2.0 + sqrt(D);
125  v = -q - u;
126 
127  * r1 = u + v - a;
128  * r1 = cbrt(u) + cbrt(v) - a;
129  * num = 1;
130  } else {
131  p = sqrt(fabs(p) / 3.0);
132  help = ( -q / ( 2.0 * p * p * p ) );
133  if ( fabs(help) > 1.0 ) {
134  help = sgn(help); // prevent rounding errors
135  }
136 
137  phi = acos(help) / 3.0;
138  double cp = cos(phi);
139  double sp = sqrt(3.0) * sin(phi);
140  * r1 = 2 * p * cp - a;
141  * r2 = -p * ( cp + sp ) - a;
142  * r3 = -p * ( cp - sp ) - a;
143 
144  // I'm getting some pretty bad accuracy, a single iteration like this would help alot
145  //* r1 -= (d + c*(*r1) + b*(*r1)*(*r1) + a*(*r1)*(*r1)*(*r1))/(c + 2*b*(*r1) + 3*a*(*r1)*(*r1));
146  //* r2 -= (d + c*(*r2) + b*(*r2)*(*r2) + a*(*r2)*(*r2)*(*r2))/(c + 2*b*(*r2) + 3*a*(*r2)*(*r2));
147  //* r3 -= (d + c*(*r3) + b*(*r3)*(*r3) + a*(*r3)*(*r3)*(*r3))/(c + 2*b*(*r3) + 3*a*(*r3)*(*r3));
148  * num = 3;
149  }
150  }
151  }
152 }
153 
154 
155 void cubic3r(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num)
156 //
157 // solves cubic equation for real roots
158 //
159 // input:
160 // a,b,c,d - coefficients of equation in form:
161 // ax^3 + bx^2 + cx + d = 0
162 //
163 // output:
164 // r1,r2,r3 - roots (only first num roots is valid)
165 // num - number of roots resolved
166 //
167 {
168  double aa, r, q, p, D, phi;
169  double help;
170 
171  if ( fabs(a) < CUBIC_ZERO ) {
172  if ( fabs(b) < CUBIC_ZERO ) {
173  if ( fabs(c) < CUBIC_ZERO ) {
174  * num = 0;
175  return;
176  } else {
177  * r1 = -d / c;
178  * num = 1;
179  return;
180  }
181  } else { //fabs(b) < CUBIC_ZERO
182  if ( ( D = c * c - 4.0 * b * d ) < 0.0 ) {
183  * num = 0;
184  return;
185  } else {
186  //*r1 = (-c + sqrt(D)) / 2.0 / b;
187  //*r2 = (-c - sqrt(D)) / 2.0 / b;
188  if ( fabs(c) < CUBIC_ZERO ) {
189  help = -d / b;
190  if ( help >= 0. ) {
191  * r1 = sqrt(help);
192  * r2 = -sqrt(help);
193  * num = 2;
194  return;
195  } else {
196  * num = 0;
197  return;
198  }
199  } else {
200  help = -( c + sgn(c) * sqrt(D) ) / 2.0;
201  * r1 = help / b;
202  * r2 = d / help;
203  * num = 2;
204  return;
205  }
206  }
207  }
208  } else {
209  aa = a;
210  a = b / aa;
211  b = c / aa;
212  c = d / aa;
213 
214  q = ( a * a - 3.0 * b ) / 9.0;
215  r = ( 2.0 * a * a * a - 9.0 * a * b + 27.0 * c ) / 54.0;
216 
217  //Hydrostatic case, in such a case q=r=0
218  if ( (fabs(q) < CUBIC_ZERO) && ( fabs(r) < CUBIC_ZERO ) ) {
219  *r1 = - a / 3.0;
220  *r2 = *r1;
221  *r3 = *r1;
222  *num = 3;
223  return;
224  }
225 
226  //D = q * q * q - r * r;
227  // if (D > 0.) {
228  // three real roots
229  help = r / sqrt(q * q * q);
230  if ( fabs(help) > 1.0 ) {
231  help = sgn(help); // prevent rounding errors
232  }
233 
234  phi = acos(help);
235  p = sqrt(q);
236 
237  * r1 = -2.0 *p *cos(phi / 3.0) - a / 3.0;
238  * r2 = -2.0 *p *cos( ( phi + 2.0 *M_PI ) / 3.0 ) - a / 3.0;
239  * r3 = -2.0 *p *cos( ( phi - 2.0 *M_PI ) / 3.0 ) - a / 3.0;
240  * num = 3;
241  /* } else {
242  *
243  * help = fabs(r) + sqrt(-D);
244  * A = -sgn(r)*cbrt(D);
245  * if (fabs(A) > CUBIC_ZERO)
246  * B = q/A;
247  * else
248  * B = 0.0;
249  *
250  *****r1 = (A+B) - a/3.0;
251  *****num = 1;
252  * }
253  */
254  return;
255  }
256 }
257 
258 
259 
260 int iperm(int val, int rank)
261 //
262 // returns iperm of val, in specific rank
263 //
264 // rank 3 : val is valid from 1,2,3 and returns for val 1,2,3 2,3,1
265 // rank 2 : -||-
266 {
267  return ( val + 1 > rank ? 1 : val + 1 );
268 }
269 
270 
271 void ls2fit(const FloatArray &x, const FloatArray &y, FloatArray &a)
272 {
273  int n = x.giveSize();
274  a.resize(3);
275  if ( n > 2 ) {
276  // Least square fitting.
277  double f1 = 0, f2 = 0, f3 = 0;
278  double sx = 0, sx2 = 0, sx3 = 0, sx4 = 0;
279  double detCi, Ci11, Ci12, Ci13, Ci22, Ci23, Ci33;
280 
281  double yi, xi, xi2, xi3, xi4;
282  for ( int i = 0; i < n; ++i ) {
283  // Calculate foot points in local coord.sys.
284  xi = x(i);
285  yi = y(i);
286  xi2 = xi * xi;
287  xi3 = xi * xi2;
288  xi4 = xi * xi3;
289 
290  // Construct the coefficient matrix.
291  sx += xi;
292  sx2 += xi2;
293  sx3 += xi3;
294  sx4 += xi4;
295 
296  // And the RHS.
297  f1 += yi;
298  f2 += xi * yi;
299  f3 += xi2 * yi;
300  }
301 
302  // Explicit inverse (avoids numerical problems)
303  Ci11 = sx2 * sx4 - sx3 * sx3;
304  Ci12 = sx2 * sx3 - sx * sx4;
305  Ci13 = sx * sx3 - sx2 * sx2;
306  Ci22 = n * sx4 - sx2 * sx2;
307  Ci23 = sx * sx2 - n * sx3;
308  Ci33 = n * sx2 - sx * sx;
309  detCi = 1 / ( n * Ci11 + sx * Ci12 + sx2 * Ci13 );
310  a(0) = ( Ci11 * f1 + Ci12 * f2 + Ci13 * f3 ) * detCi;
311  a(1) = ( Ci12 * f1 + Ci22 * f2 + Ci23 * f3 ) * detCi;
312  a(2) = ( Ci13 * f1 + Ci23 * f2 + Ci33 * f3 ) * detCi;
313  } else if ( n == 2 ) {
314  a(2) = 0;
315  a(1) = ( y(1) - y(0) ) / ( x(1) - x(0) );
316  a(0) = y(0) - a(1) * x(0);
317  } else if ( n == 1 ) {
318  a(0) = y(0);
319  a(1) = 0;
320  a(2) = 0;
321  } else {
322  a.zero();
323  }
324 }
325 
326 double signum(double i) {
327  if ( i < 0 ) {
328  return -1;
329  } else if ( i > 0 ) {
330  return 1;
331  } else {
332  return 0;
333  }
334 }
335 } // end namespace oofem
double sgn(double i)
Returns the signum of given value (if value is < 0 returns -1, otherwise returns 1) ...
Definition: mathfem.h:91
#define CUBIC_ZERO
Definition: mathfem.C:40
void cubic(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num)
Solves cubic equation for real roots.
Definition: mathfem.C:43
#define M_PI
Definition: mathfem.h:52
void cubic3r(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num)
Solves cubic equation for real roots, assuming that if cubic polynomial given then the only possibili...
Definition: mathfem.C:155
int iperm(int val, int rank)
Returns iperm of val, in specific rank.
Definition: mathfem.C:260
void ls2fit(const FloatArray &x, const FloatArray &y, FloatArray &a)
Least-square fit of 2nd degree polynomial .
Definition: mathfem.C:271
Class representing vector of real numbers.
Definition: floatarray.h:82
double cbrt(double x)
Returns the cubic root of x.
Definition: mathfem.h:109
double norm(const FloatArray &x)
Definition: floatarray.C:985
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
double signum(double i)
Returns the signum of given value (i = 0 returns 0, i < 0 returns -1, i > 0 returns 1) ...
Definition: mathfem.C:326
int giveSize() const
Returns the size of receiver.
Definition: floatarray.h:218
the oofem namespace is to define a context or scope in which all oofem names are defined.
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631

This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:29 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011