source: issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/EigenMetric.cpp@ 21759

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

#include "Metric.h"
#include "../shared/shared.h"

namespace bamg {

        /*Constructor*/
        EigenMetric::EigenMetric(const Metric& M){/*{{{*/
                /*From a metric (a11,a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/

                /*Intermediaries*/
                double a11=M.a11,a21=M.a21,a22=M.a22;
                double normM;
                double delta,b;

                /*To get the eigen values, we must solve the following equation:
                 *     | a11 - lambda    a21        |
                 * det |                            | = 0
                 *     | a21             a22-lambda |
                 *
                 * We have to solve the following polynom:
                 *  lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/

                /*Compute polynom determinant*/
                b=-a11-a22;
                delta=b*b - 4*(a11*a22-a21*a21);

                /*Compute norm of M to avoid round off errors*/
                normM=a11*a11 + a22*a22 + a21*a21;

                /*1: normM too small: eigen values = 0*/
                if(normM<1.e-30){
                        lambda1=0;
                        lambda2=0;
                        vx=1;
                        vy=0;
                }
                /*2: delta is small -> double root*/
                else if (delta < 1.e-5*normM){
                        lambda1=-b/2;
                        lambda2=-b/2;
                        vx=1;
                        vy=0;
                }
                /*3: general case -> two roots*/
                else{
                        delta     = sqrt(delta);
                        lambda1   = (-b-delta)/2.0;
                        lambda2   = (-b+delta)/2.0;

                        /*Now, one must find the eigen vectors. For that we use the following property of the inner product
                         *    <Ax,y> = <x,tAy>
                         * Here, M'(M-lambda*Id) is symmetrical, which gives:
                         *    ∀(x,y)∈R²xR² <M'x,y> = <M'y,x>
                         * And we have the following:
                         *    if y∈Ker(M'), ∀x∈R² <M'x,y> = <x,M'y> = 0
                         * We have shown that
                         *    Im(M') ⊥ Ker(M')
                         *
                         * To find the eigen vectors of M, we only have to find two vectors
                         * of the image of M' and take their perpendicular as long as they are
                         * not 0.
                         * To do that, we take the images (1,0) and (0,1):
                         *  x1 = (a11 - lambda)      x2 = a21
                         *  y1 = a21                 y2 = (a22-lambda)
                         *
                         * We take the vector that has the larger norm and take its perpendicular.*/

                        double norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21; 
                        double norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1);

                        if (norm2<norm1){
                                norm1=sqrt(norm1);
                                vx = - a21/norm1;
                                vy = (a11-lambda1)/norm1;
                        }
                        else{
                                norm2=sqrt(norm2);
                                vx = - (a22-lambda1)/norm2;
                                vy = a21/norm2;
                        }
                }

        }
        /*}}}*/
        EigenMetric::EigenMetric(double r1,double r2,const D2& vp1){/*{{{*/
                this->lambda1 = r1;
                this->lambda2 = r2;
                this->vx = vp1.x;
                this->vy = vp1.y;
        }/*}}}*/

        /*Methods*/
        void   EigenMetric::Abs(){/*{{{*/
                lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);
        }/*}}}*/
        double EigenMetric::Aniso2() const  { /*{{{*/
                return lmax()/lmin();
        }/*}}}*/
        void EigenMetric::Echo(void){/*{{{*/

                _printf_("EigenMetric:\n");
                _printf_("   lambda1: " << lambda1 << "\n");
                _printf_("   lambda2: " << lambda2 << "\n");
                _printf_("   vx: " << vx << "\n");
                _printf_("   vy: " << vy << "\n");

                return;
        }
        /*}}}*/
        double EigenMetric::hmin() const {/*{{{*/
                return sqrt(1/bamg::Max3(lambda1,lambda2,1e-30));
        }/*}}}*/
        double EigenMetric::hmax() const {/*{{{*/
                return sqrt(1/bamg::Max(bamg::Min(lambda1,lambda2),1e-30));
        }/*}}}*/
        double EigenMetric::lmax() const {/*{{{*/
                return bamg::Max3(lambda1,lambda2,1e-30);
        }/*}}}*/
        double EigenMetric::lmin() const {/*{{{*/
                return bamg::Max(bamg::Min(lambda1,lambda2),1e-30);
        }/*}}}*/
        void   EigenMetric::Min(double a) { /*{{{*/
                lambda1=bamg::Min(a,lambda1); lambda2=bamg::Min(a,lambda2) ;
        }/*}}}*/
        void   EigenMetric::Max(double a) { /*{{{*/
                //change eigen values
                lambda1=bamg::Max(a,lambda1); lambda2=bamg::Max(a,lambda2) ;
        }/*}}}*/
        void   EigenMetric::Minh(double h) {/*{{{*/
                Min(1.0/(h*h));
        }/*}}}*/
        void   EigenMetric::Maxh(double h) {/*{{{*/
                //Call Max function
                Max(1.0/(h*h));
        }/*}}}*/
        void   EigenMetric::pow(double p){/*{{{*/
                lambda1=::pow(lambda1,p);lambda2=::pow(lambda2,p);
        }/*}}}*/

}