[Getfem-commits] r5120 - in /trunk/getfem: interface/tests/matlab/ interface/tests/python/ src/getfem/

[Yves . Renard]

Author: renard
Date: Wed Nov  4 09:41:05 2015
New Revision: 5120

URL: http://svn.gna.org/viewcvs/getfem?rev=5120&view=rev
Log:
A working DG method

Modified:
    trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m
    trunk/getfem/interface/tests/python/demo_laplacian_DG.py
    trunk/getfem/src/getfem/getfem_generic_assembly.h

Modified: trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m?rev=5120&r1=5119&r2=5120&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m     (original)
+++ trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m     Wed Nov  4 
09:41:05 2015
@@ -22,15 +22,16 @@
 % Numer. Anal. vol. 39:5, pp 1749-1779, 2002.
 
 % Options for prescribing the Dirichlet condition
-dirichlet_version = 0; % 0 = simplification, 1 = with multipliers,
+dirichlet_version = 3; % 0 = simplification, 1 = with multipliers,
                        % 2 = penalization,  3 = Nitsche's method
 theta = 1;       % Nitsche's method parameter theta
 gamma0 = 0.001;  % Nitsche's method parameter gamma0 (gamma = gamma0*h)
-r = 1e8;         % Penalization parameter
+r = 1e8;         % Penalization parameter for the Dirichlet condition
 draw = true;
 quadrangles = true;
+NX = 20;
 K = 2;           % Degree of the discontinuous finite element method
-interior_penalty_factor = 1E7; % Parameter of the interior penalty term
+interior_penalty_factor = 300*NX; % Parameter of the interior penalty term
 verify_neighbour_computation = true;
 
 asize =  size(who('automatic_var654'));
@@ -38,7 +39,6 @@
 
 % trace on;
 gf_workspace('clear all');
-NX = 20;
 if (quadrangles)
   m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
 else
@@ -61,10 +61,12 @@
 border = gf_mesh_get(m,'outer faces');
 GAMMAD=1;
 gf_mesh_set(m, 'region', GAMMAD, border);
+
 % Inner edges for the interior penalty terms
 in_faces = gf_mesh_get(m,'inner faces');
 INNER_FACES=2;
 gf_mesh_set(m, 'region', INNER_FACES, in_faces);
+
 if (verify_neighbour_computation)
   TEST_FACES=3;
   adjf = gf_mesh_get(m, 'adjacent face', 43, 1);
@@ -115,13 +117,13 @@
 
 % Interior penalty terms
 gf_model_set(md, 'add initialized data', 'alpha', [interior_penalty_factor]);
-jump = '(u-Interpolate(u,neighbour_elt))';
-test_jump = '(Test_u-Interpolate(Test_u,neighbour_elt))';
-grad_mean = '((Grad_u.Normal-Interpolate(Grad_u,neighbour_elt).Normal)/2)';
-grad_test_mean = 
'((Grad_Test_u.Normal-Interpolate(Grad_Test_u,neighbour_elt).Normal)/2)';
-gf_model_set(md, 'add linear generic assembly brick', mim, 
sprintf('(%s)*(%s)', jump, grad_test_mean), INNER_FACES);
-gf_model_set(md, 'add linear generic assembly brick', mim, 
sprintf('(%s)*(%s)', test_jump, grad_mean), INNER_FACES);
-gf_model_set(md, 'add linear generic assembly brick', mim, 
sprintf('alpha*(%s)*(%s)', jump, test_jump), INNER_FACES);
+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)';
+gf_model_set(md, 'add linear generic assembly brick', mim, 
sprintf('-((%s).(%s))', grad_mean, test_jump), INNER_FACES);
+gf_model_set(md, 'add linear generic assembly brick', mim, 
sprintf('-((%s).(%s))', jump, grad_test_mean), INNER_FACES);
+gf_model_set(md, 'add linear generic assembly brick', mim, 
sprintf('-alpha*((%s).(%s))', jump, test_jump), INNER_FACES);
 
 gf_model_get(md, 'solve');
 U = gf_model_get(md, 'variable', 'u');
