In [103]:
#!/usr/bin/python
import petsc4py
import sys

petsc4py.init(sys.argv)

from petsc4py import PETSc
Print = PETSc.Sys.Print
# from MatrixOperations import *
from dolfin import *
import numpy as np
import matplotlib.pylab as plt
import os
import scipy.io
from PyTrilinos import Epetra, EpetraExt, AztecOO, ML, Amesos
from scipy2Trilinos import scipy_csr_matrix2CrsMatrix
import PETScIO as IO

m = 1000
errL2u = np.zeros((m-1,1))
errL2p = np.zeros((m-1,1))
l2uorder = np.zeros((m-1,1))
l2porder = np.zeros((m-1,1))
NN = np.zeros((m-1,1))
DoF = np.zeros((m-1,1))
Vdim = np.zeros((m-1,1))
Qdim = np.zeros((m-1,1))
Wdim = np.zeros((m-1,1))
iterations = np.zeros((m-1,1))
SolTime = np.zeros((m-1,1))
udiv = np.zeros((m-1,1))
nn = 2

dim = 2
Solving = 'Iterative'
ShowResultPlots = 'no'
ShowErrorPlots = 'no'
EigenProblem = 'no'
SavePrecond = 'no'
case = 3
parameters['linear_algebra_backend'] = 'PETSc'

xx = 4

nn = 2**(xx)
# Create mesh and define function space
nn = int(nn)
NN[xx-1] = nn

mesh = RectangleMesh(0, 0, 1, 1, nn, nn,'right')

parameters['reorder_dofs_serial'] = False
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
parameters['reorder_dofs_serial'] = False
W = V*Q
Vdim[xx-1] = V.dim()
Qdim[xx-1] = Q.dim()
Wdim[xx-1] = W.dim()
print "\n\nV:  ",Vdim[xx-1],"Q:  ",Qdim[xx-1],"W:  ",Wdim[xx-1],"\n\n"
def boundary(x, on_boundary):
    return on_boundary


if case == 1:
    u0 = Expression(("20*x[0]*pow(x[1],3)","5*pow(x[0],4)-5*pow(x[1],4)"))
    p0 = Expression("60*pow(x[0],2)*x[1]-20*pow(x[1],3)")
elif case == 2:
    u0 = Expression(("sin(pi*x[1])","sin(pi*x[0])"))
    p0 = Expression("sin(x[1]*x[0])")
elif case == 3:
    u0 =Expression(("sin(x[1])*exp(x[0])","cos(x[1])*exp(x[0])"))
    p0 = Expression("sin(x[0])*cos(x[1])")
elif case == 4:
    u0 = Expression(("sin(x[1])*exp(x[0])","cos(x[1])*exp(x[0])"))
    p0 = Expression("sin(x[0])*cos(x[1])")




bc = DirichletBC(W.sub(0),u0, boundary)
bcs = [bc]

(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
if case == 1:
    f = Expression(("120*x[0]*x[1]*(1-mu)","60*(pow(x[0],2)-pow(x[1],2))*(1-mu)"), mu = 1e0)
elif case == 2:
    f = Expression(("pi*pi*sin(pi*x[1])+x[1]*cos(x[1]*x[0])","pi*pi*sin(pi*x[0])+x[0]*cos(x[1]*x[0])"))
elif case == 3:
    f = Expression(("cos(x[0])*cos(x[1])","-sin(x[0])*sin(x[1])"))
elif case == 4:
    f = Expression(("cos(x[1])*cos(x[0])","-sin(x[1])*sin(x[0])"))




N = FacetNormal(mesh)
h = CellSize(mesh)
h_avg =avg(h)
alpha = 10.0
gamma =10.0
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg =avg(h)
d = 0
a11 = inner(grad(v), grad(u))*dx
a12 = div(v)*p*dx
a21 = div(u)*q*dx
L1  =  inner(v,f)*dx
a = a11-a12-a21

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

(p) = TrialFunction(Q)
(q) = TestFunction(Q)
p11 = inner(grad(v), grad(u))*dx
i = p*q*dx

tic()
AA, bb = assemble_system(a, L1, bcs)

A = as_backend_type(AA).mat()
print toc()
b = bb.array()
zeros = 0*b
del bb
bb = IO.arrayToVec(b)
x = IO.arrayToVec(zeros)


DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
INFO:FFC:Adjusting missing element domain to Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2).
INFO:FFC:Adjusting element degree from ? to 2
INFO:FFC:Adjusting missing element domain to Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2).
INFO:FFC:Adjusting element degree from ? to 2
INFO:FFC:Adjusting missing element domain to Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2).
INFO:FFC:Adjusting element degree from ? to 2
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.

V:   [ 2178.] Q:   [ 289.] W:   [ 2467.] 


0.0207719802856

In [105]:
P11 = assemble(p11)
bc.apply(P11)
I = assemble(i)


DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:UFL:No integrals left after transformation, returning empty form.
DEBUG:FFC:Reusing form from cache.

In [120]:
P = PETSc.Mat()
P.create()
# P.setBlockSizes(2,2)


Out[120]:
<petsc4py.PETSc.Mat at 0x47b4cb0>

In [60]:
u_is = PETSc.IS().createGeneral(range(V.dim()))
p_is = PETSc.IS().createGeneral(range(V.dim(),V.dim()+Q.dim()))

In [123]:
P.setSizes((2,2),(2,2))
# P.setBlockSizes?
# AA.view()

In [124]:
P.setType(P.Type.BAIJ)

