In [21]:

    

from dolfin import *
import numpy as np
import matplotlib.pylab as plt
from scipy.sparse import csc_matrix, spdiags
mesh = UnitSquareMesh(2,2)
parameters['reorder_dofs_serial'] = False
V = VectorFunctionSpace(mesh,"CG",2)
Q = FunctionSpace(mesh,"CG",1)
W = V*Q
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
print V.dim(), Q.dim()


    

    




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.






    




50 9

In [22]:

    

def boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(W.sub(0),Expression(("0","0")), boundary)
parameters['linear_algebra_backend'] = 'uBLAS'
BQB = assemble(inner(u,v)*dx+ div(v)*p*dx+div(u)*q*dx)
bc.apply(BQB)
BQB=BQB.sparray()


    

    




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 [23]:

    

Q = BQB[0:np.max(range(V.dim()))+1,0:np.max(range(V.dim()))+1]


    

In [24]:

    

Qdiag = Q.diagonal()


    

In [25]:

    

B = BQB[np.max(range(V.dim()))+1:np.max(range(W.dim()))+1,0:np.max(range(V.dim()))+1]
Bt = BQB[0:np.max(range(V.dim()))+1,np.max(range(V.dim()))+1:np.max(range(W.dim()))+1]


    

In [27]:

    

d = spdiags(1.0/Qdiag, 0, len(Qdiag), len(Qdiag))
B*D*BT


    

In [19]:

    

row_divide_col_reduce(Bt,Qdiag)


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-19-a24459b1f67b> in <module>()
----> 1 row_divide_col_reduce(Bt,Qdiag)

<ipython-input-18-8686e71a03d2> in row_divide_col_reduce(b, x)
      1 def row_divide_col_reduce(b, x):
      2     data = b.data.copy() / np.take(x, b.indices)
----> 3     ret = sps.csc_matrix((data, b.indices.copy(), b.indptr.copy()),
      4                          shape=b.shape)
      5     return ret.sum(axis=1)

NameError: global name 'sps' is not defined