Get function values on boundaries

Question

Hello,

I'm sure there is a simple solution to this question, but after some research and trying I can't figure out on my own how to do this...
As the title says, I want to evaluate the solution u of a PDE on (different) boundaries of my domain.
I defined something like

class Right(SubDomain):
def inside(self, x, on_boundary):
    tol = 1E-13   # tolerance for coordinate comparisons
    return on_boundary and  abs(x[0] - x0_max) < tol

class Top(SubDomain):
def inside(self, x, on_boundary):
    tol = 1E-13   # tolerance for coordinate comparisons
    return on_boundary and  abs(x[1] - x1_max) < tol

and would now like to evaluate my function (or get the maximum of it) on each of the boundaries, i.e. get the values of my function on the Top SubDomain and on the Right SubDomain.
I tried to use a MeshFunction (subdomains = MeshFunction('size_t', mesh, 2)) and mark the corresponding parts of the boundaries and iterate through the cells:

for cell_no in range(len(subdomains.array())):
  subdomain_no = subdomains.array()[cell_no]
  if (subdomain_no == 2):
        print u.vector()[cell_no]

resulting in an error:

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run

[0]PETSC ERROR: to get more information on the crash.

MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 59.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.

You may or may not see output from other processes, depending on exactly when Open MPI kills them.

Going through all coordinates of my mesh and test whether they are on the boundary would also be possible, but isn't there are simpler way to do this, using MeshFunction or FaceFuntion?

Answer 1

Here is one way to get boundary values

from dolfin import *
mesh = UnitSquareMesh(4, 4)

Right = AutoSubDomain(lambda x, on_bnd: near(x[0], 1) and on_bnd)
Top   = AutoSubDomain(lambda x, on_bnd: near(x[1], 1) and on_bnd)
# Mark a CG1 Function with different values on the two boundaries
V = FunctionSpace(mesh, 'CG', 1)
bc0 = DirichletBC(V, 1, Right)
bc1 = DirichletBC(V, 2, Top)
u = Function(V)
bc0.apply(u.vector())
bc1.apply(u.vector())

# Test for some Function v in function space V
v = interpolate(Expression("x[0]"), V)
print "Values on right = ", v.vector()[u.vector() == 1].array()
print "Values on top   = ", v.vector()[u.vector() == 2].array()

Results in:

Values on right =  [ 1.  1.  1.  1.]
Values on top   =  [ 0.    0.25  0.5   0.75  1.  ]

Comments

Thanks a lot, it works, although it seems a little bit... strange ;-) But it does what I want it to do, so I'm happy!

---

Note that the corner dof is marked by 2 with this procedure and as such you have only 4 ones, but 5 values on top. You should actually get the boundary dofs like

from numpy import where
bc0.apply(u.vector())
right_dofs = where(u.vector()==1)
bc1.apply(u.vector())
top_dofs = where(u.vector()==2)

Answer 2

Following the answer by Mikael, I changed the solution to be a bit more systematic:

from dolfin import *
from numpy import array, where

def top(x, on_boundary):
    return x[1] > 1 - DOLFIN_EPS and on_boundary

def subdomain_dofs(space, subdomain, dof):
    # Define a vector u that corresponds to boundary dofs
    # 0 -> not in the subdomain
    # 1 -> in the subdomain, 0th dimension
    # 2 -> in the subdomain, 1th dimension etc.
    bc0 = DirichletBC(space, range(1,space.num_sub_spaces()+1), subdomain)
    u = Function(space)
    bc0.apply(u.vector())

    # Create another vector and apply gather so that
    # the function works in parallel
    u_mpi =  PETScVector(mpi_comm_self())
    u.vector().gather(u_mpi, array(range(u_space.dim()), "intc"))

    # Return the dofs that correspond to input dimension
    return where(u_mpi==dof+1)[0].tolist()

mesh = UnitSquareMesh(10,10)
u_space = VectorFunctionSpace(mesh, "Lagrange", 1, 2)

u_x_top_dofs = subdomain_dofs(u_space, AutoSubDomain(top), 0)
u_y_top_dofs = subdomain_dofs(u_space, AutoSubDomain(top), 1)

print("x dofs for top boundary: "+ str(u_x_top_dofs))
print("y dofs for top boundary: "+ str(u_y_top_dofs))