galerkin.py

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

from dolfin import *
#create the mesh, finite elements and define the boundary conditions
mesh = Mesh("Square.xml")
mesh.refine()
element = FiniteElement("Lagrange", "triangle", 1)
element2 = FiniteElement("Lagrange", "triangle", 2)
class Source(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
	values[0] = 0.0

	
class InitialValue(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
	values[0] = 1.0

	
class Flux(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        values[0] = 0.0

        
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(on_boundary)
    
#define the function and set the boundary condition
v = TestFunction(element)
U = TrialFunction(element)
f = Source(element2, mesh)
g = Flux(element, mesh)
c = 0.01
T = 0.5
k = 0.01
t = k
f = Source(element2, mesh)
u0 = InitialValue(element, mesh)
v = TestFunction(element)
u = TrialFunction(element)
x0 = Vector()
x1 = Vector()
nonono = Function(element, mesh, x1)
h = MeshSize("triangle", mesh)
c = 1.0
a = (u*v+k*dot(grad(u),grad(v)))*dx
L = (u0*v+f*v)*dx
boundary = DirichletBoundary()
zero = Function(mesh, 0.0)
bc = DirichletBC(zero, mesh, boundary)
#solve the function
A = assemble(a, mesh)
out_file = File("Nonono_T.pvd")
while t < T:
    b = assemble(L, mesh)
    bc.apply(A, b, a)
    solve(A, x1, b)
    out_file << nonono
    t += k
    u0.assign(nonono)