adaptivemesh.py

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

from dolfin import *
from Rjump import *
from math import *
#############################################################################
##construct mesh
#############################################################################

# Create mesh and finite element
mesh = Mesh('Square2.xml')
element = FiniteElement("Lagrange", "triangle", 1)
element2 = FiniteElement("Lagrange", "triangle", 2)
celement = FiniteElement("Discontinuous Lagrange", "triangle", 0)

#############################################################################
##defining functions
#############################################################################
# Source term
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]))

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(on_boundary )

#############################################################################
##solving poisson equation
############################################################################

def poisson_solve():
    # Define variational problem
    v = TestFunction(element)
    u = TrialFunction(element)
    a = dot(grad(v), grad(u))*dx
    L = v*f*dx
    # Define boundary condition
    u0 = Function(mesh, 0.0)
    boundary = DirichletBoundary()
    bc = DirichletBC(u0, mesh, boundary)
    # Assemble linear system and solve using a linear solver
    A = assemble(a, mesh)
    b = assemble(L, mesh)
    bc.apply(A, b, a)
    x = Vector()
    solve(A, x, b)
    # Construct solution function
    u = Function(element, mesh, x)
    return u

########################################################################
#error control
########################################################################
def residual_calculation():
    # construct residual
    Rj = rjump(u, mesh)
    R = Rj + modulus(f)
    #construct integral form
    vc = TestFunction(celement)
    h = MeshSize("triangle", mesh)
    Lei = (h*R)*(h*R)*vc*dx
    #assemble the integral to vectors
    eix = assemble(Lei, mesh)
    #save error indicator
    sum=0
    eimf = MeshFunction("real", mesh, mesh.topology().dim())
    for c in cells(mesh):
        eimf.set(c, eix[c.index()])
        sum+=eix[c.index()]

    #file = File("ei.pvd")
    #file << eimf
    #print sqrt(sum)
    return sqrt(sum)


#############################################################################
##iteratively refine the mesh, until the error criterion satisfies
#############################################################################
#define tolerence and initial norm of R
tol = 20
R_norm = 1000

#############################################################################
##uniform refinement
#############################################################################

while 1:
    f = Source(element2, mesh)
    u=poisson_solve()
    R_norm = residual_calculation()
    if R_norm < tol:
        file = File("poisson_uniform.pvd")
        file << u
        break

    mesh.refine()
    
print 'R_norm is ',R_norm
print 'Number of vertices',mesh.numVertices()
    
print 'Uniform refinement finished!'

#R_norm is  6.40633004739
#Number of vertices 549889
#Uniform refinement finished!


#R_norm is  12.80647219
#Number of vertices 137729
#Uniform refinement finished!


############################################################################
##adaptive refinement
############################################################################
R_norm = 1000
percentage = 0.3
mesh = Mesh('Square2.xml')

while 1:
    f = Source(element2, mesh)
    u=poisson_solve()
    # construct residual
    Rj = rjump(u, mesh)
    R = Rj + modulus(f)
    
    #construct integral form
    vc = TestFunction(celement)
    h = MeshSize("triangle", mesh)
    Lei = (h*R)*(h*R)*vc*dx
       
    #assemble the integral to vectors
    eix = assemble(Lei, mesh)
    
    #load the error vector into a list and sort the list
    temp=[]
    sum=0
    #   eimf = MeshFunction("real", mesh, mesh.topology().dim())
    for c in cells(mesh):
    #     eimf.set(c, eix[c.index()])
        sum+=eix[c.index()]
        temp.append(eix[c.index()])
        
    if sqrt(sum)<tol:
        file = File("poisson_adaptive.pvd")
        file << u
        break
    
    temp.sort()
    ind=temp[int(floor(len(temp)*(1-percentage)))]
    
    cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
    cell_markers.fill(False)
    
    for c in cells(mesh):
        if(eix[c.index()]>=ind):
            cell_markers.set(c, True)
    
    mesh.refine(cell_markers)


print 'R_norm is ',sqrt(sum)
print 'Number of vertices',mesh.numVertices()
print 'adaptive refinement finished!'
#R_norm is  9.99428466847
#Number of vertices 14678
#adaptive refinement finished!


#R_norm is  18.0225840498
#Number of vertices 4616
#adaptive refinement finished!