Extract coordinates of the nodes of a subboundary

The boundary of my mesh is divided into several parts with different physical meaning. I implemented this via gmsh and

facets = df.MeshFunction('size_t',mesh,'my-mesh_facet_region.xml')

I can extract the coordinates of the whole boundary via

boundary_coordinates = facets.mesh().coordinates()

But I need the coordinates of a part of the boundary only, subdomain 1, say, and it is difficult in this case to sort these out afterwards. So, how can I get the coordinates of just the subboundary directly?

By the way, I also would like to evaluate a function u exactly at these nodes. So, how can I relate the entries of u.array()[:] with those boundary nodes?

Thank you!

Consider this minimal example:

from dolfin import *
import numpy as np

class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < 0.4 and on_boundary


mesh = UnitSquareMesh(20,20)
V = FunctionSpace(mesh, "CG", 1)
v = project(Expression("x[0]+x[1]", degree=2),V)
mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundary().mark(mf, 1)
v2d = vertex_to_dof_map(V)
dofs = []
for facet in facets(mesh):
    if mf[facet.index()] == 1:
        vertices = facet.entities(0)
        for vertex in vertices:
            dofs.append(v2d[vertex])

unique_dofs = np.array(list(set(dofs)), dtype=np.int32)
boundary_coords = V.tabulate_dof_coordinates()[unique_dofs]
for i, dof in enumerate(unique_dofs):
    print(boundary_coords[i], v.vector()[dof])
2 Likes

Thank you, dokken! This essentially solves my problem, indeed.

However, my function v is from

V2 = FunctionSpace(mesh, "CG", 2)

so

v.vector()[dof]

will not give a correct result in this case. Unfortunately,

vertex_to_dof_map(V2)

yields an error: Can only tabulate dofs on vertices.

That is not an error, that is simply from the fact that a CG-2 space has degrees of freedom on facets as well as vertices. Therefore you cannot use vertex_to_dof_map for this operation. You need a more complex operation for this:

from dolfin import *
import numpy as np

class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5 and on_boundary


mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh, "CG", 2)
v = project(Expression("x[0]+x[1]", degree=2),V)
mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundary().mark(mf, 1)

mesh.init(2,1)
dofs = []
cell_to_facets = mesh.topology()(2,1)
for cell in cells(mesh):
    facets = cell_to_facets(cell.index())
    for facet in facets:
        if mf[facet] == 1:
            dofs_ = V.dofmap().entity_closure_dofs(mesh, 1, [facet])
            for dof in dofs_:
                dofs.append(dof)
1 Like

Excellent! Thank you.

The MeshTopology() class is not adequately documented, and Its not clear what the commands mesh.init(2,1) and mesh.topology()(2,1) do. How we can get more information to use in different applicaions?

This function creates the connectivity between entities of topological dimension 2 (triangles) to entities of topological dimension 1 (lines) using indices local to a process. Each triangle will have three line segments connected, i.e.

gives you a map from each cell (local to process) to the three line segments (index local to process).

Thanks a lot for your fast answer!