Extract coordinates of the nodes of a subboundary

mesh
1

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!

2

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
3

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.

4

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
5

Excellent! Thank you.

6

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?

7

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).

8

Thanks a lot for your fast answer!