Improved performace

hdpg1d/adaptation.py

@@ -1,10 +1,14 @@
 import numpy as np
+from numpy import concatenate as cat
 import matplotlib.pyplot as plt
+import warnings
 from matplotlib import rc
 from .preprocess import shape, discretization
 
 plt.rc('text', usetex=True)
 plt.rc('font', family='serif')
+# supress the deprecation warning
+warnings.filterwarnings("ignore", ".*GUI is implemented.*")
 
 
 class hdpg1d(object):
@@ -22,8 +26,8 @@ class hdpg1d(object):
         self.kappa = coeff.diffusion
         self.coeff = coeff
         self.u = []
-        self.estError = []
-        self.trueError = []
+        self.estErrorList = []
+        self.trueErrorList = []
 
     def bc(self, case, t=None):
         # boundary condition
@@ -62,39 +66,43 @@ class hdpg1d(object):
             uNode[j - 1] = u[(j - 1) * self.numBasisFuncs]
         uNode[-1] = u[-1]
         plt.figure(1)
-        if counter in [0, 4, 9, 19]:
-            plt.plot(xSmooth, uSmooth, '-', color='C3')
-            plt.plot(self.mesh, uNode, 'C3.')
-            plt.xlabel('$x$', fontsize=17)
-            plt.ylabel('$u$', fontsize=17)
-            plt.axis([-0.05, 1.05, 0, 1.3])
-            plt.grid()
-            # plt.savefig('u_test_{}.pdf'.format(i+1))
-            # plt.show()
-            plt.draw()
-            plt.pause(1e-1)
-            plt.clf()
+        plt.plot(xSmooth, uSmooth, '-', color='C3')
+        plt.plot(self.mesh, uNode, 'C3.')
+        plt.xlabel('$x$', fontsize=17)
+        plt.ylabel('$u$', fontsize=17)
+        plt.axis([-0.05, 1.05, 0, 1.3])
+        plt.grid()
+        plt.pause(1e-1)
+        plt.clf()
 
     def meshAdapt(self, index):
         """Given the index list, adapt the mesh"""
-        in_value = np.zeros(len(index))
+        inValue = np.zeros(len(index))
         for i in np.arange(len(index)):
-            in_value[i] = (self.mesh[index[i]] +
-                           self.mesh[index[i] - 1]) / 2
-
-        self.mesh = np.sort(np.insert(self.mesh, 0, in_value))
+            inValue[i] = (self.mesh[index[i]] +
+                          self.mesh[index[i] - 1]) / 2
+        self.mesh = np.sort(np.insert(self.mesh, 0, inValue))
 
     def solveLocal(self):
         """Solve the primal problem"""
-        A, B, C, D, E, F, G, H, L, R, m = discretization(self.coeff, self.mesh)
+        if 'matLocal' in locals():
+            # if matLocal exists,
+            # only change the mesh instead of initialize again
+            matLocal.mesh = self.mesh
+        else:
+            matLocal = discretization(self.coeff, self.mesh)
+        matLocal.matGen()
         # solve
-        K = -np.concatenate((C.T, G), axis=1).dot(np.linalg.inv(
-            np.bmat([[A, -B], [B.T, D]])).dot(np.concatenate((C, E)))) + H
-        F_hat = np.array([L]).T - np.concatenate((C.T, G), axis=1).dot(np.linalg.inv(
-            np.bmat([[A, -B], [B.T, D]]))).dot(np.array([np.concatenate((R, F))]).T)
+        K = -cat((matLocal.C.T, matLocal.G), axis=1)\
+            .dot(np.linalg.inv(np.bmat([[matLocal.A, -matLocal.B], [matLocal.B.T, matLocal.D]]))
+                 .dot(cat((matLocal.C, matLocal.E)))) + matLocal.H
+        F_hat = np.array([matLocal.L]).T - cat((matLocal.C.T, matLocal.G), axis=1)\
+            .dot(np.linalg.inv(np.bmat([[matLocal.A, -matLocal.B], [matLocal.B.T, matLocal.D]])))\
+            .dot(np.array([cat((matLocal.R, matLocal.F))]).T)
         uFace = np.linalg.solve(K, F_hat)
