不规则区域上有限元方法的进一步探讨和疑惑

Galerkin码农选手

已于 2023-04-18 09:37:42 修改

阅读量476

点赞数 3

分类专栏：传统数值方法文章标签：线性代数

于 2021-01-27 15:15:17 首次发布

本文链接：https://blog.csdn.net/forrestguang/article/details/113248632

版权

传统数值方法专栏收录该内容

35 篇文章 17 订阅 ¥199.90 ¥99.00

订阅专栏

L型区域上的有限元方法

之前的传统数值方法基本上都是针对规则的矩形区域，按道理有限元方法的应用就是为不规则区域量身打造，为此本人尝试使用有限元方法来求解L型区域上的泊松方程，使用的同样是Galerkin原理。

混合边界

首先我们先尝试求解混合边界上的泊松方程，参考以下图像：
在这里插入图片描述
其中蓝色代表Dirichlet边界，橙色代表Neumann边界，有限元的基础知识本人已经在上一篇博客介绍有限元
因此我就不花篇幅介绍基函数以及试探检验函数空间的选取了，本人使用的仍然是矩形有限元基函数。下面本人重点介绍这样的L型区域在求解过程中的困难以及一些本人编程中的定义：
单元节点序号
这里面最主要的是我们需要弄清楚每个单元上面4个节点的整体序号。首先引入定义： $bounds=[[x_a,x_b,x_c],[y_a,y_b,y_c]]$ 表示L型区域的边界点，比如上面这个图中 $x_a=0,x_b=1,x_c=2;y_a=-2,y_b=0,y_c=1$ ，定义 $hx\_{}tr=[0.2,0.1]$ 表示在有限元FENET剖分中x轴剖分长度为0.2，y轴剖分长度为0.1，由此引出 $n x = [[5, 10], [10, 5]]$ 分别表示对应地方维度的单元数目。那么此时我们根据这个确定每个网格点尤其是边界上网格点的整体序号，在这里插入图片描述
比如说上面这个图第一个被我标注的点的整体序号 $2\times(nx[1,0]+nx[1,1]+1)-1$ ，第二个被我标注的点是 $nx[0,0]\times(nx[1,0]+nx[1,1]+1)+nx[1,0]$ ,例子就举这两个，其实整个计算过程和规则区域没有差别，只是需要特别注意网格点的整体序号，这个我给出代码大家慢慢体会。

import matplotlib.pyplot as plt
import numpy as np
import time


def UU(X, order,prob):#X表示(x,t)
    if prob==1:
        temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5
        if order[0]==0 and order[1]==0:
            return np.log(temp)
        if order[0]==1 and order[1]==0:#对x求偏导
            return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
        if order[0]==0 and order[1]==1:#对t求偏导
            return temp**(-1) * (20*(X[:,0]+X[:,1]) - 2*(X[:,0]-X[:,1]))
        if order[0]==2 and order[1]==0:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) ** 2 \
                   + temp**(-1) * (22)
        if order[0]==1 and order[1]==1:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
                   * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
                   + temp**(-1) * (18)
        if order[0]==0 and order[1]==2:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) ** 2 \
                   + temp**(-1) * (22)
    if prob==2:
        if order[0]==0 and order[1]==0:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   0.5*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==0:
            return (3*X[:,0]*X[:,0]-1) * \
                   0.5*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==1:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   (np.exp(2*X[:,1])-np.exp(-2*X[:,1]))
        if order[0]==2 and order[1]==0:
            return (6*X[:,0]) * \
                   0.5*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==1:
            return (3*X[:,0]*X[:,0]-1) * \
                   (np.exp(2*X[:,1])-np.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==2:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   2*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
    if prob==3:
        temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
        temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
        if order[0]==0 and order[1]==0:
            return temp1 * temp2**(-1)
        if order[0]==1 and order[1]==0:
            return (2*X[:,0]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,0])
        if order[0]==0 and order[1]==1:
            return (-2*X[:,1]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,1])
        if order[0]==2 and order[1]==0:
            return (2) * temp2**(-1) + \
                   2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
        if order[0]==1 and order[1]==1:
            return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
        if order[0]==0 and order[1]==2:
            return (-2) * temp2**(-1) + \
                   2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
    if prob == 4:
        if order == [0,0]:
            return X[:,0]*X[:,1]
        if order == [1,0]:
            return X[:,1]
        if order == [0,1]:
            return X[:,0]
        if order == [0,2]:
            return 0*X[:,0]
        if order == [2,0]:
            return 0*X[:,1]
def FF(prob,X):
    return -UU(X,[0,2],prob) - UU(X,[2,0],prob)
def Beta(X):
    #return 0*X[:,0]
    return np.exp(X[:,0])*X[:,1]**2
def NEU(X,prob,n):#定义g(x),beta > 0
    return UU(X,[1,0],prob)*n[:,0] + UU(X,[0,1],prob)*n[:,1] + Beta(X)*UU(X,[0,0],prob)
np.random.seed(1234)


