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

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

#include "../Bamgx.h"
#include "../QuadTree.h"
#include "../SetOfE4.h"

#include "../../shared/shared.h"
#include "../../include/macros.h"
#include "../../toolkits/toolkits.h"

#undef __FUNCT__ 
#define __FUNCT__ "MatVVP2x2"

namespace bamg {

        /*Constructor*/
        /*FUNCTION MatVVP2x2::MatVVP2x2(const MetricAnIso M){{{1*/
        MatVVP2x2::MatVVP2x2(const MetricAnIso M){
                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 ( n2 < 1e-30) 
                 lambda1=lambda2=0,v.x=1,v.y=0;
                else if (delta < eps*n2)
                  { 
                        lambda1=lambda2=-b/2, v.x=1,v.y=0;
                  }
                else 
                  {  //    ---  construction  de 2 vecteur dans (Im ( A - D(i) Id) ortogonal 
                        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*/

}