Z0 = mu0*c0                               def femmesh_2_dolfin_mesh(femmesh_file):
    """Convert a femmesh file to a dolfin mesh object"""
    femmesh_reader = FemmeshReader(femmesh_file)
    return MeshConverters.femmesh_reader_2_dolfin_mesh(femmesh_reader) def dolfin_ublassparse_to_scipy_csr ( A, dtype=None, imagify=False ):
    """
    convert a DOLFIN uBLASSparseMatrix to a scipy.sparse.csr_matrix()
    
    @param A: a DOLFIN uBLASSparseMatrix
    @param dtype: the numpy data type to use to store the matrix
    @param imagify: multiply the original matrix data by 1j
    """
    import scipy.sparse
    # get the sparse data from the input matrix
    (row,col,data) = A.data()   # get sparse data
    col = np.intc(col)
    row = np.intc(row)
    n = A.size(0)
    if imagify: data = data*1j
    A_sp = scipy.sparse.csr_matrix( (data,col,row), shape=(n,n), dtype=dtype)
    
    return A_sp class PointCurrentSource(CurrentSource):
    def set_position(self, position):
    def set_value(self, value):
    def get_contribution(self):
def calc_pointsource_contrib(V, source_coords, source_value): class BiCGStabSolver ( SystemSolverBase ):
    """
    An iterative Stabilised BICG solver using scipy.
    """    
    def _call_solver (self):
        """Solves the linear system (self._A)x = self._b and returns the solution vector.
        
        @return: The solution to the linear system.
        """      
        return scipy.sparse.linalg.bicgstab(self._A, self._b, M=self._M, callback=self._callback ) class GMRESSolver ( SystemSolverBase ):
    """
    An iterative GMRES solver using scipy.
    """
    def _call_solver (self):
        """Solves the linear system (self._A)x = self._b and returns the solution vector.
        
        @return: The solution to the linear system.
        """
        return scipy.sparse.linalg.gmres(self._A, self._b, M=self._M, callback=self._callback ) def calculate_residual ( A, x, b ):
    """
    Calculate the residual of the system Ax = b
    
    @param A: a matrix
    @param x: a vector
    @param b: a vector 
    """
    return np.linalg.norm( A*x - b.reshape(x.shape) ) class UMFPACKSolver ( SystemSolverBase ):
    """
    A direct UMFPACK-based solver provided by scipy
    """
    def _call_solver (self):
        """Solves the linear system (self._A)x = self._b and returns the solution vector.
        
        @return: The solution to the linear system.
        """
        scipy.sparse.linalg.use_solver(useUmfpack=True)
        return scipy.sparse.linalg.spsolve ( self._A, self._b ), 0
    
    def plot_convergence (self, x_is_time=False, show_plot=False, label=None, style='-'):
        """Output a string indicating that no convergence history is available.
        """
        print "Direct solver has no convergence history" class PyAMGSolver ( SystemSolverBase ):
    """
    A PyAMG-based iterative solver.
    """
    def _call_solver (self):
        """Solves the linear system (self._A)x = self._b and returns the solution vector.
        
        @return: The solution to the linear system.
        """
        import pyamg
        x = pyamg.solve ( self._A, self._b, verb=True, tol=1e-8, maxiter=800 )
        return x, 0
    
    def plot_convergence (self, x_is_time=False, show_plot=False, label=None, style='-'):
        """Output a string indicating that no convergence history is available.
        """
        print "PAMG solver convergence display is not yet implemented" def load_scipy_matrix_from_mat ( path, name ):
    """
    Load a scipy sparse matrix from a .mat file.
    
    @param path: the folder in which the matrix is saved
    @param name: the filename of the matrix
    """
    import scipy.io
    import scipy.sparse
    
    filename =  os.path.join ( path, name)
    
    if not os.path.exists( filename + '.mat' ):
        return None
    
    data = scipy.io.loadmat ( filename )[name]
    
    if type(data) is np.ndarray:
        matrix = data
    else:
        matrix = data.tocsr ()
         
    return matrix def save_scipy_matrix_as_mat ( path, name, matrix ):
    """
    Save a scipy sparse matrix as a .mat file.
    
    @param path: the folder in which the matrix is to be saved
    @param name: the filename of the matrix
    @param matrix: the scipy matrix to save
    """
    import scipy.io
    if not check_path ( path ):
        return False
    
    scipy.io.savemat ( os.path.join(path, name), {name: matrix }, oned_as='column' )
    
    return True
