Changeset 2928


Ignore:
Timestamp:
01/28/10 16:31:01 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added my routine of Reduction simultanee

Location:
issm/trunk/src/c/Bamgx
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Bamgx/Bamgx.cpp

    r2914 r2928  
    9696                        MetricAnIso M=Gh[i];
    9797                        MatVVP2x2 Vp(M/coef);
    98                         Vp.Maxh(hmin);
    99                         Vp.Minh(hmax);
     98                        Vp.Maxh(hmax);
     99                        Vp.Minh(hmin);
    100100                        Gh.vertices[i].m = Vp;
    101101                }

  • issm/trunk/src/c/Bamgx/Metric.h

    r2924 r2928  
    6363                MetricAnIso( Real8  a,const  MetricAnIso ma,
    6464                                        Real8  b,const  MetricAnIso mb);
    65                 int  IntersectWith(const MetricAnIso M);
     65                int  IntersectWith(const MetricAnIso M2);
     66                int  IntersectWith_new(const MetricAnIso M2);
    6667                MetricAnIso operator*(Real8 c) const {Real8 c2=c*c;return  MetricAnIso(a11*c2,a21*c2,a22*c2);}
    6768                MetricAnIso operator/(Real8 c) const {Real8 c2=1/(c*c);return  MetricAnIso(a11*c2,a21*c2,a22*c2);}
    … …  
    122123        }
    123124
    124         void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
    125         MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
     125        void SimultaneousMatrixReduction(const MetricAnIso M1,const MetricAnIso M2, MetricAnIso* M);
     126        void  SimultaneousMatrixReduction_bamg( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V);
     127        void eigen2(double* M,double lambda[2],double vp[2][2]);;
    126128
    127129        inline MetricAnIso::MetricAnIso(const MatVVP2x2 M) {
    … …  
    152154
    153155        class SaveMetricInterpole {
    154                 friend  Real8 LengthInterpole(const MetricAnIso ,const  MetricAnIso , R2 );
     156                friend Real8 LengthInterpole(const MetricAnIso ,const  MetricAnIso , R2 );
    155157                friend Real8 abscisseInterpole(const MetricAnIso ,const  MetricAnIso , R2 ,Real8 ,int );
    156158                int opt;

  • issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp

    r2917 r2928  
    3030                if ( n2 < 1e-30) lambda1=lambda2=0,v.x=1,v.y=0;
    3131
    32                 //if ???
     32                //if matrix is already diagonal and has a 0 eigen value
    3333                else if (delta < eps*n2){
    3434                        lambda1=lambda2=-b/2, v.x=1,v.y=0;

  • issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp

    r2924 r2928  
    6969        }
    7070        /*}}}*/
     71        /*FUNCTION MetricAnIso::IntersectWith_new{{{1*/
     72        int MetricAnIso::IntersectWith_new(MetricAnIso M2) {
     73                /*Get a new metric from an existing metric (M1=this)
     74                 * and a new metric given in input M2 using a
     75                 * Simultaneous Matrix Reduction
     76                 * */
     77
     78                MetricAnIso M;
     79
     80                //Compute M intersection of M1 and M2
     81                SimultaneousMatrixReduction(*this,M2,&M);
     82
     83                //check wether something has been done
     84                if(M.a11!=a11 || M.a21!=a21 || M.a22!=a22){
     85
     86                        //update this with M
     87                        a11=M.a11;
     88                        a21=M.a21;
     89                        a22=M.a22;
     90
     91                        //return 1 -> something has been done
     92                        return 1;
     93                }
     94                else{
     95
     96                        //nothing has changed, return 0
     97                        return 0;
     98                }
     99        }
     100        /*}}}1*/
    71101        /*FUNCTION MetricAnIso::IntersectWith{{{1*/
    72102        int MetricAnIso::IntersectWith(const MetricAnIso M2) {
    … …  
    81111                double l1,l2;
    82112
    83                 SimultaneousMatrixReduction(*this,M2,l1,l2,M);
     113                SimultaneousMatrixReduction_bamg(*this,M2,l1,l2,M);
    84114
    85115                R2 v1(M.x.x,M.y.x);
    … …  
    103133
    104134        /*Intermediary*/
    105         /*FUNCTION SimultaneousMatrixReduction{{{1*/
    106         void SimultaneousMatrixReduction( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
    107                 double a11=M1.a11,a21=M1.a21,a22=M1.a22;
    108                 double b11=M2.a11,b21=M2.a21,b22=M2.a22;
    109                 //  M1 v = l M2 v
    110                 // (M1 - l M2) v =0
    111                 // det (M1 - l M2) =0
    112                 // det (M1 - l M2) = a l^2 + b l + c;
    113                 // = (a11 - l * b11) * (a22 - l * b22) - (a21 - l * b21 ) ^2
    114                 const double /*c11 = a11*a11,*/ c21= a21*a21;
    115                 const double /*d11 = b11*b11,*/ d21= b21*b21;
    116                 const double a=b11*b22 - d21;
    117                 const double b=-a11*b22-a22*b11+2*a21*b21;
    118                 const double c=-c21+a11*a22;
    119                 const double bb = b*b,ac= a*c;
    120                 const double delta = bb - 4 * ac;
    121                 if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) ){
    122                         // racine double;
    123                         if (Abs(a) < 1.e-30 )
    124                          l1 = l2 = 0;
    125                         else
    126                          l1=l2=-b/(2*a);
    127                         V= D2xD2(1,0,0,1);
    128                   }
    129                 else {
    130                         const double delta2 = sqrt(delta);
    131                         l1= (-b - delta2)/(2*a);
    132                         l2= (-b + delta2)/(2*a);
    133                         // M1 v = l M2 v
    134                         //  ( (M1 - I M2) x,y)  = (x,(M1 - I M2) y) \forall y
    135                         // so Ker((M1 - I M2)) = Im((M1 - I M2))^\perp
    136                         double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
    137                         double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
    138                         double vp1x,vp1y,vp2x,vp2y;
    139 
    140                         if(s1 < s0)
    141                          s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
    142                         else
    143                          s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;
    144 
    145                         v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
    146                         s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
    147                         if(s1 < s0)
    148                          s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
    149                         else
    150                          s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
    151                         V=D2xD2(vp1x,vp2x,vp1y,vp2y);
    152                 }
    153                 return;
    154         }
    155         /*}}}1*/
    156         /*FUNCTION Intersection{{{1*/
    157         MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) {
    158                 D2xD2 M;
    159                 double l1,l2;
    160                 SimultaneousMatrixReduction(M1,M2,l1,l2,M);
    161                 R2 v0(M.x.x,M.y.x);
    162                 R2 v1(M.x.y,M.y.y);
    163                 D2xD2 M_1(M.inv());
    164                 D2xD2 D(Max(M1(v0,v0),M2(v0,v0)),0,0,Max(M1(v1,v1),M2(v1,v1)));
    165                 D2xD2 Mi(M_1.t()*D*M_1);
    166                 return MetricAnIso(Mi.x.x,0.5*(Mi.x.y+Mi.y.x),Mi.y.y);
     135        /*FUNCTION eigen2{{{1*/
     136        void eigen2(double* M,double* lambda,double vp[2][2]) {
     137                /*m=[a b]
     138                 *  [c d] */
     139
     140                double normM,detM,trM;
     141                double norm;
     142                const double a=M[2*0+0];
     143                const double b=M[2*0+1];
     144                const double c=M[2*1+0];
     145                const double d=M[2*1+1];
     146
     147                //compute norm(M)
     148                normM=pow(pow(a,2.0)+pow(b,2.0)+pow(c,2.0)+pow(d,2.0),0.5);
     149
     150                //if normM<10-30, M=zeros(2,2)
     151                if (normM<1.0e-30){
     152                        lambda[0]=0.0;
     153                        lambda[1]=0.0;
     154                        vp[0][0]=vp[0][1]=vp[1][0]=vp[1][1]=0.0;
     155                        return;
     156                }
     157
     158                //check that matrix is already diagonal if det(M)=0
     159                if(Abs(b)<(1.0e-12*normM) & Abs(c)<(1.0e-12*normM)){
     160                        lambda[0]=a;
     161                        lambda[1]=d;
     162                        vp[0][0]=1;
     163                        vp[0][1]=0;
     164                        vp[1][0]=0;
     165                        vp[1][1]=1;
     166                        return;
     167                }
     168
     169                //compute det(M)
     170                detM=a*d-c*b;
     171
     172                //not symetrical and still det<10^-30
     173                if (Abs(detM)<1.0e-30){
     174                        throw ErrorException(__FUNCT__,exprintf("Cannot find eigen values of a non invertible matrix"));
     175                }
     176
     177                //See http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/index.html
     178                trM=a+d;
     179
     180                //general case
     181                lambda[0] = 0.5*(trM + pow( trM*trM - 4.0*detM ,0.5));
     182                lambda[1] = 0.5*(trM - pow( trM*trM - 4.0*detM ,0.5));
     183
     184                if(Abs(b)>1.0e-30){ //if b!=0
     185                        norm=pow(b*b+pow(lambda[0]-a,2.0),0.5);
     186                        vp[0][0]=b/norm;
     187                        vp[1][0]=(lambda[0]-a)/norm;
     188                        norm=pow(b*b+pow(lambda[1]-a,2.0),0.5);
     189                        vp[0][1]=b/norm;
     190                        vp[1][1]=(lambda[1]-a)/norm;
     191                }
     192                else{
     193                        norm=pow(c*c+pow(lambda[0]-d,2.0),0.5);
     194                        vp[0][0]=(lambda[0]-d)/norm;
     195                        vp[1][0]=c/norm;
     196                        norm=pow(c*c+pow(lambda[1]-d,2.0),0.5);
     197                        vp[0][1]=(lambda[1]-d)/norm;
     198                        vp[1][1]=c/norm;
     199                }
    167200        }
    168201        /*}}}1*/
    … …  
    233266        }
    234267        /*}}}1*/
     268        /*FUNCTION SimultaneousMatrixReduction{{{1*/
     269        void SimultaneousMatrixReduction(MetricAnIso M1,MetricAnIso M2,MetricAnIso* pM) {
     270
     271                /*intermediary*/
     272                double a11,a21,a22;
     273                double b11,b21,b22;
     274                double mu1,mu2;
     275                double lambda1,lambda2;
     276                double detM1,detP;
     277                double N[2][2];
     278                double lambda[2];
     279                double vp[2][2];
     280                double invP[2][2];
     281                MetricAnIso M;
     282
     283                //get components of both metrics
     284                a11=M1.a11; a21=M1.a21; a22=M1.a22;
     285                b11=M2.a11; b21=M2.a21; b22=M2.a22;
     286
     287                //check diagnonal matrices (easy)
     288                if (Abs(a21)< 1.0e-20 && Abs(b21)< 1.0e-20){
     289                        M.a11=Max(a11,b11);
     290                        M.a22=Max(a22,b22);
     291                        M.a21=0;
     292                        *pM=M;
     293                        return;
     294                }
     295
     296                //Build N=M1^-1 M2 (Alauzet2003 p16)
     297                detM1=a11*a22-a21*a21;
     298                if (!detM1){
     299                        printf("a11 a21 a22 = %g %g %g\n",a11,a21,a22);
     300                        throw ErrorException(__FUNCT__,exprintf("One of the metric has a zero eigen value or cannot be inverted... (detP=%g)",detP));
     301                }
     302                N[0][0]= 1.0/detM1 * ( a22*b11-a21*b21);
     303                N[0][1]= 1.0/detM1 * ( a22*b21-a21*b22);
     304                N[1][0]= 1.0/detM1 * (-a21*b11+a11*b21);
     305                N[1][1]= 1.0/detM1 * (-a21*b21+a11*b22);
     306
     307                //now, we must find the eigen values and vectors of N
     308                eigen2(&N[0][0],lambda,vp);
     309
     310                //Now that we have P=vp, we can recompute M = M1 inter M2
     311                //
     312                //      -T [ max(lambda1, mu1)          0         ]  -1
     313                // M = P   [                                      ] P
     314                //         [        0            max(lambda2, mu2)]
     315
     316                //get det(P)
     317                detP=vp[0][0]*vp[1][1]-vp[0][1]*vp[1][0];
     318                if (!detP){
     319                        printf("vp= [%g %g ; %g %g]\n",vp[0][0],vp[0][1],vp[1][0],vp[1][1]);
     320                        throw ErrorException(__FUNCT__,exprintf("One of the metric has a zero eigen value or cannot be inverted... (detP=%g)",detM1));
     321                }
     322                invP[0][0]= 1.0/detP * ( vp[1][1]);
     323                invP[0][1]= 1.0/detP * (-vp[0][1]);
     324                invP[1][0]= 1.0/detP * (-vp[1][0]);
     325                invP[1][1]= 1.0/detP * ( vp[0][0]);
     326
     327                //Compute lambdai and mui using Rayleigh formula: lambdai = ei' M1 ei
     328                lambda1=vp[0][0]*(a11*vp[0][0]+a21*vp[1][0]) + vp[1][0]*(a21*vp[0][0]+a22*vp[1][0]);
     329                lambda2=vp[0][1]*(a11*vp[0][1]+a21*vp[1][1]) + vp[1][1]*(a21*vp[0][1]+a22*vp[1][1]);
     330                mu1    =vp[0][0]*(b11*vp[0][0]+b21*vp[1][0]) + vp[1][0]*(b21*vp[0][0]+b22*vp[1][0]);
     331                mu2    =vp[0][1]*(b11*vp[0][1]+b21*vp[1][1]) + vp[1][1]*(b21*vp[0][1]+b22*vp[1][1]);
     332
     333                //finally compute M
     334                M.a11=invP[0][0]*Max(lambda1,mu1)*invP[0][0] + invP[1][0]*Max(lambda2,mu2)*invP[1][0];
     335                M.a21=invP[0][0]*Max(lambda1,mu1)*invP[0][1] + invP[1][0]*Max(lambda2,mu2)*invP[1][1];
     336                M.a22=invP[0][1]*Max(lambda1,mu1)*invP[0][1] + invP[1][1]*Max(lambda2,mu2)*invP[1][1];
     337                *pM=M;
     338        }
     339        /*}}}1*/
     340        /*FUNCTION SimultaneousMatrixReduction_bamg{{{1*/
     341        void SimultaneousMatrixReduction_bamg( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
     342                double a11=M1.a11,a21=M1.a21,a22=M1.a22;
     343                double b11=M2.a11,b21=M2.a21,b22=M2.a22;
     344                //  M1 v = l M2 v
     345                // (M1 - l M2) v =0
     346                // det (M1 - l M2) =0
     347                // det (M1 - l M2) = a l^2 + b l + c;
     348                // = (a11 - l * b11) * (a22 - l * b22) - (a21 - l * b21 ) ^2
     349                const double /*c11 = a11*a11,*/ c21= a21*a21;
     350                const double /*d11 = b11*b11,*/ d21= b21*b21;
     351                const double a=b11*b22 - d21;
     352                const double b=-a11*b22-a22*b11+2*a21*b21;
     353                const double c=-c21+a11*a22;
     354                const double bb = b*b,ac= a*c;
     355                const double delta = bb - 4 * ac;
     356                if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) ){
     357                        // racine double;
     358                        if (Abs(a) < 1.e-30 )
     359                         l1 = l2 = 0;
     360                        else
     361                         l1=l2=-b/(2*a);
     362                        V= D2xD2(1,0,0,1);
     363                }
     364                else {
     365                        const double delta2 = sqrt(delta);
     366                        l1= (-b - delta2)/(2*a);
     367                        l2= (-b + delta2)/(2*a);
     368                        // M1 v = l M2 v
     369                        //  ( (M1 - I M2) x,y)  = (x,(M1 - I M2) y) \forall y
     370                        // so Ker((M1 - I M2)) = Im((M1 - I M2))^\perp
     371                        double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
     372                        double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
     373                        double vp1x,vp1y,vp2x,vp2y;
     374
     375                        if(s1 < s0)
     376                         s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
     377                        else
     378                         s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;
     379
     380                        v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
     381                        s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
     382                        if(s1 < s0)
     383                         s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
     384                        else
     385                         s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
     386                        V=D2xD2(vp1x,vp2x,vp1y,vp2y);
     387                }
     388                return;
     389        }
     390        /*}}}1*/
    235391        /*FUNCTION abscisseInterpole{{{1*/
    236392        Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim) {

  • issm/trunk/src/c/Bamgx/objects/Triangles.cpp

    r2924 r2928  
    20212021                                        //Get eigen values and vectors of Miv
    20222022                                        MatVVP2x2 Vp(Miv);
    2023                                         if(iv==2){
    2024                                                 printf("Hessien' = [%g %g %g]\n",dxdx_vertex[iv]*ci,dxdy_vertex[iv]*ci,dydy_vertex[iv]*ci);
    2025                                                 Vp.Echo();
    2026                                         }
    20272023
    20282024                                        //move eigen valuse to their absolute values
    … …  
    20372033                                        if(iv==2){
    20382034                                                MVp.Echo();
    2039                                                 vertices[iv].m.Echo();
    20402035                                        }
    20412036
    … …  
    20442039                                        if(iv==2){
    20452040                                                vertices[iv].m.Echo();
     2041                                                throw ErrorException(__FUNCT__,exprintf(""));
    20462042                                        }
    2047                                 }
    2048 
    2049                                 for(i=0;i<3;i++){
    2050                                         vertices[i].m.Echo();
    20512043                                }
    20522044                        }//for all fields
Note: See TracChangeset for help on using the changeset viewer.