Mercurial > fem-fenics-eugenio
view inst/example/Evolution/Evolution.m @ 239:cb6f11a09667
Move ufl to m-file in Evolution example
author | Eugenio Gianniti <eugenio.gianniti@mail.polimi.it> |
---|---|
date | Tue, 24 Jun 2014 16:55:49 +0200 |
parents | e85390ed620d |
children |
line wrap: on
line source
## Copyright (C) 2013 Marco Vassallo <gedeone-octave@users.sourceforge.net> ## ## This program is free software; you can redistribute it and/or modify it under ## the terms of the GNU General Public License as published by the Free Software ## Foundation; either version 3 of the License, or (at your option) any later ## version. ## ## This program is distributed in the hope that it will be useful, but WITHOUT ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more ## details. ## ## You should have received a copy of the GNU General Public License along with ## this program; if not, see <http://www.gnu.org/licenses/>. pkg load fem-fenics msh ufl start Evolution ufl 'element = FiniteElement("Lagrange", triangle, 1)' ufl ufl u = TrialFunction(element) ufl v = TestFunction(element) ufl ufl u0 = Coefficient(element) ufl ufl dt = Constant(triangle) ufl ufl eq = (1/dt)*(u-u0)*v*dx + "inner(grad(u), grad(v))*dx" ufl a = rhs (eq) ufl L = lhs (eq) ufl end problem = 'Evolution'; mesh = Mesh (msh2m_structured_mesh (0:0.05:1, 0:0.05:1, 1, 1:4)); V = FunctionSpace(problem, mesh); bc = DirichletBC (V, @(x,y) 1, 1:4); t = 0; dt = 0.01; T = 0.3; k = Expression ('dt', @(x,y) dt); u0 = Expression ('u0', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)); a = BilinearForm (problem, V, V, k); L = LinearForm (problem, V, k, u0); A = assemble (a, bc); # solve the problem for each time step # We need to update only the lhs while t < T t += dt; # we can pass u0 to the lhs indifferently as a fem_coeff or # as a fem_func L = LinearForm (problem, V, k, u0); b = assemble (L, bc); u = A \ b; u0 = Function ('u0', V, u); #press Q to make the plot continue plot (u0); end