3dstreaming.py

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

from dolfin import *
#create the mesh, finite elements and define the boundary conditions
mesh = Mesh("Dolfin-2.xml.gz")
scalar = FiniteElement("Lagrange", "triangle", 1)
vector = VectorElement("Lagrange", "triangle", 2)
sub_domains = MeshFunction("uint", mesh, "Subdomains.xml.gz");
velocity = Function(vector, "Velocity.xml.gz");
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)

    
#define the function and set the boundary condition
T = 5.0
k = 0.1
t = k
c = 0.0001
f = Function(scalar, mesh, 0.0)
u0 = Function(scalar, mesh, 0.0)
v = TestFunction(scalar)
u = TrialFunction(scalar)
x0 = Vector()
x1 = Vector()
u0 = Function(scalar, mesh, x0)
u1 = Function(scalar, mesh, x1)
h = 2.0 * MeshSize("triangle", mesh)
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
delta = h
boundary = DirichletBoundary()
g = Function(mesh, 1.0)
zero = Function(mesh, 0.0)
bc = DirichletBC(g, sub_domains, 1)
bc2 = DirichletBC(zero, mesh, boundary)
#solve the function
A = assemble(a, mesh)
out_file = File("Nonono_Te.pvd")
vel_file = File("Nonono_ve.pvd")
vel_file << velocity
while t < T:
    u0.assign(u1)
    b = assemble(L, mesh)
    bc.apply(A, b, a)
    bc2.apply(A, b, a)
    solve(A, x1, b)
    out_file << u1
    t += k