from dolfin import dot, cross, curl, inner, dx, ds
from sucemfem.Consts import eps0, mu0, c0, Z0
from sucemfem.Utilities.MeshIO import femmesh_2_dolfin_mesh
from sucemfem.Utilities.Converters import dolfin_ublassparse_to_scipy_csr
from sucemfem.Sources import point_source
from sucemfem.Utilities.LinalgSolvers import BiCGStabSolver, GMRESSolver, calculate_residual, UMFPACKSolver, PyAMGSolver
from sucemfem.Utilities.MatrixIO import ( load_scipy_matrix_from_mat, save_scipy_matrix_as_mat)
import dolfin
import numpy as np
import os
import sys
## Copyright (C) 2011 Stellenbosch University
##
## This file is part of SUCEM.
##
## SUCEM is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## SUCEM 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 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 
__author__ = "Evan Lezar"

"""
A simple driver to test various LA solvers.
"""



sys.path.insert(0, '../../')
del sys.path[0]

def generate_matrices ( mesh_file, order ):
    """
    Generate the matrices for a 1\lambda spherical domain with an infintesimal dipole at the origin.
    
    @param mesh_file: the full path and filename to the spherical mesh to be used (this is a femmesh file).
    @todo: change the mesh type to a gmsh mesh.
    """
    problem_data = {'l': 2.9979245800000002e-4,
                    'I': 1.0,
                    'f': 1.0e9,
                    'source_coord': np.array([0,0,0.]),
                    'order': order,
                    }
    
    lam = c0/problem_data['f']
    
    # load the mesh and scale coordinates accordingly
    mesh = femmesh_2_dolfin_mesh(mesh_file)
    mesh.init()
    mesh.coordinates()[:] *= lam
    
    
    eps_r = 1;
    mu_r = 1;
    freq = problem_data['f'];
    source_coord = problem_data['source_coord']
    source_point = dolfin.Point(*source_coord)
    source_value = np.array([0,0,1.])*problem_data['I']*problem_data['l']
    
    V = dolfin.FunctionSpace(mesh, "Nedelec 1st kind H(curl)", problem_data['order'])
    print 'DOFs: ', V.dim()

    # Define basis and bilinear form
    k_0 = 2*np.pi*freq/c0

    u = dolfin.TrialFunction(V)
    v = dolfin.TestFunction(V)

    m = eps_r*inner(v, u)*dx                # Mass form
    s = (1/mu_r)*dot(curl(v), curl(u))*dx   # Stiffness form

    n = V.cell().n
    s_0 = inner(cross(n, v), cross(n, u))*ds

    def boundary(x, on_boundary):
        return on_boundary

    # Assemble forms
    print 'assembling forms'
    M = dolfin.uBLASSparseMatrix()
    S = dolfin.uBLASSparseMatrix()
    S_0 = dolfin.uBLASSparseMatrix()
    dolfin.assemble(m, tensor=M, mesh=mesh)
    dolfin.assemble(s, tensor=S, mesh=mesh)
    dolfin.assemble(s_0, tensor=S_0, mesh=mesh)

    # Set up RHS
    b = np.zeros(M.size(0), dtype=np.complex128)
    dofnos, rhs_contrib = point_source.calc_pointsource_contrib(
        V, source_coord, source_value)


    rhs_contrib = 1j*k_0*Z0*rhs_contrib