-        u = np.linalg.inv(np.bmat([[A, -B], [B.T, D]])).dot(
-            np.array([np.concatenate((R, F))]).T - np.concatenate((C, E)).dot(uFace))
+        u = np.linalg.inv(np.bmat([[matLocal.A, -matLocal.B], [matLocal.B.T, matLocal.D]]))\
+            .dot(np.array([np.concatenate((matLocal.R, matLocal.F))]).T -
+                 cat((matLocal.C, matLocal.E)).dot(uFace))
         # self.u = u.A1
         return u.A1, uFace.A1
 
@@ -102,18 +110,24 @@ class hdpg1d(object):
         """Solve the adjoint problem"""
         # solve in the enriched space
         self.coeff.pOrder += 1
-        A, B, C, D, E, F, G, H, L, R, m = discretization(
-            self.coeff, self.mesh)
+        if 'matAdjoint' in locals():
+            matAdjoint.mesh = self.mesh
+        else:
+            matAdjoint = discretization(self.coeff, self.mesh)
+        matAdjoint.matGen()
         self.coeff.pOrder = self.coeff.pOrder - 1
         # add adjoint LHS conditions
-        F = np.zeros(len(F))
-        R[-1] = -self.bc(1)[1]
+        matAdjoint.F = np.zeros(len(matAdjoint.F))
+        matAdjoint.R[-1] = -self.bc(1)[1]
 
         # assemble global matrix LHS
-        LHS = np.bmat([[A, -B, C], [B.T, D, E], [C.T, G, H]])
+        LHS = np.bmat([[matAdjoint.A, -matAdjoint.B, matAdjoint.C],
+                       [matAdjoint.B.T, matAdjoint.D, matAdjoint.E],
+                       [matAdjoint.C.T, matAdjoint.G, matAdjoint.H]])
         # solve
-        U = np.linalg.solve(LHS.T, np.concatenate((R, F, L)))
-        return U[0:2 * len(C)], U[len(C):len(U)]
+        U = np.linalg.solve(LHS.T, np.concatenate(
+            (matAdjoint.R, matAdjoint.F, matAdjoint.L)))
+        return U[0:2 * len(matAdjoint.C)], U[len(matAdjoint.C):len(U)]
 
     def residual(self, U, hat_U, z, hat_z):
         numEle = self.numEle
@@ -273,29 +287,31 @@ class hdpg1d(object):
 
     def adaptive(self):
         tol = 1e-12
-        est_error = 10
+        estError = 10
         counter = 0
         ceilCounter = 100
-        trueError = [[], []]
-        estError = [[], []]
-        while est_error > tol or counter > ceilCounter:
+        trueErrorList = [[], []]
+        estErrorList = [[], []]
+        while estError > tol or counter > ceilCounter:
             # solve
             u, uFace = self.solveLocal()
             adjoint, adjointFace = self.solveAdjoint()
             self.u = u
