Адаптивная очистка сетки - 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])