Nonlinear problem: Fix vector to a specific value at a given mesh vertex

Question

I'm currently trying to solve a three-dimensional non-linear eigenvalue problem, for which I need to restrict my seeked function at a mesh vertex to a given value. I typically use a linear eigenvalue problem to calculate an initial guess for the nonlinear problem, which I then solve via the solve function. As a toy model in order to illustrate my problem let's consider the simple 1D Helmholtz eigenvalue problem

$$( \nabla^2 + k_n^2) u_n = 0 $$

with Neumann boundary conditions and the further constraint that the function has a particular value at a given mesh vertex $v_i$ , i.e.

$$ u_n(x) |_{x=v_i}. $$

My particular problem is that I do not know how to express the latter constraint in an ufl form.

My current non-working code looks like this:

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

parameters["linear_algebra_backend"] = "uBLAS"
parameters["reorder_dofs_serial"] = False

mesh = UnitIntervalMesh(50)

INDEX_FIXED = 21 # should match the index of a vertex in the range of mesh.coordinates()

## Define domain for point, where we want to enfore a paritcular value for the function
indexDom = MeshFunction('size_t', mesh, 0)
indexDom.set_all(0)
indexDom.set_value(INDEX_FIXED, 1)

V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

#Initial guess
u_mixed = interpolate(Expression(("cos(k*x[0])","k"), k=3*pi * 0.9), W)

duw = TrialFunction(W)
uw, kw = split(u_mixed)
uv, kv = TestFunctions(W)

# Extract current value from function
value = Constant(u_mixed.sub(0)(mesh.coordinates()[INDEX_FIXED]))

F = (inner(nabla_grad(uw), nabla_grad(uv)) - uw * uv * kw**2) * dx
# Add additional constraint to enforce the function value at the given index
# (does NOT work)
F += (uw - value) * kv * dP(1, domain_data=indexDom)

J = derivative(F, u_mixed, duw)

solve(F == 0, u_mixed, J=J) # fails accordingly

The current attempt to describe the constraint effectively does nothing and yields an empty last row in the Jacobian matrix, which explains why the Newton solver fails.

# Last row of J is empty
Jmat = uBLASSparseMatrix()
assemble(J, tensor=Jmat)
plt.spy(Jmat.array())
plt.show()

For the simple toy model, I can circumvent the issue, by setting the corresponding matrix entry by hand using a NonlinearProblem class. However, this is not elegant at all and furthermore not usable with PETSc and therefore not usable with MPI.

Thanks for any pointers in the right direction.

PS: In a prior question the problem was solved in 1D for the case that the particular vertex lies on the boundary, but this does not work in higher dimensions or for any other vertex.

PPS: So far I solved the issue by using a different constraint, namely that int u^2 dx = 1,
but this yields a full row in the Jacobian J, which can significantly slow down the solver.

Answer 1

Use the class DirichletBC. This class can be used for enforcing internal constraints too (not only boundary conditions). Make sure to use the "pointwise" method if the condition should be enforced at only one point.

Try something of the form

particular_value = ...
particular_point = "near(x[0], vi[0]) && near(x[1], vi[1])"
bc = DirichletBC(W.sub(0), particular_value, particular_point, "pointwise")
solve(F == 0, u_mixed, bcs=bc, J=J)

Comments

Thank you for the quick answer!

However, the DirichletBC class does not exactly do what I want in this case. The additional constraint in the mixed problem above is needed due to the additional dof. For this I need the constraint to be associated with the additional degree of freedom.

With Dirichlet boundary conditions, as you have suggested, I get the following nonzero pattern for the Jacobian J:

However, what I actually want is this:

I will probably just go forward and use the compute_bc_pointwise method in the Dirichlet class to write my own class which can set the appropriate entries in the matrix. However, I think that the syntax I used in the original question is actually supposed to do what I want, but the feature is simply not implemented. Shall I sent a bug report/feature request regarding this or does another method exist to achieve the same result that I have missed so far.

Thank you a lot!

---

Marie, I'm having difficulty implementing your approach at a very basic level. Sorry for the level of 'newbie' here.
If particular_point is just a string, how does FEniCS know what vi is? Also, what kind of object must particular_value be? Here is what I'm trying to work with:

mesh = UnitSquareMesh(10,10)
Vv = VectorFunctionSpace(mesh, "Lagrange", 1)
particular_value = Constant((0.,0.))
particular_point = "near(x[0], vi[0]) && near(x[1], vi[1])"
bc = DirichletBC(Vv, particular_value, particular_point, "pointwise")

---

You could replace vi in the string directly or using string replacement for instance.