[Getfem-commits] (no subject)

[Franz Chouly]

branch: devel-franz
commit 0fbe21f4342c0a97af69574e35ff699cd997b44e
Author: Franz CHOULY MCF <address@hidden>
Date:   Wed May 24 11:34:20 2017 +0200

    Demo for residual a posteriori estimator and adaptive refinement, in python 
and for Laplace equation.
---
 .../tests/python/demo_laplacian_aposteriori.py     | 80 +++++++++++++++-------
 1 file changed, 55 insertions(+), 25 deletions(-)

diff --git a/interface/tests/python/demo_laplacian_aposteriori.py 
b/interface/tests/python/demo_laplacian_aposteriori.py
index 684f1ab..56a88f8 100644
--- a/interface/tests/python/demo_laplacian_aposteriori.py
+++ b/interface/tests/python/demo_laplacian_aposteriori.py
@@ -37,7 +37,7 @@ dirichlet_coefficient = 1e10       # Penalization coefficient
 export_mesh = True                 # Draw the mesh after mesh generation or not
 
 
-# Create a simple cartesian mesh
+# Create a mesh
 mo1 = gf.MesherObject('rectangle', [0., 50.], [100., 100.])
 mo2 = gf.MesherObject('rectangle', [50., 0.], [100., 100.])
 mo3 = gf.MesherObject('union', mo1, mo2)
@@ -55,8 +55,8 @@ mesh = gf.Mesh('generate', mo, h, 3)
 #
 fb1 = mesh.outer_faces_with_direction([-1., 0.], 0.01) # Left   (Dirichlet)
 fb2 = mesh.outer_faces_with_direction([0., -1.], 0.01) # Bottom (Neumann)
-fb3 = mesh.outer_faces_in_box([-1., 49.], [101., 101.]) 
-fb4 = mesh.outer_faces_in_box([49., -1.], [101., 101.]) 
+fb3 = mesh.outer_faces_in_box([-1., 10.], [101., 101.]) 
+fb4 = mesh.outer_faces_in_box([10., -1.], [101., 101.]) 
 LEFT_BOUND=1; BOTTOM_BOUND=2; AUX_BOUND1 = 3; AUX_BOUND2 = 4; 
 mesh.set_region(  LEFT_BOUND, fb1)
 mesh.set_region(BOTTOM_BOUND, fb2)
@@ -65,11 +65,6 @@ mesh.set_region(  AUX_BOUND2, fb4)
 mesh.region_subtract(  LEFT_BOUND, AUX_BOUND2)
 mesh.region_subtract(BOTTOM_BOUND, AUX_BOUND1)
 
-if (export_mesh):
-    mesh.export_to_vtk('mesh.vtk');
-    print ('\nYou can view the mesh for instance with');
-    print ('mayavi2 -d mesh.vtk -f ExtractEdges -m Surface \n');
-
 # Create a MeshFem for u and rhs fields of dimension 1 (i.e. a scalar field)
 mfu  = gf.MeshFem(mesh, 1)
 mfP0 = gf.MeshFem(mesh, 1)
@@ -109,26 +104,61 @@ else:
   md.add_Dirichlet_condition_with_penalization(mim, 'u', dirichlet_coefficient,
                                                LEFT_BOUND)
 
-# Interior penalty terms
-# md.add_initialized_data('alpha', [interior_penalty_factor])
-# jump = "((u-Interpolate(u,neighbour_elt))*Normal)"
-# test_jump = "((Test_u-Interpolate(Test_u,neighbour_elt))*Normal)"
-# grad_mean = "((Grad_u+Interpolate(Grad_u,neighbour_elt))*0.5)"
-# grad_test_mean = "((Grad_Test_u+Interpolate(Grad_Test_u,neighbour_elt))*0.5)"
-# md.add_linear_generic_assembly_brick(mim, 
"-(({F}).({G}))-(({H}).({I}))+alpha*(({J}).({K}))".format(F=grad_mean, 
G=test_jump, H=jump, I=grad_test_mean, J=jump, K=test_jump), INNER_FACES);
+for refiter in range(5):
+
+       if (export_mesh):
+               mesh.export_to_vtk('mesh%d.vtk'%refiter);
+               print ('\nYou can view the mesh for instance with');
+               print ('mayavi2 -d mesh%d.vtk -f ExtractEdges -m Surface 
\n'%refiter);
+               
+       # Assembly of the linear system and solve.
+       md.solve()
+
+       # Main unknown
+       U = md.variable('u')
+
+       # Residual a posteriori estimator
+
+       grad_jump = '( (Grad_u-Interpolate(Grad_u,neighbour_elt)).Normal )'
+
+       bulkresidual = 'sqr(element_size*Trace(Hess_u))*Test_psi'
+       edgeresidual = '0.25*element_size*sqr(%s)*(Test_psi + 
Interpolate(Test_psi,neighbour_elt))'%grad_jump
+
+       ETA1tmp = 
gf.asm('generic',mim,1,bulkresidual,-1,md,'psi',1,mfP0,np.zeros(mfP0.nbdof()))
+       ETA1 = ETA1tmp [ ETA1tmp.size - mfP0.nbdof() : ETA1tmp.size ]
+       ETA2tmp = 
gf.asm('generic',mim,1,edgeresidual,INNER_FACES,md,'psi',1,mfP0,np.zeros(mfP0.nbdof()))
+       ETA2 = ETA2tmp [ ETA2tmp.size - mfP0.nbdof() : ETA2tmp.size ]
+       ETA  = np.sqrt ( ETA1 + ETA2 )
+
+       # Export data
+       mfu.export_to_pos('laplacian%d.pos'%refiter, U, 'Computed solution')
+       print 'You can view the solution with (for example):'
+       print 'gmsh laplacian%d.pos'%refiter
+
+       mfP0.export_to_pos('eta1_%d.pos'%refiter, ETA1, 'Bulk residual')
+       print 'You can view eta1 with (for example):'
+       print 'gmsh eta1_%d.pos'%refiter
+
+       mfP0.export_to_pos('eta2_%d.pos'%refiter, ETA2, 'Edge residual')
+       print 'You can view eta2 with (for example):'
+       print 'gmsh eta2_%d.pos'%refiter
+
+       mfu.export_to_vtk('laplacian%d.vtk'%refiter, mfu, U, 'u', mfP0, ETA1, 
'eta1', mfP0, ETA2, 'eta2')
+       print ('mayavi2 -d laplacian%d.vtk -f WarpScalar -m Surface'%refiter)
+
+       # Refine the mesh
+       dd=mfP0.basic_dof_from_cvid()
+
+       ETAElt = np.zeros(dd[0].size)
+       for i in range(dd[0].size):
+               ETAElt[i] = ETA[ dd[0][i] ] 
+
+       mesh.refine(np.where( ETAElt > 0.6*np.max(ETA) ))
+       mesh.optimize_structure()
+
 
 
-# Assembly of the linear system and solve.
-md.solve()
 
-# Main unknown
-U = md.variable('u')
 
-# Export data
-mfu.export_to_pos('laplacian.pos', U, 'Computed solution')
-print 'You can view the solution with (for example):'
-print 'gmsh laplacian.pos'
 
 
-mfu.export_to_vtk('laplacian.vtk', mfu, U, 'u')
-print ('mayavi2 -d laplacian.vtk -f WarpScalar -m Surface')