Strange output on assembly with stable dolfin-1.2.0

Question

Hi,

I have posted this question as a bug, apparently it was not. I was looking to resent changes on numbering, I would like to understand it's connection with the following test:

I am trying to assemble the following form:

element = FiniteElement("Lagrange", triangle, 1) 
u = TrialFunction(element) 
v = TestFunction(element) 
F = as_vector([1, 1]) 
a = inner(F,grad(u)) * v * dx

Then in C++ file I call:

UnitSquare mesh(1, 1); 
Test::FunctionSpace V(mesh); 
Test::BilinearForm a(V, V); 
Matrix A; 
assemble(A, a); 
info(A, true);

With dolfin-1.2.0 I get:

row 0: (0, 0) (1, 0.166667) (2, -0.166667) 
row 1: (0, 0) (1, 0.333333) (2, -0.333333) (3, 0) 
row 2: (0, 0) (1, 0.333333) (2, -0.333333) (3, 0) 
row 3: (1, 0.166667) (2, -0.166667) (3, 0)

while with dolfin-1.0.0+ I get:

row 0: (0, -0.333333) (1, 0) (2, 0) (3, 0.333333) 
row 1: (0, -0.166667) (1, 0) (3, 0.166667) 
row 2: (0, -0.166667) (2, 0) (3, 0.166667)
row 3: (0, -0.333333) (1, 0) (2, 0) (3, 0.333333)

It is not clear to me what is going on. The first matrix is singular and has zero entities and completely different than the second one.

Can someone clarify what is going on?
Thanks,
m

Answer 1

From DOLFIN 1.1.0 release, DOFs of CG spaces (and others) are being reordered to make sparsity graphs less connected. This is the difference of two matrices. Note that this should cause no trouble to user as he is supposed to avoid direct access to DOFs if possible.

Indeed, matrices should be singular because if $z(x, y)$ is solution of $$Lz=y$$ for $y\neq0$ and $L$ being differential operator associated with your form $a$, then $z(y, x)$ is also solution.

Comments

How can one get access to entity A(i,j) of the matrix that gives the same answer as dolfin-1.0?

---

As I said, user is not normally supposed to fiddle with entries of FE tensors. Nevertheless, for vertex-based space (as is CG1) you may print all the entries indexed by vertices

dof_to_vertex_map = V.dofmap().dof_to_vertex_map(mesh)
print A.array()[dof_to_vertex_map][:, dof_to_vertex_map]

In DOLFIN 1.0 you get the same output by just

print A.array()

Now I see that you use C++ interface which means you don't have Matrix::array() method available. You need to use some Matrix::get method to get values and transform their indexing using dof_to_vertex_map.