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)
In [105]:
P11 = assemble(p11)
bc.apply(P11)
I = assemble(i)
In [59]:
P = PETSc.Mat()
P.create()
# P.setBlockSizes(2,2)
Out[59]:
In [60]:
u_is = PETSc.IS().createGeneral(range(V.dim()))
p_is = PETSc.IS().createGeneral(range(V.dim(),V.dim()+Q.dim()))
In [65]:
P.setBlockSize(2)
# P.setBlockSizes?
# AA.view()
In [100]:
P.setType(P.Type.BAIJ)
In [101]:
P.setBlockSize(2)
# P.setValuesBlocked(u_is,u_is,P11)
In [102]:
P.setValuesBlocked(u_is,u_is,P11)
In [104]:
PP = BlockMatrix(2,2)
In [107]:
PP[0,0] = P11
PP[1,1] = I
In [110]:
PETScBaseMatrix.