Merge branch 'FER_11000' into amdadev

Benjamin Renard
2 parents 4933833e da266937
Showing 6 changed files with 47 additions and 112 deletions Show diff stats
CMakeLists.txt
cmake/modules/FindEigen.cmake
src/ExternLib/Tetrahedron/characteristic/CMakeLists.txt
src/ExternLib/Tetrahedron/common/tetrahedron.cc
src/ExternLib/Tetrahedron/elongation/CMakeLists.txt
src/ExternLib/Tetrahedron/planarity/CMakeLists.txt
@@ -76,12 +76,15 @@ find_package(NetCDF REQUIRED)
 find_package(DDClientCpp REQUIRED)
 find_package(Geopack REQUIRED)
 find_package(Maglib REQUIRED)
+find_package(Eigen REQUIRED)
 find_package(InternalField REQUIRED)
 find_package(Con2020 REQUIRED)
 find_package(KT17 REQUIRED)
 find_package(Cainlib REQUIRED)
 find_package(CSlim)
  
+MESSAGE( STATUS  "FER ${EIGENINCLUDE_DIR}")
+
 get_filename_component(DDCLIENTLIB_DIR ${DDCLIENTLIBRARIES} PATH)
  
 # Apply configuration
@@ -0,0 +1,16 @@
+FIND_PATH(EIGENINCLUDE_DIR Eigen
+  HINTS
+    ${USERLOCAL_ROOT}/include/eigen3/Eigen
+  NO_DEFAULT_PATH
+)
+mark_as_advanced(EIGENINCLUDE_DIR)
+
+
+
+include(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(Eigen DEFAULT_MSG EIGENINCLUDE_DIR)
+
+if(EIGEN_FOUND)
+   set(EIGENINCLUDE_DIRS ${EIGENINCLUDE_DIR})
+endif()
+
@@ -5,6 +5,7 @@ set(LIBRARY_OUTPUT_PATH ${PLUGIN_OUTPUT_PATH}/characteristic)
  
 include_directories(
 	../common/
+	${EIGENINCLUDE_DIR}
 )
  
 #Library configuration
@@ -2,13 +2,30 @@
 #include<iostream>
 #include<cmath>
 #include<algorithm>
+#include <Dense>
+#include <vector>
+
+using Eigen::MatrixXd;
+using Eigen::SelfAdjointEigenSolver;
  
 using Vector = std::vector<double>;
 using Matrix = std::vector<Vector>;
  
-#ifdef DEBUG
-void print_matrix(Matrix);
-#endif
+
+Vector getEigenvalues(const Matrix& matrix) {
+    MatrixXd A(3, 3);
+    A << matrix[0][0], matrix[0][1], matrix[0][2],
+         matrix[1][0], matrix[1][1], matrix[1][2],
+         matrix[2][0], matrix[2][1], matrix[2][2];
+
+    SelfAdjointEigenSolver<MatrixXd> solver(A);
+    Vector eigenvalues(3);
+    eigenvalues[0] = solver.eigenvalues()[0];
+    eigenvalues[1] = solver.eigenvalues()[1];
+    eigenvalues[2] = solver.eigenvalues()[2];
+
+    return eigenvalues;
+}
  
 Matrix concatenate_vectors(Vector p1, Vector p2, Vector p3, Vector p4) {
 	Matrix mat;
@@ -46,7 +63,6 @@ Matrix matrix_multiply(Matrix mat1, Matrix mat2) {
 	Matrix mat2t = transpose(mat2);
 	Matrix r;
 	int n = mat1.size();
-	int m = mat1[0].size();
 	for(int i=0;i<n;i++) {
 		// new row
 		Vector row;
@@ -59,87 +75,6 @@ Matrix matrix_multiply(Matrix mat1, Matrix mat2) {
  
 }
  
-#ifdef DEBUG
-void print_vector(Vector vect) {
-	for(int i=0;i<(int)vect.size();i++) std::cout << vect[i] << " ";
-	std::cout << std::endl;
-}
-
-void print_matrix(Matrix mat) {
-	for(int i=0;i<(int)mat.size();i++) {
-		print_vector(mat[i]);
-	}
-}
-#endif
-
-Vector vector_substract(Vector vect1, Vector vect2) {
-	Vector ans;
-	for(int i=0;i<(int)vect1.size();i++) {
-		ans.push_back(vect1[i] - vect2[i]);
-	}
-	return ans;
-}
-
-Vector vector_projection(Vector v, Vector u) {
-	Vector ans;
-	double a = scalar_product(u,v);
-	double b = scalar_product(u,u);
-	for(int i=0;i<(int)v.size();i++) {
-		ans.push_back((a/b) * u[i]);
-	}
-	return ans;
-}
-
-Vector vector_normalize(Vector vect) {
-	Vector ans;
-	double norm = scalar_product(vect, vect);
-	for(int i=0;i<(int)vect.size();i++) ans.push_back(vect[i] / sqrt(norm));
-	return ans;
-}
-
-Matrix gram_schmidt(Matrix mat) {
-	Matrix matt = transpose(mat);
-	Matrix basis;
-	Matrix Us;
-	for(int i=0;i<(int)mat.size();i++) {
-		// new basis vector
-		Vector nb = matt[i];
-		for(int j=0;j<=i-1;j++) {
-			nb = vector_substract(nb, vector_projection(matt[i], Us[j]));
-			Vector vv = vector_projection(matt[i], Us[j]);
-		}
-		Us.push_back(nb);
-		nb = vector_normalize(nb);
-		basis.push_back(vector_normalize(nb));
-	}
-	return transpose(basis);
-}
-
-bool is_upper_triangular(Matrix mat) {
-	for(int i=0;i<(int)mat.size();i++) {
-		for(int j=0;j<i;j++) {
-			if(mat[i][j] > 1e-16) return false;
-		}
-	}
-	return true;
-}
-
-Vector matrix_eigenvalues(Matrix mat) {
-	Matrix A = mat;
-	for(int k=0;k<1000000;k++) {
-		Matrix basis = gram_schmidt(A);
-		Matrix basisT = transpose(basis);
-		Matrix t1 = matrix_multiply(A, basis);
-		Matrix t2 = matrix_multiply(basisT, t1);
-		if(is_upper_triangular(t2)) break;
-
-		A = t2;
-	}
- 	Vector eig;
-	for(int i=0;i<(int)mat.size();i++) eig.push_back(A[i][i]);
-	return eig;
-}
-
 void compute_tetrahedron(Vector p1, Vector p2, Vector p3, Vector p4, double& elongation, double& planarity, double& characteristic) {
 	// make matrix by concatenating vectors
 	Matrix mat = concatenate_vectors(p1,p2,p3,p4);
@@ -149,35 +84,13 @@ void compute_tetrahedron(Vector p1, Vector p2, Vector p3, Vector p4, double&amp; elo
 	// matrix product
 	Matrix pp = matrix_multiply(matt, mat);
  
-	Vector eig = matrix_eigenvalues(pp);
-
+	// computing eigenvalues
+	Vector eig = getEigenvalues(pp);
+ 
 	// sort eigenvalues
 	sort(eig.begin(), eig.end());
  
 	planarity  = 1. - sqrt(eig[0] / eig[1]);
 	elongation = 1. - sqrt(eig[1] / eig[2]);
 	characteristic = 2 * sqrt(eig[2]);
-}
-
-#ifdef DEBUG
-Vector random_vector() {
-	Vector v;
-	for(int i=0;i<3;i++) v.push_back((double)(rand()%7));
-
-	return v;
-}
-
-int main(int argc, char ** argv) {
-	// create some random data and compute the elongation
-	for(int i=0;i<1000000;i++) {
-		// generate twelve random numbers
-		Vector p1 = random_vector();
-		Vector p2 = random_vector();
-	     	Vector p3 = random_vector();
-	       	Vector p4 = random_vector();
-		std::cout << "Elongation " << elongation(p1,p2,p3,p4) << ", planarity: " << planarity(p1,p2,p3,p4) << std::endl;
-	}
-	return 0;
-}
-#endif
-
+}
 \ No newline at end of file
@@ -4,7 +4,8 @@ PROJECT(elongation)
 set(LIBRARY_OUTPUT_PATH ${PLUGIN_OUTPUT_PATH}/elongation)
  
 include_directories(
-	../common
+	../common/
+	${EIGENINCLUDE_DIR}
 )
  
 #Library configuration
@@ -5,6 +5,7 @@ set(LIBRARY_OUTPUT_PATH ${PLUGIN_OUTPUT_PATH}/planarity)
  
 include_directories(
 	../common/
+	${EIGENINCLUDE_DIR}
 )
  
 #Library configuration