src/INTERP_KERNEL/GenMathFormulae.hxx

   1 // Copyright (C) 2007-2021  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, or (at your option) any later version.
   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 // Author : Anthony Geay (CEA/DEN)
  20
  21 #ifndef __GENMATHFORMULAE_HXX__
  22 #define __GENMATHFORMULAE_HXX__
  23
  24 #include "InterpKernelException.hxx"
  25
  26 #include <cmath>
  27
  28 namespace INTERP_KERNEL
  29 {
  30   /*!
  31    * This method computes eigenvalues of a 3x3 symmetric matrix stored with 6 values in 'matrix'. The convention chosen for 'matrix' is described here:
  32    * matrix[0]=m_xx, matrix[1]=m_yy, matrix[2]=m_zz,
  33    * matrix[3]=m_xy, matrix[4]=m_yz, matrix[5]=m_xz
  34    * This method returns the 3 eigenvalues in 'eigenVals'.
  35    *
  36    * https://en.wikipedia.org/wiki/Eigenvalue_algorithm
  37    */
  38   void computeEigenValues6(const double *matrix, double *eigenVals)
  39   {
  40     double tr=(matrix[0]+matrix[1]+matrix[2])/3.;
  41     double K[6]={matrix[0]-tr,matrix[1]-tr,matrix[2]-tr,matrix[3],matrix[4],matrix[5]};
  42     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.;
  43     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]);
  44     p/=6.;
  45     double sqp=sqrt(p);
  46     double tmp=p*sqp;
  47     double phi;
  48     if(fabs(q)<=fabs(tmp))
  49       if(tmp !=0)
  50         phi=1./3.*acos(q/tmp);
  51       else // EDF22176
  52         phi=1./3.;
  53     else
  54       phi=0.;
  55     if(phi<0.)
  56       phi+=M_PI/3.;
  57     eigenVals[0]=tr+2.*sqp*cos(phi);
  58     eigenVals[1]=tr-sqp*(cos(phi)+sqrt(3.)*sin(phi));
  59     eigenVals[2]=tr-sqp*(cos(phi)-sqrt(3.)*sin(phi));
  60   }
  61
  62   /*!
  63    * This method computes one eigenvector of a 3x3 symmetric matrix stored with 6 values in 'matrix'. The convention chosen for 'matrix' is described here:
  64    * matrix[0]=m_xx, matrix[1]=m_yy, matrix[2]=m_zz,
  65    * matrix[3]=m_xy, matrix[4]=m_yz, matrix[5]=m_xz
  66    * This method returns the eigenvector of the corresponding eigenvalue in 'eigenVal'. The returned eigenValue is normalized.
  67    */
  68   void computeEigenVectorForEigenValue6(const double *matrix, double eigenVal, double eps, double *eigenVector)
  69   {
  70     //if(fabs(eigenVal)>eps)
  71       {
  72         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};
  73         for(int i=0;i<3;i++)
  74           {
  75             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.};
  76             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];
  77             if(fabs(det)>eps)
  78               {
  79                 eigenVector[0]=(w[1]*w[5]-w[4]*w[2])/det;
  80                 eigenVector[1]=(w[2]*w[3]-w[0]*w[5])/det;
  81                 eigenVector[2]=(w[0]*w[4]-w[1]*w[3])/det;
  82                 double norm=sqrt(eigenVector[0]*eigenVector[0]+eigenVector[1]*eigenVector[1]+eigenVector[2]*eigenVector[2]);
  83                 eigenVector[0]/=norm;
  84                 eigenVector[1]/=norm;
  85                 eigenVector[2]/=norm;
  86                 return;
  87               }
  88           }
  89       }
  90       //else
  91       {
  92         eigenVector[0]=0.;
  93         eigenVector[1]=0.;
  94         eigenVector[2]=0.;
  95         return;
  96       }
  97       //throw INTERP_KERNEL::Exception("computeEigenVector : Do not succed in finding eigen vector !");
  98   }
  99 }
 100
 101 #endif