In [ ]:

    

from dolfin import *
import numpy as np
import scipy.sparse as sp
import numpy
from scipy2Trilinos import scipy_csr_matrix2CrsMatrix
from PyTrilinos import Epetra, ML, AztecOO

nn = 8
mesh = RectangleMesh(0, 0, 1, 1, nn, nn,'left')

order  = 2
Magnetic = FunctionSpace(mesh, "N1curl", order)
Lagrange = FunctionSpace(mesh, "CG", order)
L= FunctionSpace(mesh, "DG", order-1)
C = sp.csr_matrix((Magnetic.dim(),Lagrange.dim()))

parameters['linear_algebra_backend'] = 'uBLAS'

u = Function(Lagrange)
for i in range(0, Lagrange.dim() ):
    zero = np.zeros((Lagrange.dim(),1))[:,0]
    zero[i] = 1
    u.vector()[:] = zero

    uu = grad(u)
    Pv = project(uu,Magnetic)
    P = Pv.vector().array()

    index = P.nonzero()
    index = index[0]
    for x in range(0,len(index)):
        if np.abs(P[index[x]]) < 1e-3:
            P[index[x]] = 0
    print P.shape, C.shape
    pn = P.nonzero()[0]
    for j in range(0,len(pn)):
        C[pn[j],i] = P[pn[j]]
    del P


def boundary(x, on_boundary):
    return on_boundary



c = 1
(u) = TrialFunction(Magnetic)
(v) = TestFunction(Magnetic)
(p) = TrialFunction(Lagrange)
(q) = TestFunction(Lagrange)
a = dot(curl(u),curl(v))*dx + c*inner(u, v)*dx

l = (inner(grad(p),grad(q))*dx)
u0 = Expression(("sin(2*pi*x[1])*cos(2*pi*x[0])","-sin(2*pi*x[0])*cos(2*pi*x[1])"))
f = 8*pow(pi,2)*u0+c*u0
L1  = inner(v, f)*dx


bc = DirichletBC(Magnetic,u0, boundary)

# MLList = {

#     "default values" : "maxwell",
#     "max levels" : 1,
#     "output" : 10,
#     "PDE equations" : 1,
#     "increasing or decreasing" : "decreasing",
#     "aggregation: type" : "Uncoupled-MIS",
#     "aggregation: damping factor" : 1.3333,
#     "coarse: max size" : 75,
#     "aggregation: threshold" : 0.0,
#     "smoother: sweeps" : 2,
#     "smoother: damping factor" : 0.67,
#     "smoother: type" : "MLS",
#     "smoother: MLS polynomial order" : 4,
#     "smoother: pre or post" : "both",
#     "coarse: type" : "Amesos-KLU",
#     "prec type" : "MGV",
#     "print unused" : -2
# }

MLList = {
    "default values" : "maxwell",
    "max levels"                                     : 10,
    "prec type"                                        : "MGV",
    "increasing or decreasing"               : "decreasing",
    "aggregation: type"                          : "Uncoupled-MIS",
    "aggregation: damping factor"         : 4.0/3.0,
    "eigen-analysis: type"                      : "cg",
    "eigen-analysis: iterations"              : 10,
    "smoother: sweeps"                          : 1,
    "smoother: damping factor"              : 1.0,
    "smoother: pre or post"                     : "both",
    "smoother: type"                               : "Hiptmair",
    "subsmoother: type"                         : "Chebyshev",
    "subsmoother: Chebyshev alpha"    : 27.0,
    "subsmoother: node sweeps"           : 4,
    "subsmoother: edge sweeps"           : 4,
    "PDE equations" : 1,
    "coarse: type"                                   : "Amesos-MUMPS",
    "coarse: max size"                           : 128

}

comm = Epetra.PyComm()



Acurl,b = assemble_system(a,L1,bc)
Anode = assemble(l)
Acurl = Acurl.sparray()
Anode = Anode.sparray()
b = b.array()