00001 #include "mathem.h"
00002 #include "basefun.h"
00003 #include "gtopology.h"
00004
00005 #include <errno.h>
00006 #include <math.h>
00007 #include <stdlib.h>
00008 #define MATH_ZERO 1.0e-100
00009
00010 double radius (vector &x,vector &natcoord)
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 {
00024 long nne;
00025 double r;
00026 nne=x.n;
00027 vector bf(nne);
00028
00029 switch (nne){
00030 case 3:{
00031 bf[0]=natcoord[0];
00032 bf[1]=natcoord[1];
00033 bf[2]=natcoord[2];
00034 break;
00035 }
00036 case 4:{
00037 bf_lin_4_2d (bf.a,natcoord[0],natcoord[1]);
00038 break;
00039 }
00040 default:{
00041 fprintf (stderr,"\n\n unknown number of nodes is required in function radius (mathem.cpp).\n");
00042 }
00043 }
00044
00045 scprd (x,bf,r);
00046 return r;
00047 }
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 double length (vector &coorda,vector &coordb)
00058 {
00059 long i,m,n;
00060 double l;
00061
00062 m=coorda.n; n=coordb.n;
00063 if (m!=n) fprintf (stderr,"\n\n non-compatible numbers of components in function length (mathem.cpp).\n");
00064
00065 l=0.0;
00066 for (i=0;i<m;i++){
00067 l+=(coorda[i]-coordb[i])*(coorda[i]-coordb[i]);
00068 }
00069 l=sqrt (l);
00070
00071 return l;
00072 }
00073
00074 double sqr(double x)
00075 {
00076 return (x*x);
00077 }
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092 void nodal_values (double **val,vector &nx,vector &ny,vector &nz,
00093 double *lhs,long dim,long fi,long ncomp)
00094 {
00095 long i,j,jj;
00096
00097 if (dim==1){
00098 for (i=0;i<nx.n;i++){
00099 jj=fi;
00100 for (j=0;j<ncomp;j++){
00101 val[i][jj] = lhs[j*2]+lhs[j*2+1]*nx[i];
00102 jj++;
00103 }
00104 }
00105 }
00106 if (dim==2){
00107 for (i=0;i<nx.n;i++){
00108 jj=fi;
00109 for (j=0;j<ncomp;j++){
00110 val[i][jj] = lhs[j*3]+lhs[j*3+1]*nx[i]+lhs[j*3+2]*ny[i];
00111 jj++;
00112 }
00113 }
00114 }
00115 if (dim==3){
00116 for (i=0;i<nx.n;i++){
00117 jj=fi;
00118 for (j=0;j<ncomp;j++){
00119 val[i][jj] = lhs[j*4]+lhs[j*4+1]*nx[i]+lhs[j*4+2]*ny[i]+lhs[j*4+3]*nz[i];
00120 jj++;
00121 }
00122 }
00123 }
00124
00125 }
00126
00127
00128
00129
00130
00131
00132 double sgn (double i)
00133 {
00134 return (i< 0.0 ? -1.0 : 1.0);
00135 }
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148 double heaviside (double x)
00149 {
00150 if (x<0.0)
00151 return 0.0;
00152 else
00153 return 1.0;
00154
00155 }
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168 double genheaviside (double x,double eps)
00169 {
00170 double h;
00171
00172 if (x<0.0-eps)
00173 h= 0.0;
00174 if (0.0-eps <= x && x <= eps)
00175 h = (x+eps)/2.0/eps;
00176 if (eps < x)
00177 h = 1.0;
00178
00179
00180
00181 return h;
00182
00183 }
00184
00185
00186
00187
00188
00189
00190 double polynom_4 (double x,double *a)
00191 {
00192 return (a[0]*x*x*x*x + a[1]*x*x*x + a[2]*x*x + a[3]*x + a[4]);
00193 }
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 void sort_2 (long *x)
00206 {
00207 if (x[0] < x[1])
00208 return;
00209 long a = x[0]; x[0] = x[1]; x[1] = a;
00210 }
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 void sort_3 (long *x)
00222 {
00223 long a;
00224
00225 if (x[0] > x[1])
00226 {
00227 if (x[0] > x[2])
00228 {
00229 if (x[1] > x[2])
00230 {
00231 a=x[0]; x[0]=x[2]; x[2]=a;
00232 }
00233 else
00234 {
00235 a=x[0]; x[0]=x[1]; x[1]=x[2]; x[2]=a;
00236 }
00237 }
00238 else
00239 {
00240 a=x[0]; x[0]=x[1]; x[1]=a;
00241 }
00242 }
00243 else
00244 {
00245 if (x[1] > x[2])
00246 {
00247 if (x[0] > x[2])
00248 {
00249 a=x[0]; x[0]=x[2]; x[2]=x[1]; x[1]=a;
00250 }
00251 else
00252 {
00253 a=x[1]; x[1]=x[2]; x[2]=a;
00254 }
00255 }
00256 }
00257 }
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 void sort_4 (long *x)
00269 {
00270 long a;
00271
00272 if (x[0] < x[1])
00273 {
00274 if (x[1] > x[2])
00275 {
00276 if (x[1] > x[3])
00277 {
00278 a=x[1]; x[1]=x[3]; x[3]=a; sort_3(x);
00279 return;
00280 }
00281 else
00282 {
00283 a=x[1]; x[1]=x[2]; x[2]=a;
00284 }
00285 }
00286 else
00287 {
00288 if (x[2] > x[3])
00289 {
00290 a=x[2]; x[2]=x[3]; x[3]=a; sort_3(x);
00291 return;
00292 }
00293 else
00294 {
00295 return;
00296 }
00297 }
00298 }
00299 else
00300 {
00301 if (x[0] > x[2])
00302 {
00303 if (x[0] > x[3])
00304 {
00305 a=x[0]; x[0]=x[3]; x[3]=a; sort_3(x);
00306 return;
00307 }
00308 else
00309 {
00310 a=x[0]; x[0]=x[2]; x[2]=a;
00311 }
00312 }
00313 else
00314 {
00315 if (x[2] > x[3])
00316 {
00317 a=x[2]; x[2]=x[3]; x[3]=a; sort_3(x);
00318 return;
00319 }
00320 }
00321 }
00322 if (x[0]>x[1])
00323 {
00324 a=x[0]; x[0]=x[1]; x[1]=a;
00325 }
00326 }
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337 void sort_3 (double *x)
00338 {
00339 double a;
00340
00341 if (x[0] > x[1])
00342 {
00343 if (x[0] > x[2])
00344 {
00345 if (x[1] > x[2])
00346 {
00347 a=x[0]; x[0]=x[2]; x[2]=a;
00348 }
00349 else
00350 {
00351 a=x[0]; x[0]=x[1]; x[1]=x[2]; x[2]=a;
00352 }
00353 }
00354 else
00355 {
00356 a=x[0]; x[0]=x[1]; x[1]=a;
00357 }
00358 }
00359 else
00360 {
00361 if (x[1] > x[2])
00362 {
00363 if (x[0] > x[2])
00364 {
00365 a=x[0]; x[0]=x[2]; x[2]=x[1]; x[1]=a;
00366 }
00367 else
00368 {
00369 a=x[1]; x[1]=x[2]; x[2]=a;
00370 }
00371 }
00372 }
00373 }
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 long solv_polynom_2 (double a,double b,double c,double &r1,double &r2)
00389 {
00390 if (fabs(a) < MATH_ZERO)
00391 {
00392 if (fabs(b) < MATH_ZERO)
00393 {
00394 if (fabs(c) < MATH_ZERO)
00395 return -1;
00396 else
00397 return 0;
00398 }
00399 else
00400 {
00401 r1 = -c/b;
00402 return 1;
00403 }
00404 }
00405 double D = b*b - 4*a*c;
00406
00407 if (D < -MATH_ZERO)
00408 return 0;
00409
00410 if (-MATH_ZERO<D && D<MATH_ZERO)
00411 {
00412 r1 = -b/2.0/a;
00413 return 1;
00414 }
00415
00416 r1 = (-b + sqrt(D))/2.0/a;
00417 r2 = (-b - sqrt(D))/2.0/a;
00418 return 2;
00419 }
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 long solv_polynom_3 (double a, double b, double c, double d, double &r1, double &r2, double &r3)
00433 {
00434 if(fabs(a) < MATH_ZERO)
00435 return solv_polynom_2 (b,c,d,r1,r2);
00436
00437
00438 double aa, p, q, D, pq, u, v, phi;
00439 double help,pi;
00440
00441 pi=3.141592653589793238;
00442
00443 aa = a;
00444 a = b/aa;
00445 b = c/aa;
00446 c = d/aa;
00447 p = b - a * a / 3.0;
00448 q = 2.0 * a * a * a / 27.0 - a * b / 3.0 + c;
00449 pq = p * q;
00450 D = q * q / 4.0 + p * p * p / 27.0;
00451
00452 if (fabs(D) < MATH_ZERO){
00453 if (fabs(pq) < MATH_ZERO){
00454 r1 = 0.0 - a / 3.0;
00455 r2 = r1;
00456 r3 = r1;
00457 return 3;
00458 }
00459
00460 if (q < 0.0)
00461 r2 = -exp(log(-q / 2.0) / 3.0);
00462 else
00463 r2 = exp(log(q / 2.0) / 3.0);
00464
00465 r1 = -2.0 * r2 - a / 3.0;
00466 r2 -= a / 3.0;
00467 return 2;
00468 }
00469
00470 if (D > 0.0){
00471 u = -q / 2.0 + sqrt(D);
00472 v = -q / 2.0 - sqrt(D);
00473 if (u < 0.0)
00474 u = -exp(log(-u) / 3.0);
00475 else
00476 u = exp(log(u) / 3.0);
00477 if (v < 0.0)
00478 v = -exp(log(-v) / 3.0);
00479 else
00480 v = exp(log(v) / 3.0);
00481 r1 = u + v - a / 3.0;
00482 return 1;
00483 }
00484
00485 p = sqrt(fabs(p) / 3.0);
00486 help = (-q / (2.0 * p * p * p));
00487 if (fabs (help) > 1.0) help = sgn(help);
00488
00489 phi = acos(help) / 3.0;
00490 r1 = 2.0 * p * cos(phi) - a / 3.0;
00491 r2 = -2.0 * p * cos(phi - pi / 3.0) - a / 3.0;
00492 r3 = -2.0 * p * cos(phi + pi / 3.0) - a / 3.0;
00493 return 3;
00494 }
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509 long solv_polynom_4 (double *coeff,double a,double b,double acc,double *roots)
00510 {
00511 long i,j,ninflps,nint,pol,maxit,nr;
00512 double inflp[5],x,y,ya,yb;
00513
00514 maxit = 200;
00515
00516 ninflps = solv_polynom_3 (4*coeff[0],3*coeff[1],2*coeff[2],coeff[3],inflp[1],inflp[2],inflp[3]);
00517 sort_3 (inflp+1);
00518
00519 nint = 1;
00520 inflp[0] = a;
00521 for (i=0;i<ninflps;i++)
00522 if ( a < inflp[i+1] && inflp[i+1] < b ) inflp[nint++] = inflp[i+1];
00523
00524 inflp[nint] = b;
00525
00526 nr = 0;
00527 for (i=0;i<nint;i++){
00528 a = inflp[i];
00529 b = inflp[i+1];
00530 ya = polynom_4 (a,coeff);
00531 yb = polynom_4 (b,coeff);
00532
00533 if (fabs(ya) < MATH_ZERO) roots[nr++] = a;
00534 else if (fabs(yb) < MATH_ZERO) roots[nr++] = b;
00535 else if (ya*yb<0) {
00536 pol = (ya < 0.0 ? 1 : -1);
00537
00538 for (j=0;j<maxit;j++){
00539
00540 switch (j%3){
00541 case 0:{ x = (a+b)/2.0; break; }
00542 case 1:{ x = a - ya*(b-a)/(yb-ya); break; }
00543 case 2:{ x = (a - ya*(b-a)/(yb-ya))*2.0 - x; break; }
00544 }
00545
00546 y = polynom_4 (x,coeff);
00547 if (0.0 < pol*y) { b = x; yb = y; }
00548 else { a = x; ya = y; }
00549
00550 if ( 2.0*(b-a)/(b+a) < acc ){
00551 roots[nr++] = x;
00552 break;
00553 }
00554 }
00555 }
00556 else continue;
00557
00558 if (j==maxit){
00559 fprintf (stderr,"\n\n\n ERROR 1 in function solv_polynom_4, PLEASE CONTACT AUTHOR of this function(%s, line %d)\n\n\n",__FILE__,__LINE__);
00560 exit (0);
00561 }
00562
00563 if ( nr!=1 && roots[nr-2]+acc > roots[nr-1]) nr--;
00564 }
00565
00566 return nr;
00567 }
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582 long solv_1le (double a,double b,double &x)
00583 {
00584 if (fabs(a) < MATH_ZERO)
00585 {
00586 if (fabs(b) < MATH_ZERO)
00587 return -1;
00588 else
00589 return 0;
00590 }
00591
00592 x = -b/a;
00593
00594 return 1;
00595 }
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611 long solv_2le (double *a, double *b, double *c,double &x,double &y)
00612 {
00613
00614 if (MATH_ZERO < fabs(a[0]))
00615 {
00616 if (MATH_ZERO > fabs(a[1]*(a[0]+b[0])-a[0]*(a[1]+b[1])))
00617 {
00618 if (MATH_ZERO < fabs(a[1]*c[0]-a[0]*c[1]))
00619 return 0;
00620 else
00621 return -1;
00622 }
00623 }
00624 else
00625 {
00626 if (MATH_ZERO > fabs(b[1]*(a[0]+b[0])-b[0]*(a[1]+b[1])))
00627 {
00628 if (MATH_ZERO < fabs(b[1]*c[0]-b[0]*c[1]) || MATH_ZERO < fabs(b[1]*c[0]-b[0]*c[1]))
00629 return 0;
00630 else
00631 {
00632 if (MATH_ZERO > fabs(b[0]) && MATH_ZERO > fabs(b[1]) && MATH_ZERO < fabs(c[0]))
00633 return 0;
00634 else
00635 return -1;
00636 }
00637 }
00638 }
00639
00640 if (MATH_ZERO < fabs(a[0]))
00641 {
00642 y = (c[0]*a[1]-c[1]*a[0]) / (a[0]*b[1]-a[1]*b[0]);
00643 x = (-c[0]-b[0]*y) / a[0];
00644 }
00645 else
00646 {
00647 y = -c[0]/b[0];
00648 x = (-c[1]-b[1]*y) / a[1];
00649 }
00650
00651 return 1;
00652 }
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668 long solv_2nle (double *a, double *b, double *c, double *d,double *x,double *y)
00669 {
00670 long ii,jj,nr;
00671
00672 if (fabs(a[0]) < MATH_ZERO && fabs(a[1]) < MATH_ZERO)
00673 return (solv_2le (b,c,d,x[0],y[0]));
00674
00675 if (fabs(a[0]) < MATH_ZERO || fabs(a[1]) < MATH_ZERO)
00676 {
00677 if (fabs(a[1]) > MATH_ZERO)
00678 {
00679 ii=0; jj=1;
00680 }
00681 else
00682 {
00683 ii=1; jj=0;
00684 }
00685 if (fabs(b[ii]) < MATH_ZERO)
00686 {
00687 if (fabs(c[ii]) < MATH_ZERO)
00688 {
00689 if (fabs(d[ii]) < MATH_ZERO)
00690 return -1;
00691 else
00692 return 0;
00693 }
00694 else
00695 {
00696 y[0] = -d[ii]/c[ii];
00697 return solv_1le ( a[jj]*y[0]+b[jj] , c[jj]*y[0]+d[jj] ,x[0]);
00698 }
00699 }
00700 else
00701 {
00702 if (fabs(c[ii]) < MATH_ZERO)
00703 {
00704 x[0] = -d[ii]/b[ii];
00705 return solv_1le ( a[jj]*x[0]+c[jj] , b[jj]*x[0]+d[jj] ,y[0]);
00706 }
00707 else
00708 {
00709 nr = solv_polynom_2 ( -a[jj]*b[ii] , -a[jj]*d[ii] + c[ii]*b[jj]-c[jj]*b[ii] , c[ii]*d[jj]-c[jj]*d[ii] ,x[0],x[1]);
00710 if (nr>0)
00711 y[0] = (-b[ii]*x[0]-d[ii]) / c[ii];
00712 if (nr>1)
00713 y[1] = (-b[ii]*x[1]-d[ii]) / c[ii];
00714 return nr;
00715 }
00716 }
00717 }
00718
00719 nr = solv_polynom_2 ( a[0]*b[1]-a[1]*b[0] , a[0]*d[1]-a[1]*d[0] + c[0]*b[1]-c[1]*b[0] , c[0]*d[1]-c[1]*d[0] ,x[0],x[1]);
00720
00721 for (long i=0;i<nr;i++)
00722 {
00723 if (fabs(y[i] = a[0]*x[i]+c[0]) > MATH_ZERO)
00724 y[i] = (-b[0]*x[i]-d[0]) / y[i] ;
00725 else
00726 {
00727 if (fabs(b[0]*x[i]+d[0]) > MATH_ZERO)
00728 fprintf (stderr,"\n\n\n ERROR 1 in function solv_3r, PLEASE CONTACT AUTHOR of this function(%s, line %d)\n\n\n",__FILE__,__LINE__);
00729 else
00730 return -1;
00731 }
00732 }
00733 return nr;
00734 }
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 void nc_lin_3_2d (double xx,double yy,double *x,double *y,double &xi,double &eta)
00748 {
00749 double det = (x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]);
00750
00751 xi = ( (x[1]*y[2] - x[2]*y[1]) + (y[1] - y[2])*xx + (x[2] - x[1])*yy )/det;
00752 eta = ( (x[2]*y[0] - x[0]*y[2]) + (y[2] - y[0])*xx + (x[0] - x[2])*yy )/det;
00753 }
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767 void nc_quad_3_2d (double acc,double xx,double yy,double *x,double *y,double &xi,double &eta)
00768 {
00769 long i,ne,nk;
00770 double aux,coeff[5],eeta[4],accabs;
00771 double ax,bx,cx,dx,ex,fx;
00772 double ay,by,cy,dy,ey,fy;
00773 double k,l,m,n,o,p,q,r,s;
00774
00775 accabs = acc*( fabs(x[0]-x[1]) + fabs(x[1]-x[2]) + fabs(x[2]-x[0]) )/2.0;
00776
00777 ax = 2.0*x[0] + 2.0*x[2] - 4.0*x[5];
00778 bx = 2.0*x[1] + 2.0*x[2] - 4.0*x[4];
00779 cx = -x[0] - 3.0*x[2] + 4.0*x[5];
00780 dx = -x[1] - 3.0*x[2] + 4.0*x[4];
00781 ex = 4.0*x[2] + 4.0*x[3] - 4.0*x[4] - 4.0*x[5];
00782 fx = x[2] - xx;
00783
00784 ay = 2.0*y[0] + 2.0*y[2] - 4.0*y[5];
00785 by = 2.0*y[1] + 2.0*y[2] - 4.0*y[4];
00786 cy = -y[0] - 3.0*y[2] + 4.0*y[5];
00787 dy = -y[1] - 3.0*y[2] + 4.0*y[4];
00788 ey = 4.0*y[2] + 4.0*y[3] - 4.0*y[4] - 4.0*y[5];
00789 fy = y[2] - yy;
00790
00791
00792 if ( fabs(ax)+fabs(bx)+fabs(ex)+fabs(ay)+fabs(by)+fabs(ey) < accabs){
00793 nc_lin_3_2d (xx,yy,x,y,xi,eta);
00794 return;
00795 }
00796
00797 if (fabs(ax) < accabs){
00798 coeff[0] = ay*bx*bx + by*ex*ex - ey*ex*bx ;
00799 coeff[1] = 2.0*ay*bx*dx + 2.0*by*cx*ex + dy*ex*ex - ey*cx*bx - cy*ex*bx - ey*ex*dx ;
00800 coeff[2] = 2.0*ay*bx*fx + 2.0*dy*cx*ex + ay*dx*dx + by*cx*cx + fy*ex*ex - cy*cx*bx - ey*cx*dx - cy*ex*dx - ey*ex*fx ;
00801 coeff[3] = 2.0*ay*fx*dx + 2.0*fy*cx*ex + dy*cx*cx - cy*cx*dx - ey*cx*fx - cy*ex*fx ;
00802 coeff[4] = ay*fx*fx + fy*cx*cx - cy*cx*fx ;
00803
00804 ne = solv_polynom_4 (coeff,-acc,1+acc,1e-13,eeta);
00805
00806 eta = 2.0;
00807 for (i=0;i<ne;i++)
00808 {
00809 xi = -( bx*eeta[i]*eeta[i] + dx*eeta[i] + fx ) / ( cx + ex*eeta[i] );
00810 if ( -acc < xi && xi < 1+acc )
00811 {
00812 if (eta==2.0)
00813 eta = eeta[i];
00814 else
00815 fprintf (stderr,"\n\n\n error 0 in function nc_quad_3_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00816 }
00817 }
00818 if (eta==2.0)
00819 fprintf (stderr,"\n\n\n error 4 in function nc_quad_3_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00820
00821 return;
00822 }
00823
00824 q = ex*ex-4.0*ax*bx ;
00825 r = 2.0*(cx*ex - 2.0*ax*dx) ;
00826 s = cx*cx - 4.0*ax*fx ;
00827 k = 4.0*by*ax*ax - 2.0*ey*ex*ax + q*ay + ex*ex*ay ;
00828 l = 4.0*dy*ax*ax - 2.0*ey*cx*ax + r*ay + 2.0*ex*(cx*ay - ax*cy) ;
00829 m = 4.0*fy*ax*ax - 2.0*cy*cx*ax + s*ay + cx*cx*ay ;
00830 n = ( 2.0*ax*ey - ex*(ay + ay) ) * ( 2.0*ax*ey - ex*(ay + ay) ) ;
00831 o = 2.0 * ( 2.0*ax*ey - ex*(ay + ay) ) * ( 2.0*ax*cy - ay*(cx + cx) ) ;
00832 p = ( 2.0*ax*cy - ay*(cx + cx) ) * ( 2.0*ax*cy - ay*(cx + cx) ) ;
00833
00834 coeff[0] = k*k - n*q ;
00835 coeff[1] = 2.0*k*l - n*r - o*q ;
00836 coeff[2] = 2.0*m*k + l*l - n*s - o*r - p*q ;
00837 coeff[3] = 2.0*l*m - o*s - p*r ;
00838 coeff[4] = m*m - p*s ;
00839
00840 ne = solv_polynom_4 (coeff,-acc,1+acc,1e-13,eeta);
00841
00842 nk = 0;
00843 eta = 2.0;
00844 for (i=0;i<ne;i++){
00845 nk = solv_polynom_2 ( ax , cx+ex*eeta[i] , bx*eeta[i]*eeta[i]+dx*eeta[i]+fx , xi , aux );
00846
00847 if ( nk>0 && -acc < xi && xi < 1+acc )
00848 {
00849 if (eta==2.0)
00850 eta = eeta[i];
00851 else
00852 fprintf (stderr,"\n\n\n error 1 in function nc_quad_3_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00853 }
00854 if ( nk>1 && -acc < aux && aux < 1+acc )
00855 {
00856 if (eta==2.0)
00857 {
00858 xi = aux;
00859 eta = eeta[i];
00860 }
00861 else
00862 fprintf (stderr,"\n\n\n error 2 in function nc_quad_3_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00863 }
00864 }
00865 if (eta==2.0)
00866 fprintf (stderr,"\n\n\n error 3 in function nc_quad_3_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00867
00868 }
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881
00882 void nc_lin_4_2d (double acc,double xx,double yy,double *x,double *y,double &xi,double &eta)
00883 {
00884 long nroots;
00885 double ax,bx,cx,dx,ay,by,cy,dy;
00886 double a,b,c,ksi1,ksi2,accabs;
00887
00888 accabs = acc*( fabs(x[0]-x[1]) + fabs(x[1]-x[2]) + fabs(x[2]-x[3]) + fabs(x[3]-x[0]) )/2.0;
00889
00890 ax = x[0] - x[1] + x[2] - x[3];
00891 bx = x[0] - x[1] - x[2] + x[3];
00892 cx = x[0] + x[1] - x[2] - x[3];
00893 dx = x[0] + x[1] + x[2] + x[3] - 4.0 * xx;
00894
00895 ay = y[0] - y[1] + y[2] - y[3];
00896 by = y[0] - y[1] - y[2] + y[3];
00897 cy = y[0] + y[1] - y[2] - y[3];
00898 dy = y[0] + y[1] + y[2] + y[3] - 4.0 * yy;
00899
00900 a = ax * by - ay * bx;
00901 b = ax * dy - ay * dx + by * cx - bx * cy;
00902 c = cx * dy - cy * dx;
00903
00904 nroots = solv_polynom_2 (a,b,c,ksi1,ksi2);
00905
00906 if (nroots == 0) fprintf (stderr,"\n\n\n error 0 in function nc_lin_4_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00907 if (nroots == 1) xi = ksi1;
00908 if (nroots == 2)
00909 {
00910 if (( -1-acc <= ksi1 && ksi1 <= 1+acc ) && ( ksi2 < -1-acc || 1+acc < ksi2 ))
00911 xi = ksi1;
00912 else
00913 {
00914 if (( -1-acc <= ksi2 && ksi2 <= 1+acc ) && ( ksi1 < -1-acc || 1+acc < ksi1 ))
00915 xi = ksi2;
00916 else fprintf (stderr,"\n\n\n error 1 in function nc_lin_4_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00917 }
00918 }
00919
00920 if ( fabs(cx+ax*xi) < accabs )
00921 {
00922 if ( fabs(bx) < accabs )
00923 fprintf (stderr,"\n\n\n error 2 in function nc_lin_4_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00924 else
00925 {
00926 if ( fabs(ay*xi+cy) < accabs || fabs(xi+dx/bx) > acc )
00927 fprintf (stderr,"\n\n\n error 3 in function nc_lin_4_2d (%s, line %d)\n\n\n",__FILE__,__LINE__);
00928 else
00929 eta = (-dy-by*xi)/(cy+ay*xi);
00930 }
00931 }
00932 else
00933 eta = (-dx-bx*xi)/(cx+ax*xi);
00934
00935 }
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949 void nc_quad_4_2d (double acc,double xx,double yy,double *x,double *y,double &xi,double &eta)
00950 {
00951 nc_lin_4_2d (acc,xx,yy,x,y,xi,eta);
00952
00953
00954
00955
00956
00957
00958 }
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989 void interpolelem (gtopology *gt,long &nli,long **&lin,long &ntr,long **&trn,long &nqu,long **&qun,long &nte,long **&ten,long &ncu,long **&cun,long &nmp,double **&mpn,char flag)
00990 {
00991 long i;
00992 long *nodes;
00993
00994 nli = ntr = nqu = nte = ncu = nmp = 0;
00995
00996 lin = new long* [gt->ne*8];
00997 trn = new long* [gt->ne*4];
00998 qun = new long* [gt->ne*4];
00999 ten = new long* [gt->ne*1];
01000 cun = new long* [gt->ne*1];
01001 mpn = new double* [gt->ne];
01002
01003 for (i=0;i<gt->ne;i++){
01004 nodes = gt->gelements[i].nodes;
01005
01006 switch (gt->gelements[i].auxinf){
01007 case 312:{
01008 trn[ntr] = new long[3]; trn[ntr][0] = nodes[0]; trn[ntr][1] = nodes[1]; trn[ntr][2] = nodes[2]; ntr++;
01009 break;
01010 }
01011 case 412:{
01012 qun[nqu] = new long [4];
01013 if (flag == 'e') { qun[nqu][0] = nodes[0]; qun[nqu][1] = nodes[1]; qun[nqu][2] = nodes[2]; qun[nqu][3] = nodes[3]; nqu++; }
01014 if (flag == 'd') { qun[nqu][0] = nodes[0]; qun[nqu][1] = nodes[3]; qun[nqu][2] = nodes[1]; qun[nqu][3] = nodes[2]; nqu++; }
01015 break;
01016 }
01017 case 622:{
01018 trn[ntr] = new long [3]; trn[ntr][0] = nodes[0]; trn[ntr][1] = nodes[3]; trn[ntr][2] = nodes[5]; ntr++;
01019 trn[ntr] = new long [3]; trn[ntr][0] = nodes[3]; trn[ntr][1] = nodes[1]; trn[ntr][2] = nodes[4]; ntr++;
01020 trn[ntr] = new long [3]; trn[ntr][0] = nodes[4]; trn[ntr][1] = nodes[2]; trn[ntr][2] = nodes[5]; ntr++;
01021 trn[ntr] = new long [3]; trn[ntr][0] = nodes[3]; trn[ntr][1] = nodes[4]; trn[ntr][2] = nodes[5]; ntr++;
01022 if (flag == 'd'){
01023 lin[nli] = new long [2]; lin[nli][0] = nodes[0]; lin[nli][1] = nodes[3]; nli++;
01024 lin[nli] = new long [2]; lin[nli][0] = nodes[3]; lin[nli][1] = nodes[1]; nli++;
01025 lin[nli] = new long [2]; lin[nli][0] = nodes[1]; lin[nli][1] = nodes[4]; nli++;
01026 lin[nli] = new long [2]; lin[nli][0] = nodes[4]; lin[nli][1] = nodes[2]; nli++;
01027 lin[nli] = new long [2]; lin[nli][0] = nodes[2]; lin[nli][1] = nodes[5]; nli++;
01028 lin[nli] = new long [2]; lin[nli][0] = nodes[5]; lin[nli][1] = nodes[0]; nli++;
01029 }
01030 break;
01031 }
01032 case 822:{
01033 if (flag == 'e'){
01034 qun[nqu] = new long [4]; qun[nqu][0] = nodes[7]; qun[nqu][1] = nodes[0]; qun[nqu][2] = nodes[4]; qun[nqu][3] = gt->nn+nmp; nqu++;
01035 qun[nqu] = new long [4]; qun[nqu][0] = nodes[4]; qun[nqu][1] = nodes[1]; qun[nqu][2] = nodes[5]; qun[nqu][3] = gt->nn+nmp; nqu++;
01036 qun[nqu] = new long [4]; qun[nqu][0] = nodes[5]; qun[nqu][1] = nodes[2]; qun[nqu][2] = nodes[6]; qun[nqu][3] = gt->nn+nmp; nqu++;
01037 qun[nqu] = new long [4]; qun[nqu][0] = nodes[6]; qun[nqu][1] = nodes[3]; qun[nqu][2] = nodes[7]; qun[nqu][3] = gt->nn+nmp; nqu++;
01038 }
01039 if (flag == 'd'){
01040 qun[nqu] = new long [4]; qun[nqu][0] = nodes[0]; qun[nqu][1] = nodes[7]; qun[nqu][2] = nodes[4]; qun[nqu][3] = gt->nn+nmp; nqu++;
01041 qun[nqu] = new long [4]; qun[nqu][0] = nodes[1]; qun[nqu][1] = nodes[4]; qun[nqu][2] = nodes[5]; qun[nqu][3] = gt->nn+nmp; nqu++;
01042 qun[nqu] = new long [4]; qun[nqu][0] = nodes[2]; qun[nqu][1] = nodes[5]; qun[nqu][2] = nodes[6]; qun[nqu][3] = gt->nn+nmp; nqu++;
01043 qun[nqu] = new long [4]; qun[nqu][0] = nodes[3]; qun[nqu][1] = nodes[6]; qun[nqu][2] = nodes[7]; qun[nqu][3] = gt->nn+nmp; nqu++;
01044 lin[nli] = new long [2]; lin[nli][0] = nodes[0]; lin[nli][1] = nodes[4]; nli++;
01045 lin[nli] = new long [2]; lin[nli][0] = nodes[4]; lin[nli][1] = nodes[1]; nli++;
01046 lin[nli] = new long [2]; lin[nli][0] = nodes[1]; lin[nli][1] = nodes[5]; nli++;
01047 lin[nli] = new long [2]; lin[nli][0] = nodes[5]; lin[nli][1] = nodes[2]; nli++;
01048 lin[nli] = new long [2]; lin[nli][0] = nodes[2]; lin[nli][1] = nodes[6]; nli++;
01049 lin[nli] = new long [2]; lin[nli][0] = nodes[6]; lin[nli][1] = nodes[3]; nli++;
01050 lin[nli] = new long [2]; lin[nli][0] = nodes[3]; lin[nli][1] = nodes[7]; nli++;
01051 lin[nli] = new long [2]; lin[nli][0] = nodes[7]; lin[nli][1] = nodes[0]; nli++;
01052 }
01053 mpn[nmp] = new double [3];
01054 mpn[nmp][0] = (gt->gnodes[nodes[4]].x + gt->gnodes[nodes[5]].x + gt->gnodes[nodes[6]].x + gt->gnodes[nodes[7]].x) / 4.0 ;
01055 mpn[nmp][1] = (gt->gnodes[nodes[4]].y + gt->gnodes[nodes[5]].y + gt->gnodes[nodes[6]].y + gt->gnodes[nodes[7]].y) / 4.0 ;
01056 mpn[nmp][2] = i;
01057 nmp++;
01058 break;
01059 }
01060 case 413:{
01061 if (flag == 'd'){
01062 ten[nte] = new long[4]; ten[nte][0] = nodes[0]; ten[nte][1] = nodes[1]; ten[nte][2] = nodes[2]; ten[nte][3] = nodes[3]; nte++;
01063 }
01064 if (flag == 'e')
01065 fprintf (stderr,"\n\n function interpolelem does not transform tetrahedras for drawing in elixir \n\n");
01066 break;
01067 }
01068 case 813:{
01069 if (flag == 'd'){
01070 cun[ncu] = new long[8]; cun[ncu][0] = nodes[0]; cun[ncu][1] = nodes[1]; cun[ncu][2] = nodes[3]; cun[ncu][3] = nodes[2];
01071 cun[ncu][4] = nodes[4]; cun[ncu][5] = nodes[5]; cun[ncu][6] = nodes[7]; cun[ncu][7] = nodes[6]; ncu++;
01072 }
01073 if (flag == 'e')
01074 fprintf (stderr,"\n\n function interpolelem does not transform hexahedras for drawing in elixir \n\n");
01075 break;
01076 }
01077 default:{
01078 fprintf (stderr,"\n\n wrong nnedegdim in function");
01079 fprintf (stderr," interpolelem (%s, line %d)",__FILE__,__LINE__);
01080 }
01081 }
01082 }
01083 }
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096 void print_ex (gtopology *gt,const char *file, double *valnod, double *valel)
01097 {
01098 long i,nn,ne;
01099 long nli,ntr,nqu,ntet,ncub,nmp;
01100 long **linodes,**trnodes,**qunodes,**tetnodes,**cubnodes;
01101 double **mpnodes;
01102 FILE *data_ex;
01103
01104 nn = gt->nn;
01105 ne = gt->ne;
01106
01107 interpolelem (gt,nli,linodes,ntr,trnodes,nqu,qunodes,ntet,tetnodes,ncub,cubnodes,nmp,mpnodes,'e');
01108
01109 data_ex = fopen (file,"w"); if (data_ex==NULL) fprintf (stderr,"\n .ex file has not been opened.\n\n");
01110
01111
01112 fprintf (data_ex,"%10d%10d\n",7,1);
01113 fprintf (data_ex,"%10ld%10d%10ld%10ld%10d%10ld%10d%10d\n\n",nn+nmp,0,ntr,nqu,0,ntet,0,0);
01114
01115
01116 for (i=0;i<nn;i++)
01117 fprintf (data_ex,"%10ld %.6e %.6e %.6e %.6e\n",i+1,gt->gnodes[i].x,gt->gnodes[i].y,gt->gnodes[i].z,valnod[i]);
01118
01119 if (nmp)
01120 for (i=0;i<nmp;i++)
01121 fprintf (data_ex,"%10ld %.6e %.6e %.6e %.6e\n",nn+i+1,mpnodes[i][0],mpnodes[i][1],0.0,valel[(long)mpnodes[i][2]]);
01122
01123 fprintf (data_ex,"\n");
01124
01125
01126 for (i=0;i<ntr;i++)
01127 fprintf (data_ex,"%10ld%10ld%10ld%10ld\n",i+1,trnodes[i][0]+1,trnodes[i][1]+1,trnodes[i][2]+1);
01128
01129 for (i=0;i<nqu;i++)
01130 fprintf (data_ex,"%10ld%10ld%10ld%10ld%10ld\n",i+1,qunodes[i][0]+1,qunodes[i][1]+1,qunodes[i][2]+1,qunodes[i][3]+1);
01131
01132 for (i=0;i<ntet;i++)
01133 fprintf (data_ex,"%10ld%10ld%10ld%10ld%10ld\n",i+1,tetnodes[i][0]+1,tetnodes[i][1]+1,tetnodes[i][2]+1,tetnodes[i][3]+1);
01134
01135
01136 fclose (data_ex);
01137
01138 for (i=0;i<nli;i++) delete [] linodes[i];
01139 delete [] linodes;
01140 for (i=0;i<ntr;i++) delete [] trnodes[i];
01141 delete [] trnodes;
01142 for (i=0;i<nqu;i++) delete [] qunodes[i];
01143 delete [] qunodes;
01144 for (i=0;i<ntet;i++) delete [] tetnodes[i];
01145 delete [] tetnodes;
01146 for (i=0;i<ncub;i++) delete [] cubnodes[i];
01147 delete [] cubnodes;
01148 for (i=0;i<nmp;i++) delete [] mpnodes[i];
01149 delete [] mpnodes;
01150 }
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165 void print_dx (gtopology *gt, const char *file,double **valnod,double **valel,char tve, long *dimindex, char **caption,long nindex)
01166 {
01167 long i,j,k,nn,ne,dim;
01168 long nli,ntr,nqu,ntet,ncub,nmp,object,value,vali;
01169 long **linodes,**trnodes,**qunodes,**tetnodes,**cubnodes;
01170 double **mpnodes;
01171 FILE *data_dx;
01172
01173 nn = gt->nn;
01174 ne = gt->ne;
01175 dim = gt->gelements[0].auxinf%10;
01176 object = 0;
01177
01178 interpolelem (gt,nli,linodes,ntr,trnodes,nqu,qunodes,ntet,tetnodes,ncub,cubnodes,nmp,mpnodes,'d');
01179
01180 data_dx = fopen (file,"w"); if (data_dx==NULL) fprintf (stderr,"\n .dx file has not been opened.\n\n");
01181
01182
01183 fprintf (data_dx,"object %ld class array type float rank 1 shape %ld items %ld data follows\n",++object,dim,nn+nmp);
01184
01185 for (i=0;i<nn;i++){
01186 fprintf (data_dx," %20.10f",gt->gnodes[i].x);
01187 if (dim==3 || dim==2) fprintf (data_dx," %20.10f",gt->gnodes[i].y);
01188 if (dim==3) fprintf (data_dx," %20.10f",gt->gnodes[i].z);
01189 fprintf (data_dx,"\n");
01190 }
01191
01192 if (nmp)
01193 for (i=0;i<nmp;i++)
01194 fprintf (data_dx," %20.10f %20.10f\n",mpnodes[i][0],mpnodes[i][1]);
01195
01196
01197 if (ntr){
01198 fprintf (data_dx,"object %ld class array type int rank 1 shape %d items %ld data follows\n",++object,3,ntr);
01199
01200 for (i=0;i<ntr;i++){
01201 for (j=0;j<3;j++)
01202 fprintf (data_dx," %10ld",trnodes[i][j]);
01203
01204 fprintf (data_dx,"\n");
01205 }
01206
01207 fprintf (data_dx,"attribute \"element type\" string \"%s\"\n","triangles");
01208 fprintf (data_dx,"attribute \"ref\" string \"positions\"\n");
01209 }
01210
01211 if (nqu){
01212 fprintf (data_dx,"object %ld class array type int rank 1 shape %d items %ld data follows\n",++object,4,nqu);
01213
01214 for (i=0;i<nqu;i++){
01215 for (j=0;j<4;j++)
01216 fprintf (data_dx," %10ld",qunodes[i][j]);
01217
01218 fprintf (data_dx,"\n");
01219 }
01220
01221 fprintf (data_dx,"attribute \"element type\" string \"%s\"\n","quads");
01222 fprintf (data_dx,"attribute \"ref\" string \"positions\"\n");
01223 }
01224
01225 if (ntet){
01226 fprintf (data_dx,"object %ld class array type int rank 1 shape %d items %ld data follows\n",++object,4,ntet);
01227
01228 for (i=0;i<ntet;i++){
01229 for (j=0;j<4;j++)
01230 fprintf (data_dx," %10ld",tetnodes[i][j]);
01231
01232 fprintf (data_dx,"\n");
01233 }
01234
01235 fprintf (data_dx,"attribute \"element type\" string \"%s\"\n","tetrahedra");
01236 fprintf (data_dx,"attribute \"ref\" string \"positions\"\n");
01237 }
01238
01239 if (ncub){
01240 fprintf (data_dx,"object %ld class array type int rank 1 shape %d items %ld data follows\n",++object,8,ncub);
01241
01242 for (i=0;i<ncub;i++){
01243 for (j=0;j<8;j++)
01244 fprintf (data_dx," %10ld",cubnodes[i][j]);
01245
01246 fprintf (data_dx,"\n");
01247 }
01248
01249 fprintf (data_dx,"attribute \"element type\" string \"%s\"\n","cubes");
01250 fprintf (data_dx,"attribute \"ref\" string \"positions\"\n");
01251 }
01252
01253 if (nli){
01254 fprintf (data_dx,"object %ld class array type int rank 1 shape %d items %ld data follows\n",++object,2,nli);
01255
01256 for (i=0;i<nli;i++){
01257 for (j=0;j<2;j++)
01258 fprintf (data_dx," %10ld",linodes[i][j]);
01259
01260 fprintf (data_dx,"\n");
01261 }
01262
01263 fprintf (data_dx,"attribute \"element type\" string \"%s\"\n","lines");
01264 fprintf (data_dx,"attribute \"ref\" string \"positions\"\n");
01265 }
01266
01267
01268
01269 vali = 0;
01270 for (i=0;i<nindex;i++){
01271 if (dimindex[i]==1)
01272 fprintf (data_dx,"object %ld class array type float rank 0 items %ld data follows\n",++object,nn+nmp);
01273 else
01274 fprintf (data_dx,"object %ld class array type float rank 1 shape %ld items %ld data follows\n",++object,dimindex[i],nn+nmp);
01275
01276 for (j=0;j<nn;j++){
01277 for (k=0;k<dimindex[i];k++)
01278 fprintf (data_dx," %20.10f",valnod[vali+k][j]);
01279
01280 fprintf (data_dx,"\n");
01281 }
01282
01283 if (nmp)
01284 for (j=0;j<nmp;j++){
01285 for (k=0;k<dimindex[i];k++)
01286 if (tve=='n') fprintf (data_dx," %20.10f",valel[vali+k][(long)mpnodes[j][2]]);
01287 else if (tve=='t') fprintf (data_dx," %20.10f",valel[(long)mpnodes[j][2]][vali+k]);
01288 else fprintf (stderr,"\n\n\n error in print_dx in mathem.cpp \n\n\n");
01289
01290 fprintf (data_dx,"\n");
01291 }
01292
01293
01294 fprintf (data_dx,"attribute \"dep\" string \"positions\"\n");
01295 vali += dimindex[i];
01296 }
01297
01298
01299 value = 1;
01300 if (ntr) value++;
01301 if (nqu) value++;
01302 if (ntet) value++;
01303 if (ncub) value++;
01304 if (nli) value++;
01305 for (i=0;i<nindex;i++){
01306 fprintf (data_dx,"object %ld class field\n",++object);
01307 fprintf (data_dx,"component \"positions\" value 1\n");
01308 fprintf (data_dx,"component \"connections\" value 2\n");
01309 fprintf (data_dx,"component \"data\" value %ld\n",++value);
01310 if (ntr*nqu){
01311 fprintf (data_dx,"object %ld class field\n",++object);
01312 fprintf (data_dx,"component \"positions\" value 1\n");
01313 fprintf (data_dx,"component \"connections\" value 3\n");
01314 fprintf (data_dx,"component \"data\" value %ld\n",value);
01315 }
01316 }
01317
01318
01319 if (ntr*nqu){
01320 for (i=0;i<nindex;i++){
01321 fprintf (data_dx,"object %ld class compositefield\n",++object);
01322 fprintf (data_dx,"member 0 value %ld\n",++value);
01323 fprintf (data_dx,"member 1 value %ld\n",++value);
01324 }
01325 }
01326
01327
01328 if (nli){
01329 fprintf (data_dx,"object %ld class field\n",++object);
01330 fprintf (data_dx,"component \"positions\" value 1\n");
01331 if (nqu*ntr){
01332 fprintf (data_dx,"component \"connections\" value 4\n");
01333 fprintf (data_dx,"component \"data\" value 5\n");
01334 }
01335 else{
01336 fprintf (data_dx,"component \"connections\" value 3\n");
01337 fprintf (data_dx,"component \"data\" value 4\n");
01338 }
01339 }
01340
01341
01342 fprintf (data_dx,"object %ld class group\n",++object);
01343
01344 for (i=0;i<nindex;i++)
01345 fprintf (data_dx,"member \"%s\" value %ld\n",caption[i],++value);
01346
01347 if (nli)
01348 fprintf (data_dx,"member \"%s\" value %ld\n","lincon",++value);
01349
01350 fclose (data_dx);
01351
01352 for (i=0;i<nli;i++) delete [] linodes[i];
01353 delete [] linodes;
01354 for (i=0;i<ntr;i++) delete [] trnodes[i];
01355 delete [] trnodes;
01356 for (i=0;i<nqu;i++) delete [] qunodes[i];
01357 delete [] qunodes;
01358 for (i=0;i<ntet;i++) delete [] tetnodes[i];
01359 delete [] tetnodes;
01360 for (i=0;i<ncub;i++) delete [] cubnodes[i];
01361 delete [] cubnodes;
01362 for (i=0;i<nmp;i++) delete [] mpnodes[i];
01363 delete [] mpnodes;
01364 }
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374 void maxmin_3 (double *x,double &max,double &min)
01375 {
01376 if (x[0] > x[1]) if (x[0] > x[2]) if (x[1] > x[2]){ max = x[0]; min = x[2];}
01377 else { max = x[0]; min = x[1];}
01378 else { max = x[2]; min = x[1];}
01379 else if (x[1] > x[2]) if (x[0] > x[2]){ max = x[1]; min = x[2];}
01380 else { max = x[1]; min = x[0];}
01381 else { max = x[2]; min = x[0];}
01382 }
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392 void maxmin_4 (double *x,double &max,double &min)
01393 {
01394 if (x[0] > x[1]) if (x[0] > x[2]) if (x[0] > x[3]) if (x[1] < x[2]) if (x[1] < x[3]){ max = x[0]; min = x[1];}
01395 else { max = x[0]; min = x[3];}
01396 else if (x[2] < x[3]){ max = x[0]; min = x[2];}
01397 else { max = x[0]; min = x[3];}
01398 else if (x[1] < x[2]){ max = x[3]; min = x[1];}
01399 else { max = x[3]; min = x[2];}
01400 else if (x[2] > x[3]) if (x[1] < x[3]){ max = x[2]; min = x[1];}
01401 else { max = x[2]; min = x[3];}
01402 else { max = x[3]; min = x[1];}
01403
01404 else if (x[1] > x[2]) if (x[1] > x[3]) if (x[0] < x[2]) if (x[0] < x[3]){ max = x[1]; min = x[0];}
01405 else { max = x[1]; min = x[3];}
01406 else if (x[2] < x[3]){ max = x[1]; min = x[2];}
01407 else { max = x[1]; min = x[3];}
01408 else if (x[0] < x[2]){ max = x[3]; min = x[0];}
01409 else { max = x[3]; min = x[2];}
01410 else if (x[2] > x[3]) if (x[0] < x[3]){ max = x[2]; min = x[0];}
01411 else { max = x[2]; min = x[3];}
01412 else { max = x[3]; min = x[0];}
01413 }
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444 double lsm_quad(matrix &a, vector &r, vector &l, double x, double y, double x_old, double y_old, double zero, long solc)
01445 {
01446 double det_a, det_b1, det_b2, det_b3;
01447
01448 if ((a.m != 3) && (a.n != 3) && (r.n != 3) && (l.n != 3))
01449 {
01450 print_err("Wrong dimension of matrix or vectors in lest square method", __FILE__, __LINE__, __func__);
01451 abort();
01452 }
01453
01454 a[0][0] += 1.0; a[0][1] += x; a[0][2] += x*x;
01455 r[0] += y;
01456 a[1][0] += x; a[1][1] += x*x; a[1][2] += x*x*x;
01457 r[1] += x*y;
01458 a[2][0] += x*x; a[2][1] += x*x*x; a[2][2] += x*x*x*x;
01459 r[2] += x*x*y;
01460
01461 if (solc == 0)
01462 return 0.0;
01463
01464 if (solc == 2)
01465 {
01466
01467 a[0][0] -= 1.0; a[0][1] -= x_old; a[0][2] -= x_old*x_old;
01468 r[0] -= y_old;
01469 a[1][0] -= x_old; a[1][1] -= x_old*x_old; a[1][2] -= x_old*x_old*x_old;
01470 r[1] -= x_old*y_old;
01471 a[2][0] -= x_old*x_old; a[2][1] -= x_old*x_old*x_old; a[2][2] -= x_old*x_old*x_old*x_old;
01472 r[2] -= x_old*x_old*y_old;
01473 }
01474
01475
01476 det_a = a[0][0]*a[1][1]*a[2][2];
01477 det_a += a[0][1]*a[1][2]*a[2][0];
01478 det_a += a[0][2]*a[1][0]*a[2][1];
01479 det_a -= a[0][2]*a[1][1]*a[2][0];
01480 det_a -= a[0][0]*a[1][2]*a[2][1];
01481 det_a -= a[0][1]*a[1][0]*a[2][2];
01482
01483 if (fabs(det_a) < zero)
01484 {
01485 print_err("Determinant in least square method is close to zero", __FILE__, __LINE__, __func__);
01486 return 0.0;
01487 }
01488
01489 det_b1 = r[0]*a[1][1]*a[2][2];
01490 det_b1 += a[0][1]*a[1][2]*r[2];
01491 det_b1 += a[0][2]*r[1]*a[2][1];
01492 det_b1 -= a[0][2]*a[1][1]*r[2];
01493 det_b1 -= r[0]*a[1][2]*a[2][1];
01494 det_b1 -= a[0][1]*r[1]*a[2][2];
01495
01496 det_b2 = a[0][0]*r[1]*a[2][2];
01497 det_b2 += r[0]*a[1][2]*a[2][0];
01498 det_b2 += a[0][2]*a[1][0]*r[2];
01499 det_b2 -= a[0][2]*r[1]*a[2][0];
01500 det_b2 -= a[0][0]*a[1][2]*r[2];
01501 det_b2 -= r[0]*a[1][0]*a[2][2];
01502
01503 det_b3 = a[0][0]*a[1][1]*r[2];
01504 det_b3 += a[0][1]*r[1]*a[2][0];
01505 det_b3 += r[0]*a[1][0]*a[2][1];
01506 det_b3 -= r[0]*a[1][1]*a[2][0];
01507 det_b3 -= a[0][0]*r[1]*a[2][1];
01508 det_b3 -= a[0][1]*a[1][0]*r[2];
01509
01510 l[0] = det_b1/det_a;
01511 l[1] = det_b2/det_a;
01512 l[2] = det_b3/det_a;
01513
01514 return 2*l[2]*x + l[1];
01515 }
01516
01517
01518 int check_math_err()
01519 {
01520 if (errno == EDOM)
01521 {
01522 print_err("math domain error (nan) detected\n", __FILE__, __LINE__, __func__);
01523 abort();
01524 }
01525
01526 if (errno == ERANGE)
01527 {
01528 print_err("math range error (inf) detected\n", __FILE__, __LINE__, __func__);
01529 abort();
01530 }
01531 return errno;
01532 }
01533
01534
01535
01536 int check_math_errel(long eid)
01537 {
01538 if (errno == EDOM)
01539 {
01540 print_err("math domain error (nan) detected on element %ld\n", __FILE__, __LINE__, __func__, eid+1);
01541 abort();
01542 }
01543
01544 if (errno == ERANGE)
01545 {
01546 print_err("math range error (inf) detected on element %ld\n", __FILE__, __LINE__, __func__, eid+1);
01547 abort();
01548 }
01549 return errno;
01550 }
01551
01552
01553
01554 long test_math_err()
01555 {
01556 if ((errno == EDOM) || (errno == ERANGE))
01557 return 1;
01558
01559 return 0;
01560 }