自适应网格细化 - Python

Question

我目前对我的代码中不起作用的部分非常困惑，并且已经盯着它看了好几个小时。我已经创建了一些函数来使用有限元法自适应地逼近拉普拉斯方程的解，然后使用对偶加权残差估计它的误差。错误函数应该给出一个错误向量（每个元素一个错误），然后我选择最大的错误，在它们周围添加更多元素，再次求解，然后重新检查错误；但是我不知道为什么我的错误估计没有改变！

我的前 4 个函数是正确的，但我会包括它们以防有人想尝试代码：

def Poisson_Stiffness(x0):
    """Finds the Poisson equation stiffness matrix with any non uniform mesh x0"""

    x0 = np.array(x0)
    N = len(x0) - 1 # The amount of elements; x0, x1, ..., xN

    h = x0[1:] - x0[:-1]

    a = np.zeros(N+1)
    a[0] = 1 #BOUNDARY CONDITIONS
    a[1:-1] = 1/h[1:] + 1/h[:-1]
    a[-1] = 1/h[-1]
    a[N] = 1 #BOUNDARY CONDITIONS

    b = -1/h
    b[0] = 0 #BOUNDARY CONDITIONS

    c = -1/h
    c[N-1] = 0 #BOUNDARY CONDITIONS: DIRICHLET

    data = [a.tolist(), b.tolist(), c.tolist()]
    Positions = [0, 1, -1]
    Stiffness_Matrix = diags(data, Positions, (N+1,N+1))

    return Stiffness_Matrix

def NodalQuadrature(x0):
    """Finds the Nodal Quadrature Approximation of sin(pi x)"""

    x0 = np.array(x0)
    h = x0[1:] - x0[:-1]
    N = len(x0) - 1

    approx = np.zeros(len(x0))
    approx[0] = 0 #BOUNDARY CONDITIONS

    for i in range(1,N):
        approx[i] = math.sin(math.pi*x0[i])
        approx[i] = (approx[i]*h[i-1] + approx[i]*h[i])/2

    approx[N] = 0 #BOUNDARY CONDITIONS

    return approx

def Solver(x0):

    Stiff_Matrix = Poisson_Stiffness(x0)

    NodalApproximation = NodalQuadrature(x0)
    NodalApproximation[0] = 0

    U = scipy.sparse.linalg.spsolve(Stiff_Matrix, NodalApproximation)

    return U

def Dualsolution(rich_mesh,qoi_rich_node): #BOUNDARY CONDITIONS?
    """Find Z from stiffness matrix Z = K^-1 Q over richer mesh"""

    K = Poisson_Stiffness(rich_mesh)
    Q = np.zeros(len(rich_mesh))
    Q[qoi_rich_node] = 1.0

    Z = scipy.sparse.linalg.spsolve(K,Q)
    return Z

我的错误指示器函数采用近似值 Uh，用它求解的网格，发现 eta = (f - Bu)z。

def Error_Indicators(Uh,U_mesh,Z,Z_mesh,f):
    """Take in U, Interpolate to same mesh as Z then solve for eta vector"""

    u_inter = interp1d(U_mesh,Uh) #Interpolation of old mesh
    U2 = u_inter(Z_mesh) #New function u for the new mesh to use in 

    Bz = Poisson_Stiffness(Z_mesh)
    Bz = Bz.tocsr()

    eta = np.empty(len(Z_mesh))
    for i in range(len(Z_mesh)):
        for j in range(len(Z_mesh)):
            eta[i] += (f[i] - Bz[i,j]*U2[j])

    for i in range(len(Z)):
        eta[i] = eta[i]*Z[i]

    return eta

我的下一个函数似乎可以很好地使网格适应给定的错误指示器！只是不知道为什么指标似乎无论如何都保持不变？

def Mesh_Refinement(base_mesh,tolerance,refinement,z_mesh,QOI_z_mesh):
    """Solve for U on a normal mesh, Take in Z, Find error indicators, adapt. OUTPUT NEW MESH"""
    New_mesh = base_mesh
    Z = Dualsolution(z_mesh,QOI_z_mesh) #Solve dual solution only once

    f = np.empty(len(z_mesh))
    for i in range(len(z_mesh)):
        f[i] = math.sin(math.pi*z_mesh[i])

    U = Solver(New_mesh)
    eta = Error_Indicators(U,base_mesh,Z,z_mesh,f) 

    while max(abs(k) for k in eta) > tolerance:

        orderedeta = np.sort(eta) #Sort error indicators LENGTH 40
        biggest = np.flipud(orderedeta[int((1-refinement)*len(eta)):len(eta)]) 
        position = np.empty(len(biggest))
        ratio = float(len(New_mesh))/float(len(z_mesh))

        for i in range(len(biggest)):
            position[i] = eta.tolist().index(biggest[i])*ratio #GIVES WHAT NUMBER NODE TO REFINE

        refine = np.zeros(len(position))
        for i in range(len(position)):
            refine[i] = math.floor(position[i])+0.5 #AT WHAT NODE TO PUT NEW ELEMENT 5.5 ETC
        refine = np.flipud(sorted(set(refine)))

        for i in range(len(refine)):
            New_mesh = np.insert(New_mesh,refine[i]+0.5,(New_mesh[refine[i]+0.5]+New_mesh[refine[i]-0.5])/2)

        U = Solver(New_mesh)
        eta = Error_Indicators(U,New_mesh,Z,z_mesh,f)
        print eta

一个示例输入是： Mesh_Refinement(np.linspace(0,1,3),0.1,0.2,np.linspace(0,1,60),20)

我知道这里有很多代码，但我不知所措，我不知道该去哪里！

Answer 1

请考虑来自 def Error_Indicators 的这段代码：

eta = np.empty(len(Z_mesh))
for i in range(len(Z_mesh)):
    for j in range(len(Z_mesh)):
        eta[i] = (f[i] - Bz[i,j]*U2[j])

这里你覆盖了eta[i]每次j次迭代，所以内循环被证明是无用的，你可以直接转到最后一个可能的j。您是要查找 (f[i] - Bz[i,j]*U2[j]) 系列的总和吗？

eta = np.empty(len(Z_mesh))
for i in range(len(Z_mesh)):
    for j in range(len(Z_mesh)):
        eta[i] += (f[i] - Bz[i,j]*U2[j])

自适应网格细化 - Python

Adaptive mesh refinement - Python

python

math

mesh

finite-element-analysis