EigenMetric.cpp

Go to the documentation of this file.

#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);
   }/*}}}*/
 
}