Back to index

salome-med  6.5.0
GenMathFormulae.hxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #ifndef __GENMATHFORMULAE_HXX__
00021 #define __GENMATHFORMULAE_HXX__
00022 
00023 #include "InterpKernelException.hxx"
00024 
00025 #include <cmath>
00026 
00027 namespace INTERP_KERNEL
00028 {
00035   void computeEigenValues6(const double *matrix, double *eigenVals)
00036   {
00037     double tr=(matrix[0]+matrix[1]+matrix[2])/3.;
00038     double K[6]={matrix[0]-tr,matrix[1]-tr,matrix[2]-tr,matrix[3],matrix[4],matrix[5]};
00039     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.;
00040     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]);
00041     p/=6.;
00042     double sqp=sqrt(p);
00043     double tmp=p*sqp;
00044     double phi;
00045     if(fabs(q)<=fabs(tmp))
00046       phi=1./3.*acos(q/tmp);
00047     else
00048       phi=0.;
00049     if(phi<0.)
00050       phi+=M_PI/3.;
00051     eigenVals[0]=tr+2.*sqp*cos(phi);
00052     eigenVals[1]=tr-sqp*(cos(phi)+sqrt(3.)*sin(phi));
00053     eigenVals[2]=tr-sqp*(cos(phi)-sqrt(3.)*sin(phi));
00054   }
00055   
00062   void computeEigenVectorForEigenValue6(const double *matrix, double eigenVal, double eps, double *eigenVector) throw(INTERP_KERNEL::Exception)
00063   {
00064     //if(fabs(eigenVal)>eps)
00065       {
00066         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};
00067         for(int i=0;i<3;i++)
00068           {
00069             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.};
00070             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];
00071             if(fabs(det)>eps)
00072               {
00073                 eigenVector[0]=(w[1]*w[5]-w[4]*w[2])/det;
00074                 eigenVector[1]=(w[2]*w[3]-w[0]*w[5])/det;
00075                 eigenVector[2]=(w[0]*w[4]-w[1]*w[3])/det;
00076                 double norm=sqrt(eigenVector[0]*eigenVector[0]+eigenVector[1]*eigenVector[1]+eigenVector[2]*eigenVector[2]);
00077                 eigenVector[0]/=norm;
00078                 eigenVector[1]/=norm;
00079                 eigenVector[2]/=norm;
00080                 return;
00081               }
00082           }
00083       }
00084       //else
00085       {
00086         eigenVector[0]=0.;
00087         eigenVector[1]=0.;
00088         eigenVector[2]=0.;
00089         return;
00090       }
00091       //throw INTERP_KERNEL::Exception("computeEigenVector : Do not succed in finding eigen vector !");
00092   }
00093 }
00094 
00095 #endif