Box and BoxMesh Produced different Matrices

You were probably relying on a specific dof ordering, but the dof ordering is not guaranteed. DOLFIN 1.1.0 reorders the dof map by default for (a) solver efficiency and (b) to discourage users from relying on the dof ordering because it can change. As part of point (b), not making assumptions on the dof ordering dramatically increases the likelihood that user serial code will run in parallel without modification.

This is because we have reordered the dofs in 1.1. They now do not
follow vertices. A method which tabulate the map between vertices and
dofs is added to trunk and will be backported to 1.1, in a coupled of
hours.

It is tested for Functions and vector but not for Matrices. But it
should be possible to also map matrices.

map = V.dofmap().vertex_to_dof_map(mesh)
A_np = numpy.zeros((mesh.num_vertices(), mesh.num_vertices()))
B_np = A_np.copy()
A_np[map,:] = A.array()
B_np[map,:] = B.array()
A_np[:,map] = A_np[:].copy()
B_np[:,map] = B_np[:].copy()

I have no clue if the above works but I think it will once we have got
vertex_to_dof_map tested and pushed.

Johan

On 01/28/2013 10:55 PM, Oluwaseun Sharomi wrote:
> New question #220386 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/220386
>
> I have the code below that I used with FEniCs version 1.0.0
>
> from dolfin import *
> import numpy
> import scipy.io
> # Create mesh and define function space
> mesh=Box(0, 0, 0, 2, 0.7, 0.3, 20, 7, 3);
> V = FunctionSpace(mesh, 'CG', 1)
> # Define variational problem
> u = TrialFunction(V)
> v = TestFunction(V)
> chi = 1400.0
> C = 1.0
> sigmaij = as_matrix(((1.334177215189873/chi, 0, 0), (0, 0.176061776061776/chi, 0), (0, 0, 0.176061776061776/chi)))
>
> a1 = C*u*v*dx
> a2 = -inner(nabla_grad(v),sigmaij*nabla_grad(u))*dx
>
> A = assemble(a1)
> B = assemble(a2)
> A = A.array();
> B = B.array();
> msh = mesh.coordinates()[:]
> scipy.io.savemat("Nie3D.mat", {"A": A, "B": B, "msh": msh})
>
>
> and Now with the recent version, Box does not work in version 1.1.0. I changed it to BoxMesh. Now the code looks like
>
> from dolfin import *
> import numpy
> import scipy.io
> # Create mesh and define function space
> mesh=BoxMesh(0, 0, 0, 2, 0.7, 0.3, 20, 7, 3);
> V = FunctionSpace(mesh, 'CG', 1)
> # Define variational problem
> u = TrialFunction(V)
> v = TestFunction(V)
> chi = 1400.0
> C = 1.0
> sigmaij = as_matrix(((1.334177215189873/chi, 0, 0), (0, 0.176061776061776/chi, 0), (0, 0, 0.176061776061776/chi)))
>
> a1 = C*u*v*dx
> a2 = -inner(nabla_grad(v),sigmaij*nabla_grad(u))*dx
> A = assemble(a1)
> B = assemble(a2)
> A = A.array();
> B = B.array();
> msh = mesh.coordinates()[:]
> scipy.io.savemat("Nie3D.mat", {"A": A, "B": B, "msh": msh})
>
>
> The issue with my code is that Matrix A from version 1.0.0 is different from version 1.1.0.
>
> The problem is that my code does not work properly anymore.
>
>

I will give it a try and get back to you. Thank you.

I am still not getting the same Matrix I got from version 1.0.0 even by adding the code below to my code. I don't think this is ordering issue. I wish I could attach files so that you can see a spy image of the previous mesh and the new once.

map = V.dofmap().dofs()
A_np = numpy.zeros((mesh.num_vertices(), mesh.num_vertices()))
B_np = A_np.copy()
A_np[map,:] = A.array()
B_np[map,:] = B.array()
A_np[:,map] = A_np[:].copy()
B_np[:,map] = B_np[:].copy()

What type of re-ordering Algorithm is used?

1) That was not what I wrote.
2) The mapping is not pushed to trunk.
3) I will give you a notice when it is.

Johan

On 01/29/2013 06:11 AM, Oluwaseun Sharomi wrote:
> Question #220386 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/220386
>
> Status: Answered => Open
>
> Oluwaseun Sharomi is still having a problem:
> I am still not getting the same Matrix I got from version 1.0.0 even by
> adding the code below to my code. I don't think this is ordering issue.
> I wish I could attach files so that you can see a spy image of the
> previous mesh and the new once.
>
> map = V.dofmap().dofs()
> A_np = numpy.zeros((mesh.num_vertices(), mesh.num_vertices()))
> B_np = A_np.copy()
> A_np[map,:] = A.array()
> B_np[map,:] = B.array()
> A_np[:,map] = A_np[:].copy()
> B_np[:,map] = B_np[:].copy()
>

DofMap.vertex_to_dof_map is now pushed to both trunk and the 1.1.x
branch. You should now be able to try:

map = V.dofmap().vertex_to_dof_map(mesh)
A_np = numpy.zeros((mesh.num_vertices(), mesh.num_vertices()))
B_np = A_np.copy()
A_np[map,:] = A.array()
B_np[map,:] = B.array()
A_np[:,map] = A_np.copy()
B_np[:,map] = B_np.copy()

Again, not sure it works as intended for matrices, but I think so.

Johan

On 01/29/2013 08:26 AM, Johan Hake wrote:
> Question #220386 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/220386
>
> Status: Open => Answered
>
> Johan Hake proposed the following answer:
> 1) That was not what I wrote.
> 2) The mapping is not pushed to trunk.
> 3) I will give you a notice when it is.
>
> Johan
>
> On 01/29/2013 06:11 AM, Oluwaseun Sharomi wrote:
>> Question #220386 on FEniCS Project changed:
>> https://answers.launchpad.net/fenics/+question/220386
>>
>> Status: Answered => Open
>>
>> Oluwaseun Sharomi is still having a problem:
>> I am still not getting the same Matrix I got from version 1.0.0 even by
>> adding the code below to my code. I don't think this is ordering issue.
>> I wish I could attach files so that you can see a spy image of the
>> previous mesh and the new once.
>>
>> map = V.dofmap().dofs()
>> A_np = numpy.zeros((mesh.num_vertices(), mesh.num_vertices()))
>> B_np = A_np.copy()
>> A_np[map,:] = A.array()
>> B_np[map,:] = B.array()
>> A_np[:,map] = A_np[:].copy()
>> B_np[:,map] = B_np[:].copy()
>>
>

What is the easiest way to download and install trunk or branch version 1.1.x.

It works fine. Thank you.