Fix fe_approx1D_numint

.gitignore

@@ -52,3 +52,10 @@ Trash
 *.pvd
 *.vtu
 *.tar.gz
+/doc/.src/book/src/approx_fenics_by_P2.png
+/doc/.src/book/src/approx_fenics_by_P2.pdf
+/doc/.src/book/src/approx_fenics_by_P1.png
+/doc/.src/book/src/approx_fenics_by_P1.pdf
+/doc/.src/book/src/approx_fenics_by_P0.png
+/doc/.src/book/src/approx_fenics_by_P0.pdf
+/doc/.src/book/answers/ch4/

doc/.src/book/dotxt/approx.do.txt

@@ -5499,7 +5499,7 @@ To solve a system with right-hand side `b` (one-dimensional `numpy`
 array) with a sparse coefficient matrix `A`, we must use some kind of
 a sparse linear system solver. The safest choice is a method based on
 sparse Gaussian elimination. One high-quality package for
-this purpose if "UMFPACK": "https://en.wikipedia.org/wiki/UMFPACK".
+this purpose is "UMFPACK": "https://en.wikipedia.org/wiki/UMFPACK".
 It is interfaced from SciPy by
 
 !bc pycod
@@ -6014,9 +6014,9 @@ N_e = 10
 vertices, cells, dof_map = mesh_uniform(N_e, d=3, Omega=[0,1])
 phi = [basis(len(dof_map[e])-1) for e in range(N_e)]
 A, b = assemble(vertices, cells, dof_map, phi, f)
-c = np.linalg.solve(A, b)
+c = A.LUsolve(b)
 # Make very fine mesh and sample u(x) on this mesh for plotting
-x_u, u = u_glob(c, vertices, cells, dof_map,
+x_u, u, _ = u_glob(c, cells, vertices, dof_map,
                 resolution_per_element=51)
 plot(x_u, u)
 !ec
@@ -6061,7 +6061,7 @@ in the `fe_approx1D_numint` module does this for $u(x)$ and returns
 an array `x` with coordinates and an array `u` with the $u$ values:
 
 !bc pycod
-x, u = u_glob(c, vertices, cells, dof_map,
+x, u, nodes = u_glob(c, cells, vertices, dof_map,
               resolution_per_element=101)
 e = f(x) - u
 !ec
@@ -6153,7 +6153,7 @@ gives four equations to determine $C_{i,j}$ for each $i$. In mathematical
 detail,
 !bt
 \begin{align*}
-\refphi_0 (-1) &= 1,\quad \refphi_0 (1)=\refphi_0'(-1)=\refphi_i' (1)=0,\\
+\refphi_0 (-1) &= 1,\quad \refphi_0 (1)=\refphi_0'(-1)=\refphi_0' (1)=0,\\
 \refphi_1' (-1) &= 1,\quad \refphi_1 (-1)=\refphi_1(1)=\refphi_1' (1)=0,\\
 \refphi_2 (1) &= 1,\quad \refphi_2 (-1)=\refphi_2'(-1)=\refphi_2' (1)=0,\\
 \refphi_3' (1) &= 1,\quad \refphi_3 (-1)=\refphi_3'(-1)=\refphi_3 (1)=0

doc/.src/book/dotxt/varform.do.txt

@@ -1179,10 +1179,10 @@ can subtract them and use the bilinear property of $a(\cdot,\cdot)$:
 \begin{align*}
 a(B+E,v) - a(B+F, v) &= L(v) - L(v)\\
 a(B+E-(B+F),v) &= 0\\
-a(E-F),v) &= 0
+a(E-F,v) &= 0
 \end{align*}
 !et
-If $a(E-F),v) = 0$ for all $v$ in $V$, then $E-F$ must be zero everywhere
+If $a(E-F,v) = 0$ for all $v$ in $V$, then $E-F$ must be zero everywhere
 in the domain, i.e., $E=F$. Or in other words: $u=\uex$. This proves
 that the exact solution is recovered if $\uex - B$ lies in $V$., i.e.,
 can be expressed as $\sum_{j\in\If}d_j\baspsi_j$ where $\{\baspsi_j\}_{j\in\If}$
@@ -3641,7 +3641,7 @@ $B'=C\basphi_0' + D\basphi_{N_n-1}'$.
 
 Most of the terms in the linear system have already been computed
 so we concentrate on the new contribution from the boundary function.
-The integral $C\int_0^L \basphi_{0}'(x))\basphi_i'(x) \dx$, associated
+The integral $C\int_0^L \basphi_{0}'(x)\basphi_i'(x) \dx$, associated
 with the Dirichlet condition in $x=0$,  can only get
 a nonzero contribution from the first cell,
 $\Omega^{(0)}=[\xno{0},\xno{1}]$
