How to set a value in a GenericMatrix

Asked by Edwin Mai on 2013-03-18


i want to set a value in a GenericMatrix, but using i following approach doesn't work.

from dolfin import*
import numpy as numpy

n = 10
mesh = UnitSquareMesh(n, n)

V = VectorFunctionSpace(mesh, "CG", 1)
n = FacetNormal(mesh)

trace_n = TrialFunction(V)
phi = TestFunction(V)
a = inner(trace_n,phi)*ds
L = inner(n,phi)*ds
A = assemble(a)
B = assemble(L)
size = A.array().shape[0]
for i in range(0, size-1):
   if A.array()[i,:].max() == 0:
    A.array()[i,i] = 1

trace_n = Function(V)
solve(A, trace_n.vector(), B)

The command A.array()[i,i] = 1 doesn't assign the 1 to the postion (i,i) of the Matrix A.
Does it only works with A.set() command and is there an example for using it in python?

(The meaning of the code is to avoid the ident_zeros() function, because it produces
a known bug, which will be fixed in the next release of dolfin. I know that there exits a
fix by jan blechta, but trying to complie this dolfin development version with just using
cmake ., make install on an ubuntu 12.04 system causes a lot of compilation errors. So I
thought looking for an alternative might be easier.)

Thank you very much

Question information

English Edit question
DOLFIN Edit question
No assignee Edit question
Solved by:
Edwin Mai
Last query:
Last reply:
Jan Blechta (blechta) said : #1

Generally - yes, you should use supplied set() function instead to assignment to numpy array entries as these are only copies.

But - if ident_zeros() does not work on your matrix this won't work either as some diagonal entries are not preallocated.

You could use something that Garth suggested - creating custom sparsity pattern or not finalizing tensor prematurely.

Edwin Mai (edwinmai) said : #2

Ok, i will try this.

Thank you very much!