class TESET():
    def __init__(self, bounds, hx):
        self.bounds = bounds
        
        nx = [[int((bounds[0,1] - bounds[0,0])/hx[0]),
               int((bounds[0,2] - bounds[0,1])/hx[0])],
              [int((bounds[1,1] - bounds[1,0])/hx[1]),
               int((bounds[1,2] - bounds[1,1])/hx[1])]]
        self.nx = np.array(nx)
        
        self.size = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[0,1]*(self.nx[1,0] + 1)
        self.X = np.zeros([self.size,2])
        m = 0
        for i in range(self.nx[0,0] + 1):
            for j in range(self.nx[1,0] + self.nx[1,1] + 1):
                self.X[m,0] = bounds[0,0] + i*hx[0]
                self.X[m,1] = bounds[1,0] + j*hx[1]
                m = m + 1
        for i in range(self.nx[0,1]):
            for j in range(self.nx[1,0] + 1):
                self.X[m,0] = bounds[0,1] + (i + 1)*hx[0]
                self.X[m,1] = bounds[1,0] + j*hx[1]
                m = m + 1
        self.NS = self.nx[0,0] + self.nx[0,1] + 1 + self.nx[1,0] + self.nx[1,1] + self.nx[0,1] + 1 + self.nx[1,1]
        self.NX = np.zeros([self.NS,2])
        m = 0
        for i in range(self.nx[0,0] + self.nx[0,1] + 1):
            self.NX[m,0] = bounds[0,0] + i*hx[0]
            self.NX[m,1] = bounds[1,0]
            m = m + 1
        for j in range(self.nx[1,0] + self.nx[1,1]):
            self.NX[m,0] = bounds[0,0]
            self.NX[m,1] = bounds[1,0] + (j + 1)*hx[1]
            m = m + 1
        for i in range(self.nx[0,1] + 1):
            self.NX[m,0] = bounds[0,1] + i*hx[0]
            self.NX[m,1] = bounds[1,1]
            m = m + 1
        for j in range(self.nx[1,1]):
            self.NX[m,0] = bounds[0,1]
            self.NX[m,1] = bounds[1,1] + (j + 1)*hx[1]
            m = m + 1
        self.DS = self.nx[0,1] - 1 + self.nx[1,0] - 1
        self.DX = np.zeros([self.DS,2])
        m = 0
        for i in range(self.nx[0,1] - 1):
            self.DX[m,0] = bounds[0,0] + (i + 1)*hx[0]
            self.DX[m,1] = bounds[1,2]
            m = m + 1
        for j in range(self.nx[1,0] - 1):
            self.DX[m,0] = bounds[0,2]
            self.DX[m,1] = bounds[1,0] + (j + 1)*hx[1]
            m = m + 1
