Hermes2D: Example 04 (04/28/2010)

91 days ago by solin

The following text box contains the main code which you can easily modify and evaluate. First click on "Edit this" link on the upper left corner and wait for the browser response. Then click on the input text box to see an "Evaluate" link appear at its bottom. Click on this "Evaluate" link to have the code sent to our computer and interpreted. Have fun and let us know at femhub@googlegroups.com with any problems or questions.

The Code:
# Hermes2D: Example 04 (bc-dirichlet) # This example illustrates how to use non-homogeneous(nonzero) # Dirichlet boundary conditions. # # PDE: Poisson equation -Laplace u = CONST_F, where CONST_F is # a constant right-hand side. It is not difficult to see that # the function u(x,y) = (-CONST_F/4)*(x^2 + y^2) satisfies the # above PDE. Since also the Dirichlet boundary conditions # are chosen to match u(x,y), this function is the exact # solution. # # Note that since the exact solution is a quadratic polynomial, # Hermes will compute it exactly if all mesh elements are quadratic # or higher (then the exact solution lies in the finite element space). # If some elements in the mesh are linear, Hermes will only find # an approximation. # The following parameters can be changed P_INIT = 2 # Initial polynomial degree in all elements UNIFORM_REF_LEVEL = 2 # Number of initial uniform mesh refinements # Import modules from hermes2d import (Mesh, MeshView, H1Shapeset, PrecalcShapeset, H1Space, LinSystem, Solution, ScalarView, WeakForm, DummySolver) from hermes2d.examples.c04 import set_bc from hermes2d.forms import set_forms # Initialize the mesh mesh = Mesh() # Create a mesh from a list of nodes, elements, boundary and nurbs. mesh.create([ [0, -1], [1, -1], [-1, 0], [0, 0], [1, 0], [-1, 1], [0, 1], [0.707106781, 0.707106781] ], [ [0, 1, 4, 3, 0], [3, 4, 7, 0], [3, 7, 6, 0], [2, 3, 6, 5, 0] ], [ [0, 1, 1], [1, 4, 2], [3, 0, 4], [4, 7, 2], [7, 6, 2], [2, 3, 4], [6, 5, 2], [5, 2, 6] ], [ [4, 7, 45], [7, 6, 45] ]) # Initial mesh refinements for i in range(UNIFORM_REF_LEVEL): mesh.refine_all_elements() # Initialize the shapeset and the cache shapeset = H1Shapeset() pss = PrecalcShapeset(shapeset) # Create an H1 space space = H1Space(mesh, shapeset) space.set_uniform_order(P_INIT) # Set boundary conditions set_bc(space) # Enumerate degrees of freedom space.assign_dofs() # Initialize the discrete problem wf = WeakForm() set_forms(wf, -4) # Initialize the linear system and solver solver = DummySolver() sys = LinSystem(wf, solver) sys.set_spaces(space) sys.set_pss(pss) # Assemble the stiffness matrix and solve the system sys.assemble() sln = Solution() sys.solve_system(sln) # Display the solution sln.plot(filename="a.png") # Display the mesh mesh.plot(space=space, filename="b.png") # Positioning the images print """<html><table border=1><tr><td><img src="cell://a.png"></span></td><td><img src="cell://b.png" width="540" height="405"></td></tr></table></html>"""