laplace.py

In this script we solve the Laplace equation \(u_{,kk} = 0\) on a unit square domain \(Ω\) with boundary \(Γ\), subject to boundary conditions:

\[ \begin{align}\begin{aligned}u &= 0 && Γ_{\rm left}\\ ∂_n u &= 0 && Γ_{\rm bottom}\\ ∂_n u &= \cos(1) \cosh(x_1) && Γ_{\rm right}\\ u &= \cosh(1) \sin(x_0) && Γ_{\rm top} \end{aligned}\end{align} \]

This case is constructed to contain all combinations of homogenous and heterogeneous, Dirichlet and Neumann type boundary conditions, as well as to have a known exact solution:

\[u_{\rm exact} = \sin(x_0) \cosh(x_1).\]

We start by importing the necessary modules.

23from nutils import mesh, function, solver, export, cli, testing
24from nutils.expression_v2 import Namespace
25import treelog
26
27
28def main(nelems: int, etype: str, btype: str, degree: int):
29    '''
30    Laplace problem on a unit square.
31
32    .. arguments::
33
34       nelems [10]
35         Number of elements along edge.
36       etype [square]
37         Type of elements (square/triangle/mixed).
38       btype [std]
39         Type of basis function (std/spline), availability depending on the
40         selected element type.
41       degree [1]
42         Polynomial degree.
43    '''

A unit square domain is created by calling the nutils.mesh.unitsquare() mesh generator, with the number of elements along an edge as the first argument, and the type of elements (“square”, “triangle”, or “mixed”) as the second. The result is a topology object domain and a vectored valued geometry function geom.

51    domain, geom = mesh.unitsquare(nelems, etype)

To be able to write index based tensor contractions, we need to bundle all relevant functions together in a namespace. Here we add the geometry x, a scalar basis, and the solution u. The latter is formed by contracting the basis with a to-be-determined solution vector ?lhs.

58    ns = Namespace()
59    ns.x = geom
60    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
61    ns.basis = domain.basis(btype, degree=degree)
62    ns.u = function.dotarg('lhs', ns.basis)

We are now ready to implement the Laplace equation. In weak form, the solution is a scalar field \(u\) for which:

\[∀ v: ∫_Ω \frac{dv}{dx_i} \frac{du}{dx_i} - ∫_{Γ_n} v f = 0.\]

By linearity the test function \(v\) can be replaced by the basis that spans its space. The result is an integral res that evaluates to a vector matching the size of the function space.

73    res = domain.integral('∇_i(basis_n) ∇_i(u) dV' @ ns, degree=degree*2)
74    res -= domain.boundary['right'].integral('basis_n cos(1) cosh(x_1) dS' @ ns, degree=degree*2)

The Dirichlet constraints are set by finding the coefficients that minimize the error:

\[\min_u ∫_{\Gamma_d} (u - u_d)^2\]

The resulting cons array holds numerical values for all the entries of ?lhs that contribute (up to droptol) to the minimization problem. All remaining entries are set to NaN, signifying that these degrees of freedom are unconstrained.

86    sqr = domain.boundary['left'].integral('u^2 dS' @ ns, degree=degree*2)
87    sqr += domain.boundary['top'].integral('(u - cosh(1) sin(x_0))^2 dS' @ ns, degree=degree*2)
88    cons = solver.optimize('lhs', sqr, droptol=1e-15)

The unconstrained entries of ?lhs are to be determined such that the residual vector evaluates to zero in the corresponding entries. This step involves a linearization of res, resulting in a jacobian matrix and right hand side vector that are subsequently assembled and solved. The resulting lhs array matches cons in the constrained entries.

96    lhs = solver.solve_linear('lhs', res, constrain=cons)

Once all entries of ?lhs are establised, the corresponding solution can be vizualised by sampling values of ns.u along with physical coordinates ns.x, with the solution vector provided via the arguments dictionary. The sample members tri and hull provide additional inter-point information required for drawing the mesh and element outlines.