class FENET():
    def __init__(self,bounds,hx):
        self.bounds = bounds
        self.dim = 2
        nx = [[int((bounds[0,1] - bounds[0,0])/hx[0]),
               int((bounds[0,2] - bounds[0,1])/hx[0])],
              [int((bounds[1,1] - bounds[1,0])/hx[1]),
               int((bounds[1,2] - bounds[1,1])/hx[1])]]
        self.nx = np.array(nx)
        self.hx = hx
        self.gp_num = 2
        self.gp_pos = [(1 - np.sqrt(3)/3)/2,(1 + np.sqrt(3)/3)/2]
        self.Node()
        self.Unit()
    def phi(self,X,order):#[-1,1]*[-1,1]，在原点取值为1，其他网格点取值为0的基函数
        ind00 = (X[:,0] >= -1);ind01 = (X[:,0] >= 0);ind02 = (X[:,0] >= 1)
        ind10 = (X[:,1] >= -1);ind11 = (X[:,1] >= 0);ind12 = (X[:,1] >= 1)
        if order == [0,0]:
            return (ind00*~ind01*ind10*~ind11).astype('float32')*(1 + X[:,0])*(1 + X[:,1]) + \
                    (ind01*~ind02*ind10*~ind11).astype('float32')*(1 - X[:,0])*(1 + X[:,1]) + \
                    (ind00*~ind01*ind11*~ind12).astype('float32')*(1 + X[:,0])*(1 - X[:,1]) + \
                    (ind01*~ind02*ind11*~ind12).astype('float32')*(1 - X[:,0])*(1 - X[:,1])
        if order == [1,0]:
            return (ind00*~ind01*ind10*~ind11).astype('float32')*(1 + X[:,1]) + \
                    (ind01*~ind02*ind10*~ind11).astype('float32')*(-(1 + X[:,1])) + \
                    (ind00*~ind01*ind11*~ind12).astype('float32')*(1 - X[:,1]) + \
                    (ind01*~ind02*ind11*~ind12).astype('float32')*(-(1 - X[:,1]))
        if order == [0,1]:
            return (ind00*~ind01*ind10*~ind11).astype('float32')*(1 + X[:,0]) + \
                    (ind01*~ind02*ind10*~ind11).astype('float32')*(1 - X[:,0]) + \
                    (ind00*~ind01*ind11*~ind12).astype('float32')*(-(1 + X[:,0])) + \
                    (ind01*~ind02*ind11*~ind12).astype('float32')*(-(1 - X[:,0]))
    def basic(self,X,order,i):#根据网格点的存储顺序，遍历所有网格点，取基函数
        temp = (X - self.Nodes[i,:])/np.array([self.hx[0],self.hx[1]])
        if order == [0,0]:
            return self.phi(temp,order)
        if order == [1,0]:
            return self.phi(temp,order)/self.hx[0]
        if order == [0,1]:
            return self.phi(temp,order)/self.hx[1]
    def Int_basic_basic(self,i,j,u_ind):#表示第i,j个基函数的梯度的乘积，以及在第u_ind个区域的积分
        X = self.Units_Int_Points[u_ind,:,:]#[4,2]
        basic0 = np.zeros_like(X)
        basic0[:,0] = self.basic(X,[1,0],i)
        basic0[:,1] = self.basic(X,[0,1],i)
        
        basic1 = np.zeros_like(X)
        basic1[:,0] = self.basic(X,[1,0],j)
        basic1[:,1] = self.basic(X,[0,1],j)
        return ((basic0*basic1).sum(1)).mean()*self.hx[0]*self.hx[1]
    def Int_F_basic(self,i,u_ind,prob):#第i个基函数与右端项乘积，在第u_ind个单元积分
        X = self.Units_Int_Points[u_ind,:,:]
        return (FF(prob,X)*self.basic(X,[0,0],i)).mean()*self.hx[0]*self.hx[1]
    def Bou_basic_basic(self,i,j,u_ind):#第i,j个基函数在第u_ind个Neumann边界上的乘积
        X = self.Units_Bou_Points[u_ind,:,:]#[2,2]
        basic0 = self.basic(X,[0,0],i)
        
        basic1 = self.basic(X,[0,0],j)
        area = self.area[u_ind,0]
        beta = Beta(X)
        return (beta*basic0*basic1).mean()*area
    def Bou_Neu_basic(self,j,u_ind,prob):
        X = self.Units_Bou_Points[u_ind,:,:]#形状为[2,2],Neumannn边界上的高斯积分点
        area = self.area[u_ind,0]
        n = self.dir[u_ind,:,:]#形状为[2,2],Neumannn边界上的高斯积分点的法方向
        basic = self.basic(X,[0,0],j)
        g = NEU(X,prob,n)
        return (g*basic).mean()*area
    def Node(self):
        self.Nodes_size  = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[0,1]*(self.nx[1,0] + 1)
        
        self.Nodes = np.zeros([self.Nodes_size,2])
        m = 0
        for i in range(self.nx[0,0] + 1):
            for j in range(self.nx[1,0] + self.nx[1,1] + 1):
                self.Nodes[m,0] = self.bounds[0,0] + i*self.hx[0]
                self.Nodes[m,1] = self.bounds[1,0] + j*self.hx[1]
                m = m + 1
        for i in range(self.nx[0,1]):
            for j in range(self.nx[1,0] + 1):
                self.Nodes[m,0] = self.bounds[0,1] + (i + 1)*self.hx[0]
                self.Nodes[m,1] = self.bounds[1,0] + j*self.hx[1]
                m = m + 1
    def Unit(self):
        self.Units_size = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1]) + self.nx[0,1]*self.nx[1,0]
        self.Units_Nodes = np.zeros([self.Units_size,4],np.int)
        self.Units_Int_Points = np.zeros([self.Units_size,
                                          self.gp_num*self.gp_num,self.dim])
        m = 0
        for i in range(self.nx[0,0]):
            for j in range(self.nx[1,0] + self.nx[1,1]):
                self.Units_Nodes[m,0] = i*(self.nx[1,0] + self.nx[1,1] + 1) + j
                self.Units_Nodes[m,1] = i*(self.nx[1,0] + self.nx[1,1] + 1) + j + 1
                self.Units_Nodes[m,2] = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + j
                self.Units_Nodes[m,3] = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + j + 1
                n = 0
                for k in range(self.gp_num):
                    for l in range(self.gp_num):
                        self.Units_Int_Points[m,n,0] = \
                            self.bounds[0,0] + (i + self.gp_pos[k])*self.hx[0]
                        self.Units_Int_Points[m,n,1] = \
                            self.bounds[1,0] + (j + self.gp_pos[l])*self.hx[1]
                        n = n + 1
                #plt.scatter(self.Units_Int_Points[m,:,0],self.Units_Int_Points[m,:,1])
                m = m + 1
        for j in range(self.nx[1,0]):#右边第一列节点
            self.Units_Nodes[m,0] = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + j
            self.Units_Nodes[m,1] = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + j + 1
            self.Units_Nodes[m,2] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + j
            self.Units_Nodes[m,3] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + j + 1
            n = 0
            for k in range(self.gp_num):
                for l in range(self.gp_num):
                    self.Units_Int_Points[m,n,0] = \
                        self.bounds[0,1] + (0 + self.gp_pos[k])*self.hx[0]
                    self.Units_Int_Points[m,n,1] = \
                        self.bounds[1,0] + (j + self.gp_pos[l])*self.hx[1]
                    n = n + 1
            #plt.scatter(self.Units_Int_Points[m,:,0],self.Units_Int_Points[m,:,1])
            m = m + 1
        for i in range(self.nx[0,1] - 1):
            for j in range(self.nx[1,0]):
                self.Units_Nodes[m,0] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1) + j
                self.Units_Nodes[m,1] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1) + j + 1
                self.Units_Nodes[m,2] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (i + 1)*(self.nx[1,0] + 1) + j
                self.Units_Nodes[m,3] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (i + 1)*(self.nx[1,0] + 1) + j + 1
                
                n = 0
                for k in range(self.gp_num):
                    for l in range(self.gp_num):
                        self.Units_Int_Points[m,n,0] = \
                            self.bounds[0,1] + (i + 1 + self.gp_pos[k])*self.hx[0]
                        self.Units_Int_Points[m,n,1] = \
                            self.bounds[1,0] + (j + self.gp_pos[l])*self.hx[1]
                        n = n + 1
                #plt.scatter(self.Units_Int_Points[m,:,0],self.Units_Int_Points[m,:,1])
                m = m + 1
        #下面开始采集neumann边界上的高斯积分点
        self.Neu_size = self.nx[0,0] + self.nx[0,1] + self.nx[1,0] + self.nx[1,1] + self.nx[0,1] + self.nx[1,1] #Neumann边界上的线段数目
        self.Neu_Nodes = np.zeros([self.Neu_size,2],np.int)#存储每一个线段的两个端点序号
        self.Units_Bou_Points = np.zeros([self.Neu_size,self.gp_num,self.dim])#存储线段上的两个积分点坐标
        self.dir = np.zeros([self.Neu_size,self.gp_num,self.dim])#存储Neumann边界上的法方向
        self.area = np.zeros([self.Neu_size,1])#存储Neumann边界上单元的长度
        m = 0
        for i in range(self.nx[0,0]):
            self.Neu_Nodes[m,0] = i*(self.nx[1,0] + self.nx[1,1] + 1)
            self.Neu_Nodes[m,1] = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1)
            self.area[m,0] = self.hx[0]
            for k in range(self.gp_num):
                self.Units_Bou_Points[m,k,0] = self.bounds[0,0] + (i + self.gp_pos[k])*self.hx[0]
                self.Units_Bou_Points[m,k,1] = self.bounds[1,0] 
                self.dir[m,k,0] = 0.0;self.dir[m,k,1] = - 1.0
            #plt.scatter(self.Units_Bou_Points[m,:,0],self.Units_Bou_Points[m,:,1])
            m = m + 1
        for i in range(self.nx[0,1]):
            if i == 0:
                self.Neu_Nodes[m,0] = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) 
                self.Neu_Nodes[m,1] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) 
            else:
                self.Neu_Nodes[m,0] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (i - 1)*(self.nx[1,0] + 1)
                self.Neu_Nodes[m,1] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1)
                
            self.area[m,0] = self.hx[0]
            for k in range(self.gp_num):
                self.Units_Bou_Points[m,k,0] = self.bounds[0,1] + (i + self.gp_pos[k])*self.hx[0]
                self.Units_Bou_Points[m,k,1] = self.bounds[1,0] 
                self.dir[m,k,0] = 0.0;self.dir[m,k,1] = - 1.0
            #plt.scatter(self.Units_Bou_Points[m,:,0],self.Units_Bou_Points[m,:,1])
            m = m + 1
        for j in range(self.nx[1,0] + self.nx[1,1]):
            self.Neu_Nodes[m,0] = j
            self.Neu_Nodes[m,1] = j + 1
            self.area[m,0] = self.hx[1]
            for k in range(self.gp_num):
                self.Units_Bou_Points[m,k,0] = self.bounds[0,0] 
                self.Units_Bou_Points[m,k,1] = self.bounds[1,0] + (j + self.gp_pos[k])*self.hx[1]
                self.dir[m,k,0] = -1.0;self.dir[m,k,1] = 0.0
            #plt.scatter(self.Units_Bou_Points[m,:,0],self.Units_Bou_Points[m,:,1])
            m = m + 1
        for j in range(self.nx[1,1]):
            self.Neu_Nodes[m,0] = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[1,0] + j
            self.Neu_Nodes[m,1] = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[1,0] + j + 1
            
            self.area[m,0] = self.hx[1]
            for k in range(self.gp_num):
                self.Units_Bou_Points[m,k,0] = self.bounds[0,1] 
                self.Units_Bou_Points[m,k,1] = self.bounds[1,1] + (j + self.gp_pos[k])*self.hx[1]
                self.dir[m,k,0] = 1.0;self.dir[m,k,1] = 0.0
            #plt.scatter(self.Units_Bou_Points[m,:,0],self.Units_Bou_Points[m,:,1])
            m = m + 1
        for i in range(self.nx[0,1]):
            if i == 0:
                self.Neu_Nodes[m,0] = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[1,0] #--------------------
                self.Neu_Nodes[m,1] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[1,0]
            else:
                self.Neu_Nodes[m,0] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1) - 1
                self.Neu_Nodes[m,1] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (i + 1)*(self.nx[1,0] + 1) - 1
                
            self.area[m,0] = self.hx[0]
            for k in range(self.gp_num):
                self.Units_Bou_Points[m,k,0] = self.bounds[0,1] + (i + self.gp_pos[k])*self.hx[0]
                self.Units_Bou_Points[m,k,1] = self.bounds[1,1]
                self.dir[m,k,0] = 0.0;self.dir[m,k,1] = 1.0
            #plt.scatter(self.Units_Bou_Points[m,:,0],self.Units_Bou_Points[m,:,1])
            m = m + 1
    
    def matrix(self):
        A = np.zeros([self.Nodes_size,self.Nodes_size])#self.Nodes_size = (M+ 1)*(N + 1)
        for m in range(self.Units_size):#self.Units_size = M*N，第m个区域单元的4个积分点
            for k in range(4):
                ind0 = self.Units_Nodes[m,k]#self.Units_Nodes = [M*N,4],第m个区域中第k个网格点，第k个基函数编号
                for l in range(4):
                    ind1 = self.Units_Nodes[m,l]#self.Units_Nodes = [M*N,4],第m个区域中第l网格点，第k个基函数编号
                    #第m个区域上，两个基函数梯度的乘积的积分a(u,v)
                    A[ind0,ind1] += self.Int_basic_basic(ind0,ind1,m)
        for m in range(self.Neu_size):
            for k in range(2):
                ind0 = self.Neu_Nodes[m,k]
                
                for l in range(2):
                    ind1 = self.Neu_Nodes[m,l]
                    
                    A[ind0,ind1] += self.Bou_basic_basic(ind0,ind1,m)
        for i in range(1,self.nx[0,0]):
            ind = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1) - 1
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,:] = np.zeros([1,self.Nodes_size])
            A[ind,ind] = 1.0
        for j in range(1,self.nx[1,0]):
            ind = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (self.nx[0,1] - 1)*(self.nx[1,0] + 1) + j
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,:] = np.zeros([1,self.Nodes_size])
            A[ind,ind] = 1.0
        return A
    def right(self,prob):
        b = np.zeros([self.Nodes_size,1])
        for m in range(self.Units_size):
            for k in range(4):
                ind = self.Units_Nodes[m,k]#第m个单元区域内第k个网格点，第k个基函数的编号
                b[ind] += self.Int_F_basic(ind,m,prob)
        for m in range(self.Neu_size):
            for k in range(2):
                ind = self.Neu_Nodes[m,k]
                #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
                b[ind] += self.Bou_Neu_basic(ind,m,prob)
        for i in range(1,self.nx[0,0]):
            ind = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1) - 1
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'r*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        for j in range(1,self.nx[1,0]):
            ind = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (self.nx[0,1] - 1)*(self.nx[1,0] + 1) + j
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'r*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        return b
    def Uh(self,X,prob):
        uh = np.zeros(X.shape[0])
        A = self.matrix()
        b = self.right(prob)
        Nodes_V = np.linalg.solve(A,b)
        for i in range(self.Nodes_size):#self.Nodes_size = (M + 1)*(N + 1)
            # 计算数据集 关于 基函数中心 的相对位置
            uh += Nodes_V[i]*self.basic(X,[0,0],i)
        return uh.reshape(-1,1)

