Access assemble(a) matrix values and change them

Question

Suppose I have 2 PDE's, say:

a1 = inner(grad(u), grad(v))dx
b1 = 5v*dx

and

a2 = inner(25grad(u), grad(v))dx
b2 = 20vdx

Now I assemble these as follows:
A1 = assemble(a1)
L1 = assemble(b1)
A2 = assemble(a2)
L2 = assemble(b2)

Now, I want to access the values stored in matrices A1, A2 and form a new matrix A using A1 and A2. This A is finally to be used to solve PDE for u.

for ex: solve (A, u.vector(), L1)

Is this possible..?

Answer 1

Yes.
You have set, setrow, add, etc functions. See e.g. the class GenericMatrix.
You may also access the underlying matrices in e.g. PETSc or Epetra
format (depending on your backend) using the function mat().

Be aware of the sparsity structure of the matrix.

Comments

Hey Andre,
Thanks for the reply....
I have made a sample code to explain what I am trying to do... Though it is not giving any compilation error. I wonder if it is actually doing what is intended. Can you please go through it once and tell if I am right. Also, is there a way to improve it to reduce the compilation time...

from dolfin import*
import numpy

nx = ny = 20
mesh = RectangleMesh(0, 0, 0.02, 0.02, nx, ny)
#plot(mesh)
V = FunctionSpace(mesh, "CG", 1)

Tc = Constant(300)
Tp = Constant(200)

#Define Boundary Conditions
boundary = compile_subdomains('x[0] <= 0.002 && x[1] <= 0.002')
bc = DirichletBC(V, Tp, boundary)

# Initial condition
Tinitial = 'x[0] <= 0.002 && x[1] <= 0.002 ? Tp : Tc'
Temp = Expression(Tinitial, Tp = Tp, Tc = Tc)
T_1 = interpolate(Temp, V)

T = TrialFunction(V)
v = TestFunction(V)

#define Parameters
dt = 1    # time step

a1 = 3.6*1E6*T*v*dx + dt*inner(0.5*grad(T), grad(v))*dx
L1 = 3.6*1E6*T_1*v*dx
a2 = 1E6*(880 - 3.21*T_1)*T*v*dx + dt*inner(15.98*grad(T), grad(v))*dx
L2 = 1E6*(880 - 3.21*T_1)*T_1*v*dx

A1 = None
b1 = None
A2 = None
b2 = None

T = Function(V)

duration = 50         # total simulation time
t = dt           # initialise time t

A = assemble(a1)
b = assemble(L1)

while t <= duration:
    print 'time =', t
    A1 = assemble(a1)
    b1 = assemble(L1)
    A2 = assemble(a2)
    b2 = assemble(L2)
    T_1mat = T_1.vector().array()
    for i in range(0, len(T_1mat)):
        if(T_1mat[i] >= 250):
            A.array()[i] = (A1.array()[i])
            b.array()[i] = (b1.array()[i])
        else:
            A.array()[i] = (A2.array()[i])
            b.array()[i] = (b2.array()[i])
    bc.apply(A, b)
    solve(A, T.vector(), b)
    T_1.assign(T)
    t += dt 

print T.compute_vertex_values().max(), T_1.compute_vertex_values().min()
plot (T, title = 'T plot')
interactive()

Thanks again.

---

The 'array' function gives you a dense numpy array containing
copies of the values of the matrix. Changing the values
of the numpy array will not change the values in the underlying matrix.

If you want to change the values using Python, you may consider
the example in PyAMG:
http://code.google.com/p/pyamg/source/browse/branches/2.0.x/Examples/DolfinFormulation/demo.py
Though, I have not checked if data returns a copy or a view of the underlying
values. Note that this will also depend on the linear algebra backend you
use. You may also use PyTrilinos or PyPETSc to access the matrices.