Mercurial > fem-fenics-eugenio

changeset 161:7e922376b3e8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Devel folder for examples which need to be improved
author gedeone-octave <marcovass89@hotmail.it>
date Thu, 12 Sep 2013 15:23:48 +0200
parents 4d64af1cade2
children 344e9bc30b03
files devel/example/Advection-Diffusion/AD_time.m devel/example/Advection-Diffusion/Advection_Diffusion.ufl devel/example/Advection-Diffusion/Reaction.ufl devel/example/Darcy-Stokes/P2P0/DS.ufl devel/example/Darcy-Stokes/P2P0/Darcy_Stokes.m devel/example/Darcy-Stokes/P2P0/Err_eps.ufl devel/example/Darcy-Stokes/P2P0/Err_p.ufl devel/example/Darcy-Stokes/P2P0/Err_v.ufl example/Advection-Diffusion/AD_time.m example/Advection-Diffusion/Advection_Diffusion.ufl example/Advection-Diffusion/Reaction.ufl example/Darcy-Stokes/P2P0/DS.ufl example/Darcy-Stokes/P2P0/Darcy_Stokes.m example/Darcy-Stokes/P2P0/Err_eps.ufl example/Darcy-Stokes/P2P0/Err_p.ufl example/Darcy-Stokes/P2P0/Err_v.ufl
diffstat 16 files changed, 192 insertions(+), 192 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Advection-Diffusion/AD_time.m	Thu Sep 12 15:23:48 2013 +0200
@@ -0,0 +1,54 @@
+# the example is based on the one described here for the 
+# bim pkg:
+# http://wiki.octave.org/Bim_package#3D_Time_dependent_problem
+
+pkg load msh fpl fem-fenics
+
+problem = 'Advection_Diffusion';
+import_ufl_BilinearForm (problem);
+
+problem = 'Reaction';
+import_ufl_BilinearForm (problem);
+import_ufl_FunctionSpace (problem);
+
+x = linspace (0, 1, 40);
+y = z = linspace (0, 1, 20);
+msho = msh3m_structured_mesh (x, y, z, 1, 1:6);
+mshd = Mesh (msho);
+
+x = msho.p(1, :).';
+y = msho.p(2, :).';
+z = msho.p(3, :).';
+ 
+x0 = .2; sx = .1;
+y0 = .2; sy = .1;
+z0 = .8; sz = .1;
+ 
+u = exp (- ((x-x0)/(2*sx)) .^2 - ((y-y0)/(2*sy)) .^2 - ((z-z0)/(2*sz)) .^2);
+
+V  = FunctionSpace ('Reaction', mshd);
+
+f = Expression ('phi', @(x,y,z) x + y -z);
+g = Expression ('mu', @(x,y,z) 0.01);
+
+a = BilinearForm ('Advection_Diffusion', V, f, g);
+m = BilinearForm ('Reaction', V);
+A = assemble (a);
+M = assemble (m);
+
+# Lumped reaction matrix
+ML = diag (sum (M, 1));
+
+function du = f (u, t, A, M)
+  du = - M \ (A * u);
+endfunction 
+ 
+time = linspace (0, 1, 30);
+lsode_options ("integration method", "adams");
+U = lsode (@(u, t) f(u, t, A, ML), u, time);
+ 
+for ii = 1:1:numel (time)
+  name = sprintf ("u_%3.3d", ii);
+  delete ([name ".vtu"]);
+  fpl_vtk_write_field (name, msho, {U(ii,:)', 'u'}, {}, 1);
+endfor
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Advection-Diffusion/Advection_Diffusion.ufl	Thu Sep 12 15:23:48 2013 +0200
@@ -0,0 +1,9 @@
+element = FiniteElement("Lagrange", tetrahedron, 1)
+
+u = TrialFunction(element)
+v = TestFunction(element)
+
+mu = Coefficient(element)
+phi = Coefficient(element)
+
+a = mu*inner(grad(u), grad(v))*dx - u * inner (grad(phi), grad (v))*dx
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Advection-Diffusion/Reaction.ufl	Thu Sep 12 15:23:48 2013 +0200
@@ -0,0 +1,6 @@
+element = FiniteElement("Lagrange", tetrahedron, 1)
+
+u = TrialFunction(element)
+v = TestFunction(element)
+
+a = u*v*dx
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Darcy-Stokes/P2P0/DS.ufl	Thu Sep 12 15:23:48 2013 +0200
@@ -0,0 +1,17 @@
+# P2-P0 finite element
+P2 = VectorElement("CG", triangle, 2)
+P0  = FiniteElement("DG", triangle, 0)
+
+R = FiniteElement ("R", triangle, 0)
+
+W = MixedElement ([P2, P0, R])
+
+(u, p, c) = TrialFunctions(W)
+(v, q, d) = TestFunctions(W)
+
+f = Coefficient(P2)
+g = Coefficient(P0)
+ep  = Constant(triangle)
+
+a = (dot(u, v) + ep*ep*div(u)*div(v) + ep*ep*dot(rot(u), rot(v)) + div(v)*p - div(u)*q + p*d + q*c)*dx
+L = (dot (f, v) - g*q)*dx
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Darcy-Stokes/P2P0/Darcy_Stokes.m	Thu Sep 12 15:23:48 2013 +0200
@@ -0,0 +1,87 @@
+pkg load msh
+pkg load fem-fenics
+
+fem_init_env ();
+
+fem_create_all ('DS');
+fem_create_functional ('Err_v');
+fem_create_functional ('Err_eps');
+fem_create_functional ('Err_p');
+
+uex = @(x, y) [-2*pi*sin(pi*x)*sin(pi*x)*cos(pi*y)*sin(pi*y);...
+                2*pi*sin(pi*x)*cos(pi*x)*sin(pi*y)*sin(pi*y)];
+
+pex = @(x, y) sin(pi*x);
+
+gradp = @(x, y) [pi*cos(pi*x); 0];
+
+laplu = @(x, y) [-4*pi^3*(1 - 4*sin(pi*x)*sin(pi*x))*sin(pi*y)*cos(pi*y);...
+                  4*pi^3*(1 - 4*sin(pi*y)*sin(pi*y))*sin(pi*x)*cos(pi*x)];
+
+eps_vector = [1 2^-2 2^-4 2^-8 0]';
+nel_vector   = [4 8 16 32 64]';
+h_vector = 1./nel_vector;
+
+for ii = 1:numel(nel_vector)
+
+  nel = nel_vector(ii);
+  x = y = linspace (0, 1, nel + 1);
+  msho = msh2m_structured_mesh (x, y, 1, 1:4);
+  mshd = Mesh (msho);
+
+  v = Expression ('v', uex);
+  p = Expression ('p', pex);
+
+  V  = DS_FunctionSpace (mshd);
+  V0 = SubSpace (V, 0);
+
+  bc = DirichletBC (V0, @(x,y) [0, 0], 1:4);
+
+  for jj = 1:numel(eps_vector)
+
+    epsi = eps_vector(jj);
+    ff = @(x, y) uex (x, y) - epsi * epsi * laplu (x, y) - gradp (x, y);
+    f = Expression ('f', ff);
+    g = Expression ('g', @(x,y) 0);
+    ep = Expression ('ep', @(x, y) epsi);
+
+    a = DS_BilinearForm (V, ep);
+    L = DS_LinearForm (V, f, g);
+
+    [A, b] = assemble_system (a, L, bc);
+
+    sol = A \ b;
+
+    func = Function ('uh', V, sol);
+    vh = Function ('vh', func, 0);
+    ph = Function ('ph', func, 1);
+
+    M = Err_v_Functional (mshd, v, vh);
+    err_v_L2(ii,jj) = sqrt (assemble (M));
+
+    const = fem_eval (ph, [0, 0]);
+    pex = @(x, y) sin(pi*x) + const;
+    p = Expression ('p', pex);
+    M = Err_p_Functional (mshd, p, ph);
+    err_p_L2(ii,jj) = sqrt (assemble (M));
+
+    M = Err_eps_Functional (mshd, v, vh, ep);
+    err_v_eps(ii,jj) = sqrt (assemble (M));
+
+  end
+
+end
+
+err_v_L2
+err_p_L2
+err_v_eps
+pause ()
+
+rate_v_L2 = (log (err_v_L2 (1:end-1,:)) - log (err_v_L2 (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
+conv_v_L2 = mean (rate_v_L2)
+
+rate_p_L2 = (log (err_p_L2 (1:end-1,:)) - log (err_p_L2 (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
+conv_p_L2 = mean (rate_p_L2)
+
+rate_v_eps = (log (err_v_eps (1:end-1,:)) - log (err_v_eps (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
+conv_v_eps = mean (rate_v_eps)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Darcy-Stokes/P2P0/Err_eps.ufl	Thu Sep 12 15:23:48 2013 +0200
@@ -0,0 +1,7 @@
+P2 = VectorElement("CG", triangle, 2)
+
+v = Coefficient(P2)
+vh = Coefficient(P2)
+ep  = Constant(triangle)
+
+M = (dot(v-vh,v-vh) + div(v-vh)*div(v-vh) + ep*ep*div(v-vh)*div(v-vh) + ep*ep*dot(rot(v-vh), rot(v-vh)))*dx
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Darcy-Stokes/P2P0/Err_p.ufl	Thu Sep 12 15:23:48 2013 +0200
@@ -0,0 +1,6 @@
+P0  = FiniteElement("DG", triangle, 0)
+
+p = Coefficient(P0)
+ph = Coefficient(P0)
+
+M = ((p-ph)*(p-ph))*dx
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/devel/example/Darcy-Stokes/P2P0/Err_v.ufl	Thu Sep 12 15:23:48 2013 +0200
@@ -0,0 +1,6 @@
+P2 = VectorElement("CG", triangle, 2)
+
+v = Coefficient(P2)
+vh = Coefficient(P2)
+
+M = (dot((v-vh),(v-vh)))*dx
--- a/example/Advection-Diffusion/AD_time.m	Thu Sep 12 15:22:40 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,54 +0,0 @@
-# the example is based on the one described here for the 
-# bim pkg:
-# http://wiki.octave.org/Bim_package#3D_Time_dependent_problem
-
-pkg load msh fpl fem-fenics
-
-problem = 'Advection_Diffusion';
-import_ufl_BilinearForm (problem);
-
-problem = 'Reaction';
-import_ufl_BilinearForm (problem);
-import_ufl_FunctionSpace (problem);
-
-x = linspace (0, 1, 40);
-y = z = linspace (0, 1, 20);
-msho = msh3m_structured_mesh (x, y, z, 1, 1:6);
-mshd = Mesh (msho);
-
-x = msho.p(1, :).';
-y = msho.p(2, :).';
-z = msho.p(3, :).';
- 
-x0 = .2; sx = .1;
-y0 = .2; sy = .1;
-z0 = .8; sz = .1;
- 
-u = exp (- ((x-x0)/(2*sx)) .^2 - ((y-y0)/(2*sy)) .^2 - ((z-z0)/(2*sz)) .^2);
-
-V  = FunctionSpace ('Reaction', mshd);
-
-f = Expression ('phi', @(x,y,z) x + y -z);
-g = Expression ('mu', @(x,y,z) 0.01);
-
-a = BilinearForm ('Advection_Diffusion', V, f, g);
-m = BilinearForm ('Reaction', V);
-A = assemble (a);
-M = assemble (m);
-
-# Lumped reaction matrix
-ML = diag (sum (M, 1));
-
-function du = f (u, t, A, M)
-  du = - M \ (A * u);
-endfunction 
- 
-time = linspace (0, 1, 30);
-lsode_options ("integration method", "adams");
-U = lsode (@(u, t) f(u, t, A, ML), u, time);
- 
-for ii = 1:1:numel (time)
-  name = sprintf ("u_%3.3d", ii);
-  delete ([name ".vtu"]);
-  fpl_vtk_write_field (name, msho, {U(ii,:)', 'u'}, {}, 1);
-endfor
--- a/example/Advection-Diffusion/Advection_Diffusion.ufl	Thu Sep 12 15:22:40 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,9 +0,0 @@
-element = FiniteElement("Lagrange", tetrahedron, 1)
-
-u = TrialFunction(element)
-v = TestFunction(element)
-
-mu = Coefficient(element)
-phi = Coefficient(element)
-
-a = mu*inner(grad(u), grad(v))*dx - u * inner (grad(phi), grad (v))*dx
--- a/example/Advection-Diffusion/Reaction.ufl	Thu Sep 12 15:22:40 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-element = FiniteElement("Lagrange", tetrahedron, 1)
-
-u = TrialFunction(element)
-v = TestFunction(element)
-
-a = u*v*dx
--- a/example/Darcy-Stokes/P2P0/DS.ufl	Thu Sep 12 15:22:40 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,17 +0,0 @@
-# P2-P0 finite element
-P2 = VectorElement("CG", triangle, 2)
-P0  = FiniteElement("DG", triangle, 0)
-
-R = FiniteElement ("R", triangle, 0)
-
-W = MixedElement ([P2, P0, R])
-
-(u, p, c) = TrialFunctions(W)
-(v, q, d) = TestFunctions(W)
-
-f = Coefficient(P2)
-g = Coefficient(P0)
-ep  = Constant(triangle)
-
-a = (dot(u, v) + ep*ep*div(u)*div(v) + ep*ep*dot(rot(u), rot(v)) + div(v)*p - div(u)*q + p*d + q*c)*dx
-L = (dot (f, v) - g*q)*dx
--- a/example/Darcy-Stokes/P2P0/Darcy_Stokes.m	Thu Sep 12 15:22:40 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,87 +0,0 @@
-pkg load msh
-pkg load fem-fenics
-
-fem_init_env ();
-
-fem_create_all ('DS');
-fem_create_functional ('Err_v');
-fem_create_functional ('Err_eps');
-fem_create_functional ('Err_p');
-
-uex = @(x, y) [-2*pi*sin(pi*x)*sin(pi*x)*cos(pi*y)*sin(pi*y);...
-                2*pi*sin(pi*x)*cos(pi*x)*sin(pi*y)*sin(pi*y)];
-
-pex = @(x, y) sin(pi*x);
-
-gradp = @(x, y) [pi*cos(pi*x); 0];
-
-laplu = @(x, y) [-4*pi^3*(1 - 4*sin(pi*x)*sin(pi*x))*sin(pi*y)*cos(pi*y);...
-                  4*pi^3*(1 - 4*sin(pi*y)*sin(pi*y))*sin(pi*x)*cos(pi*x)];
-
-eps_vector = [1 2^-2 2^-4 2^-8 0]';
-nel_vector   = [4 8 16 32 64]';
-h_vector = 1./nel_vector;
-
-for ii = 1:numel(nel_vector)
-
-  nel = nel_vector(ii);
-  x = y = linspace (0, 1, nel + 1);
-  msho = msh2m_structured_mesh (x, y, 1, 1:4);
-  mshd = Mesh (msho);
-
-  v = Expression ('v', uex);
-  p = Expression ('p', pex);
-
-  V  = DS_FunctionSpace (mshd);
-  V0 = SubSpace (V, 0);
-
-  bc = DirichletBC (V0, @(x,y) [0, 0], 1:4);
-
-  for jj = 1:numel(eps_vector)
-
-    epsi = eps_vector(jj);
-    ff = @(x, y) uex (x, y) - epsi * epsi * laplu (x, y) - gradp (x, y);
-    f = Expression ('f', ff);
-    g = Expression ('g', @(x,y) 0);
-    ep = Expression ('ep', @(x, y) epsi);
-
-    a = DS_BilinearForm (V, ep);
-    L = DS_LinearForm (V, f, g);
-
-    [A, b] = assemble_system (a, L, bc);
-
-    sol = A \ b;
-
-    func = Function ('uh', V, sol);
-    vh = Function ('vh', func, 0);
-    ph = Function ('ph', func, 1);
-
-    M = Err_v_Functional (mshd, v, vh);
-    err_v_L2(ii,jj) = sqrt (assemble (M));
-
-    const = fem_eval (ph, [0, 0]);
-    pex = @(x, y) sin(pi*x) + const;
-    p = Expression ('p', pex);
-    M = Err_p_Functional (mshd, p, ph);
-    err_p_L2(ii,jj) = sqrt (assemble (M));
-
-    M = Err_eps_Functional (mshd, v, vh, ep);
-    err_v_eps(ii,jj) = sqrt (assemble (M));
-
-  end
-
-end
-
-err_v_L2
-err_p_L2
-err_v_eps
-pause ()
-
-rate_v_L2 = (log (err_v_L2 (1:end-1,:)) - log (err_v_L2 (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
-conv_v_L2 = mean (rate_v_L2)
-
-rate_p_L2 = (log (err_p_L2 (1:end-1,:)) - log (err_p_L2 (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
-conv_p_L2 = mean (rate_p_L2)
-
-rate_v_eps = (log (err_v_eps (1:end-1,:)) - log (err_v_eps (2:end,:)))./((log (h_vector (1:end-1)) - log (h_vector (2:end))));
-conv_v_eps = mean (rate_v_eps)
--- a/example/Darcy-Stokes/P2P0/Err_eps.ufl	Thu Sep 12 15:22:40 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-P2 = VectorElement("CG", triangle, 2)
-
-v = Coefficient(P2)
-vh = Coefficient(P2)
-ep  = Constant(triangle)
-
-M = (dot(v-vh,v-vh) + div(v-vh)*div(v-vh) + ep*ep*div(v-vh)*div(v-vh) + ep*ep*dot(rot(v-vh), rot(v-vh)))*dx
--- a/example/Darcy-Stokes/P2P0/Err_p.ufl	Thu Sep 12 15:22:40 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-P0  = FiniteElement("DG", triangle, 0)
-
-p = Coefficient(P0)
-ph = Coefficient(P0)
-
-M = ((p-ph)*(p-ph))*dx
--- a/example/Darcy-Stokes/P2P0/Err_v.ufl	Thu Sep 12 15:22:40 2013 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,6 +0,0 @@
-P2 = VectorElement("CG", triangle, 2)
-
-v = Coefficient(P2)
-vh = Coefficient(P2)
-
-M = (dot((v-vh),(v-vh)))*dx