bounds = np.array([[0.0,1.0,2.0],[0.0,1.0,2.0]])
hx_tr = [0.2,0.1]    # 网格大小
fenet = FENET(bounds,hx_tr)

#hx_te = hx_tr
hx_te = [0.05,0.05]
teset = TESET(bounds,hx_te)
X = teset.X
prob = 1

u_acc = UU(X,[0,0],prob).reshape(-1,1)

u_pred = fenet.Uh(X,prob)
#print(u_acc - u_pred)
error = max(abs(u_acc - u_pred))
print(error)

问题和疑惑
本人调试了两天代码，发现误差始终都只有 $O (h)$ ，特别是当prob=2的时候，误差太大，如果函数简单一些，比如prob=3,4误差相对较小，但是也好不到哪去，这个留给同僚自己解决。

Dirichlet边界

以下纯粹是本人调试代码的过程，为了弄清楚这个问题的所在，本人把边界全部改成了Dirichlet边界，这样误差能小一些，对于prob=4误差几乎为0，不过同样存在误差大，对不同函数差异大的问题，这里本人把代码留在这，大家体会。

import matplotlib.pyplot as plt
import numpy as np
import time


def UU(X, order,prob):#X表示(x,t)
    if prob==1:
        temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5
        if order[0]==0 and order[1]==0:
            return np.log(temp)
        if order[0]==1 and order[1]==0:#对x求偏导
            return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
        if order[0]==0 and order[1]==1:#对t求偏导
            return temp**(-1) * (20*(X[:,0]+X[:,1]) - 2*(X[:,0]-X[:,1]))
        if order[0]==2 and order[1]==0:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) ** 2 \
                   + temp**(-1) * (22)
        if order[0]==1 and order[1]==1:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
                   * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
                   + temp**(-1) * (18)
        if order[0]==0 and order[1]==2:
            return - temp**(-2) * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) ** 2 \
                   + temp**(-1) * (22)
    if prob==2:
        if order[0]==0 and order[1]==0:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   0.5*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==0:
            return (3*X[:,0]*X[:,0]-1) * \
                   0.5*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==1:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   (np.exp(2*X[:,1])-np.exp(-2*X[:,1]))
        if order[0]==2 and order[1]==0:
            return (6*X[:,0]) * \
                   0.5*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
        if order[0]==1 and order[1]==1:
            return (3*X[:,0]*X[:,0]-1) * \
                   (np.exp(2*X[:,1])-np.exp(-2*X[:,1]))
        if order[0]==0 and order[1]==2:
            return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
                   2*(np.exp(2*X[:,1])+np.exp(-2*X[:,1]))
    if prob==3:
        temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
        temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
        if order[0]==0 and order[1]==0:
            return temp1 * temp2**(-1)
        if order[0]==1 and order[1]==0:
            return (2*X[:,0]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,0])
        if order[0]==0 and order[1]==1:
            return (-2*X[:,1]) * temp2**(-1) + \
                   temp1 * (-1)*temp2**(-2) * (2*X[:,1])
        if order[0]==2 and order[1]==0:
            return (2) * temp2**(-1) + \
                   2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
        if order[0]==1 and order[1]==1:
            return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
        if order[0]==0 and order[1]==2:
            return (-2) * temp2**(-1) + \
                   2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
                   temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
                   temp1 * (-1)*temp2**(-2) * (2)
    if prob == 4:
        if order == [0,0]:
            return X[:,0]*X[:,1]
        if order == [1,0]:
            return X[:,1]
        if order == [0,1]:
            return X[:,0]
        if order == [0,2]:
            return 0*X[:,0]
        if order == [2,0]:
            return 0*X[:,1]
