Merge branch 'release/v2.0'

This is a major update that includes multiple bug fixes and performance improvement.
This version also introduces the capability to solve convection only problems.

.gitignore

@@ -1,7 +1,7 @@
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]
-
+MANIFEST
 # C extensions
 *.so
 
@@ -84,3 +84,6 @@ venv
 # Org-mode
 .org-id-locations
 *_archive
+
+# markdown preview
+.markdown*

MANIFEST

@@ -1,7 +0,0 @@
-# file GENERATED by distutils, do NOT edit
-setup.py
-hdpg1d/__init__.py
-hdpg1d/coefficients.py
-hdpg1d/discretization.py
-hdpg1d/postprocess.py
-hdpg1d/solve.py

README.md

@@ -11,12 +11,23 @@
 
 # HDPG1D<a id="sec-1" name="sec-1"></a>
 
-This is a python package that solves 1-Dimensional PDEs using hybridizable discontinuous Petrov-Galerkin discretization, a novel FEA discretization that ensures the best solution quality without extra stablization mechanism. The following image shows the solution of 1D inviscid Burger's equation using HDPG discretization. Notice that the oscillations near the shock are well controlled even without any artificial stablization mechanism.
-<a href="http://imgur.com/HrWIi4s"><img src="http://i.imgur.com/HrWIi4s.png" width="50%" height="50%" title="source: imgur.com" /></a>
+This is a python package that solves 1-Dimensional PDEs using hybridizable discontinuous Petrov-Galerkin discretization, a novel FEA discretization that ensures the best solution quality without extra stablization mechanism. The following image shows the solution of 1D inviscid Burger's equation using HDPG discretization. 
+```math
+\frac{\partial u}{\partial t} + \frac{1}{2} \frac{\partial u^2}{\partial x} = 0, \quad \text{in } \Omega \in [0,1].
+```
+Notice that the oscillations near the shock are well controlled even without any artificial stablization mechanism.
+
+<p align="center">
+<img align="centre" img src="http://i.imgur.com/HrWIi4s.png" width="50%" height="50%" title="source: imgur.com" />
+</p>
+
+The main task of this project is to build an automated CFD framework based on HPDG discretization and robust adaptive mesh refinement. The solver is currently capable of solving 1D convection-diffusion equations and is being actively developed.
 
 # Install<a id="sec-2" name="sec-2"></a>
-In terminal: 
+In the souce directory: 
 ```bash
+python setup.py sdist
+cd dist/
 pip install hdpg1d
 ```
 
@@ -25,3 +36,4 @@ In terminal, call:
 ```bash
 PGsolve
 ```
+Follow the prompts to setup the problem, visualize the solution, and check the convergence.

hdpg1d/adaptation.py

