FEniCS : How to interpolate data at vertices of (3D) cells?

Hello everyone,
I hope that one of you guys can help me because i have been stuck here for a week.
I am trying to read a gmsh file (.msh) using dolfin convert to XML and then download it with dolfin.

The thing is when i assemble the stiffness matrix and the mass matrix and try to generate an exemple (MonteCarlo) using these matrices , i have figured out that they are not the good ones .

I think that dolfin when downloaded the XML file which has already convert it from MSH , he changes the ordering of the numbering of the elements in ascending way.
I am sure that he does this, because i have tried my code on a box generated by BoxMesh of fenics and my program worked and if i check mesh.cells() and my XML file i can see that the ordering had been changed, So that’s why i am not having my good matrices

And to be honest, i don’t know how to fix this thing? Does dolfin destroys the ordering of any mesh and build another ordering (in an ascending way) ? how do you guys get the good stiffness matrix and mass one when you download a mesh of your type?

I have read the same problem founded by another person but still cannot understand how to fix that.

Thank you very much in advance

this is my code

from dolfin import *
import os
import numpy as np
from numpy import linalg
from math import exp
from math import sqrt
from math import pi
from math import log
import matplotlib.pyplot as plt
import scipy
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
import scipy.stats as stats
from scipy.optimize import minimize

alpha=1.22
beta_x=2
beta_y=0.01 #0.05
beta_z=0.01#0.05
mu=-4.25

str_os = ‘dolfin-convert ./mesh/calcul_little_specimen.msh ./calcul_little_specimen.xml’
os.system(str_os)
mesh = Mesh(“calcul_little_specimen.xml”)

V = FunctionSpace(mesh, “Lagrange”, 1)
u = TrialFunction(V)
v = TestFunction(V)
parameters[‘linear_algebra_backend’] = ‘Eigen’

S=as_matrix([[beta_x**2,0,0],[0,0,0],[0,0,0]])
kx=inner(dot(S,grad(u)),grad(v))*dx
Kx=assemble(kx)
row,col,val = as_backend_type(Kx).data()
Kx_sp=scipy.sparse.csc_matrix((val,col,row))

S=as_matrix([[0,0,0],[0,beta_y**2,0],[0,0,0]])
ky=inner(dot(S,grad(u)),grad(v))*dx
Ky=assemble(ky)
row,col,val = as_backend_type(Ky).data()
Ky_sp=scipy.sparse.csc_matrix((val,col,row))

S=as_matrix([[0,0,0],[0,0,0],[0,0,beta_z**2]])
kz=inner(dot(S,grad(u)),grad(v))*dx
Kz=assemble(kz)
row,col,val = as_backend_type(Kz).data()
Kz_sp=scipy.sparse.csc_matrix((val,col,row))

m=u *v* dx
M = assemble(m)
row,col,val = as_backend_type(M).data()
M_sp=scipy.sparse.csc_matrix((val,col,row))

from sksparse.cholmod import cholesky
factor=cholesky(M_sp,ordering_method=“natural”)
L=factor.L()

W=np.random.randn(M_sp.get_shape()[0],1)
y=L.dot(W)

#generation of an example
e=spsolve(M_sp+Kx_sp+Ky_sp+Kz_sp,alpha*y)
nu=np.ones((M_sp.get_shape()[0],1))*mu
e=e.reshape(-1,1)
e=e+nu

file_mean_kalman = File(‘results/realisation.pvd’)
func = Function(V)
func.vector().set_local(e)
func.rename(‘Kal_mean’,‘Kal_mean’)
file_mean_kalman << func
1 Like

This question appears to be a duplicate of Solution of my problem e does not match my mesh. Please don’t pose a question multiple times to try to get faster replies.

To avoid issues with ordering, you could use dolfin to create both your mass and stiffness matrices as then they will be in the same order.

Alternatively, you could get the stiffness matrix for each cell (see Element stiffness matrix), then loop through each cell.

