laplace.py¶
In this script we solve the Laplace equation \(u_{,kk} = 0\) on a unit square domain \(Ω\) with boundary \(Γ\), subject to boundary conditions:
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:
We start by importing the necessary modules.
23 | import nutils, numpy
|
The main function defines the parameter space for the script. Configurable parameters are the mesh density (in number of elements along an edge), element type (square, triangle, or mixed), type of basis function (std or spline, with availability depending on element type), and polynomial degree.
30 31 32 33 | def main(nelems: 'number of elements along edge' = 10,
etype: 'type of elements (square/triangle/mixed)' = 'square',
btype: 'type of basis function (std/spline)' = 'std',
degree: 'polynomial degree' = 1):
|
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
.
41 | domain, geom = nutils.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
.
48 49 50 51 | ns = nutils.function.Namespace()
ns.x = geom
ns.basis = domain.basis(btype, degree=degree)
ns.u = 'basis_n ?lhs_n'
|
We are now ready to implement the Laplace equation. In weak form, the solution is a scalar field \(u\) for which:
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.
62 63 | res = domain.integral('basis_n,i u_,i d:x' @ ns, degree=degree*2)
res -= domain.boundary['right'].integral('basis_n cos(1) cosh(x_1) d:x' @ ns, degree=degree*2)
|
The Dirichlet constraints are set by finding the coefficients that minimize the error:
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.
75 76 77 | sqr = domain.boundary['left'].integral('u^2 d:x' @ ns, degree=degree*2)
sqr += domain.boundary['top'].integral('(u - cosh(1) sin(x_0))^2 d:x' @ ns, degree=degree*2)
cons = nutils.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.
85 | lhs = nutils.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.
94 95 96 | bezier = domain.sample('bezier', 9)
x, u = bezier.eval(['x_i', 'u'] @ ns, lhs=lhs)
nutils.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.
101 102 103 104 | err = domain.integral('(u - sin(x_0) cosh(x_1))^2 d:x' @ ns, degree=degree*2).eval(lhs=lhs)**.5
nutils.log.user('L2 error: {:.2e}'.format(err))
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).
112 113 | 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.
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 | class test(nutils.testing.TestCase):
@nutils.testing.requires('matplotlib')
def test_default(self):
cons, lhs, err = main(nelems=4, etype='square', btype='std', degree=1)
with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
eNrbKPv1QZ3ip9sL1BgaILDYFMbaZwZj5ZnDWNfNAeWPESU=''')
with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
eNoBMgDN/7Ed9eB+IfLboCaXNKc01DQaNXM14jXyNR82ZTa+NpI2oTbPNhU3bjf7Ngo3ODd+N9c3SNEU
1g==''')
with self.subTest('L2-error'):
self.assertAlmostEqual(err, 1.63e-3, places=5)
@nutils.testing.requires('matplotlib')
def test_spline(self):
cons, lhs, err = main(nelems=4, etype='square', btype='spline', degree=2)
with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
eNqrkmN+sEfhzF0xleRbDA0wKGeCYFuaIdjK5gj2aiT2VXMAJB0VAQ==''')
with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
eNqrkmN+sEfhzF0xleRbrsauxsnGc43fGMuZJJgmmNaZ7jBlN7M08wLCDLNFZh/NlM0vmV0y+2CmZV5p
vtr8j9kfMynzEPPF5lfNAcuhGvs=''')
with self.subTest('L2-error'):
self.assertAlmostEqual(err, 8.04e-5, places=7)
@nutils.testing.requires('matplotlib')
def test_mixed(self):
cons, lhs, err = main(nelems=4, etype='mixed', btype='std', degree=2)
with self.subTest('constraints'): self.assertAlmostEqual64(cons, '''
eNorfLZF2ucJQwMC3pR7+QDG9lCquAtj71Rlu8XQIGfC0FBoiqweE1qaMTTsNsOvRtmcoSHbHL+a1UD5
q+YAxhcu1g==''')
with self.subTest('left-hand side'): self.assertAlmostEqual64(lhs, '''
eNorfLZF2ueJq7GrcYjxDJPpJstNbsq9fOBr3Gh8xWS7iYdSxd19xseMP5hImu5UZbv1xljOxM600DTW
NN/0k2mC6SPTx6Z1pnNMGc3kzdaaPjRNMbMyEzWzNOsy223mBYRRZpPNJpktMks1azM7Z7bRbIXZabNX
ZiLmH82UzS3Ns80vmj004za/ZPYHCD+Y8ZlLmVuYq5kHm9eahwDxavPF5lfNAWFyPdk=''')
with self.subTest('L2-error'):
self.assertAlmostEqual(err, 1.25e-4, places=6)
|