ex6

ex6/code/adapt.py

@@ -0,0 +1,83 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.integrate as integrate
+
+# Basis functions
+# ---- TODO: order p? ----
+def phi(k, mesh, x):
+    if x > mesh[k-1] and x < mesh[k]:
+        return (x-mesh[k-1]) / (mesh[k] - mesh[k-1])
+    elif x > mesh[k] and x < mesh[k+1]:
+        return (mesh[k+1]-x) / (mesh[k+1] - mesh[k])
+    elif x == mesh[k]:
+        return 1
+    else:
+        return 0
+def phi_prime(k, mesh, x):
+    if x > mesh[k-1] and x < mesh[k]:
+        return 1 / (mesh[k] - mesh[k-1])
+    elif x > mesh[k] and x < mesh[k+1]:
+        return -1 / (mesh[k+1] - mesh[k])
+    elif x == mesh[k]:
+        return print("oh no")
+    else:
+        return 0
+    
+# vertex flux jump between elements
+def flux_jumps(mesh, u):
+    m = len(mesh)
+
+    r = (u[2:]-u[1:-1]) / (mesh[2:]-mesh[1:-1])
+    l = (u[1:-1]-u[:-2]) / (mesh[1:-1]-mesh[:-2])
+    jumps = np.zeros(m)                      # 0 jump at bnd?
+    jumps[1:-1] = r - l
+
+    jumps[0] = jumps[1]                      # or same jump at bnd?
+    jumps[-1] = jumps[-2]                    # or same jump at bnd?
+    return jumps
+
+# h-adaptivity (refine neighboring elments, if flux jump over certain threshold)
+def adapt_h(mesh, jumps):
+    new_mesh = [mesh[0]]
+    jumps = jumps[1:-1]                      # only interior jumps (excluding bnd nodes needed for De Boor)
+    
+    # define threshold
+    threshold = 0.5 * np.max(np.abs(jumps))
+    # mark elements
+    marked_nodes = np.abs(jumps) > threshold
+    marked_el = np.zeros(len(mesh)-1, dtype=bool)
+    for i, refine in enumerate(marked_nodes):
+        if refine:
+            marked_el[i] = True
+            marked_el[i+1] = True
+
+    # make new mesh
+    for k, refine in enumerate(marked_el):
+        l = mesh[k]
+        r = mesh[k+1]
+        if refine:
+            new_mesh.extend(np.linspace(l, r, 3)[1:])
+        else:
+            new_mesh.append(r)
+
+    return np.asarray(new_mesh, dtype=float)
+
+# r-adaptivity (one iteration of De Boor's algorithm, moving mesh nodes for equidistributing mesh)
+def adapt_r(mesh, rho):
+    m = len(mesh)
+    
+    p = 0.5 * (rho[:-1] + rho[1:])                  # piecewise constant function on elements
+    
+    P = np.zeros(m)                             
+    for i in range(1,m):
+        P[i] = P[i-1] + (mesh[i]-mesh[i-1])*p[i-1]  # approx integral of p, from node 0 to node i
+    Pb = P[-1]                                      # integral of p, over whole mesh 
+
+    new_mesh = mesh.copy()
+    for j in range(1, m-1):
+        xi_j = (j)/(m-1)
+        k = np.searchsorted(P, xi_j*Pb)             # search k s.t.: P[node k-1] < xi_j*Pb < P[node k]
+        k = max(k,1)                                # if k=0
+        new_mesh[j] = mesh[k-1] + 2*(xi_j*Pb - P[k-1]) / (rho[k-1]+rho[k])    # new node  j
+    return new_mesh
+

ex6/code/task_a.py

