自适应网格的误差估计指标

Question

我一直在尝试为拉普拉斯方程的有限元法建立一个自适应网格（因此AU=F，其中A是刚度矩阵），但是在误差估计方面遇到了一些问题，因此继续使用如下网格：

[0。0.11111111 0.22222222 0.27777778 0.33333333 0.44444444 0.47222222 0.47916667 0.48090278 0.48111979 0.48133681 0.48177083 0.48263889 0.48611111 0.5 0.55555556 0.66666667 0.77777778 0.88888889 1. ]你知道吗

正如您所看到的，由于此时的误差估计值没有减少，许多节点都在同一个位置上分组重点。重点误差估计的基本思想是我们有两个网格，我们的自适应网格和一个丰富/更大的网格。我们在丰富网格上找到A，Z和F的近似值，在自适应网格上找到U，然后在精细网格上插值，并计算如下：

eta_i=（f_i-总和_j（A[i，j]U[j]））*z[i]

一个大的eta意味着元素需要改进。这是我的代码，其中我确定泊松刚度（找到A），节点求积（找到F），解算器（找到U）和对偶解（找到Z）都是正确的。我无法判断问题是在SubdivisionIndicators函数中还是在网格细化中，但是我假设是前者，因为通常问题源于错误没有减少，因此在同一个位置周围的元素堆积。你知道吗

import numpy as np
import math
import scipy
from scipy.sparse import diags
import scipy.sparse.linalg
from scipy.interpolate import interp1d
import matplotlib
import matplotlib.pyplot as plt

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): 
"""Find Z from stiffness matrix Z = K^-1 Q over richer mesh"""

K = Poisson_Stiffness(rich_mesh)
Q = np.zeros(len(rich_mesh))
if qoi in rich_mesh:

    qoi_rich_node = rich_mesh.tolist().index(qoi)
    Q[qoi_rich_node] = 1.0
else:

    nearest = find_nearest(rich_mesh,qoi)
    if nearest < qoi:

        qoi_lower_node = rich_mesh.tolist().index(nearest)
        qoi_upper_node = qoi_lower_node+1

    else:

        qoi_upper_node = rich_mesh.tolist().index(nearest)
        qoi_lower_node = qoi_upper_node-1

    Q[qoi_upper_node] = (qoi - rich_mesh[qoi_lower_node])/(rich_mesh[qoi_upper_node]-rich_mesh[qoi_lower_node])
    Q[qoi_lower_node] = (rich_mesh[qoi_upper_node] - qoi)/(rich_mesh[qoi_upper_node]-rich_mesh[qoi_lower_node])
Z = scipy.sparse.linalg.spsolve(K,Q)
return Z

def SubdivisionIndicators(richx0,basex0,qoi):
"""interpolate U, eta_i = r_i z_i = (f_i - sum(AijUj)Zi"""
U = Solver(basex0) #finds U over base mesh
U_inter = interp1d(basex0,U,kind='cubic') #interpolates U over rich mesh
U_rich = U_inter(richx0)
f = NodalQuadrature(richx0) #finds f over rich mesh
A = Poisson_Stiffness(richx0).tocsr() #finds A over rich mesh
Z = Dualsolution(richx0,qoi) #finds Z over rich mesh

eta = np.empty(len(richx0))
for i in range(len(eta)):
    eta[i] = abs((f[i] - A[i,:].dot(U_rich))*Z[i])
return eta


def Mesh_refine(base_mesh,amount_of_refinements,z_mesh,QOI):
refinedmesh = base_mesh #initialise

for i in range(amount_of_refinements):
    eta = SubdivisionIndicators(z_mesh,refinedmesh,QOI)
    maxpositionseta = np.nonzero(eta == eta.max())[0]
    maxpositionseta = np.array(maxpositionseta)
    ratio = float(len(refinedmesh))/float(len(z_mesh))

    for i in range(len(maxpositionseta)):
        print eta[maxpositionseta[i]]
        maxposition = maxpositionseta[i]*ratio
        if maxposition == int(maxposition):
            refinedmesh = np.insert(refinedmesh,maxposition,(refinedmesh[maxposition] + refinedmesh[maxposition-1])/2)
        else:
            refinedmesh = np.insert(refinedmesh,math.ceil(maxposition),(refinedmesh[math.ceil(maxposition)]+refinedmesh[math.floor(maxposition)])/2)

    print refinedmesh
return refinedmesh

网格优化(np.L空间(0,1,10),10,np.L空间（0.1100），0.5）