Hermes2D: Example 12 (03/09/2009)

142 days ago by skregmi

#! /usr/bin/env python # This example uses automatic adaptivity to solve a general second-order linear # equation with non-constant coefficients. # # PDE: -d/dx(a_11(x,y)du/dx) - d/dx(a_12(x,y)du/dy) - d/dy(a_21(x,y)du/dx) - d/dy(a_22(x,y)du/dy) # + a_1(x,y)du/dx + a_21(x,y)du/dy + a_0(x,y)u = rhs(x,y) # # Domain: arbitrary # # BC: Dirichlet for boundary marker 1: u = g_D(x,y) # Natural for any other boundary marker: (a_11(x,y)*nu_1 + a_21(x,y)*nu_2) * dudx # + (a_12(x,y)*nu_1 + s_22(x,y)*nu_2) * dudy = g_N(x,y) # Import modules from hermes2d import Mesh, MeshView, VectorView, OrderView, H1Shapeset, PrecalcShapeset, H1Space, \ WeakForm, Solution, ScalarView, LinSystem, DummySolver, RefSystem, H1OrthoHP from hermes2d.examples.c12 import set_bc, set_forms from hermes2d.examples import get_12_mesh # The following parameters can be changed: P_INIT = 1 # Initial polynomial degree of all mesh elements. THRESHOLD = 0.6 # This is a quantitative parameter of the adapt(...) function and # it has different meanings for various adaptive strategies (see below). STRATEGY = 0 # Adaptive strategy: # STRATEGY = 0 ... refine elements until sqrt(THRESHOLD) times total # error is processed. If more elements have similar errors, refine # all to keep the mesh symmetric. # STRATEGY = 1 ... refine all elements whose error is larger # than THRESHOLD times maximum element error. # STRATEGY = 2 ... refine all elements whose error is larger # than THRESHOLD. # More adaptive strategies can be created in adapt_ortho_h1.cpp. ADAPT_TYPE = 0 # Type of automatic adaptivity: # ADAPT_TYPE = 0 ... adaptive hp-FEM (default), # ADAPT_TYPE = 1 ... adaptive h-FEM, # ADAPT_TYPE = 2 ... adaptive p-FEM. ISO_ONLY = False # Isotropic refinement flag (concerns quadrilateral elements only). # ISO_ONLY = false ... anisotropic refinement of quad elements # is allowed (default), # ISO_ONLY = true ... only isotropic refinements of quad elements # are allowed. MESH_REGULARITY = -1 # Maximum allowed level of hanging nodes: # MESH_REGULARITY = -1 ... arbitrary level hangning nodes (default), # MESH_REGULARITY = 1 ... at most one-level hanging nodes, # MESH_REGULARITY = 2 ... at most two-level hanging nodes, etc. # Note that regular meshes are not supported, this is due to # their notoriously bad performance. ERR_STOP = 0.1 # Stopping criterion for adaptivity (rel. error tolerance between the # fine mesh and coarse mesh solution in percent). NDOF_STOP = 1000 # Adaptivity process stops when the number of degrees of freedom grows # over this limit. This is to prevent h-adaptivity to go on forever. # Load the mesh mesh = Mesh() mesh.load(get_12_mesh()) mesh.refine_all_elements() # Initialize the shapeset and the cache shapeset = H1Shapeset() pss = PrecalcShapeset(shapeset) # Create finite element space space = H1Space(mesh, shapeset) set_bc(space) space.set_uniform_order(P_INIT) # Enumerate basis functions space.assign_dofs() # Initialize the weak formulation wf = WeakForm(1) set_forms(wf) # Visualize solution and mesh sview = ScalarView("Coarse solution", 0, 0, 600, 1000) oview = OrderView("Polynomial orders", 1220, 0, 600, 1000) mview = MeshView("Example 12", 100, 100, 500, 500) # Matrix solver solver = DummySolver() # Adaptivity loop it = 0 ndofs = 0 done = False sln_coarse = Solution() sln_fine = Solution() while 1: print("\n---- Adaptivity step %d ---------------------------------------------\n" % (it+1)) it += 1 # Solve the coarse mesh problem ls = LinSystem(wf, solver) ls.set_spaces(space) ls.set_pss(pss) ls.assemble() ls.solve_system(sln_coarse) # View the solution and mesh sview.show(sln_coarse, lib="mayavi", filename="a%02d.png" % it, notebook=True) mview.show(mesh, space=space, lib="mpl", method="orders", filename="b%02d.png" % it, notebook=True) # Solve the fine mesh problem rs = RefSystem(ls) rs.assemble() rs.solve_system(sln_fine) # Calculate element errors and total error estimate hp = H1OrthoHP(space); err_est = hp.calc_error(sln_coarse, sln_fine) * 100 print("Error estimate: %d" % err_est) print """<html><table border=1>""" print "iter=%02d, err_est=%5.2f%%, NDOF=%d" % (it, err_est, ndofs) print """<tr><td><img src="cell://""" + ("a%02d.png" % it) + """"></td><td><img src="cell://""" + ("b%02d.png" % it) +"""" width="540" height="405"></td></tr>""" print """</table></html>""" # If err_est too large, adapt the mesh if (err_est < ERR_STOP): done = True else: hp.adapt(THRESHOLD, STRATEGY, ADAPT_TYPE)#, ISO_ONLY, MESH_REGULARITY) ndofs = space.assign_dofs() if (ndofs >= NDOF_STOP): done = True sview.show(sln_fine, lib="mayavi", filename="a.png", notebook=True) 
        