@@ -3650,7 +3650,7 @@ $\basphi_{0}'(x)\basphi_i'(x) \dx \neq 0$ only for $i=0$ and $i=1$
 (but node $i=0$ is excluded from the formulation),
 since $\basphi_{i}=0$ on the first cell if $i>1$.
 With a similar reasoning we realize that
-$D\int_0^L \basphi_{N_n-1}'(x))\basphi_i'(x) \dx$ can only get
+$D\int_0^L \basphi_{N_n-1}'(x)\basphi_i'(x) \dx$ can only get
 a nonzero contribution from the last cell.
 From the explanations of the
 calculations in ref[Section ref{fem:approx:global:linearsystem}][ in
@@ -5748,7 +5748,7 @@ and general results for elliptic PDEs.
 
 __Remark.__
 The mathematical theory of finite element methods is primarily
-developed for to stationary PDE problems of elliptic nature whose
+developed for stationary PDE problems of elliptic nature whose
 solutions are smooth. However, such problems can be solved with the
 desired accuracy by most numerical methods and pose no difficulties.
 Time-dependent problems, on the other hand, easily lead to
@@ -6445,7 +6445,7 @@ using a central difference scheme.
 Finite difference upwinding is difficult to express using finite
 elements methods, but it is closely to adding some kind of diffusion
 to the scheme.  A clever method of adding diffusion is the so-called
-*Petrov-Galerkin* method.  The Galerkin approximation consist of
+*Petrov-Galerkin* method.  The Galerkin approximation consists of
 finding $u\in V$ such that for all $v\in V$ the variational
 formulation is satisfied.
 
@@ -6612,7 +6612,7 @@ the boundary layer is not artificially smoothed.
    such that the Dirichlet values are imposed on the unknown parameters.
  * Neumann conditions, involving prescribed values of the derivative
    (or flux) of $u$, are incorporated in boundary terms arising from
-   integrating terms with second-order derivatives by part.
+   integrating terms with second-order derivatives by parts.
    Forgetting to account for the boundary terms implies the
    condition $\partial u/\partial n=0$ at parts of the boundary where
    no Dirichlet condition is set.

doc/.src/book/src/fe_approx1D_numint.py

@@ -85,17 +85,9 @@ def element_vector(f, phi, Omega_e, symbolic=True, numint=None):
                 h = Omega_e[1] - Omega_e[0]
                 detJ = h/2
                 f_func = sym.lambdify([X], f, 'mpmath')
-                # phi is function
-                integrand = lambda X: f_func(X)*phi[r](X)*detJ
-                #integrand = integrand.subs(sym.pi, np.pi)
-                # integrand may still contain symbols like sym.pi that
-                # prevents numerical evaluation...
-                try:
-                    I = mpmath.quad(integrand, [-1, 1])
-                except Exception as e:
-                    print('Could not integrate f*phi[r] numerically:')
-                    print(e)
-                    sys.exit(0)
+                phi_func = sym.lambdify([X], phi[r], 'mpmath')
+                integrand = lambda X: f_func(X)*phi_func(X)*detJ
+                I = mpmath.quad(integrand, [-1, 1])
             b_e[r] = I
     else:
         #phi = [sym.lambdify([X], phi[r]) for r in range(n)]
@@ -249,10 +241,11 @@ def approximate(f, symbolic=False, d=1, N_e=4, numint=None,
             title += ', exact integration'
         else:
             title += ', integration: %s' % numint_name
-        x_u, u, _ = u_glob(c, vertices, cells, dof_map,
+        x_u, u, _ = u_glob(c, cells, vertices, dof_map,
                            resolution_per_element=51)
         x_f = np.linspace(Omega[0], Omega[1], 10001) # mesh for f
-        import scitools.std as plt
+        import matplotlib.pyplot as plt
+        # import scitools.std as plt
         plt.plot(x_u, u, '-',
                  x_f, f(x_f), '--')
         plt.legend(['u', 'f'])
@@ -285,7 +278,12 @@ def u_glob(U, cells, vertices, dof_map, resolution_per_element=51):
         u_cell = 0  # u(x) over this cell
         for r in range(d+1):
             i = dof_map[e][r]  # global dof number
-            u_cell += U[i]*phi[r](X)
+            try:
+                u_cell += U[i]*phi[r](X)
+            except TypeError as error:
+                # if d=0, phi will not be symbolic
+                print(error, '  <--caught!')
+                u_cell += U[i]*phi[r]*np.ones_like(X)
         u_patches.append(u_cell)
         # Compute global coordinates of local nodes,
         # assuming all dofs corresponds to values at nodes
@@ -298,5 +296,6 @@ def u_glob(U, cells, vertices, dof_map, resolution_per_element=51):
     u = np.concatenate(u_patches)
     return x, u, nodes
 
+
 if __name__ == '__main__':
     pass