In [101]:
P.setBlockSize(2)
# P.setValuesBlocked(u_is,u_is,P11)

In [137]:
# P.setValueBlockedStencil(u_is,u_is,P11)
Pstencil = P11.Stencil()

In [150]:
print Pstencil


<petsc4py.PETSc._Mat_Stencil object at 0x2366690>

In [107]:
PP[0,0] = P11
PP[1,1] = I

In [166]:
P11.getSubMatrix

In [167]:
PP, Pb = assemble_system(a11+i,L1,bcs)
P = as_backend_type(PP).mat()


Found different Arguments with same counts.
Did you combine test or trial functions from different spaces?
The Arguments found are:
  v_{-1}
  v_{-2}
  v_{-1}
  v_{-2}
ERROR:UFL:Found different Arguments with same counts.
Did you combine test or trial functions from different spaces?
The Arguments found are:
  v_{-1}
  v_{-2}
  v_{-1}
  v_{-2}
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
<ipython-input-167-d6fb65ea74ec> in <module>()
----> 1 PP, Pb = assemble_system(a11+i,L1,bcs)
      2 P = as_backend_type(PP).mat()

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/dolfin/fem/assembling.pyc in assemble_system(A_form, b_form, bcs, x0, A_coefficients, b_coefficients, A_function_spaces, b_function_spaces, cell_domains, exterior_facet_domains, interior_facet_domains, reset_sparsity, add_values, finalize_tensor, keep_diagonal, A_tensor, b_tensor, backend, form_compiler_parameters)
    275     # Wrap forms
    276     A_dolfin_form = Form(A_form, A_function_spaces, A_coefficients,
--> 277                          subdomains, form_compiler_parameters)
    278     b_dolfin_form = Form(b_form, b_function_spaces, b_coefficients,
    279                          subdomains, form_compiler_parameters)

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/dolfin/fem/form.pyc in __init__(self, form, function_spaces, coefficients, subdomains, form_compiler_parameters, common_cell)
     54             self._compiled_form, module, self.form_data, prefix \
     55                                  = jit(form, form_compiler_parameters,
---> 56                                        common_cell)
     57 
     58         elif isinstance(form, ufc.form):

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/dolfin/compilemodules/jit.pyc in mpi_jit(*args, **kwargs)
     58         # Just call JIT compiler when running in serial
     59         if MPI.num_processes() == 1:
---> 60             return local_jit(*args, **kwargs)
     61 
     62         # Compile first on process 0

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/dolfin/compilemodules/jit.pyc in jit(form, form_compiler_parameters, common_cell)
    120         raise RuntimeError, "Form compiler must implement the jit function."
    121 
--> 122     return jit_compile(form, parameters=p, common_cell=common_cell)

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/ffc/jitcompiler.pyc in jit(ufl_object, parameters, common_cell)
     76         return jit_element(ufl_object, parameters)
     77     else:
---> 78         return jit_form(ufl_object, parameters, common_cell)
     79 
     80 def _auto_select_degree(elements):

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/ffc/jitcompiler.pyc in jit_form(form, parameters, common_cell)
    167     # Compute form metadata and extract preprocessed form
    168     form_data = form.compute_form_data(common_cell=common_cell,
--> 169                                        element_mapping=element_mapping)
    170 
    171     # Wrap input

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/ufl/form.pyc in compute_form_data(self, object_names, common_cell, element_mapping)
    122                                          object_names=object_names,
    123                                          common_cell=common_cell,
--> 124                                          element_mapping=element_mapping)
    125         # Always validate arguments, keeping sure that the validation works
    126         self._form_data.validate(object_names=object_names,

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/ufl/algorithms/preprocess.pyc in preprocess(form, object_names, common_cell, element_mapping)
    141     tic('extract_arguments_and_coefficients')
    142     original_arguments, original_coefficients = \
--> 143         extract_arguments_and_coefficients(form)
    144 
    145     tic('build_element_mapping')

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/ufl/algorithms/analysis.pyc in extract_arguments_and_coefficients(a)
    116 Did you combine test or trial functions from different spaces?
    117 The Arguments found are:\n%s""" % "\n".join("  %s" % f for f in arguments)
--> 118         error(msg)
    119 
    120     if len(fcounts) != len(set(fcounts.values())):

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/ufl/log.pyc in error(self, *message)
    152         "Write error message and raise an exception."
    153         self._log.error(*message)
--> 154         raise self._exception_type(self._format_raw(*message))
    155 
    156     def begin(self, *message):

UFLException: Found different Arguments with same counts.
Did you combine test or trial functions from different spaces?
The Arguments found are:
  v_{-1}
  v_{-2}
  v_{-1}
  v_{-2}

In [113]:
PP = BlockMatrix(2,2)
PP[0,0] = P11
PP[1,1] = I


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-113-4dbf51d6eda5> in <module>()
      1 PP = BlockMatrix(2,2)
----> 2 PP[0,0] = P11
      3 PP[1,1] = I

/home/mwathen/Work/FEniCS/lib/python2.7/site-packages/dolfin/cpp/la.pyc in __setitem__(self, t, m)
   4851         self._items[t] = m
   4852         i,j = t
-> 4853         self.set_block(i, j, m)
   4854 
   4855 

TypeError: in method 'BlockMatrix_set_block', argument 4 of type 'boost::shared_ptr< dolfin::GenericMatrix >'

In [169]:
Pdiag = P11.getDiagonal()

In [172]:
invDiag = 1/Pdiag

In [178]:
p =P11.diagonalScale(invDiag)

In [179]:
P11.view()