@@ -0,0 +1,184 @@
+import numpy as np
+from numpy import concatenate as cat
+from copy import copy
+import matplotlib.pyplot as plt
+import warnings
+from .preprocess import shape, discretization, boundaryCondition
+
+plt.rc('text', usetex=True)
+plt.rc('font', family='serif')
+# supress the deprecation warning
+warnings.filterwarnings("ignore", ".*GUI is implemented.*")
+
+
+class hdpg1d(object):
+    """
+    1D HDG solver
+    """
+
+    def __init__(self, coeff):
+        self.numEle = coeff.numEle
+        self.numBasisFuncs = coeff.pOrder + 1
+        self.coeff = coeff
+        self.mesh = np.linspace(0, 1, self.numEle + 1)
+        self.enrichOrder = 1
+        self.primalSoln = None
+        self.adjointSoln = None
+        self.estErrorList = [[], []]
+        self.trueErrorList = [[], []]
+
+    def separateSoln(self, soln):
+        """Separate gradState (q and u), stateFace from the given soln"""
+        gradState, stateFace = np.split(
+            soln, [len(soln) - self.numEle + 1])
+        return gradState, stateFace
+
+    def plotState(self, counter):
+        """Plot solution u with smooth higher oredr quadrature"""
+        stateSmooth = np.array([])
+        stateNode = np.zeros(self.numEle + 1)
+        xSmooth = np.array([])
+        gradState, _ = self.separateSoln(self.primalSoln)
+        halfLenState = int(len(gradState) / 2)
+        state = gradState[halfLenState:2 * halfLenState]
+        # quadrature rule
+        gorder = 10 * self.numBasisFuncs
+        xi, wi = np.polynomial.legendre.leggauss(gorder)
+        shp, shpx = shape(xi, self.numBasisFuncs)
+        for j in range(1, self.numEle + 1):
+            xSmooth = np.hstack((xSmooth, (self.mesh[(j - 1)] + self.mesh[j]) / 2 + (
+                self.mesh[j] - self.mesh[j - 1]) / 2 * xi))
+            stateSmooth = np.hstack(
+                (stateSmooth, shp.T.dot(state[(j - 1) * self.numBasisFuncs:j * self.numBasisFuncs])))
+            stateNode[j - 1] = state[(j - 1) * self.numBasisFuncs]
+        stateNode[-1] = state[-1]
+        plt.figure(1)
+        plt.plot(xSmooth, stateSmooth, '-', color='C3')
+        plt.plot(self.mesh, stateNode, 'C3.')
+        plt.xlabel('$x$', fontsize=17)
+        plt.ylabel('$u$', fontsize=17)
+        # plt.axis([-0.05, 1.05, 0, 1.3])
+        plt.grid()
+        plt.pause(5e-1)
+
+    def meshAdapt(self, index):
+        """Given the index list, adapt the mesh"""
+        inValue = np.zeros(len(index))
+        for i in np.arange(len(index)):
+            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"""
+        if 'matLocal' in locals():
+            # if matLocal exists,
+            # only change the mesh instead of initializing again
+            matLocal.mesh = self.mesh
+        else:
+            matLocal = discretization(self.coeff, self.mesh)
+        matGroup = matLocal.matGroup()
+        A, B, _, C, D, E, F, G, H, L, R = matGroup
+        # solve by exploiting the local global separation
+        K = -cat((C.T, G), axis=1)\
+            .dot(np.linalg.inv(np.bmat([[A, -B], [B.T, D]]))
+                 .dot(cat((C, E)))) + H
+        F_hat = np.array([L]).T - cat((C.T, G), axis=1)\
+            .dot(np.linalg.inv(np.bmat([[A, -B], [B.T, D]])))\
+            .dot(np.array([cat((R, F))]).T)
+        stateFace = np.linalg.solve(K, F_hat)
+        gradState = np.linalg.inv(np.bmat([[A, -B], [B.T, D]]))\
+            .dot(np.array([np.concatenate((R, F))]).T -
+                 cat((C, E)).dot(stateFace))
+        self.primalSoln = cat((gradState.A1, stateFace.A1))
+
+    def solveAdjoint(self):
+        """Solve the adjoint problem"""
+        # solve in the enriched space
+        _coeff = copy(self.coeff)
+        _coeff.pOrder = _coeff.pOrder + 1
+        if 'matAdjoint' in locals():
+            matAdjoint.mesh = self.mesh
+        else:
+            matAdjoint = discretization(_coeff, self.mesh)
+        matGroup = matAdjoint.matGroup()
+        A, B, _, C, D, E, F, G, H, L, R = matGroup
+        # add adjoint LHS conditions
+        F = np.zeros(len(F))
+        R[-1] = -boundaryCondition(1)[1]
+        # assemble global matrix LHS
+        LHS = np.bmat([[A, -B, C],
+                       [B.T, D, E],
+                       [C.T, G, H]])
+        # solve in one shoot
+        soln = np.linalg.solve(LHS.T, cat((R, F, L)))
+        self.adjointSoln = soln
+
+    def DWResidual(self):
+        if 'matResidual' in locals():
+            matResidual.mesh = self.mesh
+        else:
+            matResidual = discretization(
+                self.coeff, self.mesh, self.enrichOrder)
+        matGroup = matResidual.matGroup()
+        A, B, BonQ, C, D, E, F, G, H, L, R = matGroup
+        LHS = np.bmat([[A, -B, C],
+                       [BonQ, D, E]])
+        RHS = cat((R, F))
+        residual = np.zeros(self.numEle)
+        numEnrich = self.numBasisFuncs + self.enrichOrder
+        adjointGradState, adjointStateFace = self.separateSoln(
+            self.adjointSoln)
+        for i in np.arange(self.numEle):
+            primalResidual = (LHS.dot(self.primalSoln) - RHS).A1
+            uLength = self.numEle * numEnrich
+            stepLength = i * numEnrich
+            uDWR = primalResidual[stepLength:stepLength + numEnrich].dot(
+                (1 - adjointGradState)[stepLength:stepLength + numEnrich])
+            qDWR = primalResidual[uLength + stepLength:uLength +
+                                  stepLength + numEnrich]\
+                .dot((1 - adjointGradState)[uLength + stepLength:uLength +
+                                            stepLength + numEnrich])
+            residual[i] = uDWR + qDWR
+        # sort residual index
+        residualIndex = np.argsort(np.abs(residual))
+        # select top \theta% elements with the largest error
+        theta = 0.15
+        refineIndex = residualIndex[
+            int(self.numEle * (1 - theta)):len(residual)] + 1
+        return np.abs(np.sum(residual)), refineIndex
+
+    def adaptive(self):
+        TOL = 1e-10
+        estError = 10
+        nodeCount = 0
+        maxCount = 30
+        while estError > TOL and nodeCount < maxCount:
+            # solve
+            self.solveLocal()
+            self.solveAdjoint()
+            # plot the solution at certain counter
+            if nodeCount in [0, 4, 9, 19, maxCount]:
+                plt.clf()
+                self.plotState(nodeCount)
+            # record error
+            self.trueErrorList[0].append(self.numEle)
+            self.trueErrorList[1].append(
+                self.primalSoln[self.numEle * self.numBasisFuncs - 1])
+            estError, index = self.DWResidual()
+            self.estErrorList[0].append(self.numEle)
+            self.estErrorList[1].append(estError)
+            # adapt
+            index = index.tolist()
+            self.meshAdapt(index)
+            self.numEle = self.numEle + len(index)
+            nodeCount += 1
+            print("Iteration {}. Target function error {:.3e}.".format(
+                nodeCount, estError))
+            if nodeCount == maxCount:
+                print("Max iteration number is reached "
+                      "while the convergence criterion is not satisfied.\n"
+                      "Check the problem statement or "
+                      "raise the max iteration number, then try again.\n")
+                from .solve import runInteractive
+                runInteractive()

hdpg1d/cmd.py

@@ -1,8 +1,8 @@
 """
 command line interface
 """
-from .solve import run
+from .solve import runInteractive
 
 
 def main():
-    run()
+    runInteractive()

hdpg1d/coefficients.py

@@ -1,24 +1,33 @@
 class coefficients:
-    def __init__(self, diff, conv, flux, porder, nele):
-        self.diffusion = diff
-        self.covection = conv
-        self.flux = flux
-        self.porder = porder
-        self.nele = nele
+    def __init__(self, diff, conv, reaction, pOrder, numEle, tauPlus, tauMinus):
+        if diff == 0:
+            # set the diffusion constant to a small number
+            # to avoid division by zero error
+            diff = 1e-16
+        self.DIFFUSION = diff
+        self.CONVECTION = conv
+        self.REACTION = reaction
+        self.pOrder = pOrder
+        self.numEle = numEle
+        self.TAUPLUS = tauPlus
+        self.TAUMINUS = tauMinus
 
     @classmethod
-    def from_input(cls):
+    def fromInput(cls):
         while True:
             try:
                 print("Please provide the following coefficients.")
-                diff = float(input("Diffusion coefficient: "))
-                conv = float(input("Covection coefficient: "))
-                flux = float(input("Flux: "))
-                porder = int(input("Order of polynomials: "))
-                nele = int(input("Number of elements: "))
+                diff = float(input("Diffusion constant (float): "))
+                conv = float(input("Covection constant (float): "))
+                reaction = float(input("Reaction constant (float): "))
+                pOrder = int(input("Order of polynomials (int): "))
+                numEle = int(input("Number of elements (int): "))
+                tauPlus = float(input("Stablization parameter plus (float): "))
+                tauMinus = float(
+                    input("Stablization parameter minus (float): "))
             except ValueError:
                 print("Sorry, wrong data type.")
                 continue
             else:
                 break
-        return cls(diff, conv, flux, porder, nele)
+        return cls(diff, conv, reaction, pOrder, numEle, tauPlus, tauMinus)

hdpg1d/discretization.py

@@ -1,621 +0,0 @@
-import numpy as np
-from numpy import sin, cos, pi
-from scipy.linalg import block_diag
-import matplotlib.pyplot as plt
-from matplotlib import rc
-plt.rc('text', usetex=True)
-plt.rc('font', family='serif')
-
-
-class HDPG1d(object):
-    """
-    1d advection hdg solver outlined in 'an implicit HHDG method for
-    confusion'. Test case: /tau = 1, convection only, linear and higher order.
-    Please enter number of elements and polynomial order, i.e., HDG1d(10,2)
-    """
-    # stablization parameter
-    tau_pos = 1e-6
-    tau_neg = 1e-6
-    # advection constant
-    c = 0
-    # diffusion constant
-    kappa = 1e-6
-
-    def __init__(self, n_ele, p):
-        self.n_ele = n_ele
-        # number of the shape functions (thus porder + 1)
-        self.p = p + 1
-
-    def bc(self, case, t=None):
-        # boundary condition
-        if case == 0:
-            # advection-diffusion
-            bc = [0, 0]
-        if case == 1:
-            # simple convection
-            # bc = np.sin(2*np.pi*t)
-            # adjoint boundary
-            bc = [0, 1]
-        return bc
-
-    def shape(self, x, p):
-        """ evaluate shape functions at give locations"""
-        # coeffient matrix
-        A = np.array([np.linspace(-1, 1, p)]).T**np.arange(p)
-        C = np.linalg.inv(A).T
-        x = np.array([x]).T
-        shp = C.dot((x**np.arange(p)).T)
-        shpx = C[:, 1::1].dot((x**np.arange(p - 1) * np.arange(1, p)).T)
-        return shp, shpx
-
-    def forcing(self, x):
-        # f = np.cos(2*np.pi*x)
-        # f = 4*pi**2*sin(2*pi*x)
-        f = 1
-        return f
-
-    def mesh(self, n_ele, index, x):
-        """generate mesh"""
-        # if n_ele < 1 or n_ele > self.n_ele:
-        #   raise RuntimeError('Bad Element number')
-        in_value = np.zeros(len(index))
-        for i in np.arange(len(index)):
-            in_value[i] = (x[index[i]] + x[index[i] - 1]) / 2
-
-        x_c = np.sort(np.insert(x, 0, in_value))
-
-        x_i = np.linspace(x_c[n_ele - 1], x_c[n_ele], num=self.p)
-        dx = x_c[n_ele] - x_c[n_ele - 1]
-        return x_i, dx, x_c
-
-    def uexact(self, x):
-        # Ue = self.bc(0) + np.sin(2*pi*x)
-        Ue = (np.exp(1 / self.kappa * (x - 1)) - 1) / \
-            (np.exp(-1 / self.kappa) - 1)
-        return Ue
-
-    def matrix_gen(self, index, x):
-        n_ele = self.n_ele
-        # # space size
-        # dx = 1/n_ele
-
-        # order of polynomial shape functions
-        p = self.p
-
-        # order of gauss quadrature
-        gorder = 2 * p
-        # shape function and gauss quadrature
-        xi, wi = np.polynomial.legendre.leggauss(gorder)
-        shp, shpx = self.shape(xi, p)
-        # ---------------------------------------------------------------------
-        # advection constant
-        con = self.c
-
-        # diffusion constant
-        kappa = self.kappa
-
-        # number of nodes (solution U)
-        n_ele = self.n_ele + len(index)
-        # elemental forcing vector
-        F = np.zeros(p * n_ele)
-        for i in range(1, n_ele + 1):
-            x_i, dx_i, _ = self.mesh(i, index, x)
-            f = dx_i / 2 * \
-                shp.dot(wi * self.forcing(x[0] + 1 / 2 * (1 + xi) * dx_i))
-            F[(i - 1) * p:(i - 1) * p + p] = f
-            # elemental m (for convection only)
-            # m = dx/2*shp.dot(np.diag(wi).dot(shp.T))
-            # add boundary
-        F[0] += (con + self.tau_pos) * self.bc(0)[0]
-        F[-1] += (-con + self.tau_neg) * self.bc(0)[1]
-
-        # elemental d
-        # d = con*(-shpx.T*np.ones((gorder,p))).T.dot(np.diag(wi).dot(shp.T))
-        d = shp.dot(np.diag(wi).dot(shp.T))
-        # d[0, 0] += self.tau_pos
-        # d[-1, -1] += self.tau_neg
-
-        # elemental a
-        a = 1 / kappa * shp.dot(np.diag(wi).dot(shp.T))
-
-        # elemental b
-        b = (shpx.T * np.ones((gorder, p))).T.dot(np.diag(wi).dot(shp.T))
-
-        # elemental h
-        h = np.zeros((2, 2))
-        h[0, 0], h[-1, -1] = -con - self.tau_pos, con - self.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] = self.tau_pos, self.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 - self.tau_pos, con - self.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] = self.bc(0)[0]
-        R[-1] = -self.bc(0)[1]
-        return d, m, E, G, H, F, L, a, b, C, R
-
-    def solve_local(self, index, x):
-        """ solve the 1d advection equation wit local HDG"""
-        d, _, E, G, H, F, L, a, b, C, R = self.matrix_gen(index, x)
-        # find dx
-        dx = np.zeros(self.n_ele + len(index))
-
-        for i in range(1, self.n_ele + len(index) + 1):
-            x_i, dx_i, x_n = self.mesh(i, index, x)
-            dx[i - 1] = dx_i
-
-        # assemble global D
-        bb = np.zeros((self.p, self.p))
-        bb[0, 0] = self.tau_pos
-        bb[-1, -1] = self.tau_neg
-        D = np.repeat(dx, self.p) / 2 * block_diag(*[d] * (
-            self.n_ele + len(index))) + block_diag(*[bb] * (self.n_ele + len(index)))
-        # assemble global A
-        A = np.repeat(dx, self.p) / 2 * block_diag(*
-                                                   [a] * (self.n_ele + len(index)))
-
-        # assemble global B
-        B = block_diag(*[b] * (self.n_ele + len(index)))
-        # solve U and \lambda
-        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)
-        lamba = 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(lamba))
-        return U, lamba
-
-    def solve_adjoint(self, index, x, u, u_hat):
-        self.p = self.p + 1
-        d, _, E, G, H, F, L, a, b, C, R = self.matrix_gen(index, x)
-
-        # add boundary
-        F = np.zeros(len(F))
-        R[-1] = -self.bc(1)[1]
-
-        # find dx
-        dx = np.zeros(self.n_ele + len(index))
-        for i in range(1, self.n_ele + len(index) + 1):
-            x_i, dx_i, x_n = self.mesh(i, index, x)
-            dx[i - 1] = dx_i
-        # assemble global D
-        bb = np.zeros((self.p, self.p))
-        bb[0, 0] = self.tau_pos
-        bb[-1, -1] = self.tau_neg
-        D = np.repeat(dx, self.p) / 2 * block_diag(*[d] * (
-            self.n_ele + len(index))) + block_diag(*[bb] * (self.n_ele + len(index)))
-
-        # assemble global A
-        A = np.repeat(dx, self.p) / 2 * block_diag(*
-                                                   [a] * (self.n_ele + len(index)))
-
-        # assemble global B
-        B = block_diag(*[b] * (self.n_ele + len(index)))
-
-        # # assemble global matrix LHS
-        LHS = np.bmat([[A, -B, C], [B.T, D, E], [C.T, G, H]])
-
-        # solve U and \lambda
-        U = np.linalg.solve(LHS.T, np.concatenate((R, F, L)))
-
-        return U[0:2 * self.p * (self.n_ele + len(index))], U[2 * self.p * (self.n_ele + len(index)):len(U)]
-
-    def diffusion(self):
-        """solve 1d convection with local HDG"""
-
-        # begin and end time
-        t, T = 0, 1
-
-        # time marching step for diffusion equation
-        dt = 1e-3
-
-        d, m, E, G, H, F, L, a, b, C, R = self.matrix_gen()
-        # add time derivatives to the space derivatives (both are
-        # elmental-wise)
-        d = d + 1 / dt * m
-        # assemble global D
-        D = block_diag(*[d] * self.n_ele)
-        # assemble global A
-        A = block_diag(*[a] * self.n_ele)
-        # assemble global B
-        B = block_diag(*[b] * self.n_ele)
-
-        # initial condition
-        X = np.zeros(self.p * self.n_ele)
-        for i in range(1, self.n_ele + 1):
-            x = self.mesh(i)
-            X[(i - 1) * self.p:(i - 1) * self.p + self.p] = x
-        U = np.concatenate((pi * cos(pi * X), sin(pi * X)))
-
-        # assemble M
-        M = block_diag(*[1 / dt * m] * self.n_ele)
-
-        # time marching
-        while t < T:
-            # add boundary conditions
-            F_dynamic = F + \
-                M.dot(U[self.n_ele * self.p:2 * self.n_ele * self.p])
-
-            # assemble global matrix LHS
-            LHS = np.bmat([[A, -B, C], [B.T, D, E], [C.T, G, H]])
-
-            # solve U and \lambda
-            U = np.linalg.solve(LHS, np.concatenate((R, F_dynamic, L)))
-
-            # plot solutions
-            plt.clf()
-            plt.plot(X, U[self.n_ele * self.p:2 * self.n_ele * self.p], '-r.')
-            plt.plot(X, sin(pi * X) * np.exp(-pi**2 * t))
-            plt.ylim([0, 1])
-            plt.grid()
-            plt.pause(1e-3)
-
-        plt.close()
-        print("Diffusion equation du/dt - du^2/d^2x = 0 with u_exact ="
-              ' 6sin(pi*x)*exp(-pi^2*t).')
-        plt.figure(1)
-        plt.plot(X, U[self.n_ele * self.p:2 * self.n_ele * self.p], '-r.')
-        plt.plot(X, sin(pi * X) * np.exp(-pi**2 * T))
-        plt.xlabel('x')
-        plt.ylabel('y')
-        plt.legend(('Numberical', 'Exact'), loc='upper left')
-        plt.title('Simple Diffusion Equation Solution at t = {}'.format(T))
-        plt.grid()
-        plt.savefig('diffusion', bbox_inches='tight')
-        plt.show(block=False)
-        return U
-
-    def errorL2(self):
-        # evaluate error
-        errorL2 = 0.
-        # solve and track solving time
-        x = np.linspace(0, 1, self.n_ele + 1)
-        U, _ = self.solve_local([], x)
-        errorL2 = np.abs(U[self.p * self.n_ele - 1] - np.sqrt(self.kappa))
-        return errorL2
-
-    def conv(self):
-        p = np.arange(2, 6)
-        n_ele = 2**np.arange(1, 9)
-        legend = []
-        errorL2 = np.zeros((n_ele.size, p.size))
-        for i in range(p.size):
-            self.p = p[i]
-            for j, n in enumerate(n_ele):
-                self.n_ele = n
-                print(self.errorL2())
-                errorL2[j, i] = self.errorL2()
-                legend.append('p = {}'.format(p[i] - 1))
-        # loglog plot
-        plt.figure(3)
-        plt.loglog(n_ele, errorL2, '-o')
-        plt.legend((legend))
-        plt.xlabel('Number of nodes')
-        plt.ylabel('L2 norm(log)')
-        plt.title('Convergence rate(log)')
-        plt.grid()
-        plt.savefig('con_rate_uni.pdf', bbox_inches='tight')
-
-        # output then save the convergence rates
-        rate = np.log(errorL2[0:-1, :] /
-                      errorL2[1:errorL2.size, :]) / np.log(2)
-        for i in range(p.size):
-            print('p = {}, convergence rate = {:.4f}'.format(
-                p[i] - 1, rate[-1, i]))
-        np.savetxt('Con_rate.csv', rate, fmt='%10.4f',
-                   delimiter=',', newline='\n')
-        return n_ele, errorL2
-
-    def residual(self, U, hat_U, z, hat_z, dx, index, x_c):
-        n_ele = self.n_ele
-
-        # order of polynomial shape functions
-        p = self.p
-
-        p_l = p - 1
-        # order of gauss quadrature
-        gorder = 2 * p
-        # shape function and gauss quadrature
-        xi, wi = np.polynomial.legendre.leggauss(gorder)
-        shp, shpx = self.shape(xi, p)
-        shp_l, shpx_l = self.shape(xi, p_l)
-        # ---------------------------------------------------------------------
-        # advection constant
-        con = self.c
-
-        # diffusion constant
-        kappa = self.kappa
-
-        z_q, z_u, z_hat = np.zeros(self.p * self.n_ele), \
-            np.zeros(self.p * self.n_ele), np.zeros(self.n_ele - 1)
-
-        q, u, lamba = np.zeros(p_l * self.n_ele), \
-            np.zeros(p_l * self.n_ele), np.zeros(self.n_ele - 1)
-
-        for i in np.arange(self.p * self.n_ele):
-            z_q[i] = z[i]
-            z_u[i] = z[i + self.p * self.n_ele]
-
-        for i in np.arange(p_l * self.n_ele):
-            q[i] = U[i]
-            u[i] = U[i + p_l * self.n_ele]
-
-        for i in np.arange(self.n_ele - 1):
-            z_hat[i] = hat_z[i]
-
-        # add boundary condtions to U_hat
-        U_hat = np.zeros(self.n_ele + 1)
-        for i, x in enumerate(hat_U):
-            U_hat[i + 1] = x
-        U_hat[0] = self.bc(0)[0]
-        U_hat[-1] = self.bc(0)[1]
-
-        # L, easy in 1d
-        L = np.zeros(n_ele + 1)
-
-        # R, easy in 1d
-        RR = np.zeros(p * n_ele)
-
-        # elemental forcing vector
-        F = np.zeros(p * n_ele)
-        for i in range(1, n_ele + 1):
-            f = dx[i - 1] / 2 * \
-                shp.dot(
-                    wi * self.forcing(x_c[0] + 1 / 2 * (1 + xi) * dx[i - 1]))
-            F[(i - 1) * p:(i - 1) * p + p] = f
-
-        # elemental h
-        h = np.zeros((2, 2))
-        h[0, 0], h[-1, -1] = -con - self.tau_pos, con - self.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_l))
-        g[0, 0], g[-1, -1] = self.tau_pos, self.tau_neg
-        # mapping matrix
-        map_g_x = map_h
-        map_g_y = np.arange(p_l * n_ele, dtype=int).reshape(n_ele, p_l).T
-        # assemble global G
-        G = np.zeros((n_ele + 1, p_l * n_ele))
-        for i in range(n_ele):
-            for j in range(2):
-                m = map_g_x[j, i]
-                for k in range(p_l):
-                    n = map_g_y[k, i]
-                    G[m, n] += g[j, k]
-        G = G[1:n_ele, :]
-
-        # elemental c
-        c = np.zeros((p_l, 2))
-        c[0, 0], c[-1, -1] = -1, 1
-
-        # mapping matrix
-        map_e_x = np.arange(p_l * n_ele, dtype=int).reshape(n_ele, p_l).T
-        map_e_y = map_h
-
-        # assemble global C
-        C = np.zeros((p_l * n_ele, n_ele + 1))
-        for i in range(n_ele):
-            for j in range(p_l):
-                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)
-
-        # residual vector
-        R = np.zeros(self.n_ele)
-        for i in np.arange(self.n_ele):
-            a = dx[i] / 2 * 1 / kappa * \
-                ((shp.T).T).dot(np.diag(wi).dot(shp_l.T))
-
-            b = ((shpx.T) * np.ones((gorder, p))
-                 ).T.dot(np.diag(wi).dot(shp_l.T))
-
-            b_t = ((shpx_l.T) * np.ones((gorder, p_l))
-                   ).T.dot(np.diag(wi).dot(shp.T))
-
-            d = dx[i] / 2 * shp.dot(np.diag(wi).dot(shp_l.T))
-            d[0, 0] += self.tau_pos
-            d[-1, -1] += self.tau_neg
-
-            h = np.zeros((2, 2))
-            h[0, 0], h[-1, -1] = -con - self.tau_pos, con - self.tau_neg
-
-            g = np.zeros((2, p_l))
-            g[0, 0], g[-1, -1] = self.tau_pos, self.tau_neg
-
-            e = np.zeros((p, 2))
-            e[0, 0], e[-1, -1] = -con - self.tau_pos, con - self.tau_neg
-
-            c = np.zeros((p, 2))
-            c[0, 0], c[-1, -1] = -1, 1
-
-            m = np.zeros((2, p_l))
-            m[0, 0], m[-1, -1] = -1, 1
-            # local error
-            R[i] = (np.concatenate((a.dot(q[p_l * i:p_l * i + p_l]) + -b.dot(u[p_l * i:p_l * i + p_l]) + c.dot(U_hat[i:i + 2]),
-                                    b_t.T.dot(q[p_l * i:p_l * i + p_l]) + d.dot(u[p_l * i:p_l * i + p_l]) + e.dot(U_hat[i:i + 2]))) - np.concatenate((RR[p * i:p * i + p], F[p * i:p * i + p]))).dot(1 - np.concatenate((z_q[p * i:p * i + p], z_u[p * i:p * i + p])))
-
-        com_index = np.argsort(np.abs(R))
-        # select \theta% elements with the large error
-        theta = 0.15
-        refine_index = com_index[int(self.n_ele * (1 - theta)):len(R)]
-        self.p = self.p - 1
-
-        # global error
-        R_g = (C.T.dot(q) + G.dot(u) + H.dot(U_hat[1:-1])).dot(1 - z_hat)
-        return np.abs(np.sum(R) + np.sum(R_g)), refine_index + 1
-
-    def adaptive(self):
-        x = np.linspace(0, 1, self.n_ele + 1)
-        index = []
-        U, hat_U = self.solve_local(index, x)
-        U_adjoint, hat_adjoint = self.solve_adjoint(index, x, U, hat_U)
-        X = np.zeros(self.p * self.n_ele)
-        dx = np.zeros(self.n_ele + len(index))
-        for i in range(1, self.n_ele + 1):
-            x_i, dx_i, x_n = self.mesh(i, index, x)
-            X[(i - 1) * self.p:(i - 1) * self.p + self.p] = x_i
-            dx[i - 1] = dx_i
-
-        numAdaptive = 28
-        trueError = np.zeros((2, numAdaptive))
-        estError = np.zeros((2, numAdaptive))
-        for i in np.arange(numAdaptive):
-            est_error, index = self.residual(
-                U, hat_U, U_adjoint, hat_adjoint, dx, index, x)
-            index = index.tolist()
-            U, hat_U = self.solve_local(index, x)
-            U_adjoint, hat_adjoint = self.solve_adjoint(index, x, U, hat_U)
-            self.p = self.p - 1
-            X = np.zeros(self.p * (self.n_ele + len(index)))
-            dx = np.zeros(self.n_ele + len(index))
-            for j in range(1, self.n_ele + len(index) + 1):
-                x_i, dx_i, x_n = self.mesh(j, index, x)
-                X[(j - 1) * self.p:(j - 1) * self.p + self.p] = x_i
-                dx[j - 1] = dx_i
-            x = x_n
-            estError[0, i] = self.n_ele
-            estError[1, i] = est_error
-
-            self.n_ele = self.n_ele + len(index)
-
-            # U_1d = np.zeros(len(U))
-            # for j in np.arange(len(U)):
-            #     U_1d[j] = U[j]
-            # Unum = np.array([])
-            # Xnum = np.array([])
-            # Qnum = np.array([])
-            # for j in range(1, self.n_ele + 1):
-            #     # Gauss quadrature
-            #     gorder = 10 * self.p
-            #     xi, wi = np.polynomial.legendre.leggauss(gorder)
-            #     shp, shpx = self.shape(xi, self.p)
-            #     Xnum = np.hstack((Xnum, (X[(j - 1) * self.p + self.p - 1] + X[(j - 1) * self.p]) / 2 + (
-            #         X[(j - 1) * self.p + self.p - 1] - X[(j - 1) * self.p]) / 2 * xi))
-            #     Unum = np.hstack(
-            #         (Unum, shp.T.dot(U_1d[int(len(U) / 2) + (j - 1) * self.p:int(len(U) / 2) + j * self.p])))
-            #     Qnum = np.hstack(
-            #         (Qnum, shp.T.dot(U_1d[int((j - 1) * self.p):j * self.p])))
-            # if i in [0, 4, 9, 19]:
-            #     plt.plot(Xnum, Unum, '-', color='C3')
-            #     plt.plot(X, U[int(len(U) / 2):len(U)], '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.clf()
-            trueError[0, i] = self.n_ele
-            trueError[1, i] = U[self.n_ele * self.p - 1] - np.sqrt(self.kappa)
-            self.p = self.p + 1
-        return trueError, estError
-
-
-# test = HDPG1d(2, 2)
-
-# _, trueError, estError = test.adaptive()
-# p = np.arange(3, 4)
-# n_ele = 2**np.arange(1, 8)
-# legend = []
-# errorL2 = np.zeros((n_ele.size, p.size))
-# for i in range(p.size):
-#     test.p = p[i]
-#     for j, n in enumerate(n_ele):
-#         test.n_ele = n
-#         print(test.errorL2())
-#         errorL2[j, i] = test.errorL2()
-
-# plt.loglog(trueError[0, 0:-1], np.abs(trueError[1, 0:-1]), '-ro')
-# plt.axis([1, 250, 1e-13, 1e-2])
-# plt.loglog(n_ele, errorL2, '-o')
-# plt.loglog(estError[0, :], estError[1, :], '--', color='#1f77b4')
-# plt.xlabel('Number of elements', fontsize=17)
-# plt.ylabel('Error', fontsize=17)
-# plt.grid()
-# plt.legend(('Adaptive', 'Uniform', 'Estimator'), loc=3, fontsize=15)
-# # plt.savefig('conv{}p{}_4_test.pdf'.format(15, test.p))
-# plt.show()

hdpg1d/postprocess.py

@@ -5,14 +5,69 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 
-def convHistory(trueError, estError):
-    plt.loglog(trueError[0, 0:-1], np.abs(trueError[1, 0:-1]), '-ro')
-    plt.axis([1, 250, 1e-13, 1e-2])
-    # plt.loglog(n_ele, errorL2, '-o')
-    plt.loglog(estError[0, :], estError[1, :], '--', color='#1f77b4')
-    plt.xlabel('Number of elements', fontsize=17)
-    plt.ylabel('Error', fontsize=17)
-    plt.grid()
-    plt.legend(('Adaptive', 'Estimator'), loc=3, fontsize=15)
-    # plt.savefig('conv{}p{}_4_test.pdf'.format(15, test.p))
-    plt.show()
+class utils(object):
+    exactNumEle = 300
+    exactBasisFuncs = 5
+
+    def __init__(self, solution):
+        self.solution = solution
+        self.numEle = self.exactNumEle
+        self.exactSoln = self.solveExact()
+
+    def solveExact(self):
+        self.solution.coeff.numEle = self.exactNumEle
+        self.solution.coeff.pOrder = self.exactBasisFuncs - 1
+        self.solution.mesh = np.linspace(0, 1, self.exactNumEle + 1)
+        # approximate the exact solution for general problems
+        self.solution.solveLocal()
+        exactSoln = self.solution.separateSoln(self.solution.primalSoln)[0][
+            self.exactNumEle * self.exactBasisFuncs - 1]
+        # for the reaction diffusion test problem, we know the exact solution
+        # self.exactSol = np.sqrt(self.solution.coeff.diffusion)
+        return exactSoln
+
+    def errorL2(self):
+        errorL2 = 0.
+        numEle = self.solution.coeff.numEle
+        self.solution.numEle = numEle
+        numBasisFuncs = self.solution.coeff.pOrder + 1
+        # solve on the uniform mesh
+        self.solution.mesh = np.linspace(0, 1, numEle + 1)
+        self.solution.solveLocal()
+        gradState, _ = self.solution.separateSoln(self.solution.primalSoln)
+        errorL2 = np.abs(
+            gradState[numBasisFuncs * numEle - 1] - self.exactSoln)
+        return errorL2
+
+    def uniConv(self):
+        numBasisFuncs = np.arange(
+            self.solution.numBasisFuncs, self.solution.numBasisFuncs + 1)
+        numEle = 2**np.arange(1, 9)
+        uniError = np.zeros((numEle.size, numBasisFuncs.size))
+        for i in range(numBasisFuncs.size):
+            self.solution.coeff.pOrder = numBasisFuncs[i] - 1
+            for j, n in enumerate(numEle):
+                self.solution.coeff.numEle = n
+                uniError[j, i] = self.errorL2()
+        return numEle, uniError
+
+    def convHistory(self):
+        """Plot the uniform and adaptive convergence history"""
+        print("Please note that the error is calculated using {} "
+              "elements with polynomial order {}."
+              .format(self.exactNumEle, self.exactBasisFuncs))
+        plt.figure(2)
+        trueErrorList = self.solution.trueErrorList
+        trueErrorList[1] = np.abs(trueErrorList[1] - self.exactSoln)
+        estErrorList = self.solution.estErrorList
+        plt.loglog(trueErrorList[0],
+                   trueErrorList[1], '-ro')
+        numEle, errorL2 = self.uniConv()
+        plt.loglog(numEle, errorL2, '-o')
+        plt.loglog(estErrorList[0],
+                   estErrorList[1], '--', color='#1f77b4')
+        plt.xlabel('Number of elements', fontsize=17)
+        plt.ylabel('Error', fontsize=17)
+        plt.grid()
+        plt.legend(('Adaptive', 'Uniform', 'Estimator'), loc=3, fontsize=15)
+        plt.show()

hdpg1d/preprocess.py

@@ -0,0 +1,218 @@
+import numpy as np
+from collections import namedtuple
+from scipy.linalg import block_diag
+
+
+def shape(x, p):
+    """generate p shape functions and its first order derivative
+    at the given location x. x can be an array"""
+    A = np.array([np.linspace(-1, 1, p)]).T**np.arange(p)
+    C = np.linalg.inv(A).T
+    x = np.array([x]).T
+    shp = C.dot((x**np.arange(p)).T)
+    shpx = C[:, 1::1].dot((x**np.arange(p - 1) * np.arange(1, p)).T)
+    return shp, shpx
+
+
+def forcing(x):
+    f = 1
+    return f
+
+
+def boundaryCondition(case):
+    if case == 0:
+        # primal problem
+        bc = [0, 0]
+    if case == 1:
+        # adjoint problem
+        bc = [0, 1]
+    return bc
+
+
+class discretization(object):
+    """Given the problem statement and current mesh,
+    construct the discretization matrices"""
+
+    def __init__(self, coeff, mesh, enrich=None):
+        self.mesh = mesh
+        self.coeff = coeff
+        self.enrich = enrich
+        # the following init are for the sake of simplicity
+        self.numBasisFuncs = coeff.pOrder + 1
+        self.tau_pos = coeff.TAUPLUS
+        self.tau_neg = coeff.TAUMINUS
+        self.conv = coeff.CONVECTION
+        self.kappa = coeff.DIFFUSION
+        self.numEle = 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.numEle)
+        for i in range(1, self.numEle + 1):
+            dist[i - 1] = self.mesh[i] - self.mesh[i - 1]
+        return dist
+
+    def matGroup(self):
+        lhsMat = self.lhsGen()
+        F, L, R = lhsMat.F, lhsMat.L, lhsMat.R
+        eleMat = self.eleGen()
+        A, B, BonU, D = eleMat.A, eleMat.B, eleMat.BonU, eleMat.D
+        faceMat = self.interfaceGen()
+        C, E, G, H = faceMat.C, faceMat.E, faceMat.G, faceMat.H
+        matGroup = namedtuple(
+            'matGroup', ['A', 'B', 'BonU', 'C', 'D', 'E', 'F', 'G', 'H', 'L', 'R'])
+        return matGroup(A, B, BonU, C, D, E, F, G, H, L, R)
+
+    def lhsGen(self):
+        """Generate matrices associated with left hand side"""
+        if self.enrich is not None:
+            # when enrich is given
+            # provide matrices used in the DWR residual calculation
+            numBasisFuncs = self.numBasisFuncs + self.enrich
+            shp = self.shpEnrich
+        else:
+            numBasisFuncs = self.numBasisFuncs
+            shp = self.shp
+        # forcing vector F
+        F = np.zeros(numBasisFuncs * self.numEle)
+        for i in range(1, self.numEle + 1):
+            f = self.dist[i - 1] / 2 \
+                * shp.dot(self.wi * forcing(self.mesh[i - 1] +
+                                            1 / 2 * (1 + self.xi) *
+                                            self.dist[i - 1]))
+            F[(i - 1) * numBasisFuncs:i * numBasisFuncs] = f
+        F[0] += (self.conv + self.tau_pos) * boundaryCondition(0)[0]
+        F[-1] += (-self.conv + self.tau_neg) * boundaryCondition(0)[1]
+        # L, easy in 1d
+        L = np.zeros(self.numEle - 1)
+        # R, easy in 1d
+        R = np.zeros(numBasisFuncs * self.numEle)
+        R[0] = boundaryCondition(0)[0]
+        R[-1] = -boundaryCondition(0)[1]
+        lhsMat = namedtuple('lhsMat', ['F', 'L', 'R'])
+        return lhsMat(F, L, R)
+
+    def eleGen(self):
+        """Generate matrices associated with elements"""
+        if self.enrich is not None:
+            numBasisFuncsEnrich = self.numBasisFuncs + self.enrich
+            shpEnrich = self.shpEnrich
+            shpxEnrich = self.shpxEnrich
+            b = (self.shpx.T * np.ones((self.gqOrder, self.numBasisFuncs))
+                 ).T.dot(np.diag(self.wi).dot(shpEnrich.T))
+            # BonQ is only used in calculating DWR residual
+            BonQ = block_diag(*[b] * (self.numEle)).T
+        else:
+            numBasisFuncsEnrich = self.numBasisFuncs
+            shpEnrich = self.shp
+            shpxEnrich = self.shpx
+            BonQ = None
+        a = 1 / self.coeff.DIFFUSION * \
+            shpEnrich.dot(np.diag(self.wi).dot(self.shp.T))
+        A = np.repeat(self.dist, self.numBasisFuncs) / \
+            2 * block_diag(*[a] * (self.numEle))
+        b = (shpxEnrich.T * np.ones((self.gqOrder, numBasisFuncsEnrich))
+             ).T.dot(np.diag(self.wi).dot(self.shp.T))
+        B = block_diag(*[b] * (self.numEle))
+        d = self.coeff.REACTION * \
+            shpEnrich.dot(np.diag(self.wi).dot(self.shp.T))
+        # assemble D
+        dFace = np.zeros((numBasisFuncsEnrich, self.numBasisFuncs))
+        dFace[0, 0] = self.tau_pos
+        dFace[-1, -1] = self.tau_neg
+        dConv = -self.conv * (shpxEnrich.T * np.ones((self.gqOrder,
+                                                      numBasisFuncsEnrich)))\
+            .T.dot(np.diag(self.wi).dot(self.shp.T))
+        D = np.repeat(self.dist, self.numBasisFuncs) / 2 * \
+            block_diag(*[d] * (self.numEle)) + \
+            block_diag(*[dFace] * (self.numEle)) +\
+            block_diag(*[dConv] * (self.numEle))
+        eleMat = namedtuple('eleMat', ['A', 'B', 'BonU', 'D'])
+        return eleMat(A, B, BonQ, D)
+
+    def interfaceGen(self):
+        """Generate matrices associated with interfaces"""
+        if self.enrich is not None:
+            # when enrich is given
+            # provide matrices used in the DWR residual calculation
+            numBasisFuncs = self.numBasisFuncs + self.enrich
+        else:
+            numBasisFuncs = self.numBasisFuncs
+        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.numEle), dtype=int)
+        map_h[:, 0] = np.arange(2)
+        for i in np.arange(1, self.numEle):
+            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.numEle + 1, self.numEle + 1))
+        for i in range(self.numEle):
+            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.numEle][:, 1:self.numEle]
+
+        # elemental e
+        e = np.zeros((numBasisFuncs, 2))
+        e[0, 0], e[-1, -1] = -conv - tau_pos, conv - tau_neg
+        # mapping matrix
+        map_e_x = np.arange(numBasisFuncs * self.numEle,
+                            dtype=int).reshape(self.numEle,
+                                               numBasisFuncs).T
+        map_e_y = map_h
+        # assemble global E
+        E = np.zeros((numBasisFuncs * self.numEle, self.numEle + 1))
+        for i in range(self.numEle):
+            for j in range(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.numEle]
+
+        # elemental c
+        c = np.zeros((numBasisFuncs, 2))
+        c[0, 0], c[-1, -1] = -1, 1
+        # assemble global C
+        C = np.zeros((numBasisFuncs * self.numEle, self.numEle + 1))
+        for i in range(self.numEle):
+            for j in range(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.numEle]
+
+        # elemental g
+        g = np.zeros((2, numBasisFuncs))
+        g[0, 0], g[-1, -1] = tau_pos, tau_neg
+        # mapping matrix
+        map_g_x = map_h
+        map_g_y = np.arange(numBasisFuncs * self.numEle,
+                            dtype=int).reshape(self.numEle,
+                                               numBasisFuncs).T
+        # assemble global G
+        G = np.zeros((self.numEle + 1, numBasisFuncs * self.numEle))
+        for i in range(self.numEle):
+            for j in range(2):
+                m = map_g_x[j, i]
+                for k in range(numBasisFuncs):
+                    n = map_g_y[k, i]
+                    G[m, n] += g[j, k]
+        G = G[1:self.numEle, :]
+        faceMat = namedtuple('faceMat', ['C', 'E', 'G', 'H'])
+        return faceMat(C, E, G, H)

hdpg1d/solve.py

@@ -1,7 +1,6 @@
 from .coefficients import coefficients
-from .discretization import HDPG1d
-from .postprocess import convHistory
-import sys
+from .adaptation import hdpg1d
+from .postprocess import utils
 
 
 def queryYesNo(question, default="yes"):
@@ -17,28 +16,28 @@ def queryYesNo(question, default="yes"):
         raise ValueError("invalid default answer: '%s'" % default)
 
     while True:
-        sys.stdout.write(question + prompt)
+        print(question + prompt)
         choice = input().lower()
         if default is not None and choice == '':
             return valid[default]
         elif choice in valid:
             return valid[choice]
         else:
-            sys.stdout.write("Please respond with 'yes' or 'no' "
-                             "(or 'y' or 'n').\n")
+            print("Please respond with 'yes' or 'no' "
+                  "(or 'y' or 'n').\n")
 
 
 def getCoefficients():
     question = 'Do you want to use the default parameters?'
     isDefault = queryYesNo(question, "yes")
     if (isDefault):
-        Coeff = coefficients(1, 1, 0, 2, 2)
+        Coeff = coefficients(1e-6, 0, 1, 2, 2, 1, 1)
     else:
-        Coeff = coefficients.from_input()
+        Coeff = coefficients.fromInput()
     return Coeff
 
 
-def run():
+def menu():
     menu = {}
     menu['1.'] = "Solve with HDG."
     menu['2.'] = "Solve with HDPG."
@@ -47,15 +46,29 @@ def run():
     for key, value in sorted(menu.items()):
         print(key, value)
 
+
+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()
+
+
+def runInteractive():
+    menu()
     selection = input("Please Select: ")
-    if selection == '1':
-        hdgCoeff = getCoefficients()
-        hdgSolution = HDPG1d(hdgCoeff.nele, hdgCoeff.porder)
-        trueError, estError = hdgSolution.adaptive()
-        convHistory(trueError, estError)
-    elif selection == '2':
-        print("In development...")
-    elif selection == '3':
-        print("Bye.")
-    else:
-        print("Unknown Option Selected!")
+    while True:
+        if selection == '1':
+            hdgSolve()
+            break
+        elif selection == '2':
+            print("In development...")
+        elif selection == '3':
+            print("Bye.")
+            break
+        else:
+            print("Unknown Option Selected!")
+            continue

requirements.txt

@@ -0,0 +1,3 @@
+          numpy>=1.10
+          matplotlib>=2.0
+          scipy>=0.19

setup.py

@@ -7,7 +7,7 @@ import os
 
 
 setup(name='hdpg1d',
-      version='1.2',
+      version='2.0',
       description='An 1D finite element solver using hybridizable discontinuous\
       Petrov-Galerkin method',
       author='Keyi Ni',