Apply a DirichletBC to a single node

mesh
1

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)
1 Like
2

on_boundary will always be false if you use method="pointwise". If you define Gamma as

class Gamma(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], lenght, DOLFIN_EPS) and near(x[1], height/2, DOLFIN_EPS)

it will work

1 Like
3

Hi,

if i understand your question correctly you want to apply a normal Direchlet bc to your mesh. As far as i can tell you don’t have to specifically create a class for that, but it should not matter.

def cb(x, on_boundary):
     return on_boundary and near(x[0], lenght, tol) and near(x[1], height/2, tol)

bc = DirichletBC(ME.sub(0), Constant(1.0), cb)
bcs = [bc]

The only problem i can really think of is that the tolerance is to small and it never “finds” the point you want to apply the boundary condition to.
I sadly don’t have any knowledge on the problem you are solving, but some time ago i had problems, that were solved by applying a Neumann condition.

Best regards

4

Thanks a lot! It seems that was the problem, as now the results seem realistic but I noticed on my sub_domains output that the point is still not being marked. Why would be the reason there?

5

That would be because according to your definition of sub_domains the mesh function will contain edges, i.e. the entity one dimension less than the cells of your 2d mesh.
And since your Gamma().inside() is only satisfied for one point of the mesh, no edge will be marked, because two points would be needed for that.

6

How could I mark the point? Following the reasoning I tried sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 2) . While I don’t get any errors while executing my code, I get a visualization error on Paraview due to Cell Connectivity.

7

What would be your expected output of such a visualization?

8

Just the discrete nodes of the mesh with the respective integer marks.

9

I can’t think of a direct solution. However, maybe the following hack will suffice for your needs (I didn’t want to use your code since your mesh has so many points. I’m sure you can adapt this.).

from dolfin import * 

mesh = UnitSquareMesh(4,4)
V = FunctionSpace(mesh, FiniteElement("Lagrange", mesh.ufl_cell(), 1))

class PointOnBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0.5) and near(x[1], 0.)

subd = MeshFunction("size_t", mesh, 0)

ptb = PointOnBoundary()
ptb.mark(subd,1)


u = Function(V)
dtv = dof_to_vertex_map(V)
u.vector().set_local(subd.array()[dtv])
File("points.pvd") << u

Essentially we define a function on a linear nodal function space and directly set the dof values to the values of the mesh function.

For the visualization I used paraview, chose the representation “points”, set the point size to 15 and checked the box “render point as spheres” in the properties tab. The result looks like this:

meshfun2
meshfun2.png603×745 19.5 KB

(I loaded the file twice and displayed the second dataset as a wireframe with solid color)

1 Like
10

That is exactly what I meant. Thanks a lot for your help!

11

Klunkean gave you the solution to the problem you posed. I’d like suggest, however, that perhaps the problem you posed isn’t the problem you want to solve. For a 2D problem, the value of an H^1 function at a point isn’t well defined. The practical implications of this are that: (a) your problem may not converge with mesh refinement, because the problem defined at the continuous level (with a point Dirichlet BC) is undefined. And (b) the stiffness matrix will have a pretty high condition number. As an alternative, you might want to replace your point condition with an integral condition. An example is given in the neumann-poisson demo. The discrete problem will be better behaved because it’s a discretization of a well-posed continuous problem. Good luck!

12

Precisely. You can’t solve what is unsolvable.