DOLFIN

MeshFunction and SubMesh

Asked by Simone Pezzuto

Given a marked mesh (region and facets), I'd like to create submesh with the corresponding marked elements.

Example (quite self-explanatory):

mesh = Mesh("vessel.xml.gz")
dfun = MeshFunctionSizet(mesh, "vessel_physical_region.xml.gz")
bfun = MeshFunctionSizet(mesh, "vessel_facet_region.xml.gz")

mesh_fluid = SubMesh(mesh, dfun, 1)
mesh_struct = SubMesh(mesh, dfun, 2)

bfun_fluid = give_me_a_name(bfun, mesh_fluid)
bfun_struct = give_me_a_name(bfun, mesh_struct)

The former problem was to apply assemble a problem on the submesh keeping the correct marked elements.
Of course I cannot use directly bfun, because facets ids are different.

Is there an equivalent to "give_me_a_name" in Dolfin?

Here my work around (not really efficient):

def give_me_a_name(bfun, submesh) :

  # creates the submesh bfun
  bfun_submesh = MeshFunctionSizet(submesh, bfun.dim())

  # submesh to mesh mappings
  cfun = submesh.data().mesh_function("parent_cell_indices")
  pfun = submesh.data().mesh_function("parent_vertex_indices")

  # iterate over the cells
  for cell in cells(submesh) :

    # extracts the facets and the corresponding one of the full mesh
    facets_sub = cell.entities(2)
    facets = Cell(mesh, cfun.array()[cell.index()]).entities(2)

    for facet_sub in facets_sub :
      # mapped vertices
      vert_mapped = pfun.array()[Facet(submesh, facet_sub).entities(0)]
      # find the corresponding facet on the full mesh
      for facet in facets :
        # vertices
        vert = Facet(mesh, facet).entities(0)
        # intersection
        common = set(vert).intersection(set(vert_mapped))
        if len(common) == 3 :
          bfun_submesh.array()[facet_sub] = bfun.array()[facet]
          break

  return bfun_submesh

Sub-question: is there a way to know the type ("size_t", "uint", ...) of a MeshFunction?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
Last query:
Last reply:
Revision history for this message
Best Johan Hake (johan-hake) said : 		#1

With trunk I think this should be possible. However it works only for
markers stored as MeshValueCollections in the mesh. Try:

 mesh = Mesh("vessel.xml.gz")
 dfun = MeshFunctionSizet(mesh, "vessel_physical_region.xml.gz")
 bfun = MeshFunctionSizet(mesh, "vessel_facet_region.xml.gz")
 mesh.domains().markers(2).assign(bfun)

 mesh_fluid = SubMesh(mesh, dfun, 1)
 mesh_struct = SubMesh(mesh, dfun, 2)

The two sub mesh should now include correct facet markers. If you want
the cell markers to propagate too, you can do the assignment for that
function too, but it might not make sense as all cells have the same
marker anyway.

You can now assemble forms over these meshes without extracting the
facet function. But if you want to access them you can do:

  mesh_fluid = mesh_fluid.domains().facet_domains()
  mesh_struct = mesh_struct.domains().facet_domains()

Disclaimer:
The code above is not tested. It might be some small errors.

Johan

On 03/02/2013 09:01 PM, Simone Pezzuto wrote:
> New question #223248 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/223248
>
> Given a marked mesh (region and facets), I'd like to create submesh with the corresponding marked elements.
>
> Example (quite self-explanatory):
>
> mesh = Mesh("vessel.xml.gz")
> dfun = MeshFunctionSizet(mesh, "vessel_physical_region.xml.gz")
> bfun = MeshFunctionSizet(mesh, "vessel_facet_region.xml.gz")
>
> mesh_fluid = SubMesh(mesh, dfun, 1)
> mesh_struct = SubMesh(mesh, dfun, 2)
>
> bfun_fluid = give_me_a_name(bfun, mesh_fluid)
> bfun_struct = give_me_a_name(bfun, mesh_struct)
>
> The former problem was to apply assemble a problem on the submesh keeping the correct marked elements.
> Of course I cannot use directly bfun, because facets ids are different.
>
> Is there an equivalent to "give_me_a_name" in Dolfin?
>
> Here my work around (not really efficient):
>
> def give_me_a_name(bfun, submesh) :
>
> # creates the submesh bfun
> bfun_submesh = MeshFunctionSizet(submesh, bfun.dim())
>
> # submesh to mesh mappings
> cfun = submesh.data().mesh_function("parent_cell_indices")
> pfun = submesh.data().mesh_function("parent_vertex_indices")
>
> # iterate over the cells
> for cell in cells(submesh) :
>
> # extracts the facets and the corresponding one of the full mesh
> facets_sub = cell.entities(2)
> facets = Cell(mesh, cfun.array()[cell.index()]).entities(2)
>
> for facet_sub in facets_sub :
> # mapped vertices
> vert_mapped = pfun.array()[Facet(submesh, facet_sub).entities(0)]
> # find the corresponding facet on the full mesh
> for facet in facets :
> # vertices
> vert = Facet(mesh, facet).entities(0)
> # intersection
> common = set(vert).intersection(set(vert_mapped))
> if len(common) == 3 :
> bfun_submesh.array()[facet_sub] = bfun.array()[facet]
> break
>
> return bfun_submesh
>
>
> Sub-question: is there a way to know the type ("size_t", "uint", ...) of a MeshFunction?
>
>

