Obtain coordinates defined by mesh function

Question

Hello,

I use a Meshfunction to define various interfaces inside a mesh (labeled from 1 to 4) :

my_facets_meshfunction = MeshFunction("size_t", mesh, "my_facet_region.xml")

(where geometry in .xml file is defined and generated with gmsh)

I can use it to define boundary conditions, for example if label 1 corresponds to an exterior set of facets on which I want to impose a displacement :

V = FunctionSpace(mesh, "Lagrange", 1)
bc = DirichletBC(V, Constant(0.0), my_facets_meshfunction, 1)

Now, for some post-processing, I want to obtain the list of coordinates along a given set of my_facets_meshfunction (for instance defined with label 4, so I can follow the evolution along an interface), has someone an idea of how I can do it ?

Many thanks in advance !

Answer 1

Maybe use a combination of the SubsetIterator and MeshEntityIterator classes is a good option for you.

Comments

Thanks for your answer, it works with :

It_facet = SubsetIterator(exterior_facets_meshfunction,4)
for c in It_facet:
    print(c.midpoint().x())

I understand c is a MeshEntity, is there any difference with the object returned for example by:

It_mesh = vertices(mesh)

on which I can iterate with :

    for c in It_mesh:
        print(c.point().x())

In the first case, I use midpoint() and in the other point()... ?

---

In both cases the type of the object returned is the same, but with meanings completely different.

In the first case c.midpoint.x() is the x-coordinate of the midpoint of the edge (supposing that your entities have dimension 1), because you are iterating over facets. whereas in the second case you are iterating over all vertices of the mesh, hence, c.point().x() returns the x-coordinate of each vertex.

---

Hi, you could also use the information that the bc object has computed. In its construction it has already iterated over the facets and got the dofs. Now it only remains to "convert" them to vertices. Consider

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, 'CG', 1)
bc = DirichletBC(V, Constant(0), 'near(x[0], 0) || near(x[0], 1)')

gdim = mesh.geometry().dim()
d2v = dof_to_vertex_map(V)
for dof in bc.get_boundary_values():
    vertex = Vertex(mesh, d2v[dof])
    print [vertex.x(i) for i in range(gdim)]

---

Hi,

I have two remarks/questions on your comments :

Using SubsetIterator I can't access to c.point(), what I get when I try is :

'MeshEntity' object has no attribute 'point'

Is there a reason for this ?

Is it possible to iterate directly over the dofs of a meshfunction's subset (corresponding to a given label), without constructing the bc object ?

Many thanks for your help !

---

The error message occurs because c is not a vertex (is a MeshEntitty of geometrical dimension equal to 1). Try this:

It_facet = SubsetIterator(exterior_facets_meshfunction,4)
for f in It_facet:
  for v in vertices(f):
      print v.point().x(), v.point().y()

(maybe an statement for not repeat the output of some vertices is necesary in the loop)