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