src/INTERP_KERNEL/GenMathFormulae.hxx

Go to the documentation of this file.

// Copyright (C) 2007-2011  CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//

#ifndef __GENMATHFORMULAE_HXX__
#define __GENMATHFORMULAE_HXX__

#include "InterpKernelException.hxx"

#include <cmath>

namespace INTERP_KERNEL
{
  void computeEigenValues6(const double *matrix, double *eigenVals)
  {
    double tr=(matrix[0]+matrix[1]+matrix[2])/3.;
    double K[6]={matrix[0]-tr,matrix[1]-tr,matrix[2]-tr,matrix[3],matrix[4],matrix[5]};
    double q=(K[0]*K[1]*K[2]+2.*K[4]*K[5]*K[3]-K[0]*K[4]*K[4]-K[2]*K[3]*K[3]-K[1]*K[5]*K[5])/2.;
    double p=K[0]*K[0]+K[1]*K[1]+K[2]*K[2]+2*(K[3]*K[3]+K[4]*K[4]+K[5]*K[5]);
    p/=6.;
    double sqp=sqrt(p);
    double tmp=p*sqp;
    double phi;
    if(fabs(q)<=fabs(tmp))
      phi=1./3.*acos(q/tmp);
    else
      phi=0.;
    if(phi<0.)
      phi+=M_PI/3.;
    eigenVals[0]=tr+2.*sqp*cos(phi);
    eigenVals[1]=tr-sqp*(cos(phi)+sqrt(3.)*sin(phi));
    eigenVals[2]=tr-sqp*(cos(phi)-sqrt(3.)*sin(phi));
  }
  
  void computeEigenVectorForEigenValue6(const double *matrix, double eigenVal, double eps, double *eigenVector) throw(INTERP_KERNEL::Exception)
  {
    //if(fabs(eigenVal)>eps)
      {
        const double m9[9]={matrix[0]-eigenVal,matrix[3],matrix[5],matrix[3],matrix[1]-eigenVal,matrix[4],matrix[5],matrix[4],matrix[2]-eigenVal};
        for(int i=0;i<3;i++)
          {
            double w[9]={m9[0+3*i],m9[1+3*i],m9[2+3*i],m9[0+(3*(i+1))%6],m9[1+(3*(i+1))%6],m9[2+(3*(i+1))%6],1.,1.,1.};
            double det=w[0]*w[4]*w[8]+w[1]*w[5]*w[6]+w[2]*w[3]*w[7]-w[0]*w[5]*w[7]-w[1]*w[3]*w[8]-w[2]*w[4]*w[6];
            if(fabs(det)>eps)
              {
                eigenVector[0]=(w[1]*w[5]-w[4]*w[2])/det;
                eigenVector[1]=(w[2]*w[3]-w[0]*w[5])/det;
                eigenVector[2]=(w[0]*w[4]-w[1]*w[3])/det;
                double norm=sqrt(eigenVector[0]*eigenVector[0]+eigenVector[1]*eigenVector[1]+eigenVector[2]*eigenVector[2]);
                eigenVector[0]/=norm;
                eigenVector[1]/=norm;
                eigenVector[2]/=norm;
                return;
              }
          }
      }
      //else
      {
        eigenVector[0]=0.;
        eigenVector[1]=0.;
        eigenVector[2]=0.;
        return;
      }
      //throw INTERP_KERNEL::Exception("computeEigenVector : Do not succed in finding eigen vector !");
  }
}

#endif