Mark nodes from "facet_region.xml"

Question

I have my geometry from Gmsh.
I have mark all phisycal lines in Gmsh and I could export all meshs from fenics.

But I couldnt mark nodes from Gmsh.

Can I mark nodes from xml facet?

mesh = Mesh("Square4.xml")
subdomains = MeshFunction("size_t", mesh, "Square4_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "Square4_facet_region.xml")

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

I would like mark nodes from boundaries.

Anyone have any ideia?

Answer 1

You can iterate over marked facets and then indentify the nodes in each facet or you can apply a boundary condition on the marked facets and then map the dofs to the node indices.

Try something like:

import numpy as np
from  dolfin import *

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

boundaries = FacetFunction("size_t", mesh)
boundary().mark(boundaries, 1)

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)

## First approach
visited = []
nodes1 = []
for f in facets(mesh):
  if boundaries[f.index()] == 1:
    for v in f.entities(0):
      if v not in visited:
        nodes1.append(int(v))
        visited.append(int(v))

## Second approach    
# Identify boundary dofs
u = Function(V)
bc = DirichletBC(V, Constant(1.0), boundaries, 1)
bc.apply(u.vector())
dofs_bc = list(np.where(u.vector()[:] == 1.0))

# Identify boundary nodes
d2v = dof_to_vertex_map(V)
nodes2 = d2v[dofs_bc]

## Testing
test1 = VertexFunction("size_t", mesh)
test1.array()[nodes1] = 1
test2 = VertexFunction("size_t", mesh)
test2.array()[nodes2] = 1

plot(test1, title="1")
plot(test2, title="2")
interactive()

Comments

Thanks Hernan!