In this script we solve the Laplace problem on a unit square that has the bottom-right quadrant removed (a.k.a. an L-shaped domain) with Dirichlet boundary conditions matching the harmonic function

$\sqrt[3]{x^2 + y^2} \cos\left(\tfrac23 \arctan\frac{y+x}{y-x}\right),$

shifted by 0.5 such that the origin coincides with the middle of the unit square. This variation of a well known benchmark problem is known to converge suboptimally under uniform refinement due to a singular gradient in the reentrant corner. This script demonstrates that optimal convergence can be restored by using adaptive refinement.

15from nutils import mesh, function, solver, util, export, cli, testing
16from nutils.expression_v2 import Namespace
17import numpy
18import treelog
19
20
21def main(etype: str, btype: str, degree: int, nrefine: int):
22    '''
23    Adaptively refined Laplace problem on an L-shaped domain.
24
25    .. arguments::
26
27       etype [square]
28         Type of elements (square/triangle/mixed).
29       btype [h-std]
30         Type of basis function (h/th-std/spline), with availability depending on
31         the configured element type.
32       degree [2]
33         Polynomial degree
34       nrefine [5]
35         Number of refinement steps to perform.
36    '''
37
38    domain, geom = mesh.unitsquare(2, etype)
39
40    x, y = geom - .5
41    exact = (x**2 + y**2)**(1/3) * function.cos(function.arctan2(y+x, y-x) * (2/3))
42    domain = domain.trim(exact-1e-15, maxrefine=0)
43    linreg = util.linear_regressor()
44
45    with treelog.iter.fraction('level', range(nrefine+1)) as lrange:
46        for irefine in lrange:
47
48            if irefine:
49                refdom = domain.refined
50                ns.refbasis = refdom.basis(btype, degree=degree)
51                indicator = refdom.integral('∇_k(refbasis_n) ∇_k(u) dV' @ ns, degree=degree*2).eval(lhs=lhs)
52                indicator -= refdom.boundary.integral('refbasis_n ∇_k(u) n_k dS' @ ns, degree=degree*2).eval(lhs=lhs)
53                supp = ns.refbasis.get_support(indicator**2 > numpy.mean(indicator**2))
54                domain = domain.refined_by(refdom.transforms[supp])
55
56            ns = Namespace()
57            ns.x = geom
58            ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
59            ns.basis = domain.basis(btype, degree=degree)
60            ns.u = function.dotarg('lhs', ns.basis)
61            ns.du = ns.u - exact
62
63            sqr = domain.boundary['trimmed'].integral('u^2 dS' @ ns, degree=degree*2)
64            cons = solver.optimize('lhs', sqr, droptol=1e-15)
65
66            sqr = domain.boundary.integral('du^2 dS' @ ns, degree=7)
67            cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
68
69            res = domain.integral('∇_k(basis_n) ∇_k(u) dV' @ ns, degree=degree*2)
70            lhs = solver.solve_linear('lhs', res, constrain=cons)
71
72            ndofs = len(ns.basis)
73            error = function.sqrt(domain.integral(['du du dV', '∇_k(du) ∇_k(du) dV'] @ ns, degree=7)).eval(lhs=lhs)
74            rate, offset = linreg.add(numpy.log(len(ns.basis)), numpy.log(error))
75            treelog.user('ndofs: {ndofs}, L2 error: {error[0]:.2e} ({rate[0]:.2f}), H1 error: {error[1]:.2e} ({rate[1]:.2f})'.format(ndofs=len(ns.basis), error=error, rate=rate))
76
77            bezier = domain.sample('bezier', 9)
78            x, u, du = bezier.eval(['x_i', 'u', 'du'] @ ns, lhs=lhs)
79            export.triplot('sol.png', x, u, tri=bezier.tri, hull=bezier.hull)
80            export.triplot('err.png', x, du, tri=bezier.tri, hull=bezier.hull)
81
82    return ndofs, error, lhs


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 perform four refinement steps with quadratic basis functions starting from a triangle mesh run python3 adaptivity.py etype=triangle degree=2 nrefine=4 (view log).

91if __name__ == '__main__':
92    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.

101class test(testing.TestCase):
102
104        ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='square', degree=2)
105        with self.subTest('degrees of freedom'):
106            self.assertEqual(ndofs, 149)
107        with self.subTest('L2-error'):
108            self.assertAlmostEqual(error[0], 0.00065, places=5)
109        with self.subTest('H1-error'):
110            self.assertAlmostEqual(error[1], 0.03461, places=5)
111        with self.subTest('left-hand side'):
112            self.assertAlmostEqual64(lhs, '''
113                eNo1j6FrQmEUxT8RBi4KllVfMsl3z/nK4zEmLC6bhsKCw2gSw5IPFsymGbZiWnr+By8Ii7Yhsk3BMtC4
114                Z9sJ223ncs85vzvmM9+Yhix8hDIjtnkdHqQSdDDDj1Qajr5qPXN/07MZ2vI4V7UOIvmdO/oEZY45xYDn
115                oR7ikLHAHVpcs2A1TLhChDO+MOeWt5xjYzm6fOQrGxxiZPeoMGaf37hCyU72hB0u6PglPcQcKxRI/KUd
116                7AYLvMPpsqGkCTPumzWf+qV92kKevjK36ozDP/FSnh1iteWiqWuf+oMaKuyKaC1i52rKPokiF2WLA/20
117                bya+ZCPbWKRPpvgFaedebw==''')
118
120        ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='triangle', degree=2)
121        with self.subTest('degrees of freedom'):
122            self.assertEqual(ndofs, 98)
123        with self.subTest('L2-error'):
124            self.assertAlmostEqual(error[0], 0.00138, places=5)
125        with self.subTest('H1-error'):
126            self.assertAlmostEqual(error[1], 0.05324, places=5)
127        with self.subTest('left-hand side'):
128            self.assertAlmostEqual64(lhs, '''
129                eNprMV1oesqU2VTO1Nbko6myWbhpq+kckwST90avjRgYzptYm+YYMwBBk3GQWavZb1NXs2+mm83um1WY
130                bQbyXYEiQWbKZjNM7wJVzjBlYICoPW8CMiXH+LXRR9NwoPkg82xN5IB2MZu2mGabSBnnAbGscYEJj3GV
131                YQAQg/TVGfaA7RI0BsErRjeNeowDgDQPmF9gkmciaJxtArGjzrAKCGWNpYAQAL0kOBE=''')
132
133    def test_mixed_linear(self):
134        ndofs, error, lhs = main(nrefine=2, btype='h-std', etype='mixed', degree=1)
135        with self.subTest('degrees of freedom'):
136            self.assertEqual(ndofs, 34)
137        with self.subTest('L2-error'):
138            self.assertAlmostEqual(error[0], 0.00450, places=5)
139        with self.subTest('H1-error'):
140            self.assertAlmostEqual(error[1], 0.11683, places=5)
141        with self.subTest('left-hand side'):
142            self.assertAlmostEqual64(lhs, '''
143                eNprMT1u6mQyxUTRzMCUAQhazL6b3jNrMYPxp5iA5FtMD+lcMgDxHa4aXzS+6HDV+fKO85cMnC8zMBzS
144                AQDBThbY''')