Help for the conflict results from Dolfin assembling

Question

Hi all,

I am using Dolfin to calculate test function's assembling boundary integration result. Because of further programming, I planned to build a class to store it. However, when I debug my class and the assembling result, which is "test_mesh.fgamma()" in the code below, with the result directly got from dolfin, which is "f_gamma", they are not the same.

I am new to python and programming. So appreciated to your help and suggestions.

import numpy as np
import math
from dolfin import *
class SquareShape:
## Build the fluid dynamic model
# Initialization
def __init__(self,nx,ny):
    self.nx = nx
    self.ny = ny
    self.Ne = Constant(1/0.6)
    self.mesh = UnitSquareMesh(nx, ny)
    # Define function spaces
    self.T = FunctionSpace(self.mesh, "CG", 1)
    self.te = TrialFunction(self.T)
    self.de = TestFunction(self.T)
def fgamma(self):
    def heater_boundary(x,on_boundary): # also is the inlet boundary
        return on_boundary and near(0, x[0]) and (x[1] > 0.5)
    boundary_parts = MeshFunction("uint", self.mesh, self.mesh.topology().dim()-1)
    class heater_domain(SubDomain):
        def inside(self,x,on_boundary):
            return on_boundary and heater_boundary(x,on_boundary)
    Gamma_h = heater_domain()
    Gamma_h.mark(boundary_parts,0)
    f_gamma_express = self.Ne*self.de*ds(0)
    f_gamma_vector = assemble(f_gamma_express, exterior_facet_domains=boundary_parts)
    return (f_gamma_vector.array())

nx = 3
ny = 3
mesh = UnitSquareMesh(nx, ny)
Lend = 0
Rend = 1
T = FunctionSpace(mesh, "CG", 1)
# Create mesh function over special facets
boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)
# Define boundary domains
def heater_boundary(x,on_boundary): return on_boundary and near(Lend, x[0]) and (x[1] > 0.5)
class heater_domain(SubDomain):
         def inside(self,x,on_boundary): return on_boundary and heater_boundary(x,on_boundary)
Gamma_h = heater_domain()
Gamma_h.mark(boundary_parts,0)
Ne = Constant(1/0.6)
te = TrialFunction(T)
de = TestFunction(T)
f_gamma_express = Ne*de*ds(0)
f_gamma_vector = assemble(f_gamma_express, exterior_facet_domains=boundary_parts)
f_gamma = f_gamma_vector.array()
# compare two results
test_mesh = SquareShape(nx,ny)
print np.abs(f_gamma - test_mesh.fgamma()).max()