def FF(prob,X):
    return -UU(X,[0,2],prob) - UU(X,[2,0],prob)
def Beta(X):
    return 0*X[:,0]
    #return np.exp(X[:,0])*X[:,1]**2
def NEU(X,prob,n):#定义g(x),beta > 0
    return UU(X,[1,0],prob)*n[:,0] + UU(X,[0,1],prob)*n[:,1] + Beta(X)*UU(X,[0,0],prob)
np.random.seed(1234)


class TESET():
    def __init__(self, bounds, hx):
        self.bounds = bounds
        
        nx = [[int((bounds[0,1] - bounds[0,0])/hx[0]),
               int((bounds[0,2] - bounds[0,1])/hx[0])],
              [int((bounds[1,1] - bounds[1,0])/hx[1]),
               int((bounds[1,2] - bounds[1,1])/hx[1])]]
        self.nx = np.array(nx)
        
        self.size = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[0,1]*(self.nx[1,0] + 1)
        self.X = np.zeros([self.size,2])
        m = 0
        for i in range(self.nx[0,0] + 1):
            for j in range(self.nx[1,0] + self.nx[1,1] + 1):
                self.X[m,0] = bounds[0,0] + i*hx[0]
                self.X[m,1] = bounds[1,0] + j*hx[1]
                m = m + 1
        for i in range(self.nx[0,1]):
            for j in range(self.nx[1,0] + 1):
                self.X[m,0] = bounds[0,1] + (i + 1)*hx[0]
                self.X[m,1] = bounds[1,0] + j*hx[1]
                m = m + 1
