-
Notifications
You must be signed in to change notification settings - Fork 76
/
ex34.py
48 lines (35 loc) · 1001 Bytes
/
ex34.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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
r"""Euler-Bernoulli beam.
This example solves the Euler-Bernoulli beam equation
.. math::
(EI u'')'' = 1 \quad \text{in $[0,1]$},
with the boundary conditions
:math:`u(0)=u'(0) = 0` and using cubic Hermite elements.
The analytical solution gives :math:`u(1)=1/8`.
"""
from skfem import *
m = MeshLine().refined(3)
e = ElementLineHermite()
basis = InteriorBasis(m, e)
@BilinearForm
def bilinf(u, v, w):
from skfem.helpers import dd, ddot
return ddot(dd(u), dd(v))
@LinearForm
def linf(v, w):
return 1.0 * v
A = asm(bilinf, basis)
f = asm(linf, basis)
D = basis.find_dofs({
'left': m.facets_satisfying(lambda x: x[0] == 0),
})
x = solve(*condense(A, f, D=D))
# compare to analytical solution
err = max(x[basis.nodal_dofs[0]]) - 1. / 8.
print(err)
if __name__ == '__main__':
from os.path import splitext
from sys import argv
name = splitext(argv[0])[0]
from skfem.visuals.matplotlib import *
plot(basis, x, Nrefs=3)
savefig(f'{name}_solution.png')