FEniCS：如何在（3D）单元格的顶点插入数据？

计算科学 芬尼克斯
2021-12-20 15:06:28

我正在尝试获得插值函数f（在 3D 中）在单元格的所有顶点。我提取单元格的所有顶点，然后将值分配给每个顶点：如果它在半径范围内R，然后我分配值，例如 3.91。如果它在球体之外，那么我将值赋值为 0。我让它在没有错误消息的情况下运行，但是当我在不是顶点的点处计算函数值时，它不会给我 3.91 或 0。我在做什么这里有什么问题吗？

这是我的代码的一部分：

提取所有单元格的顶点，然后导出这些点并使用另一个软件为每个点分配值（例如 3.91 表示球内的点，0 表示球外的点）

coor = mesh.coordinates()
numpy.savetxt('meshforE.txt',coor)

我得到所有顶点的值，然后将其插入到函数 f

qvalues2 = numpy.loadtxt('qdata.txt')
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.vector()[:] = qvalues2

然后我读(xp, yp)了点z=0平面并评估这些点的功能

with open('xpdata.txt') as g:
    xp = g.readlines()
print "xp[1]=", xp[1]
with open('ypdata.txt') as h:
    yp = h.readlines()
print "yp[2]=", yp[2]

for i in range(len(xp)):
    g_in[i] = f(xp[i],yp[i],0.0)

现在当我绘制gi，它看起来不像它应该是的，即圆圈内的常数（3.91）R=50, 外为 0。任何帮助，将不胜感激。

4个回答

定义一个自定义Expression然后将其插入到您的函数空间中应该可以解决问题。退房

from dolfin import *                                                        

mesh = UnitSquareMesh(20, 20)                                               

class CharCircle(Expression):                                               
    def eval(self, value, x):                                               
        xm0 = x[0] - 0.5                                                    
        xm1 = x[1] - 0.5                                                    
        if xm0*xm0 + xm1*xm1 < 0.4**2:                                      
            value[0] = 1.0                                                  
        else:                                                               
            value[0] = 0.0                                                  

V = FunctionSpace(mesh, 'CG', 1)                                            
u = Function(V)                                                             
u.interpolate(CharCircle())                                                 

plot(u)                                                                     
interactive()

我假设 qvalues2 在顶点有一些计算值。您不能直接将这些分配给您的自由度向量，因为自由度不遵循顶点编号。

但是，您可以尝试：

vertex_to_dof = V.dofmap().vertex_to_dof_map(mesh)
f.vector()[:] = qvalues2[vertex_to_dof]

您正在使用分段线性函数 ( V = FunctionSpace(mesh, "CG", 1))，因此如果您在不是顶点的点上评估函数，您将获得线性插值。你没有说你正在使用的网格，但如果它不是曲线的，我不会假设那个圆上的点与元素边界上的点重合（插值是恒定的）。

如果您在所有顶点上都有函数的信息，则具有连续 galerkin 度 1 的函数空间的表示仅在这些顶点处是精确的。在任何其他点评估它会给出点之间的线性插值。这就是为什么您在值之间观察到 3.91 的不同值。

如果您希望球内的任何评估完全等于分配的常数，则需要每个单元格上的函数信息。当所有顶点都在球体内时，也许通过计算一个细胞在球体内。然后你可以用不连续的 galerking 度 0 来表示函数。任何点的任何评估都比预期值。

其它你可能感兴趣的问题
上一篇U D U⊤UDU⊤LAPACK/Eigen 中的分解例程？ 下一篇评估多项式符号极限的最快方法