Revision history for this message
Simone Pezzuto (junki-gnu) said : 		#2

Thanks Johan, it works like a charm. It's worth mentioning that you may assign also dfun to 3d markers:

mesh.domains().markers(3).assign(dfun)

mesh_fluid = SubMesh(mesh, 1)
mesh_struct = SubMesh(mesh, 2)

*Remark*
the code doesn't work in parallel, because SubMesh is not yet supported.

 Simone

Revision history for this message
Simone Pezzuto (junki-gnu) said : 		#3

Thanks Johan Hake, that solved my question.

Revision history for this message
Johan Hake (johan-hake) said : 		#4

On 03/03/2013 09:30 PM, Simone Pezzuto wrote:
> Question #223248 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/223248
>
> Status: Answered => Solved
>
> Simone Pezzuto confirmed that the question is solved:
> Thanks Johan, it works like a charm. It's worth mentioning that you may
> assign also dfun to 3d markers:
>
> mesh.domains().markers(3).assign(dfun)
>
> mesh_fluid = SubMesh(mesh, 1)
> mesh_struct = SubMesh(mesh, 2)

Yes, I added that functionality too ;)

> *Remark*
> the code doesn't work in parallel, because SubMesh is not yet supported.

Yes not sure how this can be done properly. There are some FIXMEs in the
SubMesh code which give some hints. This might help parallelize the
different submesh, but then you need to make sure the dofs are correctly
partitioned. This might or might not be the case. My knowledge is too
limited here.

Johan

Revision history for this message
Simone Pezzuto (junki-gnu) said : 		#5

A possibility is to associate each SubMesh to a specific processor.
This would be nice for domain decomposition methods, but the number of processor must be a multiple of the submeshes.

In some some, it's like using the submeshes a coloring for the partitions.

 Simone

Revision history for this message
Johan Hake (johan-hake) said : 		#6

On 03/04/2013 10:11 AM, Simone Pezzuto wrote:
> Question #223248 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/223248
>
> Simone Pezzuto posted a new comment:
> A possibility is to associate each SubMesh to a specific processor.

I do not think this would scale for most usage of SubMesh. One could be
to repartion the mesh. That would, however, involve a lot of
communication and you would loose connection to the original mesh.

> This would be nice for domain decomposition methods, but the number of processor must be a multiple of the submeshes.

Which it very seldom is I guess?

> In some some, it's like using the submeshes a coloring for the
> partitions.

What do you mean with that?

Johan

> Simone
>

Revision history for this message
Simone Pezzuto (junki-gnu) said : 		#7

2013/3/4 Johan Hake <email address hidden>

> > This would be nice for domain decomposition methods, but the number of
> processor must be a multiple of the submeshes.
>
> Which it very seldom is I guess?
>

True, but a segregated approach to fluid-structure interaction is in fact a
domain decomposition method (with no overlap.) It could be useful to
associate
a submesh to a specific process (or a set of processes).

Maybe we can reformulate the condition as "the number of processes must be
at least as much as the number of submeshes". This is more reasonable.

>
> > In some some, it's like using the submeshes a coloring for the
> > partitions.
>
> What do you mean with that?
>

I was thinking to the fact that a partition into submeshes can be in fact a
partition of the whole mesh. Or, it's a constraint on the forthcoming
partition,
in the sense that the submeshes must be the union of partitions (i.e. the
partitioner cannot select vertices or elements belonging to different
submeshes).

>
> Johan
>
> > Simone
> >
>
> --
> You received this question notification because you asked the question.
>