Mercurial > fem-fenics-eugenio
changeset 112:c15e5434a512
New examples: resolution of the Darcy-Stokes flow. We follow the example
described in Mardal-Tai-Winther, 2002, pp 1610
* Darcy_Stokes.m : script file for the solution of the problem
* DS.ufl: definition of the problem using P2-P0 element
* Erroreps.ufl: compute velocity error in the energy norm
* Errorp.ufl: compute L2 pressure error
* Errorv.ufl: compute L2 velocity error
| author | gedeone-octave <marcovass89@hotmail.it> |
|---|---|
| date | Sat, 24 Aug 2013 15:39:52 +0200 |
| parents | 011ad0c352f7 |
| children | 065ae2bff9ef |
| files | example/Darcy-Stokes/P2P0/DS.ufl example/Darcy-Stokes/P2P0/Darcy_Stokes.m example/Darcy-Stokes/P2P0/Erroreps.ufl example/Darcy-Stokes/P2P0/Errorp.ufl example/Darcy-Stokes/P2P0/Errorv.ufl |
| diffstat | 5 files changed, 123 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/example/Darcy-Stokes/P2P0/DS.ufl Sat Aug 24 15:39:52 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/example/Darcy-Stokes/P2P0/Darcy_Stokes.m Sat Aug 24 15:39:52 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/example/Darcy-Stokes/P2P0/Erroreps.ufl Sat Aug 24 15:39:52 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