class FENET():
    def __init__(self,bounds,hx):
        self.bounds = bounds
        self.dim = 2
        nx = [[int((bounds[0,1] - bounds[0,0])/hx[0]),
               int((bounds[0,2] - bounds[0,1])/hx[0])],
              [int((bounds[1,1] - bounds[1,0])/hx[1]),
               int((bounds[1,2] - bounds[1,1])/hx[1])]]
        self.nx = np.array(nx)
        self.hx = hx
        self.gp_num = 2
        self.gp_pos = [(1 - np.sqrt(3)/3)/2,(1 + np.sqrt(3)/3)/2]
        self.Node()
        self.Unit()
    def phi(self,X,order):#[-1,1]*[-1,1]，在原点取值为1，其他网格点取值为0的基函数
        ind00 = (X[:,0] >= -1);ind01 = (X[:,0] >= 0);ind02 = (X[:,0] >= 1)
        ind10 = (X[:,1] >= -1);ind11 = (X[:,1] >= 0);ind12 = (X[:,1] >= 1)
        if order == [0,0]:
            return (ind00*~ind01*ind10*~ind11).astype('float32')*(1 + X[:,0])*(1 + X[:,1]) + \
                    (ind01*~ind02*ind10*~ind11).astype('float32')*(1 - X[:,0])*(1 + X[:,1]) + \
                    (ind00*~ind01*ind11*~ind12).astype('float32')*(1 + X[:,0])*(1 - X[:,1]) + \
                    (ind01*~ind02*ind11*~ind12).astype('float32')*(1 - X[:,0])*(1 - X[:,1])
        if order == [1,0]:
            return (ind00*~ind01*ind10*~ind11).astype('float32')*(1 + X[:,1]) + \
                    (ind01*~ind02*ind10*~ind11).astype('float32')*(-(1 + X[:,1])) + \
                    (ind00*~ind01*ind11*~ind12).astype('float32')*(1 - X[:,1]) + \
                    (ind01*~ind02*ind11*~ind12).astype('float32')*(-(1 - X[:,1]))
        if order == [0,1]:
            return (ind00*~ind01*ind10*~ind11).astype('float32')*(1 + X[:,0]) + \
                    (ind01*~ind02*ind10*~ind11).astype('float32')*(1 - X[:,0]) + \
                    (ind00*~ind01*ind11*~ind12).astype('float32')*(-(1 + X[:,0])) + \
                    (ind01*~ind02*ind11*~ind12).astype('float32')*(-(1 - X[:,0]))
    def basic(self,X,order,i):#根据网格点的存储顺序，遍历所有网格点，取基函数
        temp = (X - self.Nodes[i,:])/np.array([self.hx[0],self.hx[1]])
        if order == [0,0]:
            return self.phi(temp,order)
        if order == [1,0]:
            return self.phi(temp,order)/self.hx[0]
        if order == [0,1]:
            return self.phi(temp,order)/self.hx[1]
    def Int_basic_basic(self,i,j,u_ind):#表示第i,j个基函数的梯度的乘积，以及在第u_ind个区域的积分
        X = self.Units_Int_Points[u_ind,:,:]#[4,2]
        basic0 = np.zeros_like(X)
        basic0[:,0] = self.basic(X,[1,0],i)
        basic0[:,1] = self.basic(X,[0,1],i)
        
        basic1 = np.zeros_like(X)
        basic1[:,0] = self.basic(X,[1,0],j)
        basic1[:,1] = self.basic(X,[0,1],j)
        return ((basic0*basic1).sum(1)).mean()*self.hx[0]*self.hx[1]
    def Int_F_basic(self,i,u_ind,prob):#第i个基函数与右端项乘积，在第u_ind个单元积分
        X = self.Units_Int_Points[u_ind,:,:]
        return (FF(prob,X)*self.basic(X,[0,0],i)).mean()*self.hx[0]*self.hx[1]
    
    def Node(self):
        self.Nodes_size  = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[0,1]*(self.nx[1,0] + 1)
        
        self.Nodes = np.zeros([self.Nodes_size,2])
        m = 0
        for i in range(self.nx[0,0] + 1):
            for j in range(self.nx[1,0] + self.nx[1,1] + 1):
                self.Nodes[m,0] = self.bounds[0,0] + i*self.hx[0]
                self.Nodes[m,1] = self.bounds[1,0] + j*self.hx[1]
                m = m + 1
        for i in range(self.nx[0,1]):
            for j in range(self.nx[1,0] + 1):
                self.Nodes[m,0] = self.bounds[0,1] + (i + 1)*self.hx[0]
                self.Nodes[m,1] = self.bounds[1,0] + j*self.hx[1]
                m = m + 1
    def Unit(self):
        self.Units_size = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1]) + self.nx[0,1]*self.nx[1,0]
        self.Units_Nodes = np.zeros([self.Units_size,4],np.int)
        self.Units_Int_Points = np.zeros([self.Units_size,
                                          self.gp_num*self.gp_num,self.dim])
        m = 0
        for i in range(self.nx[0,0]):
            for j in range(self.nx[1,0] + self.nx[1,1]):
                self.Units_Nodes[m,0] = i*(self.nx[1,0] + self.nx[1,1] + 1) + j
                self.Units_Nodes[m,1] = i*(self.nx[1,0] + self.nx[1,1] + 1) + j + 1
                self.Units_Nodes[m,2] = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + j
                self.Units_Nodes[m,3] = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + j + 1
                n = 0
                for k in range(self.gp_num):
                    for l in range(self.gp_num):
                        self.Units_Int_Points[m,n,0] = \
                            self.bounds[0,0] + (i + self.gp_pos[k])*self.hx[0]
                        self.Units_Int_Points[m,n,1] = \
                            self.bounds[1,0] + (j + self.gp_pos[l])*self.hx[1]
                        n = n + 1
                #plt.scatter(self.Units_Int_Points[m,:,0],self.Units_Int_Points[m,:,1])
                m = m + 1
        for j in range(self.nx[1,0]):#右边第一列节点
            self.Units_Nodes[m,0] = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + j
            self.Units_Nodes[m,1] = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + j + 1
            self.Units_Nodes[m,2] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + j
            self.Units_Nodes[m,3] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + j + 1
            n = 0
            for k in range(self.gp_num):
                for l in range(self.gp_num):
                    self.Units_Int_Points[m,n,0] = \
                        self.bounds[0,1] + (0 + self.gp_pos[k])*self.hx[0]
                    self.Units_Int_Points[m,n,1] = \
                        self.bounds[1,0] + (j + self.gp_pos[l])*self.hx[1]
                    n = n + 1
            #plt.scatter(self.Units_Int_Points[m,:,0],self.Units_Int_Points[m,:,1])
            m = m + 1
        for i in range(self.nx[0,1] - 1):
            for j in range(self.nx[1,0]):
                self.Units_Nodes[m,0] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1) + j
                self.Units_Nodes[m,1] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1) + j + 1
                self.Units_Nodes[m,2] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (i + 1)*(self.nx[1,0] + 1) + j
                self.Units_Nodes[m,3] = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (i + 1)*(self.nx[1,0] + 1) + j + 1
                
                n = 0
                for k in range(self.gp_num):
                    for l in range(self.gp_num):
                        self.Units_Int_Points[m,n,0] = \
                            self.bounds[0,1] + (i + 1 + self.gp_pos[k])*self.hx[0]
                        self.Units_Int_Points[m,n,1] = \
                            self.bounds[1,0] + (j + self.gp_pos[l])*self.hx[1]
                        n = n + 1
                #plt.scatter(self.Units_Int_Points[m,:,0],self.Units_Int_Points[m,:,1])
                m = m + 1
        
    
    def matrix(self):
        A = np.zeros([self.Nodes_size,self.Nodes_size])#self.Nodes_size = (M+ 1)*(N + 1)
        for m in range(self.Units_size):#self.Units_size = M*N，第m个区域单元的4个积分点
            for k in range(4):
                ind0 = self.Units_Nodes[m,k]#self.Units_Nodes = [M*N,4],第m个区域中第k个网格点，第k个基函数编号
                for l in range(4):
                    ind1 = self.Units_Nodes[m,l]#self.Units_Nodes = [M*N,4],第m个区域中第l网格点，第k个基函数编号
                    #第m个区域上，两个基函数梯度的乘积的积分a(u,v)
                    A[ind0,ind1] += self.Int_basic_basic(ind0,ind1,m)
        for j in range(self.nx[1,0] + self.nx[1,1] + 1):
            ind = j
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,:] = np.zeros([1,self.Nodes_size])
            A[ind,ind] = 1.0
        for i in range(self.nx[0,0] + 1):
            ind = i*(self.nx[1,0] + self.nx[1,1] + 1)
            A[ind,:] = np.zeros([1,self.Nodes_size])
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,ind] = 1.0
        for i in range(self.nx[0,1]):
            ind = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1)
            A[ind,:] = np.zeros([1,self.Nodes_size])
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,ind] = 1.0
        for j in range(self.nx[1,1] + 1):
            ind = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[1,0] + j
            A[ind,:] = np.zeros([1,self.Nodes_size])
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,ind] = 1.0
        for i in range(self.nx[0,1] + 1):
            ind = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1) - 1
            A[ind,:] = np.zeros([1,self.Nodes_size])
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,ind] = 1.0
        for i in range(1,self.nx[0,0]):
            ind = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1) - 1
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,:] = np.zeros([1,self.Nodes_size])
            A[ind,ind] = 1.0
        for j in range(self.nx[1,0] + 1):
            ind = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (self.nx[0,1] - 1)*(self.nx[1,0] + 1) + j
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            A[ind,:] = np.zeros([1,self.Nodes_size])
            A[ind,ind] = 1.0
        
        return A
    def right(self,prob):
        b = np.zeros([self.Nodes_size,1])
        for m in range(self.Units_size):
            for k in range(4):
                ind = self.Units_Nodes[m,k]#第m个单元区域内第k个网格点，第k个基函数的编号
                b[ind] += self.Int_F_basic(ind,m,prob)
        
        for j in range(self.nx[1,0] + self.nx[1,1] + 1):
            ind = j
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        for i in range(self.nx[0,0] + 1):
            ind = i*(self.nx[1,0] + self.nx[1,1] + 1)
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        for i in range(self.nx[0,1]):
            ind = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1)
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        for j in range(self.nx[1,1] + 1):
            ind = self.nx[0,0]*(self.nx[1,0] + self.nx[1,1] + 1) + self.nx[1,0] + j
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        for i in range(self.nx[0,1] + 1):
            ind = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + i*(self.nx[1,0] + 1) - 1
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        for i in range(1,self.nx[0,0]):
            ind = (i + 1)*(self.nx[1,0] + self.nx[1,1] + 1) - 1
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        for j in range(self.nx[1,0] + 1):
            ind = (self.nx[0,0] + 1)*(self.nx[1,0] + self.nx[1,1] + 1) + (self.nx[0,1] - 1)*(self.nx[1,0] + 1) + j
            #plt.plot(self.Nodes[ind,0],self.Nodes[ind,1],'b*')
            b[ind,0] = UU(self.Nodes[ind:ind + 1,:],[0,0],prob)
        return b
    def Uh(self,X,prob):
        uh = np.zeros(X.shape[0])
        A = self.matrix()
        b = self.right(prob)
        Nodes_V = np.linalg.solve(A,b)
        for i in range(self.Nodes_size):#self.Nodes_size = (M + 1)*(N + 1)
            # 计算数据集 关于 基函数中心 的相对位置
            uh += Nodes_V[i]*self.basic(X,[0,0],i)
        return uh.reshape(-1,1)