WARNING: Output truncated!  

----------------------------------------------
  This is Hermes2D - a C++ library for rapid 
prototyping of adaptive FEM and hp-FEM solvers
     developed by the hp-FEM group at UNR
    and distributed under the GPL license.
   For more details visit http://hpfem.org/.
----------------------------------------------

---- Adaptivity step 1 ---------------------------------------------

Creating matrix sparse structure...
  (dofs: 3, nnz: 7, size: 0.0 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 0 sec)
Linearizer: 1588 verts, 3014 tris in 0.01 sec
Creating matrix sparse structure...
  (dofs: 63, nnz: 759, size: 0.0 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 0.04 sec)
Error estimate: 69

iter=01, err_est=69.65%, NDOF=0




I refined elements: 3

---- Adaptivity step 2 ---------------------------------------------

Creating matrix sparse structure...
  (dofs: 20, nnz: 244, size: 0.0 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 0.01 sec)
Linearizer: 3545 verts, 6779 tris in 0.01 sec
Creating matrix sparse structure...
  (dofs: 175, nnz: 4371, size: 0.1 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 0.21 sec)
Error estimate: 34

iter=02, err_est=34.33%, NDOF=20




I refined elements: 3

---- Adaptivity step 3 ---------------------------------------------

Creating matrix sparse structure...
  (dofs: 38, nnz: 544, size: 0.0 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 0.04 sec)
Linearizer: 5406 verts, 10487 tris in 0.01 sec
Creating matrix sparse structure...
  (dofs: 289, nnz: 8347, size: 0.1 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 0.38 sec)
Error estimate: 9

iter=03, err_est= 9.25%, NDOF=38





...

Assembling stiffness matrix...
  (stages: 1, time: 0.44 sec)
Linearizer: 5259 verts, 10204 tris in 0.01 sec
Creating matrix sparse structure...
  (dofs: 949, nnz: 66679, size: 0.8 MB, time: 0.01 sec)
Assembling stiffness matrix...
  (stages: 1, time: 2.73 sec)
Error estimate: 0

iter=07, err_est= 0.53%, NDOF=168




I refined elements: 4

---- Adaptivity step 8 ---------------------------------------------

Creating matrix sparse structure...
  (dofs: 273, nnz: 16747, size: 0.2 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 0.75 sec)
Linearizer: 5297 verts, 10275 tris in 0.01 sec
Creating matrix sparse structure...
  (dofs: 1469, nnz: 124997, size: 1.5 MB, time: 0.02 sec)
Assembling stiffness matrix...
  (stages: 1, time: 4.93 sec)
Error estimate: 0

iter=08, err_est= 0.24%, NDOF=273




I refined elements: 3

---- Adaptivity step 9 ---------------------------------------------

Creating matrix sparse structure...
  (dofs: 315, nnz: 21649, size: 0.3 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 1.21 sec)
Linearizer: 5274 verts, 10231 tris in 0.01 sec
Creating matrix sparse structure...
  (dofs: 1659, nnz: 155919, size: 1.8 MB, time: 0.02 sec)
Assembling stiffness matrix...
  (stages: 1, time: 6.51 sec)
Error estimate: 0

iter=09, err_est= 0.16%, NDOF=315




I refined elements: 3

---- Adaptivity step 10
---------------------------------------------

Creating matrix sparse structure...
  (dofs: 363, nnz: 28933, size: 0.3 MB, time: 0 sec)
Assembling stiffness matrix...
  (stages: 1, time: 1.24 sec)
Linearizer: 5279 verts, 10240 tris in 0.01 sec
Creating matrix sparse structure...
  (dofs: 1877, nnz: 201447, size: 2.3 MB, time: 0.03 sec)
Assembling stiffness matrix...