# Hermes2D: Example 08 (system)
# This example explains how to create two spaces over a mesh and use them
# to solve a simple problem of linear elasticity. At the end, VonMises
# filter is used to visualize the stress.
#
# PDE: Lame equations of linear elasticity
#
# BC: du_1/dn = f_0 on Gamma_3 and du_1/dn = 0 on Gamma_2, Gamma_4, Gamma_5
# du_2/dn = f_1 on Gamma_3 and du_2/dn = 0 on Gamma_2, Gamma_4, Gamma_5
# u_1 = 0 and u_2 = 0 on Gamma_1
# The following parameters can be changed
P_INIT = 8
# Import modules
from hermes2d import Mesh, MeshView, H1Shapeset, PrecalcShapeset, H1Space, \
LinSystem, WeakForm, DummySolver, Solution, ScalarView, VonMisesFilter
from hermes2d.examples.c08 import set_bc, set_forms
from hermes2d.examples import get_sample_mesh
# Initialize mesh
mesh = Mesh()
# Create a mesh from a list of nodes, elements, boundary and nurbs.
mesh.create([
[0.1, 0],
[0.07071067809999999, 0.07071067809999999],
[0, 0.1],
[0.1707106781, 0],
[0.1707106781, 0.07071067809999999],
[0.1707106781, 0.1707106781],
[0.07071067809999999, 0.1707106781],
[0, 0.1707106781],
[1, 0],
[1, 0.07071067809999999],
[1, 0.1707106781],
[1, 1],
[0.1707106781, 1],
[0.07071067809999999, 1],
[0, 1],
[-0.9, 1],
[-0.9, 0.1707106781],
[-0.9, 0.1]
], [
[4, 1, 0, 0],
[6, 2, 1, 0],
[3, 4, 0, 0],
[6, 7, 2, 0],
[1, 4, 5, 6, 0],
[3, 8, 9, 4, 0],
[4, 9, 10, 5, 0],
[5, 10, 11, 12, 0],
[6, 5, 12, 13, 0],
[7, 6, 13, 14, 0],
[16, 7, 14, 15, 0],
[17, 2, 7, 16, 0]
], [
[1, 0, 5],
[2, 1, 5],
[0, 3, 1],
[3, 8, 1],
[8, 9, 2],
[9, 10, 2],
[10, 11, 2],
[11, 12, 3],
[12, 13, 3],
[13, 14, 3],
[14, 15, 3],
[15, 16, 4],
[17, 2, 5],
[16, 17, 4]
], [
[1, 0, -45],
[2, 1, -45]
])
# Initialize the shapeset and the cache
shapeset = H1Shapeset()
pss = PrecalcShapeset(shapeset)
# Create the x displacement space
xdisp = H1Space(mesh, shapeset)
xdisp.set_uniform_order(P_INIT)
# Create the y displacement space
ydisp = H1Space(mesh, shapeset)
ydisp.set_uniform_order(P_INIT)
# Set boundary conditions
set_bc(xdisp, ydisp)
# Enumerate degrees of freedom
ndof = xdisp.assign_dofs(0)
ndof += ydisp.assign_dofs(ndof)
# Initialize the weak formulation
wf = WeakForm(2)
set_forms(wf)
# Initialize the linear system and solver
solver = DummySolver()
sys = LinSystem(wf, solver)
sys.set_spaces(xdisp, ydisp)
sys.set_pss(pss)
# Assemble the stiffness matrix and solve the system
xsln = Solution()
ysln = Solution()
sys.assemble()
sys.solve_system(xsln, ysln, lib="scipy")
# Visualize the solution
view = ScalarView("Von Mises stress [Pa]", 50, 50, 1200, 600)
E = float(200e9)
nu = 0.3
l = (E * nu) / ((1 + nu) * (1 - 2*nu))
mu = E / (2*(1 + nu))
stress = VonMisesFilter(xsln, ysln, mu, l)
view.show(stress, lib="mayavi", filename="a.png", notebook=True)
# Display the mesh
mesh.plot(space=xdisp, filename="b.png")
# Positioning images
print """<html><table border=1><tr><td><img src="cell://a.png" width="540" height="405"></td><td><img src="cell://b.png" width="640" height="405"></td></tr></table></html>"""
