Hi everyone,

I’m working on a modified cahn-hilliard demo, where I want to apply a given concentration as a Dirichlet boundary condition but only on a single node. On the documentation of DirichletBC says the for pointloads one should use the “pointwise” method, but somehow it won’t be applied. I marked my mesh to see what the problem was but the point is not being marked. I think I’m using the SubDomain class wrong because my InitialCondition class (UserExpression subclass) does assign the initial conditions at the node the way I want to. It’d would be really great if someone could help me find out what’s wrong.

Cheers

My working code:

```
from fenics import *
class InitialConditions(UserExpression):
def eval(self, values, x):
if near(x[0], lenght, DOLFIN_EPS) and near(x[1], height/2, DOLFIN_EPS):
values[0] = 1.0
else:
values[0] = 0.0
values[1] = 0.0
def value_shape(self):
return (2,)
class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L, bcs):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bcs = bcs
def F(self, b, x):
assemble(self.L, tensor=b)
for bc in self.bcs:
bc.apply(b, x)
def J(self, A, x):
assemble(self.a, tensor=A)
for bc in self.bcs:
bc.apply(A)
class Gamma(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], lenght, DOLFIN_EPS) and near(x[1], height/2, DOLFIN_EPS)
# Output file
file = File("Results_2D/output.pvd", "compressed")
# Model parameters
M = 1.0e00
lmbda = 1.0e02 # surface parameter
dt = 5.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
# Geometry parameters and mesh parameters
lenght = 6. # [mm]
height = 30. # [mm]
Nelem = 100 # [-]
# Form compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# Create mesh and build function space
mesh = RectangleMesh(Point(0., 0.), Point(lenght, height), Nelem, Nelem*5, diagonal='right')
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
boundary = Gamma()
boundary.mark(sub_domains, 1)
File("subdomains2D.pvd") << sub_domains
# 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
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)
File("Results/InitCond2D.pvd", "compressed") << u0.split()[0]
# Compute the chemical potential df/dc
c = variable(c)
f = 100*c**2*(1-c)**2
dfdc = diff(f, 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*M*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# Boundary conditions
bc = DirichletBC(ME.sub(0), Constant(1.0), Gamma(), method='pointwise')
bcs = [bc]
# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L, bcs)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6
# Step in time
t = 0.0
Nsteps = 50
T = Nsteps*dt
file << (u.split()[0], t)
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u.split()[0], t)
```