# Path: sucemfem/ProblemConfigurations/EMVectorWaveEigenproblem.py
# class EigenProblem(EMProblem):
#     FormCombiner = CombineForms        
# 
# Path: sucemfem/ProblemConfigurations/EMVectorWaveEigenproblem.py
# class DefaultEigenSolver(object):
#     def set_eigenproblem(self, eigenproblem):
#         """Sets initialised instance of EigenProblem to solve"""
#         self.eigenproblem = eigenproblem
# 
#     def set_sigma(self, sigma):
#         """Spectrum shift (sigma) to apply to k^2"""
#         self.sigma = sigma
# 
#     def solve_problem(self, nev, ncv=None):
#         """Solve problem and return the eigenvalues and eigenvectors
#         
#         @param nev: Number of eigenpairs to compute
#         @keyword ncv: Number of Arnoldi basisvectors to use. 
#             (default: 2*nev+1).
#         
#         @rtype: (C{numpy.array}, C{numpy.array})
#         @return: (eigs_w, eigs_v) -- A tupple consisting of the eigenvalues and eigenvectors of the eigensystem.
#             The eigenvalues are returned as an array of n_eig values corresponding to k^2 for the mode's resonant wavenumber.
#             The eigenvectors are returned as a 2D array of shape (n_eig, problem_dim), with row i corresponding to the 
#             modal distributions associated with the i-th eigenvalue. 
#         """
#         M = self.eigenproblem.system_matrices['M']
#         S = self.eigenproblem.system_matrices['S']
#         solve_mat = S - self.sigma*M
#         lu = dolfin.LUSolver(solve_mat)
#         lu.parameters["reuse_factorization"] = True
#         lu.parameters["report"] = False
#         bb = dolfin.Vector(M.size(0))
#         xx = dolfin.Vector(M.size(0))        
#         def sigma_solve(b):
#             bb[:] = b
#             lu.solve(xx, bb)                                    
#             return xx[:]
#         M_matvec = lambda x: M*x                    
#               
#         #speigs in ARPACK has been removed in scipy 0.9/0.10
#         #eigs is now in scipy.sparse.linalg.arpack
#         #therefore force M_matvec to have "shape" and "dtype" attributes to comply with "eigs" requirements
#         #also perform the sigma_solve function for spectrum shift to sigma
#         class RM:
#             def __init__(self, M):
#                 self.M = M
#                 self.shape = (M.size(0), M.size(1))
#                 self.dtype = np.dtype('d')                              
#             def matvec(self, x):                                
#                 return sigma_solve(self.M*x)
#                                
# #        eigs_w, eigs_v = speigs.ARPACK_gen_eigs(
# #            M_matvec, sigma_solve, M.size(0), self.sigma, nev, ncv=ncv)
#                                             
#         eigs_w, eigs_v= arpack.eigs(RM(M), k=nev, sigma=self.sigma, which='LM', ncv=ncv)
#       
#         return eigs_w, eigs_v.T
# 
# Path: sucemfem/Consts.py
# Z0 = mu0*c0                              


import sys
import numpy as N
import os
import dolfin as dol
from sucemfem.ProblemConfigurations.EMVectorWaveEigenproblem import EigenProblem
from sucemfem.ProblemConfigurations.EMVectorWaveEigenproblem import DefaultEigenSolver
from sucemfem.Consts import c0
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with SUCEM. If not, see <http://www.gnu.org/licenses/>. 
##
## Contact: cemagga@gmail.com 
# Authors:
# Neilen Marais <nmarais@gmail.com>
from __future__ import division

sys.path.append('../../')




# Define mesh
# mesh = dol.UnitCube(1,1,1)
# mesh.coordinates()[:] *= [cdims.a,cdims.b,cdims.c]
#mesh_file = 'lee_mittra92_fig6b.xml'
#mesh_file = 'lee_mittra92_fig6c.xml'
mesh_file = '../../examples/albani_bernardi74/mesh/albani_bernardi74_fig2VII.xml'
materials_mesh_file = "%s_physical_region%s" % (os.path.splitext(mesh_file))
mesh = dol.Mesh(mesh_file)
material_mesh_func = dol.MeshFunction('uint', mesh, materials_mesh_file)
materials = {1000:dict(eps_r=16),
             1001:dict(eps_r=1)}
order = 3
sigma = 1.5

# Set up eigen problem

