一千萬個為什麽

搜索

FEniCS:如何在(3D)單元的頂點插入數據?

我試圖在單元格的所有頂點獲得插值函數$ 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

然後我在$ z = 0 $ plane上讀取了(xp,yp)點並評估這些點的功能

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)

現在當我繪制$ g_i $時,它看起來不像它應該是什麽,即在圓圈裏面的常數(3.91)$ R = 50 $,0在外面。 任何幫助,將不勝感激。

最佳答案

定義自定義 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()

轉載註明原文: FEniCS:如何在(3D)單元的頂點插入數據?