@@ -0,0 +1,104 @@
+import numpy as np
+from adapt import *
+
+# PDE:
+#  -u''(x) = f(x)       x in (-1,1)
+#    u(-1) = -arctan(p)
+#    u'(1) = p / (p^2 + 1)
+# 
+# weak form
+# int u'v' dx = int f(x) * v(x) dx + p/(p^2+1) * v(1)
+# 
+
+# rhs
+def f(x):
+    return 2 * p**3 * x / (p**2 * x**2 + 1)**2
+
+# Stiffness and Load
+def K_loc(k, mesh):
+    K_loc = np.zeros((2,2))
+    K_loc[0,0] = integrate.quad(lambda x: phi_prime(k-1, mesh, x)**2, mesh[k-1], mesh[k])[0]
+    K_loc[1,0] = integrate.quad(lambda x: phi_prime(k-1, mesh, x)*phi_prime(k, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[0,1] = integrate.quad(lambda x: phi_prime(k, mesh, x)*phi_prime(k-1, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[1,1] = integrate.quad(lambda x: phi_prime(k, mesh, x)**2, mesh[k-1], mesh[k])[0]
+    return K_loc
+def F_loc(k, mesh):
+    F_loc = np.zeros(2)
+    F_loc[0] = integrate.quad(lambda x: f(x) * phi(k-1, mesh, x), mesh[k-1], mesh[k])[0]
+    F_loc[1] = integrate.quad(lambda x: f(x) * phi(k, mesh, x), mesh[k-1], mesh[k])[0]
+    return F_loc
+
+# Assembling
+def Assemble(mesh):
+    m = len(mesh)
+    n = m-1
+
+    K = np.zeros((m,m))
+    F = np.zeros(m)
+    for k in range(1,m):
+        K[k-1:k+1,k-1:k+1] += K_loc(k, mesh)
+        F[k-1:k+1] += F_loc(k, mesh)
+
+    # Boundary conditions
+    # Dirichlet: u(-1) = -arctan(p)
+    K[0,:] = 0
+    K[0,0] = 1
+    F[0]   = -np.arctan(p)
+    # Neumann: u'(1) = p / (p^2 + 1)
+    F[-1] += p/(p**2+1) * phi(n, mesh, 1)
+    return K,F
+
+def plotting(mesh, u, comment):
+    exact_x = np.linspace(-1,1,1000)
+    exact = np.arctan(p*exact_x)
+    plt.plot(exact_x, exact, "--", linewidth=1, color="red", label="exact")
+    plt.title(f"p = {p} | n = {n} | {comment}")
+    plt.xlabel("x")
+    plt.ylabel("u(x)")
+    plt.plot(mesh, u, "-o", label="u_h")
+    plt.xticks(mesh, labels=[])
+    plt.legend()
+    plt.grid(True)
+    plt.tight_layout()
+    plt.savefig("task_a.png", dpi=300)
+    plt.show()
+    return 0
+
+##############################################################################
+
+# parameters
+p_list = [5,10,20,100]
+for p in p_list:
+
+    # h-adaptivity
+    mesh = np.array([-1.0, -0.2, 0, 0.7, 1.0])            # with 0 as node
+    # mesh = np.array([-1.0,-0.131,0.372,1.0])              # without 0 as node
+
+    # r-adaptivity
+    # mesh = np.linspace(-1, 1, 11)                         # with 0 as node
+    # mesh = np.linspace(-1, 1, 10)                         # without 0 as node
+
+    m = len(mesh)
+    n = m-1
+
+    K, F = Assemble(mesh)                                 # assemble
+    u = np.linalg.solve(K, F)                             # solve
+    jumps = flux_jumps(mesh, u)                           # flux jumps
+    plotting(mesh, u, "before adapting")                  # plotting
+
+    iterations = 6
+    for it in range(iterations):
+        mesh = adapt_h(mesh, jumps)                       # h-adaptivity
+        # mesh = adapt_r(mesh, np.abs(jumps))               # r-adaptivity (positive density (jumps)!)
+        m = len(mesh)
+        n = m-1
+
+        K, F = Assemble(mesh)                             # assemble
+        u = np.linalg.solve(K, F)                         # solve
+        jumps = flux_jumps(mesh, u)                       # flux jumps
+        # plotting(mesh, u, f"iteration {it+1}")            # plotting each iteration
+
+    # print(jumps)
+    plotting(mesh, u, f"after {iterations} iterations")   # plotting
+
+

ex6/code/task_b.py

@@ -0,0 +1,100 @@
+import numpy as np
+from adapt import *
+
+# PDE:
+#  -(lambda(x)u'(x))' = 0         x in (0,1)
+#                u(0) = 0
+#                u(1) = 1
+#           lambda(x) = |  1     x in (0,1/sqrt(2))
+#                       | 10     x in (1/sqrt(2),1)   
+
+# rhs
+def f(x):
+    return 0
+
+def lam(x):
+    if x >= 0 and x <= 1/np.sqrt(2):
+        return 1
+    elif x <= 1 and x > 1/np.sqrt(2):
+        return 10
+    else:
+        return print("lambda undefined")
+
+# Stiffness and Load
+def K_loc(k, mesh):
+    K_loc = np.zeros((2,2))
+    K_loc[0,0] = integrate.quad(lambda x: lam(x)*phi_prime(k-1, mesh, x)**2, mesh[k-1], mesh[k])[0]
+    K_loc[1,0] = integrate.quad(lambda x: lam(x)*phi_prime(k-1, mesh, x)*phi_prime(k, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[0,1] = integrate.quad(lambda x: lam(x)*phi_prime(k, mesh, x)*phi_prime(k-1, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[1,1] = integrate.quad(lambda x: lam(x)*phi_prime(k, mesh, x)**2, mesh[k-1], mesh[k])[0]
+    return K_loc
+def F_loc(k, mesh):
+    F_loc = np.zeros(2)
+    F_loc[0] = integrate.quad(lambda x: f(x) * phi(k-1, mesh, x), mesh[k-1], mesh[k])[0]
+    F_loc[1] = integrate.quad(lambda x: f(x) * phi(k, mesh, x), mesh[k-1], mesh[k])[0]
+    return F_loc
+
+# Assembling
+def Assemble(mesh):
+    m = len(mesh)
+    n = m-1
+
+    K = np.zeros((m,m))
+    F = np.zeros(m)
+    for k in range(1,m):
+        K[k-1:k+1,k-1:k+1] += K_loc(k, mesh)
+        F[k-1:k+1] += F_loc(k, mesh)
+
+    # Boundary conditions
+    # Dirichlet: u(0) = 0
+    K[0,:] = 0
+    K[0,0] = 1
+    F[0]   = 0
+    # Dirichlet: u(1) = 1
+    K[-1,:] = 0
+    K[-1,-1] = 1
+    F[-1]   = 1
+    return K,F
+
+def plotting(mesh, u, comment):
+    exact_x = [0, 1/np.sqrt(2), 1]
+    exact = [0, 10/(np.sqrt(2)+9), 1]
+    plt.plot(exact_x, exact, "--", linewidth=1, color="red", label="exact")
+    plt.title(f"n = {n} | {comment}")
+    plt.xlabel("x")
+    plt.ylabel("u(x)")
+    plt.plot(mesh, u, "-o", label="u_h")
+    plt.xticks(mesh, labels=[])
+    plt.legend()
+    plt.grid(True)
+    plt.tight_layout()
+    plt.savefig("task_b.png", dpi=300)
+    plt.show()
+    return 0
+
+##############################################################################
+
+mesh = np.linspace(0, 1, 10)
+m = len(mesh)
+n = m-1
+
+K, F = Assemble(mesh)                                 # assemble
+u = np.linalg.solve(K, F)                             # solve
+jumps = flux_jumps(mesh, u)                           # flux jumps
+plotting(mesh, u, "before adapting")                  # plotting
+
+iterations = 3
+for it in range(iterations):
+    # mesh = adapt_h(mesh, jumps)                       # h-adaptivity
+    mesh = adapt_r(mesh, np.abs(jumps))               # r-adaptivity (positive density (jumps)!)
+    m = len(mesh)
+    n = m-1
+
+    K, F = Assemble(mesh)                             # assemble
+    u = np.linalg.solve(K, F)                         # solve
+    jumps = flux_jumps(mesh, u)                       # flux jumps
+    # plotting(mesh, u, f"iteration {it+1}")            # plotting each iteration
+
+# print(jumps)
+plotting(mesh, u, f"after {iterations} iterations")   # plotting
+

ex6/code/task_c.py

@@ -0,0 +1,100 @@
+import numpy as np
+from adapt import *
+
+# Peclet problem:
+#  -u''(x) + pu'(x) = 0    x in (0,1)
+#              u(0) = 0
+#              u(1) = 1
+
+# rhs
+def f(x):
+    return 0
+
+# Stiffness and Load
+def K_loc(k, mesh):
+    K_loc = np.zeros((2,2))
+    K_loc[0,0] = integrate.quad(lambda x: phi_prime(k-1, mesh, x)**2, mesh[k-1], mesh[k])[0]
+    K_loc[1,0] = integrate.quad(lambda x: phi_prime(k-1, mesh, x)*phi_prime(k, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[0,1] = integrate.quad(lambda x: phi_prime(k, mesh, x)*phi_prime(k-1, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[1,1] = integrate.quad(lambda x: phi_prime(k, mesh, x)**2, mesh[k-1], mesh[k])[0]
+
+    K_loc[0,0] += p*integrate.quad(lambda x: phi_prime(k-1, mesh, x) * phi(k-1, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[1,0] += p*integrate.quad(lambda x: phi_prime(k-1, mesh, x) * phi(k, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[0,1] += p*integrate.quad(lambda x: phi_prime(k, mesh, x) * phi(k-1, mesh, x), mesh[k-1], mesh[k])[0]
+    K_loc[1,1] += p*integrate.quad(lambda x: phi_prime(k, mesh, x) * phi(k, mesh, x), mesh[k-1], mesh[k])[0]
+    return K_loc
+def F_loc(k, mesh):
+    F_loc = np.zeros(2)
+    F_loc[0] = integrate.quad(lambda x: f(x) * phi(k-1, mesh, x), mesh[k-1], mesh[k])[0]
+    F_loc[1] = integrate.quad(lambda x: f(x) * phi(k, mesh, x), mesh[k-1], mesh[k])[0]
+    return F_loc
+
+# Assembling
+def Assemble(mesh):
+    m = len(mesh)
+    n = m-1
+
+    K = np.zeros((m,m))
+    F = np.zeros(m)
+    for k in range(1,m):
+        K[k-1:k+1,k-1:k+1] += K_loc(k, mesh)
+        F[k-1:k+1] += F_loc(k, mesh)
+
+    # Boundary conditions
+    # Dirichlet: u(0) = 0
+    K[0,:] = 0
+    K[0,0] = 1
+    F[0]   = 0   
+    # Dirichlet: u(1) = 1
+    K[-1,:] = 0
+    K[-1,-1] = 1
+    F[-1]   = 1
+    return K,F
+
+def plotting(mesh, u, comment):
+    exact_x = np.linspace(0,1,1000)
+    exact = (np.exp(p*exact_x)-1)/(np.exp(p)-1)
+    plt.plot(exact_x, exact, "--", linewidth=1, color="red", label="exact")
+    plt.title(f"p = {p} | n = {n} | {comment}")
+    plt.xlabel("x")
+    plt.ylabel("u(x)")
+    plt.plot(mesh, u, "-o", label="u_h")
+    plt.xticks(mesh, labels=[])
+    plt.legend()
+    plt.grid(True)
+    plt.tight_layout()
+    plt.savefig("task_c.png", dpi=300)
+    plt.show()
+    return 0
+
+##############################################################################
+
+# parameters
+p_list = [-70, 70]
+for p in p_list:
+    mesh = np.linspace(0, 1, 30)
+    m = len(mesh)
+    n = m-1
+
+    K, F = Assemble(mesh)                                 # assemble
+    u = np.linalg.solve(K, F)                             # solve
+    jumps = flux_jumps(mesh, u)                           # flux jumps
+    plotting(mesh, u, "before adapting")                  # plotting
+
+    # NOTE: mesh very different every iteration for adapt_r()
+    iterations = 21
+    for it in range(iterations):
+        # mesh = adapt_h(mesh, jumps)                       # h-adaptivity
+        mesh = adapt_r(mesh, np.abs(jumps))               # r-adaptivity (positive density (jumps)!)
+        m = len(mesh)
+        n = m-1
+
+        K, F = Assemble(mesh)                             # assemble
+        u = np.linalg.solve(K, F)                         # solve
+        jumps = flux_jumps(mesh, u)                       # flux jumps
+        # plotting(mesh, u, f"iteration {it+1}")            # plotting each iteration
+
+    # print(jumps)
+    plotting(mesh, u, f"after {iterations} iterations")   # plotting
+
+