计算科学 - Dolfin convert：如何在（3D）单元格的顶点插入数据？ - 吾爱随笔录

我希望你们中的一个人能帮助我，因为我在这里被困了一个星期。我正在尝试使用 dolfin convert to XML 读取 gmsh 文件 (.msh)，然后使用 dolfin 下载它。

问题是当我组装刚度矩阵和质量矩阵并尝试使用这些矩阵生成示例（蒙特卡洛）时，我发现它们不是好的矩阵。

我认为 dolfin 在下载已经从 MSH 转换的 XML 文件时，他以升序方式更改了元素编号的顺序。我确信他会这样做，因为我已经在 fenics 的 BoxMesh 生成的盒子上尝试了我的代码，并且我的程序正常工作，如果我检查 mesh.cells() 和我的 XML 文件，我可以看到排序已更改，所以这就是为什么我没有我的好矩阵

老实说，我不知道如何解决这个问题？dolfin 是否会破坏任何网格的排序并建立另一个排序（以升序方式）？当您下载您的类型的网格时，你们如何获得良好的刚度矩阵和质量矩阵？

我读过另一个人发现的同样问题，但仍然不明白如何解决这个问题。

非常感谢您提前

这是我的代码

from dolfin import *

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=uvdx

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