非常规网格上的离散拉普拉斯算子 (python)

Discrete Laplacian on a non-regular mesh (python)

我已经为非常规网格(使用 scipy.spatial.Delaunay 函数创建)编写了 laplacien 函数。 我没有错误但结果不正确:特征向量正确但特征值太高(绝对值)。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.spatial

def rect_drum(L,H,U):
    vals = []
    val = 0
    k = 1
    l = 1
    while val >= -U:
        while val >= -U:
            val = -np.pi**2*((k/L)**2+(l/H)**2)
            if val >= -U:
                vals.append(val)
            l += 1
        l = 1
        k += 1
        val = -np.pi**2*((k/L)**2+(l/H)**2)
    return np.array(vals)


def count_vp(tab,U):
    #count the n eigenvalues ​​greater than equal to -U in the array tab
    return tab[tab>=-U]

def in_curve(f,fargs,shape,a):
    points = [] # the points inside the curve
    for j in range(shape[0]):
        for i in range(shape[1]):
            if f(i*a,j*a,*fargs) < 0:
                points.append([i*a,j*a])

    return np.array(points)

def triang(points,a,f,fargs,bord):
    tri_points = points.copy()
    tri_points[:,1] *= np.sqrt(3)
    tri_points2 = np.vstack((points,bord))
    tri_points2[:,1] *= np.sqrt(3)
    tri_points2[:,0] += a/2
    tri_points2[:,1] += np.sqrt(3)/2*a
    fin = np.vstack((tri_points,tri_points2))
    i = 0
    eps = 0.01
    while i < len(fin):
        if f(fin[i,0]+eps,fin[i,1]+eps,*fargs) > 0:
            fin = np.delete(fin,i,0)
            i -= 1
        i += 1
    return np.vstack((fin,bord)),len(fin),len(bord)


def tri_ang(points,ind,p0):
    # sort the points in trigonometric order
    vec=np.arctan2((points-p0)[:,1],(points-p0)[:,0])
    values = []
    dtype = [('val',float),('n',int)]
    
    for i in range(len(vec)):
        values.append((vec[i],i))
    values = np.sort(np.array(values,dtype),order='val') 
    new_points = []
    new_ind = []
    for tup in values:
        new_points.append(points[tup[1]])
        new_ind.append(ind[tup[1]])
    return np.array(new_points),np.array(new_ind)




def M(points,tri,Nint):
    
    indptr,ind = tri.vertex_neighbor_vertices
    W = np.zeros((Nint,Nint)) # cotangents matrix
    A = np.zeros((Nint,1)) # surfaces vertex array for each point i (A[i])
    
    for i in range(Nint):
        tot = 0
        nhb_ind = ind[indptr[i]:indptr[i+1]] # indices of the points close to the point of index k
        nhb = points[nhb_ind] # their coordinates
        nhb,nhb_ind = tri_ang(nhb,nhb_ind,points[i]) #the coordinates (nhb) and  (nhb_ind) of each neighbor of i
        
        for j in range(len(nhb_ind)):
            vec = nhb[j]-points[i] # a vector connecting the point to his neighbor of index 0
            vec_av = nhb[j-1]-points[i] # another vector but with the Vosin from before
            if j+1 >= len(nhb_ind):
                k = 0
            else:
                k = j+1
            vec_ap = nhb[k]-points[i] # another vector but with the next neighbor
            
            # another vector but with the next neighbor
            A[i] += 0.5/3*np.linalg.norm(np.cross(vec,vec_av))
            
            if nhb_ind[j] < Nint:
                # we use the vector and scalar product to calculate the cotangents: A.B/||AxB||
                cotan_alpha = np.dot(vec_av,vec_av-vec)/np.linalg.norm(np.cross(vec_av,vec_av-vec))
                cotan_beta = np.dot(vec_ap,vec_ap-vec)/np.linalg.norm(np.cross(vec_ap,vec_ap-vec))
                # Wij value :
                W[i,nhb_ind[j]] = -0.5*(cotan_alpha+cotan_beta)
                
                tot += cotan_alpha+cotan_beta
    
        W[i,i] = -0.5*tot # diagonal values
    
    return (1/A)*W


def rect(x,y,L,H,x0=0,y0=0):
    if 0<x-x0<L and 0<y-y0<H:
        return -1
    else:
        return 1

