src/INTERP_KERNEL/GenMathFormulae.hxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #ifndef __GENMATHFORMULAE_HXX__
  21 #define __GENMATHFORMULAE_HXX__
  22
  23 #include "InterpKernelException.hxx"
  24
  25 #include <cmath>
  26
  27 namespace INTERP_KERNEL
  28 {
  29   /*!
  30    * This method computes eigenvalues of a 3x3 symetric matrix stored with 6 values in 'matrix'. The convension chosen for 'matrix' is described here:
  31    * matrix[0]=m_xx, matrix[1]=m_yy, matrix[2]=m_zz,
  32    * matrix[3]=m_xy, matrix[4]=m_yz, matrix[5]=m_xz
  33    * This method returns the 3 eigenvalues in 'eigenVals'.
  34    */
  35   void computeEigenValues6(const double *matrix, double *eigenVals)
  36   {
  37     double tr=(matrix[0]+matrix[1]+matrix[2])/3.;
  38     double K[6]={matrix[0]-tr,matrix[1]-tr,matrix[2]-tr,matrix[3],matrix[4],matrix[5]};
  39     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.;
  40     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]);
  41     p/=6.;
  42     double sqp=sqrt(p);
  43     double tmp=p*sqp;
  44     double phi;
  45     if(fabs(q)<=fabs(tmp))
  46       phi=1./3.*acos(q/tmp);
  47     else
  48       phi=0.;
  49     if(phi<0.)
  50       phi+=M_PI/3.;
  51     eigenVals[0]=tr+2.*sqp*cos(phi);
  52     eigenVals[1]=tr-sqp*(cos(phi)+sqrt(3.)*sin(phi));
  53     eigenVals[2]=tr-sqp*(cos(phi)-sqrt(3.)*sin(phi));
  54   }
  55
  56   /*!
  57    * This method computes one eigenvector of a 3x3 symetric matrix stored with 6 values in 'matrix'. The convension chosen for 'matrix' is described here:
  58    * matrix[0]=m_xx, matrix[1]=m_yy, matrix[2]=m_zz,
  59    * matrix[3]=m_xy, matrix[4]=m_yz, matrix[5]=m_xz
  60    * This method returns the eigenvector of the corresponding eigenvalue in 'eigenVal'. The returned eigenValue is normalized.
  61    */
  62   void computeEigenVectorForEigenValue6(const double *matrix, double eigenVal, double eps, double *eigenVector) throw(INTERP_KERNEL::Exception)
  63   {
  64     //if(fabs(eigenVal)>eps)
  65       {
  66         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};
  67         for(int i=0;i<3;i++)
  68           {
  69             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.};
  70             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];
  71             if(fabs(det)>eps)
  72               {
  73                 eigenVector[0]=(w[1]*w[5]-w[4]*w[2])/det;
  74                 eigenVector[1]=(w[2]*w[3]-w[0]*w[5])/det;
  75                 eigenVector[2]=(w[0]*w[4]-w[1]*w[3])/det;
  76                 double norm=sqrt(eigenVector[0]*eigenVector[0]+eigenVector[1]*eigenVector[1]+eigenVector[2]*eigenVector[2]);
  77                 eigenVector[0]/=norm;
  78                 eigenVector[1]/=norm;
  79                 eigenVector[2]/=norm;
  80                 return;
  81               }
  82           }
  83       }
  84       //else
  85       {
  86         eigenVector[0]=0.;
  87         eigenVector[1]=0.;
  88         eigenVector[2]=0.;
  89         return;
  90       }
  91       //throw INTERP_KERNEL::Exception("computeEigenVector : Do not succed in finding eigen vector !");
  92   }
  93 }
  94
  95 #endif