105    bezier = domain.sample('bezier', 9)
106    x, u = bezier.eval(['x_i', 'u'] @ ns, lhs=lhs)
107    export.triplot('solution.png', x, u, tri=bezier.tri, hull=bezier.hull)

To confirm that our computation is correct, we use our knowledge of the analytical solution to evaluate the L2-error of the discrete result.

112    err = domain.integral('(u - sin(x_0) cosh(x_1))^2 dV' @ ns, degree=degree*2).eval(lhs=lhs)**.5
113    treelog.user('L2 error: {:.2e}'.format(err))
114
115    return cons, lhs, err

If the script is executed (as opposed to imported), nutils.cli.run() calls the main function with arguments provided from the command line. For example, to keep with the default arguments simply run python3 laplace.py (view log). To select mixed elements and quadratic basis functions add python3 laplace.py etype=mixed degree=2 (view log).

124if __name__ == '__main__':
125    cli.run(main)

Once a simulation is developed and tested, it is good practice to save a few strategic return values for regression testing. The nutils.testing module, which builds on the standard unittest framework, facilitates this by providing nutils.testing.TestCase.assertAlmostEqual64() for the embedding of desired results as compressed base64 data.

134class test(testing.TestCase):
135
136    def test_default(self):
137        cons, lhs, err = main(nelems=4, etype='square', btype='std', degree=1)
138        with self.subTest('constraints'):
139            self.assertAlmostEqual64(cons, '''
140                eNrbKPv1QZ3ip9sL1BgaILDYFMbaZwZj5ZnDWNfNAeWPESU=''')
141        with self.subTest('left-hand side'):
142            self.assertAlmostEqual64(lhs, '''
143                eNoBMgDN/7Ed9eB+IfLboCaXNKc01DQaNXM14jXyNR82ZTa+NpI2oTbPNhU3bjf7Ngo3ODd+N9c3SNEU
144                1g==''')
145        with self.subTest('L2-error'):
146            self.assertAlmostEqual(err, 1.63e-3, places=5)
147
148    def test_spline(self):
149        cons, lhs, err = main(nelems=4, etype='square', btype='spline', degree=2)
150        with self.subTest('constraints'):
151            self.assertAlmostEqual64(cons, '''
152                eNqrkmN+sEfhzF0xleRbDA0wKGeCYFuaIdjK5gj2aiT2VXMAJB0VAQ==''')
153        with self.subTest('left-hand side'):
154            self.assertAlmostEqual64(lhs, '''
155                eNqrkmN+sEfhzF0xleRbrsauxsnGc43fGMuZJJgmmNaZ7jBlN7M08wLCDLNFZh/NlM0vmV0y+2CmZV5p
156                vtr8j9kfMynzEPPF5lfNAcuhGvs=''')
157        with self.subTest('L2-error'):
158            self.assertAlmostEqual(err, 8.04e-5, places=7)
159
160    def test_mixed(self):
161        cons, lhs, err = main(nelems=4, etype='mixed', btype='std', degree=2)
162        with self.subTest('constraints'):
163            self.assertAlmostEqual64(cons, '''
164                eNorfLZF2ucJQwMC3pR7+QDG9lCquAtj71Rlu8XQIGfC0FBoiqweE1qaMTTsNsOvRtmcoSHbHL+a1UD5
165                q+YAxhcu1g==''')
166        with self.subTest('left-hand side'):
167            self.assertAlmostEqual64(lhs, '''
168                eNorfLZF2ueJq7GrcYjxDJPpJstNbsq9fOBr3Gh8xWS7iYdSxd19xseMP5hImu5UZbv1xljOxM600DTW
169                NN/0k2mC6SPTx6Z1pnNMGc3kzdaaPjRNMbMyEzWzNOsy223mBYRRZpPNJpktMks1azM7Z7bRbIXZabNX
170                ZiLmH82UzS3Ns80vmj004za/ZPYHCD+Y8ZlLmVuYq5kHm9eahwDxavPF5lfNAWFyPdk=''')
171        with self.subTest('L2-error'):
172            self.assertAlmostEqual(err, 1.25e-4, places=6)