OOFEM 3.0
Loading...
Searching...
No Matches
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 - 2025 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
38namespace oofem {
39// measure dependent constant
40#define CUBIC_ZERO 1.0e-100
41
42
43void 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
155void 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 * r2 = -2.0 *p *cos( ( phi + 2.0 *M_PI ) / 3.0 ) - a / 3.0;
240 * r3 = -2.0 *p *cos( ( phi - 2.0 *M_PI ) / 3.0 ) - a / 3.0;
241 * num = 3;
242 /* } else {
243 *
244 * help = fabs(r) + sqrt(-D);
245 * A = -sgn(r)*cbrt(D);
246 * if (fabs(A) > CUBIC_ZERO)
247 * B = q/A;
248 * else
249 * B = 0.0;
250 *
251 *****r1 = (A+B) - a/3.0;
252 *****num = 1;
253 * }
254 */
255 return;
256 }
257}
258
259
260
261int iperm(int val, int rank)
262//
263// returns iperm of val, in specific rank
264//
265// rank 3 : val is valid from 1,2,3 and returns for val 1,2,3 2,3,1
266// rank 2 : -||-
267{
268 return ( val + 1 > rank ? 1 : val + 1 );
269}
270
271
272void ls2fit(const FloatArray &x, const FloatArray &y, FloatArray &a)
273{
274 int n = x.giveSize();
275 a.resize(3);
276 if ( n > 2 ) {
277 // Least square fitting.
278 double f1 = 0, f2 = 0, f3 = 0;
279 double sx = 0, sx2 = 0, sx3 = 0, sx4 = 0;
280 double detCi, Ci11, Ci12, Ci13, Ci22, Ci23, Ci33;
281
282 double yi, xi, xi2, xi3, xi4;
283 for ( int i = 0; i < n; ++i ) {
284 // Calculate foot points in local coord.sys.
285 xi = x(i);
286 yi = y(i);
287 xi2 = xi * xi;
288 xi3 = xi * xi2;
289 xi4 = xi * xi3;
290
291 // Construct the coefficient matrix.
292 sx += xi;
293 sx2 += xi2;
294 sx3 += xi3;
295 sx4 += xi4;
296
297 // And the RHS.
298 f1 += yi;
299 f2 += xi * yi;
300 f3 += xi2 * yi;
301 }
302
303 // Explicit inverse (avoids numerical problems)
304 Ci11 = sx2 * sx4 - sx3 * sx3;
305 Ci12 = sx2 * sx3 - sx * sx4;
306 Ci13 = sx * sx3 - sx2 * sx2;
307 Ci22 = n * sx4 - sx2 * sx2;
308 Ci23 = sx * sx2 - n * sx3;
309 Ci33 = n * sx2 - sx * sx;
310 detCi = 1 / ( n * Ci11 + sx * Ci12 + sx2 * Ci13 );
311 a(0) = ( Ci11 * f1 + Ci12 * f2 + Ci13 * f3 ) * detCi;
312 a(1) = ( Ci12 * f1 + Ci22 * f2 + Ci23 * f3 ) * detCi;
313 a(2) = ( Ci13 * f1 + Ci23 * f2 + Ci33 * f3 ) * detCi;
314 } else if ( n == 2 ) {
315 a(2) = 0;
316 a(1) = ( y(1) - y(0) ) / ( x(1) - x(0) );
317 a(0) = y(0) - a(1) * x(0);
318 } else if ( n == 1 ) {
319 a(0) = y(0);
320 a(1) = 0;
321 a(2) = 0;
322 } else {
323 a.zero();
324 }
325}
326
327double signum(double i) {
328 if ( i < 0 ) {
329 return -1;
330 } else if ( i > 0 ) {
331 return 1;
332 } else {
333 return 0;
334 }
335}
336} // end namespace oofem
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::giveSize
Index giveSize() const
Returns the size of receiver.
Definition floatarray.h:261
oofem::FloatArray::zero
void zero()
Zeroes all coefficients of receiver.
Definition floatarray.C:683
CUBIC_ZERO
#define CUBIC_ZERO
Definition mathfem.C:40
M_PI
#define M_PI
Definition mathfem.h:52
oofem::norm
double norm(const FloatArray &x)
Definition floatarray.C:1026
oofem::cbrt
double cbrt(double x)
Returns the cubic root of x.
Definition mathfem.h:122
oofem::cubic3r
void cubic3r(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num)
Definition mathfem.C:155
oofem::signum
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:327
oofem::cubic
void cubic(double a, double b, double c, double d, double *r1, double *r2, double *r3, int *num)
Definition mathfem.C:43
oofem::sgn
double sgn(double i)
Returns the signum of given value (if value is < 0 returns -1, otherwise returns 1).
Definition mathfem.h:104
oofem::iperm
int iperm(int val, int rank)
Definition mathfem.C:261
oofem::ls2fit
void ls2fit(const FloatArray &x, const FloatArray &y, FloatArray &a)
Definition mathfem.C:272
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011