streamline.py

利用有限元算法计算各类微分方程

# Modified by Anders Logg, Johan Jansson 2008

__author__ = "Kristian B. Oelgaard ([e-mail])"
__date__ = "2007-11-14 -- 2008-01-04"
__copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
__license__  = "GNU LGPL Version 2.1"

from dolfin import *

# Load mesh and create finite elements
mesh = Mesh("Dolfin-2.xml.gz")
scalar = FiniteElement("Lagrange", "triangle", 1)
vector = VectorElement("Lagrange", "triangle", 2)

# Load subdomains and velocity
sub_domains = MeshFunction("uint", mesh, "Subdomains.xml.gz");
velocity = Function(vector, "Velocity.xml.gz");

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(on_boundary and
                    x[0] < 1.0 - DOLFIN_EPS and
                    x[0] > 0.0 + DOLFIN_EPS and
                    x[1] > 0.0 + DOLFIN_EPS and
                    x[1] < 1.0 - DOLFIN_EPS)


# Initialise source function and previous solution function
f = Function(scalar, mesh, 0.0)
u0 = Function(scalar, mesh, 0.0)

# Parameters
T = 5.0
k = 0.1
t = k
c = 0.0001

# Test and trial functions
v = TestFunction(scalar)
u = TrialFunction(scalar)

# Functions
x0 = Vector()
x1 = Vector()
u0 = Function(scalar, mesh, x0)
u1 = Function(scalar, mesh, x1)
h = 2.0 * MeshSize("triangle", mesh)
#h = MeshSize("triangle", mesh)
#h = 0.0

# Variational problem for cG(1)dG(0) for convection-diffusion
#standard Galerkin method
a = v*u*dx + k*(v*dot(velocity, grad(u))*dx + c*dot(grad(v), grad(u))*dx)
L = v*u0*dx - k*v*f*dx

#streamline diffusion
delta = h   #norm(veloctiy)=1
#a = v*u*dx + k*(dot(velocity, grad(u))*(v+delta*dot(velocity,grad(v)))*dx + c*dot(grad(v), grad(u))*dx)
#L = v*u0*dx - k*v*f*dx


# Set up boundary condition
boundary = DirichletBoundary()
g = Function(mesh, 1.0)
zero = Function(mesh, 0.0)
bc = DirichletBC(g, sub_domains, 1)
bc2 = DirichletBC(zero, mesh, boundary)

# Assemble matrix
A = assemble(a, mesh)

# Output file
out_file = File("temperature.pvd")
vel_file = File("velocity.pvd")

vel_file << velocity

# Time-stepping
while t < T:
    # Copy solution from previous interval
    u0.assign(u1)
    # Assemble vector and apply boundary conditions
    b = assemble(L, mesh)
    bc.apply(A, b, a)
    bc2.apply(A, b, a)
    # Solve the linear system
    solve(A, x1, b)
    # Save the solution to file
    out_file << u1
    # Move to next interval
    t += k