from hermes2d import (Mesh, MeshView, H1Shapeset, PrecalcShapeset, H1Space,
WeakForm, Solution, DummySolver, LinSystem, ScalarView, RefSystem,
H1OrthoHP, set_verbose)
from hermes2d.examples.c22 import set_bc, set_forms
from enthought.mayavi import mlab
set_verbose(False)
def calc(p_init = 1, threshold=0.3, strategy=0, h_only=True, err_tol=1, \
interactive_plot=False, show_mesh=True, show_graph=False, nb=True):
print "Version 11/02/2009"
# creating a mesh
mesh = Mesh()
mesh.create([
[0, 0],
[1, 0],
[1, 1],
[0, 1],
], [
[2, 3, 0, 1, 0],
], [
[0, 1, 1],
[1, 2, 1],
[2, 3, 1],
[3, 0, 1],
], [])
# initial mesh refinements (more types available)
mesh.refine_all_elements()
# defining shapeset
shapeset = H1Shapeset()
pss = PrecalcShapeset(shapeset)
# defining finite element space
space = H1Space(mesh, shapeset)
set_bc(space)
space.set_uniform_order(p_init)
# choosing weak formulation (to be improved soon)
wf = WeakForm(1)
set_forms(wf)
# basic and fine mesh solutions
sln = Solution()
rsln = Solution()
solver = DummySolver()
# visualization
view = ScalarView("Solution")
mview = MeshView("Mesh")
iter = 0
graph = []
print "Calculating..."
# adaptivity loop
while 1:
space.assign_dofs()
sys = LinSystem(wf, solver)
sys.set_spaces(space)
sys.set_pss(pss)
sys.assemble()
sys.solve_system(sln)
dofs = sys.get_matrix().shape[0]
if interactive_plot:
view.show(sln, lib="mayavi", notebook=nb, filename="a%02d.png" % iter)
if show_mesh:
mview.show(mesh, lib="mpl", method="orders", notebook=nb, filename="b%02d.png" % iter, space=space)
rsys = RefSystem(sys)
rsys.assemble()
rsys.solve_system(rsln)
hp = H1OrthoHP(space)
err_est = hp.calc_error(sln, rsln)*100
print """<html><table border=1>"""
print "iter=%02d, err_est=%5.2f%%, NDOF=%d" % (iter, err_est, dofs)
print """<tr><td><img src="cell://""" + ("a%02d.png" % iter) + """"></td><td><img src="cell://""" + ("b%02d.png" % iter) +"""" width="540" height="405"></td></tr>"""
print """</table></html>"""
graph.append([dofs, err_est])
if err_est < err_tol:
break
hp.adapt(threshold, strategy, h_only)
iter += 1
if not interactive_plot:
view.show(sln, lib="mayavi", notebook=nb)
if show_graph:
from numpy import array
graph = array(graph)
import pylab
pylab.clf()
pylab.plot(graph[:, 0], graph[:, 1], "ko", label="error estimate")
pylab.plot(graph[:, 0], graph[:, 1], "k-")
pylab.title("Error Convergence for the Inner Layer Problem")
pylab.legend()
pylab.xlabel("Degrees of Freedom")
pylab.ylabel("Error [%]")
pylab.yscale("log")
pylab.grid()
pylab.savefig("graph.png")