# Path: sucemfem/ProblemConfigurations/EMVectorWaveEigenproblem.py
# class EigenProblem(EMProblem):
#     FormCombiner = CombineForms        
# 
# Path: sucemfem/ProblemConfigurations/EMVectorWaveEigenproblem.py
# class DefaultEigenSolver(object):
#     def set_eigenproblem(self, eigenproblem):
#         """Sets initialised instance of EigenProblem to solve"""
#         self.eigenproblem = eigenproblem
# 
#     def set_sigma(self, sigma):
#         """Spectrum shift (sigma) to apply to k^2"""
#         self.sigma = sigma
# 
#     def solve_problem(self, nev, ncv=None):
#         """Solve problem and return the eigenvalues and eigenvectors
#         
#         @param nev: Number of eigenpairs to compute
#         @keyword ncv: Number of Arnoldi basisvectors to use. 
#             (default: 2*nev+1).
#         
#         @rtype: (C{numpy.array}, C{numpy.array})
#         @return: (eigs_w, eigs_v) -- A tupple consisting of the eigenvalues and eigenvectors of the eigensystem.
#             The eigenvalues are returned as an array of n_eig values corresponding to k^2 for the mode's resonant wavenumber.
#             The eigenvectors are returned as a 2D array of shape (n_eig, problem_dim), with row i corresponding to the 
#             modal distributions associated with the i-th eigenvalue. 
#         """
#         M = self.eigenproblem.system_matrices['M']
#         S = self.eigenproblem.system_matrices['S']
#         solve_mat = S - self.sigma*M
#         lu = dolfin.LUSolver(solve_mat)
#         lu.parameters["reuse_factorization"] = True
#         lu.parameters["report"] = False
#         bb = dolfin.Vector(M.size(0))
#         xx = dolfin.Vector(M.size(0))        
#         def sigma_solve(b):
#             bb[:] = b
#             lu.solve(xx, bb)                                    
#             return xx[:]
#         M_matvec = lambda x: M*x                    
#               
#         #speigs in ARPACK has been removed in scipy 0.9/0.10
#         #eigs is now in scipy.sparse.linalg.arpack
#         #therefore force M_matvec to have "shape" and "dtype" attributes to comply with "eigs" requirements
#         #also perform the sigma_solve function for spectrum shift to sigma
#         class RM:
#             def __init__(self, M):
#                 self.M = M
#                 self.shape = (M.size(0), M.size(1))
#                 self.dtype = np.dtype('d')                              
#             def matvec(self, x):                                
#                 return sigma_solve(self.M*x)
#                                
# #        eigs_w, eigs_v = speigs.ARPACK_gen_eigs(
# #            M_matvec, sigma_solve, M.size(0), self.sigma, nev, ncv=ncv)
#                                             
#         eigs_w, eigs_v= arpack.eigs(RM(M), k=nev, sigma=self.sigma, which='LM', ncv=ncv)
#       
#         return eigs_w, eigs_v.T
# 
# Path: sucemfem/Consts.py
# Z0 = mu0*c0                              


# Path: sandbox/eig/eig_driver.py
import sys
import numpy as N
import os
import dolfin as dol
from sucemfem.ProblemConfigurations.EMVectorWaveEigenproblem import EigenProblem
from sucemfem.ProblemConfigurations.EMVectorWaveEigenproblem import DefaultEigenSolver
from sucemfem.Consts import c0

## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with SUCEM. If not, see <http://www.gnu.org/licenses/>. 
##
## Contact: cemagga@gmail.com 
# Authors:
# Neilen Marais <nmarais@gmail.com>
from __future__ import division

sys.path.append('../../')




# Define mesh
# mesh = dol.UnitCube(1,1,1)
# mesh.coordinates()[:] *= [cdims.a,cdims.b,cdims.c]
#mesh_file = 'lee_mittra92_fig6b.xml'
#mesh_file = 'lee_mittra92_fig6c.xml'
mesh_file = '../../examples/albani_bernardi74/mesh/albani_bernardi74_fig2VII.xml'
materials_mesh_file = "%s_physical_region%s" % (os.path.splitext(mesh_file))
mesh = dol.Mesh(mesh_file)
material_mesh_func = dol.MeshFunction('uint', mesh, materials_mesh_file)
materials = {1000:dict(eps_r=16),
             1001:dict(eps_r=1)}
order = 3
sigma = 1.5

# Set up eigen problem

ep = EigenProblem()