Not getting the desired coordinates of the submesh

General
1

I am creating a sub mesh by marking the subdomains. Following code I am using for this

mesh = UnitSquareMesh(2, 2)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] < DOLFIN_EPS) and on_boundary
   
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
left = Left()
sub_domains.set_all(0)
left.mark(sub_domains, 1)
sub_mesh = SubMesh(mesh, sub_domains, 1)
# coordinates of the sub mesh
sub_mesh.coordinates()

I am getting the following coordinates (0,0), (0,0.5), (0.5, 0.5) while in the original mesh, this subdomain has the mesh coordinates (0,0), (0,0.5),(0,1) .
Why is this happening? Is there any problem with my code?

2

SubMesh is not made for meshes of co-dimension 1 (i.e. facets).
You can check this by print(submesh.topology().dim()).
To create a submesh of facets, either use BoundaryMesh or MeshView.

Personally, I would advice you to move to DOLFINx, as the support for these kind of operations have been vastly improved (on main branch).

3

Thanks, Now it is woking fine. Can we solve an ordinary differential equation on this boundary mesh?

4

There are several example of this.
In legacy fenics you can use pointintegral solver. In dolfinx I would go a slightly different route, and use external software such as gotranx by @Henrik_Nicolay_Finsb (GitHub - finsberg/gotranx: Next generation ODE translator) and use it with external operator by @Andrash GitHub - a-latyshev/dolfinx-external-operator: Extension of DOLFINx implementing the concept of external operator

1 Like