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

[Yves . Renard]

Author: renard
Date: Tue Nov  3 11:03:44 2015
New Revision: 5118

URL: http://svn.gna.org/viewcvs/getfem?rev=5118&view=rev
Log:
some corrections and supplementary tests

Modified:
    trunk/getfem/interface/src/gf_mesh_get.cc
    trunk/getfem/interface/tests/matlab/demo_laplacian.m
    trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m
    trunk/getfem/interface/tests/python/demo_laplacian_DG.py
    trunk/getfem/src/getfem_generic_assembly.cc

Modified: trunk/getfem/interface/src/gf_mesh_get.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_mesh_get.cc?rev=5118&r1=5117&r2=5118&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_mesh_get.cc   (original)
+++ trunk/getfem/interface/src/gf_mesh_get.cc   Tue Nov  3 11:03:44 2015
@@ -819,9 +819,25 @@
        outer_faces(*pmesh, in, out, "box");
        );
 
-    
+    /address@hidden CVFIDs = ('adjacent face', @int cvid, @int fid)
+    Return convex face of the neighbour element if it exists.
+    If the convex have more than one neighbour
+    relativley to the face ``f`` (think to bar elements in 3D for instance),
+    return the first face found. @*/
+    sub_command
+      ("adjacent face", 2, 2, 0, 1,
+       check_empty_mesh(pmesh);
+       size_type cv = in.pop().to_convex_number(*pmesh);
+       short_type f = 
in.pop().to_face_number(pmesh->structure_of_convex(cv)->nb_faces());
+       bgeot::convex_face cvf = pmesh->adjacent_face(cv, f);
+       getfem::mesh_region flst;
+       if (cvf.cv != size_type(-1))
+        flst.add(cvf.cv,cvf.f);
+       out.pop().from_mesh_region(flst);
+       );
+
     /address@hidden CVFIDs = ('faces from cvid'[, @ivec CVIDs][, 'merge'])
-    Return a list of convexes faces from a list of convex #id.
+    Return a list of convex faces from a list of convex #id.
 
     `CVFIDs` is a two-rows matrix, the first row lists convex #ids,
     and the second lists face numbers (local number in the convex).

Modified: trunk/getfem/interface/tests/matlab/demo_laplacian.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_laplacian.m?rev=5118&r1=5117&r2=5118&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_laplacian.m        (original)
+++ trunk/getfem/interface/tests/matlab/demo_laplacian.m        Tue Nov  3 
11:03:44 2015
@@ -53,7 +53,7 @@
 border = gf_mesh_get(m,'outer faces');
 % Mark it as boundary GAMMAD=1
 GAMMAD=1;
-gf_mesh_set(m, 'boundary', GAMMAD, border);
+gf_mesh_set(m, 'region', GAMMAD, border);
 if (draw)
   gf_plot_mesh(m, 'regions', [GAMMAD]); % the boundary edges appear in red
   pause(1);

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=5118&r1=5117&r2=5118&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m     (original)
+++ trunk/getfem/interface/tests/matlab/demo_laplacian_DG.m     Tue Nov  3 
11:03:44 2015
@@ -16,7 +16,10 @@
 % Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
 
 % A Poisson problem solved with a Discontinuous Galerkin method (or
-% interior penalty method).
+% interior penalty method). See for instance
+% "Unified analysis of discontinuous Galerkin methods for elliptic
+% problems", D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini, SIAM J.
+% Numer. Anal. vol. 39:5, pp 1749-1779, 2002.
 
 % Options for prescribing the Dirichlet condition
 dirichlet_version = 0; % 0 = simplification, 1 = with multipliers,
@@ -28,6 +31,7 @@
 quadrangles = true;
 K = 2;           % Degree of the discontinuous finite element method
 interior_penalty_factor = 1E7; % Parameter of the interior penalty term
+verify_neighbour_computation = true;
 
 asize =  size(who('automatic_var654'));
 if (asize(1)) draw = false; end;
@@ -56,11 +60,24 @@
 % Detect the border of the mesh
 border = gf_mesh_get(m,'outer faces');
 GAMMAD=1;
-gf_mesh_set(m, 'boundary', GAMMAD, border);
+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, 'boundary', INNER_FACES, in_faces);
+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);
+  if (size(adjf,2) == 0)
+    error('Adjacent face not found');
+  end
+  gf_mesh_set(m, 'region', TEST_FACES, [[43; 1], adjf]);
+  if (draw)
+    figure(2); 
+    gf_plot_mesh(m, 'regions', [TEST_FACES], 'convexes', 'on');
+    figure(1)
+  end
+end
 
 % Mesh plot
 if (draw)
@@ -121,8 +138,19 @@
 
 disp(sprintf('H1 norm of error: %g', err));
 
