Hello, everyone
I am trying to get \chi_1 in the problem below with Neumann Boundary Condition. And I wish the Neumann Boundary Condition is defined on the edges, where \chi_1[0] = -1, \chi_1[0] = 1, \chi_1[1] = -1, \chi_1[1] = 1.
So I wrote the following code:
from dolfin import *
from mshr import *
mesh = RectangleMesh.create([Point(-1,-1), Point(1, 1)],[10, 10],CellType.Type.quadrilateral)
tol = -1E-1
subdomains = MeshFunction("size_t", mesh, 2)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return True if x[1] <= tol else False
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return True if x[1] >= tol else False
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(subdomains, 0)
subdomain_1.mark(subdomains, 1)
V0 = FunctionSpace(mesh,"DG",0)
k = Function(V0)
kvalues = [2, 8] # values of k in the two subdomains
for cell in range(len(subdomains.array())):
subdomain_cell = subdomains.array()[cell]
k.vector()[cell] = kvalues[subdomain_cell]
class LeftBoundary(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0], -1)
class RightBoundary(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[0], 1)
class TopBoundary(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1], 1)
class BottomBoundary(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and near(x[1], -1)
left_bc = LeftBoundary()
right_bc = RightBoundary()
top_bc = TopBoundary()
bottom_bc = BottomBoundary()
Ve = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
Re = FiniteElement('Real', mesh.ufl_cell(), 0)
X = FunctionSpace(mesh, MixedElement([Ve, Re]))
dx = Measure('dx', subdomain_data=subdomains)
n = FacetNormal(mesh)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries)
bottom_bc.mark(boundaries,1)
right_bc.mark(boundaries,2)
top_bc.mark(boundaries,3)
left_bc.mark(boundaries,4)
u, m = TrialFunctions(X)
v, r = TestFunctions(X)
F1 = inner(nabla_grad(u) + Constant((1,0)), k*nabla_grad(v))*dx + u*r*dx + m*v*dx
- inner(Constant((1,0)),n('-'))*v*dS(1) - inner(Constant((1,0)),n('-'))*v*dS(2)
- inner(Constant((1,0)),n('-'))*v*dS(3) - inner(Constant((1,0)),n('-'))*v*dS(4)
a1, l1 = lhs(F1), rhs(F1)
x = Function(X)
solve(a1 == l1, x)
uh_x, mh_x = x.split()
File("uh_x.pvd") << uh_x
But the result seems like unreasonable. I am wondering what is the problem in this code. And I also tried to print the code below but got 0.
print(assemble(inner(Constant((1,0)),n('-'))*dS(1) + inner(Constant((1,0)),n('-'))*dS(2)
+ inner(Constant((1,0)),n('-'))*dS(3) + inner(Constant((1,0)),n('-'))*dS(4)
+ Constant(0)*dx(domain=mesh, subdomain_data=subdomains)))
Is there someone who might be able to help me with this? I would appreciate it.