From 0ad4f3c026473a4cc4873329bc273c95e74662e5 Mon Sep 17 00:00:00 2001
From: ageay <ageay>
Date: Mon, 4 Oct 2010 14:44:47 +0000
Subject: [PATCH] *** empty log message ***

---
 src/INTERP_KERNEL/GenMathFormulae.hxx | 95 +++++++++++++++++++++++++++
 src/INTERP_KERNEL/Makefile.am         |  1 +
 2 files changed, 96 insertions(+)
 create mode 100644 src/INTERP_KERNEL/GenMathFormulae.hxx

diff --git a/src/INTERP_KERNEL/GenMathFormulae.hxx b/src/INTERP_KERNEL/GenMathFormulae.hxx
new file mode 100644
index 000000000..92de33f7c
--- /dev/null
+++ b/src/INTERP_KERNEL/GenMathFormulae.hxx
@@ -0,0 +1,95 @@
+//  Copyright (C) 2007-2010  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
+{
+  /*!
+   * This method computes eigenvalues of a 3x3 symetric matrix stored with 6 values in 'matrix'. The convension chosen for 'matrix' is described here:
+   * matrix[0]=m_xx, matrix[1]=m_yy, matrix[2]=m_zz,
+   * matrix[3]=m_xy, matrix[4]=m_yz, matrix[5]=m_xz
+   * This method returns the 3 eigenvalues in 'eigenVals'.
+   */
+  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));
+  }
+  
+  /*!
+   * This method computes one eigenvector of a 3x3 symetric matrix stored with 6 values in 'matrix'. The convension chosen for 'matrix' is described here:
+   * matrix[0]=m_xx, matrix[1]=m_yy, matrix[2]=m_zz,
+   * matrix[3]=m_xy, matrix[4]=m_yz, matrix[5]=m_xz
+   * This method returns the eigenvector of the corresponding eigenvalue in 'eigenVal'. The returned eigenValue is normalized.
+   */
+  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
diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am
index 67e780fc0..afbe7285d 100644
--- a/src/INTERP_KERNEL/Makefile.am
+++ b/src/INTERP_KERNEL/Makefile.am
@@ -141,6 +141,7 @@ VTKNormalizedUnstructuredMesh.hxx	\
 VTKNormalizedUnstructuredMesh.txx	\
 VectorUtils.hxx                         \
 VolSurfFormulae.hxx                     \
+GenMathFormulae.hxx                     \
 VolSurfUser.hxx                         \
 VolSurfUser.txx                         \
 CurveIntersector.hxx                    \
-- 
2.39.2