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

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() )
    bottom_label = 2
    bottom_bd = bottom()
    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!


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

1 Answer

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 
Please login to add an answer/comment or follow this question.

Similar posts:
Search »