source: issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp@ 3332

#include <cstdio>
#include <cstring>
#include <cmath>
#include <ctime>

#include "Metric.h"

namespace bamg {

        /*Constructor*/
        /*FUNCTION MatVVP2x2::MatVVP2x2(const Metric M){{{1*/
        MatVVP2x2::MatVVP2x2(const Metric M){
                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MatVVP2x2)*/

                double a11=M.a11,a21=M.a21,a22=M.a22;
                const double eps = 1.e-5;
                double c11 = a11*a11, c22 = a22*a22, c21= a21*a21;
                double b=-a11-a22,c=-c21+a11*a22;
                double delta = b*b - 4 * c ;
                double n2=(c11+c22+c21);

                //if norm(M)<10^30 -> M=zeros(2,2)
                if ( n2 < 1e-30) lambda1=lambda2=0,v.x=1,v.y=0;

                //if matrix is already diagonal and has a 0 eigen value
                else if (delta < eps*n2){ 
                        lambda1=lambda2=-b/2, v.x=1,v.y=0;
                }

                //general case: construction of 2 eigen vectors
                else {
                        delta     = sqrt(delta);
                        lambda1   = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0;
                        double v0 = a11-lambda1, v1 = a21,v2 = a22 - lambda1;
                        double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
                        if(s1 < s0)
                         s0=sqrt(s0),v.x=v1/s0,v.y=-v0/s0;
                        else
                         s1=sqrt(s1),v.x=v2/s1,v.y=-v1/s1;
                };
        }
        /*}}}1*/

        /*Methods*/
        /*FUNCTION MatVVP2x2::Echo {{{1*/
        void MatVVP2x2::Echo(void){

                printf("MatVVP2x2:\n");
                printf("   lambda1: %g\n",lambda1);
                printf("   lambda2: %g\n",lambda2);
                printf("   v.x: %g\n",v.x);
                printf("   v.y: %g\n",v.y);

                return;
        }
        /*}}}*/

}