-if (err > 2E-4)
-   error('Laplacian test: error to big');
+if (verify_neighbour_computation)
+  A=gf_asm('generic', mim, 1, 'u*Test_u*(Normal.Normal)', TEST_FACES, md);
+  B=gf_asm('generic', mim, 1, 
'-Interpolate(u,neighbour_elt)*Interpolate(Test_u,neighbour_elt)*(Interpolate(Normal,neighbour_elt).Normal)',
 TEST_FACES, md);
+  err_v = norm(A-B);
+  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)
+    error('Test on neighbour element computation: error to big');
+  end
 end
 
+if (err > 2E-4)
+  error('Laplacian test: error to big');
+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=5118&r1=5117&r2=5118&view=diff
==============================================================================
--- trunk/getfem/interface/tests/python/demo_laplacian_DG.py    (original)
+++ trunk/getfem/interface/tests/python/demo_laplacian_DG.py    Tue Nov  3 
11:03:44 2015
@@ -1,8 +1,7 @@
 #!/usr/bin/env python
-# -*- coding: UTF8 -*-
 # Python GetFEM++ interface
 #
-# Copyright (C) 2004-2015 Yves Renard, Julien Pommier.
+# Copyright (C) 2015-2015 Yves Renard
 #
 # This file is a part of GetFEM++
 #
@@ -24,18 +23,25 @@
   This program is used to check that python-getfem is working. This is
   also a good example of use of GetFEM++.
 
-  $Id: demo_laplacian.py 4429 2013-10-01 13:15:15Z renard $
+  Poisson problem solved with a Discontinuous Galerkin method (or
+  interior penalty method). See for instance
+  "Unified analysis of discontinuous Galerkin methods for elliptic
+  problems", D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini, SIAM J.
+  Numer. Anal. vol. 39:5, pp 1749-1779, 2002.
+
+  $Id: demo_laplacian_DG.py 4429 2013-10-01 13:15:15Z renard $
 """
 # Import basic modules
 import getfem as gf
 import numpy as np
 
 ## Parameters
-NX = 100                           # Mesh parameter.
+NX = 20                            # Mesh parameter.
 Dirichlet_with_multipliers = True  # Dirichlet condition with multipliers
                                    # or penalization
 dirichlet_coefficient = 1e10       # Penalization coefficient
-interior_penalty_factor = 1e6      # Parameter of the interior penalty term
+interior_penalty_factor = 1e7      # Parameter of the interior penalty term
+verify_neighbour_computation = True;
 
 
 # Create a simple cartesian mesh
@@ -73,6 +79,14 @@
 INNER_FACES=4
 m.set_region(INNER_FACES, in_faces)
 
+if (verify_neighbour_computation):
+  TEST_FACES=5
+  adjf = m.adjacent_face(42, 0);
+  if (len(adjf) != 2):
+    print ('No adjacent edge found, change the element number')
+    exit(1)
+  m.set_region(TEST_FACES, np.array([[42,adjf[0][0]], [0,adjf[1][0]]]));
+  
 
 # Interpolate the exact solution (Assuming mfu is a Lagrange fem)
 Ue = mfu.eval('y*(y-1)*x*(x-1)+x*x*x*x*x')
@@ -125,8 +139,13 @@
 
 # Interior penalty term
 md.add_initialized_data('alpha', [interior_penalty_factor])
-md.add_linear_generic_assembly_brick(mim, 
'alpha*(u-Interpolate(u,neighbour_elt))*Test_u - 
alpha*(u-Interpolate(u,neighbour_elt))*Interpolate(Test_u,neighbour_elt)', 
INNER_FACES)
-
+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);
 
 gf.memstats()
 # md.listvar()
@@ -148,6 +167,17 @@
 print 'You can view the solution with (for example):'
 print 'gmsh laplacian.pos'
 
+if (verify_neighbour_computation):
+  A=gf.asm('generic', mim, 1, 'u*Test_u*(Normal.Normal)', TEST_FACES, md)
+  B=gf.asm('generic', mim, 1, 
'-Interpolate(u,neighbour_elt)*Interpolate(Test_u,neighbour_elt)*(Interpolate(Normal,neighbour_elt).Normal)',
 TEST_FACES, md)
+  err_v = np.linalg.norm(A-B)
+  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):
+    print 'Test on neighbour element computation: error to big: ', err_v
+    exit(1)
+  
 if (H1error > 1e-3):
     print 'Error too large !'
     exit(1)

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5118&r1=5117&r2=5118&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Tue Nov  3 11:03:44 2015
@@ -681,7 +681,6 @@
     }
 
     void clear(void) { 
-      std::string ga_tree_to_string(const ga_tree &tree);
       if (root) clear_node_rec(root); root = current_node = 0;
     }