Trouble applying dim=0 mesh function boundary definition using DirichletBC

Question

Hey guys

For different reasons, I need to generate the mesh using my own program. I write two Dolfin XML files, one for the mesh itself:

<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://www.fenicsproject.org">
  <mesh celltype="tetrahedron" dim="3">
    <vertices size="428868">
      <vertex index="0" x="0.0" y="0.0" z="413.54"/>
      <vertex index="1" x="2200.0" y="0.0" z="272.73"/>
      ...

And one for the boundaries (boundary.xml, gzipped boundary.xml.gz). I'm marking them on every vertex using a 0-dimensional mesh_function:

<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
  <mesh_function type="uint" dim="0" size="428868">
    <entity index="0" value="1"/>
    <entity index="1" value="1"/>
    ...

Now I would like to apply DirichletBC on every DOF where my boundary mesh function is set ==2. I tried it like this:

subdomains = MeshFunction('size_t', mesh, boundary.xml.gz')
bc = DirichletBC(V, Constant(0), subdomains, 2)

But I can clearly see in my solution, that the DirichletBC is not being applied. I also tried pointwise:

subdomains = MeshFunction('size_t', mesh, 'boundary.xml.gz')
bc = DirichletBC(V, zero, subdomains, 2, 'pointwise')

But then I get this error:

*** Error:   Unable to computing Dirichlet boundary values, pointwise search.
*** Reason:  A SubDomain is required for pointwise search.
*** Where:   This error was encountered inside DirichletBC.cpp.

Do I have to create a FacetFunction using my mesh function for use in DirichletBC? How would I do that? I know that the boundary meshes provided by gmsh are of dimension (mesh-1), hence Facets.

There could be other ways marking my boundary tough. I have an irregular 3D domain decomposed into tetraedra. I need to mark all boundary points except the boundary with the lowest z coordinates. Since I know which points are the lowest when generating the mesh, I thought it would be best to mark them and import this definition into FEniCS.

Related Questions that I read but still cannot put together my puzzle:
- http://fenicsproject.org/qa/8780/what-is-the-best-way-to-set-a-subdomain-to-a-given-value
- http://fenicsproject.org/qa/6267/dirichlet-boundary-condition-element-from-imported-xml-mesh
- http://fenicsproject.org/qa/2986/how-to-define-boundary-condition-for-mesh-generated-by-gmsh
- http://fenicsproject.org/qa/8918/dolfin-xml-subdomain-format

Answer 1

This answer contains the code of the answer proposed by maartent. DirichletBC needs a MeshFunction defined on the facets (in this case triangles). Hence one can be created, using the MeshFunction defined on the vertices.

from dolfin import *

mesh = Mesh('cosmo2.xml.gz')
vertexfunc = MeshFunction('size_t', mesh, 'cosmo2_vertexboundary.xml.gz')

# create facetfunction
facetfunction = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
facetfunction.set_all(0)

for f in facets(mesh):
    # check if every vertex on this facet has the same boundary marker; if yes store it in the facetfunction
    last = -1
    isboundary = True
    for v in vertices(f):
        val = vertexfunc[v]
        if last==-1:
            # first iteration
            last = val
        else:
            # any other iteration
            if last != val:
                isboundary = False
                break

    if isboundary:
        facetfunction[f] = last

File('cosmo2_boundary.xml.gz') << facetfunction

Comments

I also tried defining the facet MeshFunction only on the BoundaryMesh; but the numbering of the vertex MeshFunction is of course not compatible with the numbering on the BoundaryMesh; hence I abandoned the idea. I think there could be a way tough; but I'm not sure if DirichletBC accepts subdomains that are defined on a BoundaryMesh.

Answer 2

Assuming I understood your question correctly, my plan of attack would be

create a new FacetFunction, initialize to zero
iterate over all facets, for each facet
    get its vertices and their value from your 0-dim meshfunction
    if these values are all equal (let's say to n)
        then mark this facet as 'n'

You can then use this new FacetFunction in your DirichletBC. The first link you mentioned should get you on your way.

Comments

Thanks for your solution. I will definitely go that way then. But it seems to me that it's more like a makeshift. I'd like to go parallel with this code, and loops like that won't be distributed to different cores automatically I guess (that would amaze me). Shouldn't there be a better solution? Do I have to create and import my boundary markers differently?

The problem when thinking about creating the BoundaryMesh is that I don't know the Facet numbering FEniCS uses when creating my mesh...

How could I create a FacetFunction using my own program that can be imported directly into FEniCS?

---

It isn't the most elegant solution and there might be others, but I am not aware of them yet.

If I am not mistaken you can write FacetFunction's to a file, and import it again when needed. I couldn't find this in the documentation, though, so you'll have to try it yourself or ask someone else.

---

Yes you can create this MeshFunction on the Facet and save it to reuse it. Actually that is a good idea, since my mesh is (still) static. Thanks for your answer and hint!