diff --git a/.gitignore b/.gitignore index d6560bb..f0f4fc6 100644 --- a/.gitignore +++ b/.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/ diff --git a/doc/.src/book/dotxt/approx.do.txt b/doc/.src/book/dotxt/approx.do.txt index 530c1b0..20c6612 100644 --- a/doc/.src/book/dotxt/approx.do.txt +++ b/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 diff --git a/doc/.src/book/dotxt/varform.do.txt b/doc/.src/book/dotxt/varform.do.txt index d5c518f..727065e 100644 --- a/doc/.src/book/dotxt/varform.do.txt +++ b/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. diff --git a/doc/.src/book/src/fe_approx1D_numint.py b/doc/.src/book/src/fe_approx1D_numint.py index 43b349e..ed58833 100644 --- a/doc/.src/book/src/fe_approx1D_numint.py +++ b/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