-            self.plotU(counter)
-            trueError[0].append(self.numEle)
-            trueError[1].append(np.abs(
+            # plot the solution at certain counter
+            if counter in [0, 4, 9, 19]:
+                self.plotU(counter)
+            # record error
+            trueErrorList[0].append(self.numEle)
+            trueErrorList[1].append(np.abs(
                 u[self.numEle * self.numBasisFuncs - 1] - np.sqrt(self.kappa)))
-            est_error, index = self.residual(
-                u, uFace, adjoint, adjointFace)
-            estError[0].append(self.numEle)
-            estError[1].append(est_error)
+            estError, index = self.residual(u, uFace, adjoint, adjointFace)
+            estErrorList[0].append(self.numEle)
+            estErrorList[1].append(estError)
             # adapt
             index = index.tolist()
             self.meshAdapt(index)
             self.numEle = self.numEle + len(index)
             counter += 1
         # save error
-        self.trueError = trueError
-        self.estError = estError
+        self.trueErrorList = trueErrorList
+        self.estErrorList = estErrorList

hdpg1d/postprocess.py

@@ -41,14 +41,14 @@ class utils(object):
     def convHistory(self):
         """Plot the uniform and adaptive convergence history"""
         plt.figure(2)
-        trueError = self.solution.trueError
-        estError = self.solution.estError
-        plt.loglog(trueError[0],
-                   trueError[1], '-ro')
+        trueErrorList = self.solution.trueErrorList
+        estErrorList = self.solution.estErrorList
+        plt.loglog(trueErrorList[0],
+                   trueErrorList[1], '-ro')
         numEle, errorL2 = self.uniConv()
         plt.loglog(numEle, errorL2, '-o')
-        plt.loglog(estError[0],
-                   estError[1], '--', color='#1f77b4')
+        plt.loglog(estErrorList[0],
+                   estErrorList[1], '--', color='#1f77b4')
         plt.xlabel('Number of elements', fontsize=17)
         plt.ylabel('Error', fontsize=17)
         plt.grid()

hdpg1d/preprocess.py

@@ -14,8 +14,6 @@ def shape(x, p):
 
 
 def forcing(x):
-    # f = np.cos(2*np.pi*x)
-    # f = 4*pi**2*sin(2*pi*x)
     f = 1
     return f
 
@@ -33,120 +31,144 @@ def bc(case, t=None):
     return bc
 
 
-def discretization(coeff, mesh):
-    """Given the problem statement, construct the discretization matrice"""
-    # p is the number of basis functions
-    p = coeff.pOrder + 1
-    tau_pos = coeff.tauPlus
-    tau_neg = coeff.tauMinus
-
-    # order of gauss quadrature
-    gorder = 2 * p
-    # shape function and gauss quadrature
-    xi, wi = np.polynomial.legendre.leggauss(gorder)
-    shp, shpx = shape(xi, p)
-
-    con = coeff.convection
-    kappa = coeff.diffusion
-
-    n_ele = len(mesh) - 1
-    dist = np.zeros(n_ele)
-    # elemental forcing vector
-    F = np.zeros(p * n_ele)
-    for i in range(1, n_ele + 1):
-        dist[i - 1] = mesh[i] - mesh[i - 1]
-        f = dist[i - 1] / 2 * shp.dot(
-            wi * forcing(mesh[i - 1] + 1 / 2 * (1 + xi) * dist[i - 1]))
-        F[(i - 1) * p:(i - 1) * p + p] = f
-    F[0] += (con + tau_pos) * bc(0)[0]
-    F[-1] += (-con + tau_neg) * bc(0)[1]
-
-    d = shp.dot(np.diag(wi).dot(shp.T))
-
-    # assemble global D
-    d_face = np.zeros((p, p))
-    d_face[0, 0] = tau_pos
-    d_face[-1, -1] = tau_neg
-    D = np.repeat(dist, p) / 2 * block_diag(*
-                                            [d] * (n_ele)) + block_diag(*[d_face] * (n_ele))
-
-    a = 1 / kappa * shp.dot(np.diag(wi).dot(shp.T))
-    A = np.repeat(dist, p) / 2 * block_diag(*[a] * (n_ele))
-    b = (shpx.T * np.ones((gorder, p))).T.dot(np.diag(wi).dot(shp.T))
-    B = block_diag(*[b] * (n_ele))
-
-    # elemental h
-    h = np.zeros((2, 2))
-    h[0, 0], h[-1, -1] = -con - tau_pos, con - tau_neg
-    # mappinng matrix
-    map_h = np.zeros((2, n_ele), dtype=int)
-    map_h[:, 0] = np.arange(2)
-    for i in np.arange(1, n_ele):
-        map_h[:, i] = np.arange(
-            map_h[2 - 1, i - 1], map_h[2 - 1, i - 1] + 2)
-    # assemble H and eliminate boundaries
-    H = np.zeros((n_ele + 1, n_ele + 1))
-    for i in range(n_ele):
-        for j in range(2):
-            m = map_h[j, i]
-            for k in range(2):
-                n = map_h[k, i]
-                H[m, n] += h[j, k]
-    H = H[1:n_ele][:, 1:n_ele]
-
-    # elemental g
-    g = np.zeros((2, p))
-    g[0, 0], g[-1, -1] = tau_pos, tau_neg
-    # mapping matrix
-    map_g_x = map_h
-    map_g_y = np.arange(p * n_ele, dtype=int).reshape(n_ele, p).T
-    # assemble global G
-    G = np.zeros((n_ele + 1, p * n_ele))
-    for i in range(n_ele):
-        for j in range(2):
-            m = map_g_x[j, i]
-            for k in range(p):
-                n = map_g_y[k, i]
-                G[m, n] += g[j, k]
-    G = G[1:n_ele, :]
-
-    # elemental e
-    e = np.zeros((p, 2))
-    e[0, 0], e[-1, -1] = -con - tau_pos, con - tau_neg
-    # mapping matrix
-    map_e_x = np.arange(p * n_ele, dtype=int).reshape(n_ele, p).T
-    map_e_y = map_h
-    # assemble global E
-    E = np.zeros((p * n_ele, n_ele + 1))
-    for i in range(n_ele):
-        for j in range(p):
-            m = map_e_x[j, i]
-            for k in range(2):
-                n = map_e_y[k, i]
-                E[m, n] += e[j, k]
-    E = E[:, 1:n_ele]
-
-    # L, easy in 1d
-    L = np.zeros(n_ele - 1)
-
-    # elemental c
-    c = np.zeros((p, 2))
-    c[0, 0], c[-1, -1] = -1, 1
-
-    # assemble global C
-    C = np.zeros((p * n_ele, n_ele + 1))
-    for i in range(n_ele):
-        for j in range(p):
-            m = map_e_x[j, i]
-            for k in range(2):
-                n = map_e_y[k, i]
-                C[m, n] += c[j, k]
-    C = C[:, 1:n_ele]
-    # L, easy in 1d
-    L = np.zeros(n_ele - 1)
-
-    # R, easy in 1d
-    R = np.zeros(p * n_ele)
-    R[0] = bc(0)[0]
-    R[-1] = -bc(0)[1]
-    return A, B, C, D, E, F, G, H, L, R, m
+class discretization(object):
+    """Given the problem statement, construct the discretization matrices"""
+
+    def __init__(self, coeff, mesh, enrich=None):
+        self.mesh = mesh
+        self.coeff = coeff
+        self.enrich = enrich
+        # p is the number of basis functions
+        self.numBasisFuncs = coeff.pOrder + 1
+        self.tau_pos = coeff.tauPlus
+        self.tau_neg = coeff.tauMinus
+        self.conv = coeff.convection
+        self.kappa = coeff.diffusion
+        self.n_ele = len(mesh) - 1
+        self.dist = self.distGen()
+        # shape function and gauss quadrature
+        self.gqOrder = 5 * self.numBasisFuncs
+        self.xi, self.wi = np.polynomial.legendre.leggauss(self.gqOrder)
+        self.shp, self.shpx = shape(self.xi, self.numBasisFuncs)
+        # enrich the space if the enrich argument is given
+        if enrich is not None:
+            self.shpEnrich, self.shpxEnrich = shape(
+                self.xi, self.numBasisFuncs + enrich)
+
+    def distGen(self):
+        dist = np.zeros(self.n_ele)
+        for i in range(1, self.n_ele + 1):
+            dist[i - 1] = self.mesh[i] - self.mesh[i - 1]
+        return dist
+
+    def matGen(self):
+        self.lhsGen()
+        self.eleGen()
+        self.interfaceGen()
+
+    def lhsGen(self):
+        """Generate matrices associated with left hand side"""
+        # elemental forcing vector F
+        F = np.zeros(self.numBasisFuncs * self.n_ele)
+        for i in range(1, self.n_ele + 1):
+            f = self.dist[i - 1] / 2 * self.shp.dot(
+                self.wi * forcing(self.mesh[i - 1] + 1 / 2 * (1 + self.xi) * self.dist[i - 1]))
+            F[(i - 1) * self.numBasisFuncs:i * self.numBasisFuncs] = f
+        F[0] += (self.conv + self.tau_pos) * bc(0)[0]
+        F[-1] += (-self.conv + self.tau_neg) * bc(0)[1]
+        # L, easy in 1d
+        L = np.zeros(self.n_ele - 1)
+        # R, easy in 1d
+        R = np.zeros(self.numBasisFuncs * self.n_ele)
+        R[0] = bc(0)[0]
+        R[-1] = -bc(0)[1]
+        self.F, self.L, self.R = F, L, R
+
+    def eleGen(self):
+        """Generate matrices associated with elements"""
+        a = 1 / self.kappa * self.shp.dot(np.diag(self.wi).dot(self.shp.T))
+        A = np.repeat(self.dist, self.numBasisFuncs) / \
+            2 * block_diag(*[a] * (self.n_ele))
+        b = (self.shpx.T * np.ones((self.gqOrder, self.numBasisFuncs))
+             ).T.dot(np.diag(self.wi).dot(self.shp.T))
+        B = block_diag(*[b] * (self.n_ele))
+        d = self.shp.dot(np.diag(self.wi).dot(self.shp.T))
+        # assemble global D
+        d_face = np.zeros((self.numBasisFuncs, self.numBasisFuncs))
+        d_face[0, 0] = self.tau_pos
+        d_face[-1, -1] = self.tau_neg
+        D = np.repeat(self.dist, self.numBasisFuncs) / 2 * \
+            block_diag(*[d] * (self.n_ele)) + \
+            block_diag(*[d_face] * (self.n_ele))
+        self.A, self.B, self.D = A, B, D
+
+    def interfaceGen(self):
+        """Generate matrices associated with interfaces"""
+        tau_pos, tau_neg, conv = self.tau_pos, self.tau_neg, self.conv
+        # elemental h
+        h = np.zeros((2, 2))
+        h[0, 0], h[-1, -1] = -conv - tau_pos, conv - tau_neg
+        # mappinng matrix
+        map_h = np.zeros((2, self.n_ele), dtype=int)
+        map_h[:, 0] = np.arange(2)
+        for i in np.arange(1, self.n_ele):
+            map_h[:, i] = np.arange(
+                map_h[2 - 1, i - 1], map_h[2 - 1, i - 1] + 2)
+        # assemble H and eliminate boundaries
+        H = np.zeros((self.n_ele + 1, self.n_ele + 1))
+        for i in range(self.n_ele):
+            for j in range(2):
+                m = map_h[j, i]
+                for k in range(2):
+                    n = map_h[k, i]
+                    H[m, n] += h[j, k]
+        H = H[1:self.n_ele][:, 1:self.n_ele]
+
+        # elemental e
+        e = np.zeros((self.numBasisFuncs, 2))
+        e[0, 0], e[-1, -1] = -conv - tau_pos, conv - tau_neg
+        # mapping matrix
+        map_e_x = np.arange(self.numBasisFuncs * self.n_ele,
+                            dtype=int).reshape(self.n_ele, self.numBasisFuncs).T
+        map_e_y = map_h
+        # assemble global E
+        E = np.zeros((self.numBasisFuncs * self.n_ele, self.n_ele + 1))
+        for i in range(self.n_ele):
+            for j in range(self.numBasisFuncs):
+                m = map_e_x[j, i]
+                for k in range(2):
+                    n = map_e_y[k, i]
+                    E[m, n] += e[j, k]
+        E = E[:, 1:self.n_ele]
+
+        # elemental c
+        c = np.zeros((self.numBasisFuncs, 2))
+        c[0, 0], c[-1, -1] = -1, 1
+        # assemble global C
+        C = np.zeros((self.numBasisFuncs * self.n_ele, self.n_ele + 1))
+        for i in range(self.n_ele):
+            for j in range(self.numBasisFuncs):
+                m = map_e_x[j, i]
+                for k in range(2):
+                    n = map_e_y[k, i]
+                    C[m, n] += c[j, k]
+        C = C[:, 1:self.n_ele]
+
+        # elemental g
+        g = np.zeros((2, self.numBasisFuncs))
+        g[0, 0], g[-1, -1] = tau_pos, tau_neg
+        # mapping matrix
+        map_g_x = map_h
+        map_g_y = np.arange(self.numBasisFuncs * self.n_ele,
+                            dtype=int).reshape(self.n_ele, self.numBasisFuncs).T
+        # assemble global G
+        G = np.zeros((self.n_ele + 1, self.numBasisFuncs * self.n_ele))
+        for i in range(self.n_ele):
+            for j in range(2):
+                m = map_g_x[j, i]
+                for k in range(self.numBasisFuncs):
+                    n = map_g_y[k, i]
+                    G[m, n] += g[j, k]
+        G = G[1:self.n_ele, :]
+
+        self.C, self.E, self.G, self.H = C, E, G, H

hdpg1d/solve.py

@@ -34,7 +34,7 @@ def getCoefficients():
     if (isDefault):
         Coeff = coefficients(1e-6, 0, 0, 2, 2, 1e-6, 1e-6)
     else:
-        Coeff = coefficients.from_input()
+        Coeff = coefficients.fromInput()
     return Coeff
 
 
@@ -50,9 +50,11 @@ def menu():
 
 def hdgSolve():
     hdgCoeff = getCoefficients()
+    print("Solving...")
     hdgSolution = hdpg1d(hdgCoeff)
     # solve the problem adaptively and plot convergence history
     hdgSolution.adaptive()
+    print("Problem solved. Please check the convergence plot.")
     utils(hdgSolution).convHistory()