[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r5118 - in /trunk/getfem: interface/src/ interface/test
From: |
Yves . Renard |
Subject: |
[Getfem-commits] r5118 - in /trunk/getfem: interface/src/ interface/tests/matlab/ interface/tests/python/ src/ |
Date: |
Tue, 03 Nov 2015 10:03:45 -0000 |
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;
}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r5118 - in /trunk/getfem: interface/src/ interface/tests/matlab/ interface/tests/python/ src/,
Yves . Renard <=