If you still have difficulties, please reduce your code to a minimal failing example (see Read before posting: How do I get my question answered?) so that we can more easily run it and work out what the issue is. It’s also helpful if you wrap your code in three backticks (```) so that it is formatted as code.

3 Likes

Thank you for your response.

I have just modified the title of the problem so it can be more clear, i did that because i couldn’t modify the other one.

So, the thing is I cannot reduce the code this is it and it’s not too hard just calculate stiffness matrix in all directions and the mass matrix and then calculate e.

The thing is when i plot e on the mesh, i can see that my field is not ok, since i have not more correlation on direction x because i made beta_x>beta_y>beta_z

So my question is : i don’t Know if i am representing well e on my mesh? I don’t understand where is the trick and how can i correct it ?

import os
import numpy as np
from numpy import linalg
import scipy
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
import scipy.stats as stats
from scipy.optimize import minimize

alpha=1.22
beta_x=2
beta_y=0.01
beta_z=0.01
mu=-4.25
str_os = ‘dolfin-convert ./mesh/calcul_little_specimen.msh ./calcul_little_specimen.xml’
os.system(str_os)
mesh = Mesh(“calcul_little_specimen.xml”)
V = FunctionSpace(mesh, “Lagrange”, 1)
u = TrialFunction(V)
v = TestFunction(V)
parameters[‘linear_algebra_backend’] = ‘Eigen’

S=as_matrix([[beta_x**2,0,0],[0,0,0],[0,0,0]])
kx=inner(dot(S,grad(u)),grad(v))*dx
Kx=assemble(kx)
row,col,val = as_backend_type(Kx).data()
Kx_sp=scipy.sparse.csc_matrix((val,col,row))

S=as_matrix([[0,0,0],[0,beta_y**2,0],[0,0,0]])
ky=inner(dot(S,grad(u)),grad(v))*dx
Ky=assemble(ky)
row,col,val = as_backend_type(Ky).data()
Ky_sp=scipy.sparse.csc_matrix((val,col,row))

S=as_matrix([[0,0,0],[0,0,0],[0,0,beta_z**2]])
kz=inner(dot(S,grad(u)),grad(v))*dx
Kz=assemble(kz)
row,col,val = as_backend_type(Kz).data()
Kz_sp=scipy.sparse.csc_matrix((val,col,row))

m=u *v* dx
M = assemble(m)
row,col,val = as_backend_type(M).data()
M_sp=scipy.sparse.csc_matrix((val,col,row))

from sksparse.cholmod import cholesky
factor=cholesky(M_sp,ordering_method=“natural”)
L=factor.L()

W=np.random.randn(M_sp.get_shape()[0],1)
y=L.dot(W)

#generation of an example
e=spsolve(M_sp+Kx_sp+Ky_sp+Kz_sp,alpha*y)
nu=np.ones((M_sp.get_shape()[0],1))*mu
e=e.reshape(-1,1)
e=e+nu

file_mean_kalman = File(‘results/realisation.pvd’)
func = Function(V)
func.vector().set_local(e)
func.rename(‘Kal_mean’,‘Kal_mean’)
file_mean_kalman << func

First thing with reducing your code is to use a built in mesh.
As Matthew said, you can Get the mass matrix per element using a local solver.
This can then be put in a DG function of appropriate order, and visualized using XDMFFile.write_checkpoint .

Note that when you use PVD, the corresponding field is interpolated into a CG-1 field, Which is not suitable for mass matrices.

Hello,

@dokken and @mscroggs Thank you very much for your response.

Below is the code with the built in mesh and it works perfectly fine, i am using beta_x > beta_y=beta_z , so as you can see if you run the code that the field is correlated in direction x but on both direction y and z the gradient of the field is discontinuous. (this is what i want to see on the same mesh using gmsh but i don’t see it )

So I think i didn’t clear enough my problem, so i am stuck now because if I want to reproduce the same algorithm but just using a gmsh mesh ( I need to use it because i have relatively complex geometry to test on) as the code in my previous comment ( I mean really changing just the mesh) in this case i don’t see the field anymore, even when I use the same geometry(It means BoxMesh but reproduced with gmsh it means an msh file)

That’s why , i don’t really understand why the code is doing this , i mean i feel like when plotting e on my gmsh mesh the values do not correspond to the vertex in an order way or something like that) ? I cannot understand that in any way and i can not explain it. I hope i was clear enough and i hope i find a solution and an explication for my problem.

PS :
1.Does Boxmesh give me tetrahedron elements of the mesh It means 3D elements?
2.The gmsh i am using is a 3D mesh of the tetrahedron elements with four nodes (c3d8)
3. gmsh version is 3.0.6
4.pythhon 3.8.5

Thank you all.
Best regards.

import os
import numpy as np
from numpy import linalg
import scipy
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
import scipy.stats as stats
from scipy.optimize import minimize

alpha=1.22
beta_x=2
beta_y=0.01
beta_z=0.01
mu=-4.25

P0=Point(0.,0.,0.)
P1=Point(16.,0.6,0.6)
nx=100
ny=24
nz=24
mesh = BoxMesh(P0,P1,nx,ny,nz)

V = FunctionSpace(mesh, “Lagrange”, 1)
u = TrialFunction(V)
v = TestFunction(V)
parameters[‘linear_algebra_backend’] = ‘Eigen’

S=as_matrix([[beta_x**2,0,0],[0,0,0],[0,0,0]])
kx=inner(dot(S,grad(u)),grad(v))*dx
Kx=assemble(kx)
row,col,val = as_backend_type(Kx).data()
Kx_sp=scipy.sparse.csc_matrix((val,col,row))

S=as_matrix([[0,0,0],[0,beta_y**2,0],[0,0,0]])
ky=inner(dot(S,grad(u)),grad(v))*dx
Ky=assemble(ky)
row,col,val = as_backend_type(Ky).data()
Ky_sp=scipy.sparse.csc_matrix((val,col,row))

S=as_matrix([[0,0,0],[0,0,0],[0,0,beta_z**2]])
kz=inner(dot(S,grad(u)),grad(v))*dx
Kz=assemble(kz)
row,col,val = as_backend_type(Kz).data()
Kz_sp=scipy.sparse.csc_matrix((val,col,row))

m=u *v* dx
M = assemble(m)
row,col,val = as_backend_type(M).data()
M_sp=scipy.sparse.csc_matrix((val,col,row))

from sksparse.cholmod import cholesky
factor=cholesky(M_sp,ordering_method=“natural”)
L=factor.L()

W=np.random.randn(M_sp.get_shape()[0],1)
y=L.dot(W)

#generation of an example
e=spsolve(M_sp+Kx_sp+Ky_sp+Kz_sp,alpha*y)
nu=np.ones((M_sp.get_shape()[0],1))*mu
e=e.reshape(-1,1)
e=e+nu

file_mean_kalman = File(‘results/realisation.pvd’)
func = Function(V)
func.vector().set_local(e)
func.rename(‘Kal_mean’,‘Kal_mean’)
file_mean_kalman << func

here is a code a code where i reduce the domain, i will post the msh mesh also.
Why just by changing the type of mesh, it means changing boxMesh with a dolfin converted msh mesh i do not get the same result.

you can activate alternatevily a type of mesh and see the difference, is e maybe not applied on vertices of mesh correctly when i use a gmsh ?

from dolfin import *
import os
import numpy as np
from numpy import linalg
from math import exp
from math import sqrt
from math import pi
from math import log
import matplotlib.pyplot as plt
import scipy
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
import scipy.stats as stats
from scipy.optimize import minimize
import vtk
from vtk import *



# Ce code est fait pour la détermination de la matrice de masse ( en thermique) gloable
     #mais les matrices de rigidité suivant x , suivant y et suivant z

#nx = 1
#ny = 1
     
alpha=1.22
beta_x=1
beta_y=0.01 #0.05
beta_z=0.01 #0.05
mu=-4.25   
#mesh = UnitSquareMesh(nx, ny)

os.system('dolfin-convert mesh.msh mesh.xml')
mesh = Mesh("mesh.xml")
#mesh = BoxMesh(Point(0,0, 0),Point(4,0.6, 0.6), 60, 10, 10)


V = FunctionSpace(mesh, "Lagrange", 1)

u = TrialFunction(V)
v = TestFunction(V)
#Explicitly cast the matrix to the backend type
parameters['linear_algebra_backend'] = 'Eigen'
# compute the stiffness/mass matrix associated with the bilinear form a
#Stifness matrix along x
S=as_matrix([[beta_x**2,0,0],[0,0,0],[0,0,0]])
kx=inner(dot(S,grad(u)),grad(v))*dx
Kx=assemble(kx)
row,col,val = as_backend_type(Kx).data()
Kx_sp=scipy.sparse.csc_matrix((val,col,row)) 

#Stifness matrix along y
S=as_matrix([[0,0,0],[0,beta_y**2,0],[0,0,0]])
ky=inner(dot(S,grad(u)),grad(v))*dx
Ky=assemble(ky)
row,col,val = as_backend_type(Ky).data()
Ky_sp=scipy.sparse.csc_matrix((val,col,row)) 



#Stifness matrix along z
S=as_matrix([[0,0,0],[0,0,0],[0,0,beta_z**2]])
kz=inner(dot(S,grad(u)),grad(v))*dx
Kz=assemble(kz)
row,col,val = as_backend_type(Kz).data()
Kz_sp=scipy.sparse.csc_matrix((val,col,row)) 


#The mass matrix
m=u*v*dx
M = assemble(m)
row,col,val = as_backend_type(M).data()
M_sp=scipy.sparse.csc_matrix((val,col,row)) 

from sksparse.cholmod import cholesky
factor=cholesky(M_sp,ordering_method="natural")
L=factor.L()
# Wiener
W=np.random.randn(M_sp.get_shape()[0],1)
y=L.dot(W)

#generation of an example
e=spsolve(M_sp+Kx_sp+Ky_sp+Kz_sp,alpha*y)
nu=np.ones((M_sp.get_shape()[0],1))*mu
e=e.reshape(-1,1)
e=e+nu



file_mean_kalman = File('results/realisation_cube.pvd')
func = Function(V)
func.vector().set_local(e)
func.rename('Kal_mean','Kal_mean')
file_mean_kalman << func 
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
86
1 4 0.6 0
2 4 0.6 0.6
3 4 0 0
4 4 0 0.6
5 0 0.6 0
6 0 0.6 0.6
7 0 0 0
8 0 0 0.6
9 0 0.296058 0.6
10 0.323332 0.6 0.6
11 0.728028 0.6 0.6
12 1.14718 0.6 0.6
13 1.56633 0.6 0.6
14 1.98548 0.6 0.6
15 2.40463 0.6 0.6
16 2.82378 0.6 0.6
17 3.24293 0.6 0.6
18 3.66208 0.6 0.6
19 4 0.296058 0.6
20 0.323332 0 0.6
21 0.728028 0 0.6
22 1.14718 0 0.6
23 1.56633 0 0.6
24 1.98548 0 0.6
25 2.40463 0 0.6
26 2.82378 0 0.6
27 3.24293 0 0.6
28 3.66208 0 0.6
29 0 0.296058 0
30 0.323332 0.6 0
31 0.728028 0.6 0
32 1.14718 0.6 0
33 1.56633 0.6 0
34 1.98548 0.6 0
35 2.40463 0.6 0
36 2.82378 0.6 0
37 3.24293 0.6 0
38 3.66208 0.6 0
39 4 0.285033 0
40 0.323332 0 0
41 0.728028 0 0
42 1.14718 0 0
43 1.56633 0 0
44 1.98548 0 0
45 2.40463 0 0
46 2.82378 0 0
47 3.24293 0 0
48 3.66208 0 0
49 4 0.6 0.297824
50 0 0.6 0.297824
51 4 0 0.286798
52 0 0 0.297824
53 0.509599 0.299931 0.6
54 3.4698 0.299959 0.6
55 0.93594 0.299997 0.6
56 1.35644 0.3 0.6
57 1.77587 0.3 0.6
58 2.19516 0.3 0.6
59 2.61481 0.3 0.6
60 3.03653 0.299998 0.6
61 0.507507 0.298635 0
62 3.18864 0.3 0
63 0.935582 0.29994 0
64 1.35649 0.299997 0
65 1.77651 0.3 0
66 2.19859 0.3 0
67 2.6337 0.3 0
68 0.509598 0.6 0.299962
69 0.935937 0.6 0.299998
70 1.35643 0.6 0.3
71 1.7758 0.6 0.3
72 2.19481 0.6 0.3
73 2.61299 0.6 0.299998
74 3.0272 0.6 0.299959
75 3.4232 0.6 0.299326
76 3.76133 0.6 0.294362
77 0.509601 0 0.299962
78 3.46983 0 0.299862
79 0.93594 0 0.299998
80 1.35644 0 0.3
81 1.77587 0 0.3
82 2.19516 0 0.3
83 2.61481 0 0.3
84 3.03653 0 0.299994
85 4 0.294919 0.295164
86 1.14541 0.299998 0.385121
$EndNodes
$Elements
176
1 4 2 1 2 64 56 86 80
2 4 2 1 2 12 56 86 70
3 4 2 1 2 29 52 7 40
4 4 2 1 2 86 55 69 79
5 4 2 1 2 57 81 80 64
6 4 2 1 2 57 23 80 81
7 4 2 1 2 57 56 80 23
8 4 2 1 2 69 68 31 63
9 4 2 1 2 45 44 82 66
10 4 2 1 2 71 33 34 65
11 4 2 1 2 19 76 85 54
12 4 2 1 2 68 61 31 63
13 4 2 1 2 44 65 82 66
14 4 2 1 2 9 53 68 10
15 4 2 1 2 62 84 60 67
16 4 2 1 2 80 43 42 64
17 4 2 1 2 9 50 6 10
18 4 2 1 2 23 57 24 81
19 4 2 1 2 63 77 79 53
20 4 2 1 2 9 61 68 53
21 4 2 1 2 75 62 48 78
22 4 2 1 2 63 77 53 61
23 4 2 1 2 78 48 47 62
24 4 2 1 2 12 11 69 55
25 4 2 1 2 65 82 66 58
26 4 2 1 2 80 43 64 81
27 4 2 1 2 58 81 82 24
28 4 2 1 2 22 80 86 56
29 4 2 1 2 84 26 60 83
30 4 2 1 2 18 76 19 54
31 4 2 1 2 53 68 69 63
32 4 2 1 2 49 85 19 76
33 4 2 1 2 53 11 69 68
34 4 2 1 2 53 55 69 11
35 4 2 1 2 53 63 69 79
36 4 2 1 2 77 53 61 9
37 4 2 1 2 9 52 20 8
38 4 2 1 2 57 64 80 56
39 4 2 1 2 64 70 86 56
40 4 2 1 2 60 54 62 75
41 4 2 1 2 57 65 58 81
42 4 2 1 2 53 55 79 69
43 4 2 1 2 62 74 67 60
44 4 2 1 2 71 33 64 70
45 4 2 1 2 12 55 86 56
46 4 2 1 2 64 32 69 70
47 4 2 1 2 70 71 57 64
48 4 2 1 2 64 65 81 43
49 4 2 1 2 9 52 50 61
50 4 2 1 2 81 65 44 43
51 4 2 1 2 39 1 76 85
52 4 2 1 2 29 52 40 61
53 4 2 1 2 58 72 14 15
54 4 2 1 2 62 46 47 84
55 4 2 1 2 75 54 18 17
56 4 2 1 2 62 75 38 37
57 4 2 1 2 67 46 62 84
58 4 2 1 2 22 56 86 55
59 4 2 1 2 72 15 59 73
60 4 2 1 2 60 59 16 74
61 4 2 1 2 36 73 67 74
62 4 2 1 2 72 58 59 15
63 4 2 1 2 58 81 24 57
64 4 2 1 2 73 59 67 74
65 4 2 1 2 31 61 68 30
66 4 2 1 2 53 10 11 68
67 4 2 1 2 56 80 23 22
68 4 2 1 2 36 74 67 62
69 4 2 1 2 64 63 69 32
70 4 2 1 2 40 52 77 61
71 4 2 1 2 55 21 53 79
72 4 2 1 2 62 54 78 75
73 4 2 1 2 53 21 77 79
74 4 2 1 2 71 14 58 72
75 4 2 1 2 16 59 73 74
76 4 2 1 2 28 85 78 54
77 4 2 1 2 67 83 45 46
78 4 2 1 2 77 52 9 61
79 4 2 1 2 82 83 66 59
80 4 2 1 2 39 38 76 1
81 4 2 1 2 51 48 78 85
82 4 2 1 2 58 59 25 82
83 4 2 1 2 39 85 76 48
84 4 2 1 2 26 83 25 59
85 4 2 1 2 86 79 69 64
86 4 2 1 2 60 17 75 74
87 4 2 1 2 36 74 62 37
88 4 2 1 2 65 72 58 66
89 4 2 1 2 45 83 67 66
90 4 2 1 2 78 54 85 76
91 4 2 1 2 64 79 69 63
92 4 2 1 2 9 53 20 77
93 4 2 1 2 75 38 48 62
94 4 2 1 2 78 62 47 84
95 4 2 1 2 34 65 66 72
96 4 2 1 2 62 74 75 37
97 4 2 1 2 64 80 86 79
98 4 2 1 2 62 74 60 75
99 4 2 1 2 82 66 58 59
100 4 2 1 2 50 61 29 30
101 4 2 1 2 51 39 3 48
102 4 2 1 2 58 82 25 24
103 4 2 1 2 72 66 59 58
104 4 2 1 2 86 64 69 70
105 4 2 1 2 60 54 27 84
106 4 2 1 2 62 54 60 84
107 4 2 1 2 16 17 60 74
108 4 2 1 2 78 54 62 84
109 4 2 1 2 53 21 20 77
110 4 2 1 2 39 48 76 38
111 4 2 1 2 82 65 44 81
112 4 2 1 2 71 57 58 14
113 4 2 1 2 57 65 64 71
114 4 2 1 2 28 19 4 85
115 4 2 1 2 60 59 67 83
116 4 2 1 2 58 65 82 81
117 4 2 1 2 51 28 4 85
118 4 2 1 2 16 15 73 59
119 4 2 1 2 22 79 86 80
120 4 2 1 2 69 63 31 32
121 4 2 1 2 5 50 29 30
122 4 2 1 2 33 32 64 70
123 4 2 1 2 49 1 85 76
124 4 2 1 2 66 35 34 72
125 4 2 1 2 22 21 55 79
126 4 2 1 2 42 79 80 64
127 4 2 1 2 42 79 64 63
128 4 2 1 2 71 65 58 57
129 4 2 1 2 57 81 64 65
130 4 2 1 2 12 69 86 55
131 4 2 1 2 20 52 9 77
132 4 2 1 2 67 73 36 35
133 4 2 1 2 27 26 60 84
134 4 2 1 2 82 83 45 66
135 4 2 1 2 67 46 84 83
136 4 2 1 2 51 85 78 28
137 4 2 1 2 67 66 73 35
138 4 2 1 2 27 54 78 84
139 4 2 1 2 25 83 82 59
140 4 2 1 2 71 65 64 33
141 4 2 1 2 72 73 59 66
142 4 2 1 2 12 70 86 69
143 4 2 1 2 40 61 77 41
144 4 2 1 2 68 61 50 30
145 4 2 1 2 63 68 61 53
146 4 2 1 2 79 63 42 41
147 4 2 1 2 9 50 68 61
148 4 2 1 2 71 65 34 72
149 4 2 1 2 2 19 18 76
150 4 2 1 2 61 63 77 41
151 4 2 1 2 51 39 48 85
152 4 2 1 2 75 76 48 38
153 4 2 1 2 9 10 68 50
154 4 2 1 2 18 54 75 76
155 4 2 1 2 48 76 78 85
156 4 2 1 2 60 54 75 17
157 4 2 1 2 27 28 78 54
158 4 2 1 2 19 85 28 54
159 4 2 1 2 50 52 29 61
160 4 2 1 2 75 54 78 76
161 4 2 1 2 49 19 2 76
162 4 2 1 2 71 72 58 65
163 4 2 1 2 75 78 48 76
164 4 2 1 2 59 83 66 67
165 4 2 1 2 67 59 60 74
166 4 2 1 2 66 72 73 35
167 4 2 1 2 59 73 67 66
168 4 2 1 2 12 13 56 70
169 4 2 1 2 70 13 57 71
170 4 2 1 2 22 55 86 79
171 4 2 1 2 70 56 57 13
172 4 2 1 2 26 59 60 83
173 4 2 1 2 70 64 57 56
174 4 2 1 2 77 63 79 41
175 4 2 1 2 67 84 60 83
176 4 2 1 2 71 57 14 13
$EndElements
'''

My main problem with this example is the following:

  1. You are importing several used variables (such as vtk and maths imports), as well as having a wildcard import of vtk.
  2. As most people do not have sksparse on their system, or (CHOLMOD), it requires installation of several additional dependencies (and most people don’t know what these models do and assume of the input).

As you can see if you are solving any problem with these matrices, they do not change the ordering of the output:

from dolfin import *
import os     

os.system('dolfin-convert mesh.msh mesh.xml')
mesh = Mesh("mesh.xml")
#mesh = BoxMesh(Point(0,0, 0),Point(4,0.6, 0.6), 15, 3, 3)


V = FunctionSpace(mesh, "Lagrange", 1)


u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
M = inner(u,v)*dx
L = inner(x[0]*x[1],v)*dx
uh = Function(V)
solve(M==L, uh)

File("u_out.pvd") << uh