bounds = np.array([[0.0,1.0,2.0],[0.0,1.0,2.0]])
hx_tr = [0.2,0.1]    # 网格大小
fenet = FENET(bounds,hx_tr)

#hx_te = hx_tr
hx_te = [0.05,0.05]
teset = TESET(bounds,hx_te)
X = teset.X
prob = 4

u_acc = UU(X,[0,0],prob).reshape(-1,1)

u_pred = fenet.Uh(X,prob)
#print(u_acc - u_pred)
error = max(abs(u_acc - u_pred))
print(error)

有限元的探讨到此为止，接下本人目标投向多重网格方法。

Galerkin码农选手

关注

3
点赞
踩
2

收藏

觉得还不错? 一键收藏
打赏
0
评论
不规则区域上有限元方法的进一步探讨和疑惑

L型区域上的有限元方法之前的传统数值方法基本上都是针对规则的矩形区域，按道理有限元方法的应用就是为不规则区域量身打造，为此本人尝试使用有限元方法来求解L型区域上的泊松方程，使用的同样是Galerkin原理。混合边界首先我们先尝试求解混合边界上的泊松方程，参考以下图像：其中蓝色代表Dirichlet边界，橙色代表Neumann边界，有限元的基础知识本人已经在上一篇博客介绍有限元因此我就不花篇幅介绍基函数以及试探检验函数空间的选取了，本人使用的仍然是矩形有限元基函数。下面本人重点介绍这样的L型区域在
复制链接

扫一扫