def rect_rim(L,H,a,x0=0,y0=0):
    tab1 = np.arange(x0,L+x0,a)[:,np.newaxis]
    h = np.hstack((tab1,H*np.ones((len(tab1),1))+y0))
    b = np.hstack((tab1,np.zeros((len(tab1),1))+y0))
    tab2 = np.arange(y0+a,H+y0,a)[:,np.newaxis]
    g = np.hstack((np.zeros((len(tab2),1))+x0,tab2))
    d = np.hstack((L*np.ones((len(tab2),1))+x0,tab2))
    hp = np.array([[L+x0,H+y0]])
    bp = np.array([[L+x0,0]])
    return np.vstack((h,b,g,d,hp,bp))
    



# sample with a square 1*1

L = 1
H = 1

dl = 0.05
sol = in_curve(rect,[L,H],(100,100),dl)
sol_tri,Nint,Nbord = triang(sol,dl,rect,[L,H],rect_rim(L,H,dl))

# plt.plot(sol_tri[:,0],sol_tri[:,1],linestyle="",marker="+",label="tri")
# plt.plot(sol[:,0],sol[:,1],linestyle="",marker="x")
# plt.legend()
# plt.show()

# triangulation
tri = scipy.spatial.Delaunay(sol_tri)
# plt.triplot(sol_tri[:,0],sol_tri[:,1],tri.simplices)
# plt.show()

M = M(sol_tri,tri,Nint)

valp,vecp = np.linalg.eig(M) # eigenvalues and eigenvectors
vecp = np.real(vecp)

# comparison with the exact solution:
T = 1000
U = np.arange(0,T,1)
NUsim = np.array([len(count_vp(valp,u)) for u in U])
NU = np.array([len(rect_drum(L,H,u)) for u in U])
plt.plot(U,NUsim,label='simulation')
plt.plot(U,NU,label='exacts')
plt.legend()
plt.show()


# 3D plot of an eigenvector
vecp_tot = np.vstack((vecp,np.zeros((Nbord,Nint))))
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(sol_tri[:,0],sol_tri[:,1],vecp_tot[:,0],triangles=tri.simplices)
plt.show()

拉普拉斯算子是名为“M”的函数。

“in_curve 函数”return 由 f(x,y,*fargs) < 0(样本中的正方形)定义的曲线内的点。

“三角”函数 return 点添加点(三角形网格)。该函数对曲线的边缘使用另一个函数(为了最精确),在示例中它是“rect_rim”函数。

我使用了 https://en.wikipedia.org/wiki/Discrete_Laplace_operator 中给出的公式(“网格拉普拉斯算子”)。

我已经解决了我的问题:这是一个标志和一个边缘问题。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.spatial

def rect_drum(L,H,U):
    vals = []
    val = 0
    k = 1
    l = 1
    while val >= -U:
        while val >= -U:
            val = -np.pi**2*((k/L)**2+(l/H)**2)
            if val >= -U:
                vals.append(val)
            l += 1
        l = 1
        k += 1
        val = -np.pi**2*((k/L)**2+(l/H)**2)
    return np.array(vals)


def count_vp(tab,U):
    #count the n eigenvalues ​​greater than equal to -U in the array tab
    return tab[tab>=-U]

def in_curve(f,fargs,shape,a):
    points = [] # the points inside the curve
    for j in range(shape[0]):
        for i in range(shape[1]):
            if f(i*a,j*a,*fargs) < 0:
                points.append([i*a,j*a])

    return np.array(points)

def triang(points,a,f,fargs,bord):
    tri_points = points.copy()
    tri_points[:,1] *= np.sqrt(3)
    tri_points2 = np.vstack((points,bord))
    tri_points2[:,1] *= np.sqrt(3)
    tri_points2[:,0] += a/2
    tri_points2[:,1] += np.sqrt(3)/2*a
    fin = np.vstack((tri_points,tri_points2))
    i = 0
    eps = 0.01
    while i < len(fin):
        if f(fin[i,0]+eps,fin[i,1]+eps,*fargs) > 0:
            fin = np.delete(fin,i,0)
            i -= 1
        i += 1
    return np.vstack((fin,bord)),len(fin),len(bord)


