### Solving PDEs on mesh boundaries - Better way to calculate gradient along mesh boundary?

I am developing some code for fluid structure interaction using ALE methods. This requires the ability to solve a PDE along the boundary of a mesh - in other words, calculating the gradient along the boundary. Below, I am trying to solve a simple 2nd order PDE (d^2u/dx^2) = 0 along the bottom boundary of a unit square mesh (the x-axis). (The solution is linear). My initial approach was to subdivide the mesh into the "bottom" and the "the_rest". I place 2 point-wise Dirichlet boundary conditions on the bottom left and right corners [u(0,0)=0 and u(1,0) = 1], and set the remainder of the mesh (inside and boundary) = 0. Then I dot the gradient of the test and trial functions with the vector {1,0} to get the required directional derivative.

```
from fenics import *
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
if __name__ == "__main__":
ndiv = 50
mesh = UnitSquareMesh(ndiv,ndiv)
poly_order = 1
V = FunctionSpace(mesh, 'Lagrange', poly_order)
# Define subdomains
# the bottom of the unit square, want to solve a pde along this boundary
class bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs( x[1]) < DOLFIN_EPS
# The mesh boundary minus the bottom
class the_rest_bd(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and not abs( x[1]) < DOLFIN_EPS
# The inner mesh minus the bottom
class the_rest_inner(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] -1) > DOLFIN_EPS and abs(x[0]) > DOLFIN_EPS and \
abs(x[1] -1) > DOLFIN_EPS and abs(x[1]) > DOLFIN_EPS
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domain = MeshFunction("size_t", mesh, mesh.topology().dim() )
boundaries.set_all(0)
sub_domain.set_all(0)
bottom_label = 2
bottom_bd = bottom()
bottom_bd.mark(boundaries,bottom_label)
the_rest_label = 1
the_rest_bd = the_rest_bd()
the_rest_inner = the_rest_inner()
the_rest_bd.mark(boundaries, the_rest_label)
the_rest_inner.mark(sub_domain, the_rest_label)
ds = ds(subdomain_data=boundaries)
#%%
# Define the dirichlet boundary conditions (bottom two corners)
bc_left = DirichletBC(V, Constant(0), 'fabs(x[0]) < DOLFIN_EPS && fabs(x[1]) < DOLFIN_EPS', method = 'pointwise')
bc_right = DirichletBC(V, Constant(1), 'fabs(x[0]-1) < DOLFIN_EPS && fabs(x[1]) < DOLFIN_EPS', method = 'pointwise')
# the rest of the domain will be set to 0
bc_the_rest_bd = DirichletBC(V, Constant(0), the_rest_bd)
bc_the_rest_inner = DirichletBC(V, Constant(0), the_rest_inner)
bcs = [bc_left, bc_right, bc_the_rest_bd, bc_the_rest_inner]
# Define the trial and test function
u = TrialFunction(V)
v = TestFunction(V)
# Define the forcing function
f = Expression("0", degree=1)
# projection vector (dot with grad to get directional derivative)
n = Expression(('1','0' ), degree=1)
# Define the bilinear form
a = dot(dot(grad(u),n), dot(grad(v),n))*ds(2)
# Define the forcing
F = f*v*ds(2)
# Compute solution
uh = Function(V)
solve(a == F, uh, bcs)
```

However, this approach results in the following error when calling solve. Anyways, as I am relatively new to fenics and constantly find myself underestimating its capabilities, I am thinking there must be a much better/"cleaner" way of taking a gradient of a function along a boundary. I am running the latest fenics on spyder/anaconda python 3.6. Thanks!

```
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to set given (local) rows to identity matrix.
*** Reason: some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where: This error was encountered inside PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset:
*** -------------------------------------------------------------------------
```

### 1 Answer

You can't solve the linear system $\textbf{A}u=b$**A**`u`=`b` directly because the matrix in the LHS has zero values in its diagonal. You can fix this problem following the suggestion in the error message:

```
# Discrete solution
uh = Function(V)
# Assemble system
A = assemble(a, keep_diagonal=True)
b = assemble(F)
# Apply boundary conditions
[bc.apply(A, b) for bc in bcs]
# Solve system
solve(A, uh.vector(), b)
```

I don't know if this is the better approach but I would have used it too.