Merge pull request KratosMultiphysics#813 from KratosMultiphysics/fea…

…ture-adding-eigensolver-factory-correcting-eigen-solvers

Correcting and improving eigen solvers. Adding a eigen_solver_factory

applications/ExternalSolversApplication/custom_python/add_linear_solvers_to_python.cpp

@@ -68,12 +68,6 @@ void  AddLinearSolversToPython()
 
     using namespace boost::python;
 
-    class_<TLinearSolverType<std::complex<double>>, TLinearSolverType<std::complex<double>>::Pointer, boost::noncopyable>(
-        "ComplexLinearSolver").def(self_ns::str(self));
-    class_<TDirectSolverType<std::complex<double>>,
-           TDirectSolverType<std::complex<double>>::Pointer,
-           bases<TLinearSolverType<std::complex<double>>>,
-           boost::noncopyable>("ComplexDirectSolver").def(self_ns::str(self));
 
     //***************************************************************************
     //linear solvers

applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_eigensolver.py

@@ -72,19 +72,8 @@ def _create_linear_solver(self):
         This overrides the base class method and replaces the usual linear solver
         with an eigenvalue problem solver.
         """
-        if self.eigensolver_settings["solver_type"].GetString() == "FEAST":
-            feast_system_solver_settings = self.eigensolver_settings["linear_solver_settings"]
-            if feast_system_solver_settings["solver_type"].GetString() == "skyline_lu":
-                # default built-in feast system solver
-                linear_solver = ExternalSolversApplication.FEASTSolver(self.eigensolver_settings)
-            elif feast_system_solver_settings["solver_type"].GetString() == "pastix":
-                feast_system_solver = ExternalSolversApplication.PastixComplexSolver(feast_system_solver_settings)
-                linear_solver = ExternalSolversApplication.FEASTSolver(self.eigensolver_settings, feast_system_solver)
-            else:
-                raise Exception("Unsupported FEAST system solver_type: " + feast_system_solver_settings["solver_type"].GetString())
-        else:
-            raise Exception("Unsupported eigensolver solver_type: " + self.eigensolver_settings["solver_type"].GetString())
-        return linear_solver
+        import eigen_solver_factory
+        return eigen_solver_factory.ConstructSolver(self.eigensolver_settings)
 
     def _create_mechanical_solver(self):
         eigen_scheme = self.get_solution_scheme() # The scheme defines the matrices of the eigenvalue problem.

applications/StructuralMechanicsApplication/tests/eigen_test/Eigen_3D3N_Thin_Circle_test_parameters.json

@@ -31,7 +31,7 @@
             "lambda_max": 4.0e5,
             "search_dimension": 15,
             "linear_solver_settings":{
-                "solver_type": "skyline_lu"
+                "solver_type": "complex_skyline_lu_solver"
             }
 	},
         "problem_domain_sub_model_part_list" : ["Parts_plate"],

applications/StructuralMechanicsApplication/tests/eigen_test/Eigen_Q4_Thick_2x2_Plate_test_parameters.json

@@ -29,7 +29,7 @@
         "lambda_max": 1.0e4,
         "search_dimension": 9,
         "linear_solver_settings":{
-            "solver_type": "skyline_lu"
+            "solver_type": "complex_skyline_lu_solver"
         }
 	},
         "problem_domain_sub_model_part_list" : ["Parts_plate"],

applications/StructuralMechanicsApplication/tests/eigen_test/Eigen_TL_3D8N_Cube_test_parameters.json

@@ -29,7 +29,7 @@
         "lambda_max": 2.5e3,
         "search_dimension": 6,
         "linear_solver_settings":{
-            "solver_type": "skyline_lu"
+            "solver_type": "complex_skyline_lu_solver"
         }
 	},
         "problem_domain_sub_model_part_list" : ["Parts_Solid"],

kratos/includes/ublas_complex_interface.h

@@ -0,0 +1,79 @@
+//    |  /           | 
+//    ' /   __| _` | __|  _ \   __| 
+//    . \  |   (   | |   (   |\__ `
+//   _|\_\_|  \__,_|\__|\___/ ____/ 
+//                   Multi-Physics  
+//
+//  License:		 BSD License 
+//					 Kratos default license: kratos/license.txt
+//
+//  Main authors:    Michael Andre
+//
+
+#if !defined(KRATOS_UBLAS_COMPLEX_INTERFACE_H_INCLUDED )
+#define  KRATOS_UBLAS_COMPLEX_INTERFACE_H_INCLUDED
+
+// System includes
+#include <string>
+#include <iostream>
+
+// External includes
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+#include <boost/numeric/ublas/vector_sparse.hpp>
+#include <boost/numeric/ublas/vector_expression.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_sparse.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
+#include <boost/numeric/ublas/hermitian.hpp>
+#include <boost/numeric/ublas/banded.hpp>
+#include <boost/numeric/ublas/triangular.hpp>
+#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/operation.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+#include <boost/numeric/ublas/operation_sparse.hpp>
+
+// Project includes
+#include "includes/define.h"
+
+namespace Kratos
+{
+
+///@name Type Definitions
+///@{
+
+using namespace boost::numeric::ublas;
+
+typedef boost::numeric::ublas::vector<std::complex<double>> ComplexVector;
+typedef unit_vector<std::complex<double>> ComplexUnitVector;
+typedef zero_vector<std::complex<double>> ComplexZeroVector;
+typedef scalar_vector<std::complex<double>> ComplexScalarVector;
+typedef mapped_vector<std::complex<double>> ComplexSparseVector;
+typedef compressed_vector<std::complex<double>> ComplexCompressedVector;
+typedef coordinate_vector<std::complex<double>> ComplexCoordinateVector;
+typedef vector_range<ComplexVector> ComplexVectorRange;
+typedef vector_slice<ComplexVector> ComplexVectorSlice;
+typedef matrix<std::complex<double>> ComplexMatrix;
+typedef identity_matrix<std::complex<double>> ComplexIdentityMatrix;
+typedef zero_matrix<std::complex<double>> ComplexZeroMatrix;
+typedef scalar_matrix<std::complex<double>> ComplexScalarMatrix;
+typedef triangular_matrix<std::complex<double>> ComplexTriangularMatrix;
+typedef symmetric_matrix<std::complex<double>> ComplexSymmetricMatrix;
+typedef hermitian_matrix<std::complex<double>> ComplexHermitianMatrix;
+typedef banded_matrix<std::complex<double>> ComplexBandedMatrix;
+typedef mapped_matrix<std::complex<double>> ComplexSparseMatrix;
+typedef compressed_matrix<std::complex<double>> ComplexCompressedMatrix;
+typedef coordinate_matrix<std::complex<double>> ComplexCoordinateMatrix;
+typedef matrix_row<ComplexMatrix> ComplexMatrixRow;
+typedef matrix_column<ComplexMatrix> ComplexMatrixColumn;
+typedef matrix_vector_range<ComplexMatrix> ComplexMatrixVectorRange;
+typedef matrix_vector_slice<ComplexMatrix> ComplexMatrixVectorSlice;
+typedef matrix_range<ComplexMatrix> ComplexMatrixRange;
+typedef matrix_slice<ComplexMatrix> ComplexMatrixSlice;
+
+///@}
+
+}  // namespace Kratos.
+
+#endif // KRATOS_UBLAS_COMPLEX_INTERFACE_H_INCLUDED  defined 

kratos/linear_solvers/linear_solver.h

@@ -205,13 +205,18 @@ class LinearSolver
         return false;
     }
 
-    /** Eigenvalue and eigenvector solve method for derived eigensolvers */
+    /** Eigenvalue and eigenvector solve method for derived eigensolvers 
+     * @param K: The stiffness matrix
+     * @param M: The mass matrix
+     * @param Eigenvalues: The vector containing the eigen values
+     * @param Eigenvectors: The matrix containing the eigen vectors
+     */
     virtual  void Solve(SparseMatrixType& K,
                         SparseMatrixType& M,
                         DenseVectorType& Eigenvalues,
                         DenseMatrixType& Eigenvectors)
     {}
-
+    
     /** Some solvers may require a minimum degree of knowledge of the structure of the matrix. To make an example
      * when solving a mixed u-p problem, it is important to identify the row associated to v and p.
      * another example is the automatic prescription of rotation null-space for smoothed-aggregation solvers