def tri_ang(points,ind,p0):
    # sort the points in trigonometric order
    vec=np.arctan2((points-p0)[:,1],(points-p0)[:,0])
    values = []
    dtype = [('val',float),('n',int)]
    
    for i in range(len(vec)):
        values.append((vec[i],i))
    values = np.sort(np.array(values,dtype),order='val') 
    new_points = []
    new_ind = []
    for tup in values:
        new_points.append(points[tup[1]])
        new_ind.append(ind[tup[1]])
    return np.array(new_points),np.array(new_ind)


def Laplacian(points,tri,Nint):
    
    indptr,ind = tri.vertex_neighbor_vertices
    W = np.zeros((Nint,Nint)) # cotangents matrix
    A = np.zeros((Nint,1)) # surfacesvertex aray of point i (A[i])
    
    for i in range(Nint):
        tot = 0
        nhb_ind = ind[indptr[i]:indptr[i+1]] # indices of the points close to the point of index k
        nhb = points[nhb_ind] # their coordinates
        nhb,nhb_ind = tri_ang(nhb,nhb_ind,points[i]) #the coordinates (nhb) and  (nhb_ind) of each neighbor of i
        
        for j in range(len(nhb_ind)):
            vec = nhb[j]-points[i] # a vector connecting the point to his neighbor of index 0
            vec_av = nhb[j-1]-points[i] # another vector but with the Vosin from before
            if j+1 >= len(nhb_ind):
                k = 0
            else:
                k = j+1
            vec_ap = nhb[k]-points[i] # another vector but with the next neighbor
            
            # we use the cross product to calculate the areas of the triangles: ||AxB||/2:
            A[i] += 0.5/3*np.linalg.norm(np.cross(vec,vec_av))
            
            # we use the cross product and scalar product to calculate the cotangents: A.B/||AxB||
            cotan_alpha = np.dot(vec_av,vec_av-vec)/np.linalg.norm(np.cross(vec_av,vec_av-vec))
            cotan_beta = np.dot(vec_ap,vec_ap-vec)/np.linalg.norm(np.cross(vec_ap,vec_ap-vec))
            
            tot += cotan_alpha+cotan_beta
            
            if nhb_ind[j] < Nint:
                W[i,nhb_ind[j]] = 0.5*(cotan_alpha+cotan_beta)
        
        W[i,i] = -0.5*tot # diagonal values
    
    return (1/A)*W


def rect(x,y,L,H,x0=0,y0=0):
    if 0<x-x0<L and 0<y-y0<H:
        return -1
    else:
        return 1

def rect_rim(L,H,a,x0=0,y0=0):
    tab1 = np.arange(x0,L+x0,a)[:,np.newaxis]
    h = np.hstack((tab1,H*np.ones((len(tab1),1))+y0))
    b = np.hstack((tab1,np.zeros((len(tab1),1))+y0))
    tab2 = np.arange(y0+a,H+y0,a)[:,np.newaxis]
    g = np.hstack((np.zeros((len(tab2),1))+x0,tab2))
    d = np.hstack((L*np.ones((len(tab2),1))+x0,tab2))
    hp = np.array([[L+x0,H+y0]])
    bp = np.array([[L+x0,0]])
    return np.vstack((h,b,g,d,hp,bp))
    



# sample with a square 1*1

L = 1
H = 1

dl = 0.04
sol = in_curve(rect,[L,H],(100,100),dl)
sol_tri,Nint,Nbord = triang(sol,dl,rect,[L,H],rect_rim(L,H,dl))

# triangulation
tri = scipy.spatial.Delaunay(sol_tri)

M = Laplacian(sol_tri,tri,Nint)

valp,vecp = np.linalg.eig(M) # eigenvalues and eigenvectors
vecp = np.real(vecp)

# comparison with the exact solution:
T = 1000
U = np.arange(0,T,1)
NUsim = np.array([len(count_vp(valp,u)) for u in U])
NU = np.array([len(rect_drum(L,H,u)) for u in U])
plt.plot(U,NUsim,label='simulation')
plt.plot(U,NU,label='exacts')
plt.legend()
plt.show()


# 3D plot of an eigenvector
mode = 0 # change this for an another mode
vecp_tot = np.vstack((vecp,np.zeros((Nbord,Nint))))
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(sol_tri[:,0],sol_tri[:,1],vecp_tot[:,mode],triangles=tri.simplices)
plt.show()

备注:

1- 最高特征值是假的:这是离散化的结果。

2- 如果 dl 太小,我们有错误的特征向量和特征值(在 valp 的顶部和 vecp 的第一个向量),这可能是由于网格划分的质量。