Cahn Hillard boundary condition multiple boundary conditions?

General
1

I want to implement different boundary condition (in this case nueman boundary condition with different value on each boundary) .I have used multiple subdomain method but it is not working properly despite the code is running successfully

python
import random
from dolfin import *

class Left(SubDomain):
    def inside(self, x, on_boundary):
        tol= 1E-10
        return near(x[0], -200.0,tol)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        tol= 1E-10
        return near(x[0], 200.0,tol)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        tol= 1E-10
        return near(x[1], -200.0,tol)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        tol= 1E-10
        return near(x[1], 200.0,tol)


left = Left()
top = Top()
right = Right()
bottom = Bottom()




        
# model parameters

g_L=Constant(-0.5)
g_R=Constant(0.5)
g_B=Constant(0.5)
g_T=Constant(-0.6)
lmbda,L,H,r = 2,400.0,400.0,40.0
k,M,theta = 1.29,10,0.5
a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10 = 0,0,8.072789087,-81.24549382,408.0297321,-1244.129167,2444.046270,-3120.635139,2506.663551,-1151.003178,230.2006355
w = 0.1
dt = 0.5

mesh = RectangleMesh(Point(-L/2,-H/2),Point(L/2,H/2),200,200)
facet = MeshFunction("bool", mesh, 0)
facet.set_all(False)       

class Inside(SubDomain):
       def inside(self,x,on_boundary):
              return near(sqrt(pow((x[0]-0),2)+pow((x[1]-0),2)),40,5)#sqrt(pow((x[0]-400/2),2)+pow((x[1]-400/2),2)) <=(20+10)
Inside().mark(facet, True)
mesh = refine(mesh, facet)




boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)

left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)


# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)
        
    
tol1 = 1e-14
c_i = Expression(("sqrt(pow((x[0]-0),2)+pow((x[1]-0),2)) <= (r+tol1) ? 1 : 0.0065","0"),degree =0,L=L,H=H,r=r,tol1 = tol1)

# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True



                 

P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)

# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

# Create intial conditions and interpolate
c_init= c_i
u.interpolate(c_init)
u0.interpolate(c_init)

# Compute the chemical potential df/dc
c = variable(c)
f1 = w*(a0 +a1*c+a2*c**2+a3*c**3+a4*c**4+a5*c**5+a6*c**6+a7*c**7+a8*c**8+a9*c**9+a10*c**10)
dfdc = diff(f1, c)

# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu

#

# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(M*grad(mu_mid), grad(q))*dx +g_L*q*ds(1) + g_T*q*ds(2) + g_R*q*ds(3) + g_B*q*ds(4)
L1 = mu*v*dx - dfdc*v*dx - dot(lmbda*grad(c), grad(v))*dx
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
file = File("sol10/pfpr10.pvd", "compressed")

# Step in time
t = 0.0
T = 50*dt
while (t < T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    file << (u.split()[0], t)