@@ -145,7 +147,7 @@
   A=gf_asm('generic', mim, 1, '(Grad_u.Normal)*(Grad_Test_u.Normal)', 
TEST_FACES, md);
   B=gf_asm('generic', mim, 1, 
'(Interpolate(Grad_u,neighbour_elt).Normal)*(Interpolate(Grad_Test_u,neighbour_elt).Normal)',
 TEST_FACES, md);
   err_v = err_v + norm(A-B);
-  if (err_v > 1E-14)
+  if (err_v > 1E-13)
     error('Test on neighbour element computation: error to big');
   end
 end

Modified: trunk/getfem/interface/tests/python/demo_laplacian_DG.py
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/python/demo_laplacian_DG.py?rev=5120&r1=5119&r2=5120&view=diff
==============================================================================
--- trunk/getfem/interface/tests/python/demo_laplacian_DG.py    (original)
+++ trunk/getfem/interface/tests/python/demo_laplacian_DG.py    Wed Nov  4 
09:41:05 2015
@@ -40,7 +40,7 @@
 Dirichlet_with_multipliers = True  # Dirichlet condition with multipliers
                                    # or penalization
 dirichlet_coefficient = 1e10       # Penalization coefficient
-interior_penalty_factor = 1e7      # Parameter of the interior penalty term
+interior_penalty_factor = 1e2*NX   # Parameter of the interior penalty term
 verify_neighbour_computation = True;
 
 
@@ -137,15 +137,15 @@
                                                DIRICHLET_BOUNDARY_NUM2,
                                                'DirichletData')
 
-# Interior penalty term
+# Interior penalty terms
 md.add_initialized_data('alpha', [interior_penalty_factor])
-jump = "(u-Interpolate(u,neighbour_elt))";
-test_jump = "(Test_u-Interpolate(Test_u,neighbour_elt))";
-grad_mean = "((Grad_u.Normal-Interpolate(Grad_u,neighbour_elt).Normal)/2)";
-grad_test_mean = 
"((Grad_Test_u.Normal-Interpolate(Grad_Test_u,neighbour_elt).Normal)/2)";
-# md.add_linear_generic_assembly_brick(mim, "({F})*({G})".format(F=jump, 
G=grad_test_mean), INNER_FACES);
-# md.add_linear_generic_assembly_brick(mim, "({F})*({G})".format(F=test_jump, 
G=grad_mean), INNER_FACES);
-md.add_linear_generic_assembly_brick(mim, "alpha*({F})*({G})".format(F=jump, 
G=test_jump), INNER_FACES);
+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}))".format(F=grad_mean, 
G=test_jump), INNER_FACES);
+md.add_linear_generic_assembly_brick(mim, "-(({F}).({G}))".format(F=jump, 
G=grad_test_mean), INNER_FACES);
+md.add_linear_generic_assembly_brick(mim, "alpha*(({F}).({G}))".format(F=jump, 
G=test_jump), INNER_FACES);
 
 gf.memstats()
 # md.listvar()
@@ -174,7 +174,7 @@
   A=gf.asm('generic', mim, 1, '(Grad_u.Normal)*(Grad_Test_u.Normal)', 
TEST_FACES, md)
   B=gf.asm('generic', mim, 1, 
'(Interpolate(Grad_u,neighbour_elt).Normal)*(Interpolate(Grad_Test_u,neighbour_elt).Normal)',
 TEST_FACES, md)
   err_v = err_v + np.linalg.norm(A-B)
-  if (err_v > 1E-14):
+  if (err_v > 1E-13):
     print 'Test on neighbour element computation: error to big: ', err_v
     exit(1)
   

Modified: trunk/getfem/src/getfem/getfem_generic_assembly.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_generic_assembly.h?rev=5120&r1=5119&r2=5120&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_generic_assembly.h   (original)
+++ trunk/getfem/src/getfem/getfem_generic_assembly.h   Wed Nov  4 09:41:05 2015
@@ -566,8 +566,6 @@
   struct ga_instruction_set;
 
   class ga_function {
-    // gerer les arguments et leur modif éventuelle ..
-    // Il faut pouvoir deriver les fonctions et obtenir leur gradients
     mutable ga_workspace local_workspace;
     std::string expr;
     mutable ga_instruction_set *gis;