adaptivity.pyΒΆ

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.

15
import nutils, numpy

The main function defines the parameter space for the script. Configurable parameters are the element type (square, triangle, or mixed), type of basis function (std or spline, with availability depending on element type), polynomial degree, and the number of refinement steps to perform before quitting (by default the script will run forever).

23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def main(etype: 'type of elements (square/triangle/mixed)' = 'square',
         btype: 'type of basis function (h/th-std/spline)' = 'h-std',
         degree: 'polynomial degree' = 2,
         nrefine: 'number of refinement steps (-1 for unlimited)' = -1):

  domain, geom = nutils.mesh.unitsquare(2, etype)

  x, y = geom - .5
  exact = (x**2 + y**2)**(1/3) * nutils.function.cos(nutils.function.arctan2(y+x, y-x) * (2/3))
  domain = domain.trim(exact-1e-15, maxrefine=0)
  linreg = nutils.util.linear_regressor()

  for irefine in nutils.log.count('level'):

    ns = nutils.function.Namespace()
    ns.x = geom
    ns.basis = domain.basis(btype, degree=degree)
    ns.u = 'basis_n ?lhs_n'
    ns.du = ns.u - exact

    sqr = domain.boundary['trimmed'].integral('u^2 d:x' @ ns, degree=degree*2)
    cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)

    sqr = domain.boundary.integral('du^2 d:x' @ ns, degree=7)
    cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)

    res = domain.integral('basis_n,k u_,k d:x' @ ns, degree=degree*2)
    lhs = nutils.solver.solve_linear('lhs', res, constrain=cons)

    ndofs = len(ns.basis)
    error = domain.integral('<du^2, du_,k du_,k>_i d:x' @ ns, degree=7).eval(lhs=lhs)**.5
    rate, offset = linreg.add(numpy.log(len(ns.basis)), numpy.log(error))
    nutils.log.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))

    bezier = domain.sample('bezier', 9)
    x, u, du = bezier.eval(['x_i', 'u', 'du'] @ ns, lhs=lhs)
    nutils.export.triplot('sol.png', x, u, tri=bezier.tri, hull=bezier.hull)
    nutils.export.triplot('err.png', x, du, tri=bezier.tri, hull=bezier.hull)

    if irefine == nrefine:
      break

    refdom = domain.refined
    ns.refbasis = refdom.basis(btype, degree=degree)
    indicator = refdom.integral('refbasis_n,k u_,k d:x' @ ns, degree=degree*2).eval(lhs=lhs)
    indicator -= refdom.boundary.integral('refbasis_n u_,k n_k d:x' @ ns, degree=degree*2).eval(lhs=lhs)
    supp = ns.refbasis.get_support(indicator**2 > numpy.mean(indicator**2))

    domain = domain.refined_by(ns.refbasis.transforms[supp])

  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).

81
82
if __name__ == '__main__':
  nutils.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.

 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
class test(nutils.testing.TestCase):

  @nutils.testing.requires('matplotlib')
  def test_square_quadratic(self):
    ndofs, error, lhs = main(nrefine=2, etype='square', degree=2)
    with self.subTest('degrees of freedom'):
      self.assertEqual(ndofs, 149)
    with self.subTest('L2-error'):
      self.assertAlmostEqual(error[0], 0.00065, places=5)
    with self.subTest('H1-error'):
      self.assertAlmostEqual(error[1], 0.03461, places=5)
    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
      eNo1j6FrQmEUxT8RBi4KllVfMsl3z/nK4zEmLC6bhsKCw2gSw5IPFsymGbZiWnr+By8Ii7Yhsk3BMtC4
      Z9sJ223ncs85vzvmM9+Yhix8hDIjtnkdHqQSdDDDj1Qajr5qPXN/07MZ2vI4V7UOIvmdO/oEZY45xYDn
      oR7ikLHAHVpcs2A1TLhChDO+MOeWt5xjYzm6fOQrGxxiZPeoMGaf37hCyU72hB0u6PglPcQcKxRI/KUd
      7AYLvMPpsqGkCTPumzWf+qV92kKevjK36ozDP/FSnh1iteWiqWuf+oMaKuyKaC1i52rKPokiF2WLA/20
      bya+ZCPbWKRPpvgFaedebw==''')

  @nutils.testing.requires('matplotlib')
  def test_triangle_quadratic(self):
    ndofs, error, lhs = main(nrefine=2, etype='triangle', degree=2)
    with self.subTest('degrees of freedom'):
      self.assertEqual(ndofs, 98)
    with self.subTest('L2-error'):
      self.assertAlmostEqual(error[0], 0.00138, places=5)
    with self.subTest('H1-error'):
      self.assertAlmostEqual(error[1], 0.05324, places=5)
    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
      eNprMV1oesqU2VTO1Nbko6myWbhpq+kckwST90avjRgYzptYm+YYMwBBk3GQWavZb1NXs2+mm83um1WY
      bQbyXYEiQWbKZjNM7wJVzjBlYICoPW8CMiXH+LXRR9NwoPkg82xN5IB2MZu2mGabSBnnAbGscYEJj3GV
      YQAQg/TVGfaA7RI0BsErRjeNeowDgDQPmF9gkmciaJxtArGjzrAKCGWNpYAQAL0kOBE=''')

  @nutils.testing.requires('matplotlib')
  def test_mixed_linear(self):
    ndofs, error, lhs = main(nrefine=2, etype='mixed', degree=1)
    with self.subTest('degrees of freedom'):
      self.assertEqual(ndofs, 34)
    with self.subTest('L2-error'):
      self.assertAlmostEqual(error[0], 0.00450, places=5)
    with self.subTest('H1-error'):
      self.assertAlmostEqual(error[1], 0.11683, places=5)
    with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
      eNprMT1u6mQyxUTRzMCUAQhazL6b3jNrMYPxp5iA5FtMD+lcMgDxHa4aXzS+6HDV+fKO85cMnC8zMBzS
      AQDBThbY''')