Extract the nodal values of a function on the boundary of the mesh

General
1

Hello,
I am using legacy Fenics. Given a function defined on a mesh, I would like to extract its values only on the points which belong to the mesh boundary. I have looked online for this, all I could find is either

  • answers how to define a class which defines a boundary which is then used in DirichletBC to impose the BC on that boundary.
  • answers for Dolfinx, see also this
  • Answers which say how to get the function value on a point

Here is what I would like:

f = Function( Q )
#this gives the vector of nodal values
v = f.vector().get_local()

How may I prune from v the values which lie on Q.mesh()?

Thank you

2

See for instance: Question on Dolfin poisson equation for source measurement on borders - #26 by dokken

3

Thank you. This solution works only for functions in space FunctionSpace(mesh, 'CG', 1), not for functions in CG2:

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(4, 4)

V = FunctionSpace(mesh, 'CG', 2)
u = Function(V)
u.interpolate(Expression('x[0]+ 2*x[1]', degree=1))
vertex_function = MeshFunction("size_t", mesh, 0)

class BoundaryMarker(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary


BoundaryMarker().mark(vertex_function, 1)
boundary_vertices = np.asarray(vertex_function.where_equal(1))

v_to_d = vertex_to_dof_map(V)
dofs = v_to_d[boundary_vertices]

x = V.tabulate_dof_coordinates()
for dof in dofs:
    print(x[dof], u.vector()[dof])

$python3 solve.py
Traceback (most recent call last):

File "solve.py", line 21, in <module>

v_to_d = vertex_to_dof_map(V)

RuntimeError:
[...]
*** Error: Unable to tabulate dof to vertex map.
*** Reason: Can only tabulate dofs on vertices.
*** Where: This error was encountered inside DofMap.cpp.
[...]

Do you know how to make this work also for functions in other spaces ? I found this answer but then I don’t know how to do the equivalent of vertex_to_dof_map with the approach suggested in there.

Thank you

4

Following your link, I found a solution by using a dummy space of polynomials with one degree, here it is. This function returns the coordinates of the boundary points of a mesh, then the nodal values of a function (even if this function is defined on a polynomial space with degree larger than one) can be obtained by evaluating the function at these points.

def boundary_points(mesh):
    # create a dummy function space of degree 1 which will be used only to extract the boundary points
    Q_dummy = FunctionSpace( mesh, 'CG', 1 )

    vertex_to_degree_of_freedom_map = vertex_to_dof_map( Q_dummy )

    vertex_function = MeshFunction( "size_t", mesh, 0 )

    vertex_function.set_all( 0 )
    BoundaryMarker().mark( vertex_function, 1 )

    boundary_vertices = np.asarray( vertex_function.where_equal( 1 ) )

    degrees_of_freedom = vertex_to_degree_of_freedom_map[boundary_vertices]

    x = Q_dummy.tabulate_dof_coordinates()
    x = x[degrees_of_freedom]
    
    return x