Removed ex12.ipynb. Renamed the examples.

README.md

@@ -8,7 +8,7 @@ Easy to use finite element assemblers and the related tools.
 ## Installation
 
 The most recent release can be installed using `pip install scikit-fem`.
-Otherwise you can just clone this repository.
+For more recent developments, you can just clone this repository.
 
 ## Feature highlights
 
@@ -18,7 +18,6 @@ Otherwise you can just clone this repository.
   H(curl)-conforming elements
 - No complex dependencies: Most of the functionality require only SciPy and
   NumPy
-- Native Python: No need to compile anything 
 
 ## Usage
 
@@ -56,4 +55,4 @@ Please see the directory *examples* for more instructions.
 ## Contributors
 
 - Tom Gustafsson (Maintainer)
-- [Geordie McBain](https://github.com/gdmcbain) (Bug reports, examples)
+- [Geordie McBain](https://github.com/gdmcbain)

examples/ex12.py

@@ -0,0 +1,56 @@
+"""
+Author: gdmcbain
+
+In this example 'pygmsh' is used to generate a disk, replacing the default
+square of MeshTri() in ex1.py.
+
+A basic postprocessing step in finite element analysis is evaluating linear
+forms over the solution. For the boundary value problem of ex1.py, the integral
+of the solution (normalized by the area) is the 'Boussinesq k-factor' [1]; for
+the square it's roughly 0.03514, for the circle 1/π/8 ≐ 0.03979. Linear forms
+are easily evaluated in skfem using the 1-D arrays assembled using the
+@linear_form decorator. In ex1.py, the linear form required for simple
+integration happens to be the same one used on the right-hand side of the
+differential equation, so it's already to hand.
+"""
+
+from skfem import *
+
+import numpy as np
+
+from pygmsh import generate_mesh
+from pygmsh.built_in import Geometry
+
+geom = Geometry()
+geom.add_physical_surface(geom.add_circle([0.] * 3, 1., .5**3).plane_surface,
+                          'disk')
+points, cells = generate_mesh(geom)[:2]
+m = MeshTri(points[:, :2].T, cells['triangle'].T)
+
+e = ElementTriP1()
+map = MappingAffine(m)
+basis = InteriorBasis(m, e, map, 2)
+
+@bilinear_form
+def laplace(u, du, v, dv, w):
+    return du[0]*dv[0] + du[1]*dv[1]
+
+@linear_form
+def load(v, dv, w):
+    return 1.0*v
+
+A = asm(laplace, basis)
+b = asm(load, basis)
+
+I = m.interior_nodes()
+
+x = 0*b
+x[I] = solve(*condense(A, b, I=I))
+
+area = b @ np.ones_like(x)
+k = b @ x / area**2
+print('area = {:.4f} (exact = {:.4f})'.format(area, np.pi))
+print('k = {:.5f} (exact = 1/8/pi = {:.5f})'.format(k, 1/np.pi/8))
+
+m.plot3(x)
+m.show()

examples/ex13.py

@@ -1,56 +1,77 @@
 """
-Author: gdmcbain
-
-In this example 'pygmsh' is used to generate a disk, replacing the default
-square of MeshTri() in ex1.py.
-
-A basic postprocessing step in finite element analysis is evaluating linear
-forms over the solution. For the boundary value problem of ex1.py, the integral
-of the solution (normalized by the area) is the 'Boussinesq k-factor' [1]; for
-the square it's roughly 0.03514, for the circle 1/π/8 ≐ 0.03979. Linear forms
-are easily evaluated in skfem using the 1-D arrays assembled using the
-@linear_form decorator. In ex1.py, the linear form required for simple
-integration happens to be the same one used on the right-hand side of the
-differential equation, so it's already to hand.
+Author: gdmcbain.
+
+Here's another extension of examples/ex1.py, still solving the Laplace
+equation but now with mixed boundary conditions, two parts isopotential
+(charged and earthed) and the rest insulated. The isopotential parts are
+tagged during the construction of the geometry in pygmsh, as introduced in
+ex13.py.
+
+The example is ∇²u = 0 in Ω = {(x, y) : 1 < x² + y² < 4, 0 < θ < π/2},
+where tan θ = y/x, with u = 0 on y = 0 and u = 1 on x = 0. Although these
+boundaries would be simple enough to identify using the coordinates and
+skfem.assembly.find_dofs as in ex3.py, the present technique generalizes to
+more complicated shapes.
+
+The exact solution is u = 2 θ / π. The field strength is |∇ u|² = 4/π² (x² + y²)
+so the conductance (for unit potential difference and conductivity) is
+‖∇ u‖² = 2 ln 2 / π.
 """
 
 from skfem import *
+from skfem.models.poisson import laplace, mass, unit_load
 
 import numpy as np
 
 from pygmsh import generate_mesh
 from pygmsh.built_in import Geometry
 
 geom = Geometry()
-geom.add_physical_surface(geom.add_circle([0.] * 3, 1., .5**3).plane_surface,
-                          'disk')
-points, cells = generate_mesh(geom)[:2]
-m = MeshTri(points[:, :2].T, cells['triangle'].T)
-
-e = ElementTriP1()
-map = MappingAffine(m)
-basis = InteriorBasis(m, e, map, 2)
+points = []
+lines = []
+radii = [1., 2.]
+lcar = .1
+points.append(geom.add_point([0.] * 3, lcar))  # centre
+for x in radii:
+    points.append(geom.add_point([x, 0., 0.], lcar))
+for y in reversed(radii):
+    points.append(geom.add_point([0., y, 0.], lcar))
+lines.append(geom.add_line(*points[1:3]))
+geom.add_physical_line(lines[-1], 'ground')
+lines.append(geom.add_circle_arc(points[2], points[0], points[3]))
+lines.append(geom.add_line(points[3], points[4]))
+geom.add_physical_line(lines[-1], 'positive')
+lines.append(geom.add_circle_arc(points[4], points[0], points[1]))
+geom.add_physical_surface(
+    geom.add_plane_surface(geom.add_line_loop(lines)), 'domain')
 
-@bilinear_form
-def laplace(u, du, v, dv, w):
-    return du[0]*dv[0] + du[1]*dv[1]
+pts, cells, _, cell_data, field_data = generate_mesh(
+    geom, prune_vertices=False)
 
-@linear_form
-def load(v, dv, w):
-    return 1.0*v
+mesh = MeshTri(pts[:, :2].T, cells['triangle'].T)
+boundaries = {bc:
+              np.unique(cells['line'][cell_data['line']['gmsh:physical'] ==
+                                      field_data[bc][0]])
+              for bc in field_data if field_data[bc][1] == 1}
 
+elements = ElementTriP1()
+basis = InteriorBasis(mesh, elements, MappingAffine(mesh), 2)
 A = asm(laplace, basis)
-b = asm(load, basis)
+b = asm(unit_load, basis)
 
-I = m.interior_nodes()
+dofs = np.setdiff1d(np.arange(0, mesh.p.shape[1]),
+                    np.union1d(boundaries['positive'],
+                               boundaries['ground']))
 
-x = 0*b
-x[I] = solve(*condense(A, b, I=I))
+u = 0.*b
+u[boundaries['positive']] = 1.
+u[dofs] = solve(*condense(A, 0.*b, u, dofs))
 
-area = b @ np.ones_like(x)
-k = b @ x / area**2
-print('area = {:.4f} (exact = {:.4f})'.format(area, np.pi))
-print('k = {:.5f} (exact = 1/8/pi = {:.5f})'.format(k, 1/np.pi/8))
+u_exact = 2 * np.arctan2(mesh.p[1, :], mesh.p[0, :]) / np.pi
+u_error = u - u_exact
+print('L2 error =', np.sqrt(u_error @ asm(mass, basis) @ u_error))
+print('conductance = {:.4f} (exact = 2 ln 2 / pi = {:.4f})'.format(
+    u @ A @ u, 2 * np.log(2) / np.pi))
 
-m.plot3(x)
-m.show()
+mesh.plot3(u)
+mesh.show()

examples/ex14.py

@@ -1,77 +1,42 @@
 """
-Author: gdmcbain.
+Author: gdmcbain
 
-Here's another extension of examples/ex1.py, still solving the Laplace
-equation but now with mixed boundary conditions, two parts isopotential
-(charged and earthed) and the rest insulated. The isopotential parts are
-tagged during the construction of the geometry in pygmsh, as introduced in
-ex13.py.
+Another simple modification of examples/ex1.py, this time showing how to pass
+impose inhomogeneous Dirichlet conditions with condense. 
+The forcing term is suppressed for simplicity. The boundary values
+are set as the real part x**2 - y**2 of an analytic complex function z**2 which
+is harmonic and so that is the exact solution through the domain.
 
-The example is ∇²u = 0 in Ω = {(x, y) : 1 < x² + y² < 4, 0 < θ < π/2},
-where tan θ = y/x, with u = 0 on y = 0 and u = 1 on x = 0. Although these
-boundaries would be simple enough to identify using the coordinates and
-skfem.assembly.find_dofs as in ex3.py, the present technique generalizes to
-more complicated shapes.
-
-The exact solution is u = 2 θ / π. The field strength is |∇ u|² = 4/π² (x² + y²)
-so the conductance (for unit potential difference and conductivity) is
-‖∇ u‖² = 2 ln 2 / π.
+This is checked quantitatively by computing the integral of the squared
+magnitude of the gradient, by evaluating the quadratic form associated with the
+laplacian at the solution; the exact value is 8/3.
 """
 
 from skfem import *
-from skfem.models.poisson import laplace, mass, unit_load
-
-import numpy as np
 
-from pygmsh import generate_mesh
-from pygmsh.built_in import Geometry
+m = MeshTri()
+m.refine(4)
 
-geom = Geometry()
-points = []
-lines = []
-radii = [1., 2.]
-lcar = .1
-points.append(geom.add_point([0.] * 3, lcar))  # centre
-for x in radii:
-    points.append(geom.add_point([x, 0., 0.], lcar))
-for y in reversed(radii):
-    points.append(geom.add_point([0., y, 0.], lcar))
-lines.append(geom.add_line(*points[1:3]))
-geom.add_physical_line(lines[-1], 'ground')
-lines.append(geom.add_circle_arc(points[2], points[0], points[3]))
-lines.append(geom.add_line(points[3], points[4]))
-geom.add_physical_line(lines[-1], 'positive')
-lines.append(geom.add_circle_arc(points[4], points[0], points[1]))
-geom.add_physical_surface(
-    geom.add_plane_surface(geom.add_line_loop(lines)), 'domain')
+e = ElementTriP1()
+map = MappingAffine(m)
+basis = InteriorBasis(m, e, map, 2)
 
-pts, cells, _, cell_data, field_data = generate_mesh(
-    geom, prune_vertices=False)
+@bilinear_form
+def laplace(u, du, v, dv, w):
+    return du[0]*dv[0] + du[1]*dv[1]
 
-mesh = MeshTri(pts[:, :2].T, cells['triangle'].T)
-boundaries = {bc:
-              np.unique(cells['line'][cell_data['line']['gmsh:physical'] ==
-                                      field_data[bc][0]])
-              for bc in field_data if field_data[bc][1] == 1}
+@linear_form
+def load(v, dv, w):
+    return 1.0*v
 
-elements = ElementTriP1()
-basis = InteriorBasis(mesh, elements, MappingAffine(mesh), 2)
 A = asm(laplace, basis)
-b = asm(unit_load, basis)
-
-dofs = np.setdiff1d(np.arange(0, mesh.p.shape[1]),
-                    np.union1d(boundaries['positive'],
-                               boundaries['ground']))
+b = asm(load, basis)
 
-u = 0.*b
-u[boundaries['positive']] = 1.
-u[dofs] = solve(*condense(A, 0.*b, u, dofs))
+I = m.interior_nodes()
 
-u_exact = 2 * np.arctan2(mesh.p[1, :], mesh.p[0, :]) / np.pi
-u_error = u - u_exact
-print('L2 error =', np.sqrt(u_error @ asm(mass, basis) @ u_error))
-print('conductance = {:.4f} (exact = 2 ln 2 / pi = {:.4f})'.format(
-    u @ A @ u, 2 * np.log(2) / np.pi))
+u = (([1., 1.j] @ m.p) ** 2).real          # x**2 - y**2
+u[I] = solve(*condense(A, 0.*b, u, I))
 
-mesh.plot3(u)
-mesh.show()
+print('||grad u||**2 = {:.4f} (exact = 8/3 = {:.4f})'.format(u @ A @ u, 8/3))
+m.plot3(u)
+m.show()

skfem/assembly.py

@@ -452,12 +452,12 @@ def refinterp(self, interp, Nrefs=1):
 
 def asm(kernel, ubasis, vbasis=None, w=None, nthreads=1, assemble=True):
     """
-    Assembly using a kernel function.
+    Assembly finite element matrices.
 
     Parameters
     ----------
     kernel : function handle
-        See Examples.
+        See InteriorBasis or FacetBasis.
     ubasis : InteriorBasis
     vbasis : (OPTIONAL) InteriorBasis
     w : (OPTIONAL) ndarray