Fem-fenics: Difference between revisions

Jump to navigation Jump to search

Fem-fenics (edit)

Revision as of 16:30, 8 September 2013

6,224 bytes added ,  8 September 2013
no edit summary
No edit summary
No edit summary
Line 601: Line 601:
</div>
</div>
<div style="clear:both"></div>
<div style="clear:both"></div>
=== Incompressible Navier-Stokes equations ===
<div style="width: 100%;">
  <div style="float:left; width: 48%">
{{Code|Define HyperElasticity problem with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px">
pkg load fem-fenics msh
import_ufl_Problem ("TentativeVelocity");
import_ufl_Problem ("VelocityUpdate");
import_ufl_Problem ("PressureUpdate");
# We can either load the mesh from the file as in Dolfin but we can also use the msh pkg to generate the L-shape domain
name = [tmpnam ".geo"];
fid = fopen (name, "w");
fputs (fid,"Point (1)  = {0, 0, 0, 0.1};\n");
fputs (fid,"Point (2)  = {1, 0, 0, 0.1};\n");
fputs (fid,"Point (3)  = {1, 0.5, 0, 0.1};\n");
fputs (fid,"Point (4)  = {0.5, 0.5, 0, 0.1};\n");
fputs (fid,"Point (5) = {0.5, 1, 0, 0.1};\n");
fputs (fid,"Point (6) = {0, 1, 0,0.1};\n");
fputs (fid,"Line (1)  = {5, 6};\n");
fputs (fid,"Line (2) = {2, 3};\n");
fputs (fid,"Line(3) = {6,1,2};\n");
fputs (fid,"Line(4) = {5,4,3};\n");
fputs (fid,"Line Loop(7) = {3,2,-4,1};\n");
fputs (fid,"Plane Surface(8) = {7};\n");
fclose (fid);
msho = msh2m_gmsh (canonicalize_file_name (name)(1:end-4),"scale", 1,"clscale", .2);
unlink (canonicalize_file_name (name));
mesh = Mesh (msho);
# Define function spaces (P2-P1). From ufl file
#  V = VectorElement("CG", triangle, 2)
#  Q = FiniteElement("CG", triangle, 1)
V = FunctionSpace ('VelocityUpdate', mesh);
Q = FunctionSpace ('PressureUpdate', mesh);
# Define trial and test functions. From ufl file
#  u = TrialFunction(V)
#  p = TrialFunction(Q)
#  v = TestFunction(V)
#  q = TestFunction(Q)
# Set parameter values. From ufl file
#  nu = 0.01
dt = 0.01;
T = 3.;
# Define boundary conditions
noslip = DirichletBC (V, @(x,y) [0; 0], [3, 4]);
outflow = DirichletBC (Q, @(x,y) 0, 2);
# Create functions
u0 = Expression ('u0', @(x,y) [0; 0]);
# Define coefficients
k = Constant ('k', dt);
f = Constant ('f', [0; 0]);
# Tentative velocity step. From ufl file
#  eq = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
#      nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
a1 = BilinearForm ('TentativeVelocity', V, k);
# Pressure update. From ufl file
#  a = inner(grad(p), grad(q))*dx
#  L = -(1/k)*div(u1)*q*dx
a2 = BilinearForm ('PressureUpdate', Q);
# Velocity update
#  a = inner(u, v)*dx
#  L = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
a3 = BilinearForm ('VelocityUpdate', V);
# Assemble matrices
A1 = assemble (a1, noslip);
A3 = assemble (a3, noslip);
# Time-stepping
t = dt; i = 0;
while t < T
  # Update pressure boundary condition
  inflow = DirichletBC (Q, @(x,y) sin(3.0*t), 1);
  # Compute tentative velocity step
  "Computing tentative velocity"
  L1 = LinearForm ('TentativeVelocity', V, k, u0, f);
  b1 = assemble (L1, noslip);
  utmp = A1 \ b1;
  u1 = Function ('u1', V, utmp);
  # Pressure correction
  "Computing pressure correction"
  L2 = LinearForm ('PressureUpdate', Q, u1, k);
  [A2, b2] = assemble_system (a2, L2, inflow, outflow);
  ptmp = A2 \ b2;
  p1 = Function ('p1', Q, ptmp);
  # Velocity correction
  "Computing velocity correction"
  L3 = LinearForm ('VelocityUpdate', V, k, u1, p1);
  b3 = assemble (L3, noslip);
  ut = A3 \ b3;
  u1 = Function ('u0', V, ut);
  # Plot solution
  plot (p1);
  plot (u1);
  # Save to file
  save (p1, sprintf ("p_%3.3d", ++i));
  save (u1, sprintf ("u_%3.3d", i));
  # Move to next time step
  u0 = u1;
  t += dt
end
</syntaxhighlight>}}
  </div>
  <div style="float:right; width: 48%">
{{Code|Define NS problem with fenics python|<syntaxhighlight lang="python" style="font-size:13px">
from dolfin import *
# Load mesh from file
mesh = Mesh("lshape.xml.gz")
# Define function spaces (P2-P1)
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
# Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)
# Set parameter values
dt = 0.01
T = 3
nu = 0.01
# Define time-dependent pressure boundary condition
p_in = Expression("sin(3.0*t)", t=0.0)
# Define boundary conditions
noslip  = DirichletBC(V, (0, 0),
                      "on_boundary && \
                      (x[0] < DOLFIN_EPS | x[1] < DOLFIN_EPS | \
                      (x[0] > 0.5 - DOLFIN_EPS && x[1] > 0.5 - DOLFIN_EPS))")
inflow  = DirichletBC(Q, p_in, "x[1] > 1.0 - DOLFIN_EPS")
outflow = DirichletBC(Q, 0, "x[0] > 1.0 - DOLFIN_EPS")
bcu = [noslip]
bcp = [inflow, outflow]
# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)
# Define coefficients
k = Constant(dt)
f = Constant((0, 0))
# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
    nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
# Pressure update
a2 = inner(grad(p), grad(q))*dx
L2 = -(1/k)*div(u1)*q*dx
# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"
# Create files for storing solution
ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")
# Time-stepping
t = dt
while t < T + DOLFIN_EPS:
    # Update pressure boundary condition
    p_in.t = t
    # Compute tentative velocity step
    begin("Computing tentative velocity")
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, u1.vector(), b1, "gmres", "default")
    end()
    # Pressure correction
    begin("Computing pressure correction")
    b2 = assemble(L2)
    [bc.apply(A2, b2) for bc in bcp]
    solve(A2, p1.vector(), b2, "gmres", prec)
    end()
    # Velocity correction
    begin("Computing velocity correction")
    b3 = assemble(L3)
    [bc.apply(A3, b3) for bc in bcu]
    solve(A3, u1.vector(), b3, "gmres", "default")
    end()
    # Plot solution
    plot(p1, title="Pressure", rescale=True)
    plot(u1, title="Velocity", rescale=True)
    # Save to file
    ufile << u1
    pfile << p1
    # Move to next time step
    u0.assign(u1)
    t += dt
    print "t =", t
# Hold plot
interactive()
</syntaxhighlight>}}
  </div>
</div>
<div style="clear:both"></div>
[[Category:OctaveForge]]
[[Category:OctaveForge]]
[[Category:Packages]]
[[Category:Packages]]
Gedeone
60

edits

Retrieved from "https://wiki.octave.org/Special:MobileDiff/3137"

Navigation menu