[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Franz Chouly |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Wed, 24 May 2017 05:34:47 -0400 (EDT) |
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')