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

378
views
2
6 months ago by

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

# 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
***
***
*** 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:
*** -------------------------------------------------------------------------

3
6 months ago by

You can't solve the linear system  $\textbf{A}u=b$Au=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.

Thanks! Marking as answered because this works, but if I find a more elegant solution, I will post it here.
written 6 months ago by Alexander Niewiarowski