基于node和cell进行加密，离散得到错误矩阵

[zhf-0]

常用的网格生成方法是

n = 10
box = [0, 1, 0, 1, 0, 1]
mesh = TetrahedronMesh.from_box(box, nx=n, ny=n, nz=n)

但是我们需要对网格进行加密（个人实现的程序，不是fealpy中提供的加密程序），再基于加密后的node和cell生成新的网格对象

new_node, new_cell = refine(node, cell)
new_mesh = TetrahedronMesh(new_node, new_cell)

但是基于new_mesh离散得到的矩阵，与理论上的矩阵相差很大。

我检查了网格加密代码，并进行了初步测试和可试化分析，没有找到明显错误。想请教一下，在生成的new_node和new_cell中，是否需要保证 new_cell每行中，节点编号的某种顺序？是不是因为new_cell中节点编号的顺序不对，导致生成的矩阵有问题？

[weihuayi]

@zhf-0 单元节点的顺序要保证满足右手法则。另外，你的加密方式是什么？

[zhf-0]

我们参考的加密算法是 Jurgen Bey 的《Simplicial grid refinement: on Freudenthal’s algorithm and the optimal number of congruence classes》，更容易构造粗细两种网格之间的插值矩阵，但是在实现的过程中没有考虑加密后单元节点编号的顺序。您说的右手法则具体含义是什么？有参考资料吗？

[weihuayi]

给定一个四面体单元 $T = (x_0, x_1, x_2, x_3)$, 满足右手法则，向量 $x_0x_1$, $x_0x_2$ 和 $x_0x_3$ 的混合积是一个正数

[zhf-0]

好的，谢谢魏老师。

[zhf-0]

魏老师，我在运行fealpy提供的网格加密算法过程中遇到了问题：

• 运行fealpy/example/fem/poisson_tetrahedronmesh_lfem_dirichletbc_3d.py 程序，网格加密后的误差收敛阶是正常的
• 运行fealpy/example/fem/MaxwellWithFNSpace_3d.py程序，误差收敛阶也是正常的，只不过网格不是加密得到的，而是每次重新生成网格。输出的误差与收敛阶如下所示

The 0-th computation:
19
0.22128886098480263
The 1-th computation:
98
0.11498432454669744
The 2-th computation:
604
0.058020118786895426
The 3-th computation:
4184
0.029074342209268616
The 4-th computation:
31024
0.014545042817131228
ratio
1.9245132922003885
1.981800915800057
1.9955780381644947
1.998917608893162

现在我在 fealpy/example/fem/MaxwellWithFNSpace_3d.py的基础上，通过加密网格的方式得到新的网格，部分程序如下所示

pde = PDE()
maxit = 5
errorType = ['$|| E - E_h||_{\Omega,0}$']
errorMatrix = np.zeros(maxit)
NDof = np.zeros(maxit, dtype=np.int_)
mesh = pde.init_mesh(2**0)

for i in range(maxit):
    print("The {}-th computation:".format(i))

    space = FirstNedelecFiniteElementSpace3d(mesh)

    gdof = space.dof.number_of_global_dofs()
    NDof[i] = gdof
    print(gdof)

    bc = DirichletBC(space, pde.dirichlet) 

    M = space.mass_matrix()
    A = space.curl_matrix()
    b = space.source_vector(pde.source)
    B = A-M 

    Eh = space.function()
    #B, b = bc.apply(B, b, Eh)
    isDDof = space.set_dirichlet_bc(pde.dirichlet, Eh)
    b[isDDof] = Eh[isDDof]

    bdIdx = np.zeros(B.shape[0], dtype=np.int_)
    bdIdx[isDDof] = 1
    Tbd = spdiags(bdIdx, 0, B.shape[0], B.shape[0])
    T = spdiags(1-bdIdx, 0, B.shape[0], B.shape[0])
    B = T@B + Tbd

    Eh[:] = spsolve(B, b)
    # 计算误差
    errorMatrix[i] = space.integralalg.error(pde.solution, Eh)
    print(errorMatrix[i])

    if i < maxit-1:
        mesh.uniform_refine()

print('ratio')
for i in range(1,maxit):
    print(errorMatrix[i-1]/errorMatrix[i])

但是基于以上方式离散求解，误差收敛阶出错了，输出内容如下所示

The 0-th computation:
19
0.22128886098480263
The 1-th computation:
98
0.6295276457587462
The 2-th computation:
604
0.6560324946586827
The 3-th computation:
4184
0.6821608320315334
The 4-th computation:
31024
0.664079404764903
ratio
0.3515157157523899
0.9595982682020555
0.9616976874866323
1.0272278091097127

我更新了fealpy到最新版本，为什么在poisson问题中，加密网格可以得到正确的结果，但是在Maxwell问题中，加密相同网格结果出错了？

[weihuayi]

@zhf-0 我已经安排同学来检查测试了。

[zhf-0]

好的，谢谢魏老师

[weihuayi]

@zhf-0 程序已经更新了，你测试一下吧

[zhf-0]

测试的收敛阶是正常的，谢谢魏老师

[zhf-0]

魏老师，您好。我想问一下，fealpy程序中四面体网格加密算法有参考资料吗？我们想基于粗细两套网格构造插值矩阵，但是两套网格的自由度编号不好对应，比如粗网格中1个单元对应细网格内哪8个单元？粗网格中的1条长边对应细网格内哪2条短边？

[weihuayi]

@zhf-0 你是说 bisect 自适应吗？

[zhf-0]

我想问的是fealpy/mesh/tetrahedron_mesh.py文件中，class TetrahedronMesh中的uniform_refine()函数用的加密算法。这个uniform_refine() 函数用的是bisect自适应加密吗？

[weihuayi]

@zhf-0 我在 TetrahedronMesh.uniform_refine 增加了一个参数 retrunim=False，设置成 True 以后，可以返回一致加密过程中的线性元和分片常数元的插值矩阵列表，你可以测试一下

[zhf-0]

好的， 谢谢！