Hermes2D: Example 05 (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 05 (bc-neumann) # This example shows how to define Neumann boundary conditions. In addition, # you will see how a Filter is used to visualize gradient of the solution # # PDE: Poisson equation -Laplace u = f, where f = CONST_F # # BC: u = 0 on Gamma_4 (edges meeting at the re-entrant corner) # du/dn = CONST_GAMMA_1 on Gamma_1 (bottom edge) # du/dn = CONST_GAMMA_2 on Gamma_2 (top edge, circular arc, and right-most edge) # du/dn = CONST_GAMMA_3 on Gamma_3 (left-most edge) # # You can play with the parameters below. For most choices of the four constants, # the solution has a singular (infinite) gradient at the re-entrant corner. # Therefore we visualize not only the solution but also its gradient. # The following parameters can be changed P_INIT = 4 # Initial polynomial degree in all elements CORNER_REF_LEVEL = 12 # Number of mesh refinements towards the re-entrant corner # Import modules from hermes2d import Mesh, MeshView, H1Shapeset, PrecalcShapeset, H1Space, \ LinSystem, Solution, ScalarView, WeakForm, DummySolver from hermes2d.examples.c05 import set_bc, set_forms from hermes2d.examples.c05 import set_forms as set_forms_surf 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, 3] ], [ [4, 7, 45], [7, 6, 45] ]) # Initial mesh refinements mesh.refine_towards_vertex(3, CORNER_REF_LEVEL) # 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(1) set_forms(wf, -1) set_forms_surf(wf) # 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 sln = Solution() sys.assemble() sys.solve_system(sln) # Show the solution sln.plot(filename="a.png") # Show 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>"""