uniform refinement.py

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

from dolfin import *
from Rjump import *
from math import *
#create mesh, define the finite element, functions and boundary condition
mesh = Mesh('Square2.xml')
element = FiniteElement("Lagrange", "triangle", 1)
element2 = FiniteElement("Lagrange", "triangle", 2)
celement = FiniteElement("Discontinuous Lagrange", "triangle", 0)

class Source(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        values[0] = 4.0/DOLFIN_PI*200*200*(1-200*(x[0]*x[0]+x[1]*x[1]))* exp(-200*(x[0]*x[0]+x[1]*x[1]))


class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(on_boundary )

#create functions for solving the poisson equation and controlling the error
def SolvingPoisson():
    v = TestFunction(element)
    u = TrialFunction(element)
    a = dot(grad(v), grad(u))*dx
    L = v*f*dx
    u0 = Function(mesh, 0.0)
    boundary = DirichletBoundary()
    bc = DirichletBC(u0, mesh, boundary)
    A = assemble(a, mesh)
    b = assemble(L, mesh)
    bc.apply(A, b, a)
    x = Vector()
    solve(A, x, b)
    u = Function(element, mesh, x)
    return u


def ControllingError():
    Rju = rjump(u, mesh)
    R = Rju + modulus(f)
    vc = TestFunction(celement)
    h = MeshSize("triangle", mesh)
    L = (h*R)*(h*R)*vc*dx
    N_eix = assemble(L, mesh)
    sum=0
    N_eimf = MeshFunction("real", mesh, mesh.topology().dim())
    for c in cells(mesh):
        N_eimf.set(c, N_eix[c.index()])
        sum+=N_eix[c.index()]
    file = File("Nonono.pvd")
    file << N_eimf
    print sqrt(sum)
    return sqrt(sum)

#solve the problem
tol = 15
R_norm = 1000
while 1:
    f = Source(element2, mesh)
    u = SolvingPoisson ()
    R_norm = ControllingError ()
    
    if R_norm < tol:
        file = File("Nonono_uniform.pvd")
        file << u
        break
    
    mesh.refine()

    
print 'R_norm is equal to ',R_norm
print 'Number of vertices',mesh.numVertices()