How to define boundary conditions on an irregular geometry?

[sergeyg]

I would like to solve a Laplace equation on a 2D domain built as a polygon. I would like to set up a combination of Dirichlet boundary conditions on different portions of the boundary.
How do I select a desired portion of a boundary (i.e, a portion of a polygon)?

To be more specific here is an example I was playing with but I know that I did not correctly specify the condition for the top boundary.
How could I specify a boundary conditions for the top boundary?

define a domain

domain_vertices = [
Point(0, 0),
Point(1, 0),
Point(1, 1),
Point(0.75, 1.1),
Point(0.7, 1.2),
Point(0.6, 1.2),
Point(0.5, 1.15),
Point(0.2, 1.1),
Point(0, 1),
Point(0, 0),
]

domain = Polygon(domain_vertices)
mesh = generate_mesh(domain,20)
V = FunctionSpace(mesh, ‘P’, 1)

Define boundary condition

def boundary_top(x, on_boundary):
return on_boundary and near(x[1], 1) # this is wrong

def boundary_bottom(x, on_boundary):
return on_boundary and near(x[1], 0)

u_top = Constant(2.0)
u_bottom = Constant(1.0)

bc_top = DirichletBC(V,u_top,boundary_top)
bc_bottom = DirichletBC(V,u_bottom,boundary_bottom)

bcs = [bc_top, bc_bottom]

Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0) # no source for the Laplace equation
a = dot(grad(u), grad(v))dx
L = fv*dx

Compute solution

u = Function(V)
solve(a == L, u, bcs)

Plot solution and mesh

plot(u,title=“solution”)
plot(mesh,title=“mesh”)

plt.show()

[dokken]

For more complicated geometries , its common to use external software, like gmsh/pygmsh. With these softwares, you can set physical markers for cells snd facets, Which can be loaded into dolfin as described here:

Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

This can be done using meshio and writing the explicit data to xdmf file. import meshio msh = meshio.read("test.msh") meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]})) meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]}, cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}})) from dolfin import * mesh = Mesh() with XDMFFile("mesh.xdmf") …

[sergeyg]

Thanks for a reply. I can try external software if necessary, but is there still a way of marking parts of a boundary inside fenics? I already have a mesh build in fenics.

[dokken]

Yes, if you have some Notion of Which parts of the boundary should be marked, you could follow the strategy outlined in this tutorial. In for the mesh above, i guess that would be x[0]>1e-14 and x[0]<1-eps and x[1]>1 and on_boundary

[sergeyg]

Thanks again, I have read that tutorial. In that tutorial the irregular boundary corresponds to a single domain. My problem is different. My earlier example above was not a very good example. My actual problem has a more irregular geometry where facets do not align with the x and y axis. For example, the domain below with the number next to the boundary corresponding to different subdomains each of which would have a different boundary condition.

Is there a way identifying boundaries in such a case inside fenics?

[mkondratyev85]

Hello @dokken. What if I have a mesh saved in a file with some markers set on some cells? Say 1 and 0 on different cells. I can load the mesh from file in my fenics python script, but I don’t understand how can I use marker values to set boundary conditions. All functions that set bc that i’ve found in the examples respect only coorindates of the vertex in space. Thank you.

[dokken]

Then you need to load a MeshFunction. As shown in Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio you load the mesh function mf And use that In the DirichletBC, As shown in the Stokes demo

[mkondratyev85]

Thank you for your reply. Quoting from your example:

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

I have a file mesh.xdmf where I have several fields, and among one of them I have a cell field named ‘composition’ with some cells have value 0 and others have value 1. Can I load a MeshFunction from this single mesh? Something like this:

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mf = cpp.mesh.MeshFunctionSizet(mesh, 'composition_field')
#                                                     ^
#                                                     |
#                                                     |
#          some how set just the field name from the mesh

or do I have to create another mesh (i.e. mf.xdmf from your example) with some external package. If so, how should I name the data field with 0s and 1s in this mesh?

[BillS]

You can use functions that describe your boundary curves (lines) and use the fact that the boundary DOF locations must lie on (or near) the curve within the endpoint constraints for each curve segment. For example:

class BorderBC(SubDomain):
    def inside(self, x, on_boundary):
        r = (y1 - y0) * (x[0] - x0) / (x1 - x0) + y0 - x[1] 
        return on_boundary and (near(r, 0.0, tol) and (x[0] >= x0) and (x[0] <= x1))

x0,y0 is the start point of your line segment and x1, y1 is the end point. If the point sits on the line, r will vanish (or be very small). The last two conditions of the test line set the limits of the line segment. You can use x[1] instead of x[0] if it is more appropriate. Do this for each line segment of your boundary.

[dokken]

Please provide a minimal example with an error message if the suggested code does not work. But you Need to save boundary data as a separate file for dolfin to read it in. If you are reading in cell data, you can store it alongside the mesh

[mkondratyev85]

That’s the snipped I’m trying to use:

mesh = Mesh()
with XDMFFile("/tmp/mesh.xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 2) 
with XDMFFile("/tmp/mf.xdmf") as infile:
    infile.read(mvc, "bc")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

V = FunctionSpace(mesh, 'P', 1)

bc_intrusion = DirichletBC(V, Constant(100), mf, 2)

And then I get the following error:

Traceback (most recent call last):
  File "ft01_poisson.py", line 31, in <module>
    bc_intrusion = DirichletBC(V, Constant(100), mf, 2)
  File "/usr/lib/python3/dist-packages/dolfin/fem/dirichletbc.py", line 131, in __init__
    super().__init__(*args)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  User MeshFunction is not a facet MeshFunction (dimension is wrong).
*** Where:   This error was encountered inside DirichletBC.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

I could upload the meshes if it is required.

[dokken]

As the error message suggest, your mesh function does not consist of boundary data, bit cell data. This is due to:
mvc = MeshValueCollection("size_t", mesh, 2)
Which for a 2D mesh should be:
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
Which is 1, and not 2 in the case of a 2D mesh

[mkondratyev85]

But my mesh is 3D. I was naively using this line from the example without understaing of the meaning of function arguments.

[mkondratyev85]

When I use this line instead of mine:

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)

I still get the same error. Is this the correct line if the mesh is in 3D?

[dokken]

It should be correct, so could you supply the mesh and boundary data file?

[mkondratyev85]

File meshes.zip is in the following link:

easyupload.io

[dokken]

So inspecting your boundary data in Paraview, I see that this is not boundary data, but cell data (The elements are tetrahedrons). Boundary data has to be based on a FacetMesh (made up by triangles). May I ask what kind of mesh generator you used to generate your mesh?

[mkondratyev85]

That is just me lacking understanding of basic terminology. I create my meshes using pyvista package in python. So, some of my cells have value of 2 and I want to set initial condition to these cells and then solve Poisson equation. In pyvista I can assign values to vertexes and to cells.

[dokken]

Hi again,
So you can convert your cell data to boundary data with the following functions.
By replacing the cellfunction i made in the code with the one you are loading into dolfin, and the set criteria with the corresponding value from your mesh file and you should be good to go

from dolfin import *

mesh = UnitCubeMesh(4,4,4)
mesh.init()
class top_half(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0]<=0.5)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
top_half().mark(domains,1)


boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
cell_to_facet = mesh.topology()(mesh.topology().dim(), mesh.topology().dim()-1)
for cell in cells(mesh):
    for face in cell_to_facet(cell.index()):
        facet = Facet(mesh, face)
        if facet.exterior():
            if domains.array()[cell.index()]==1:
                boundaries.set_value(facet.index(), 1)

XDMFFile("bds.xdmf").write(boundaries)

[mkondratyev85]

Still doesn’t work. I’m not sure if I understood correctly what do I have to change in your code, but here is the whole script I used:

from dolfin import *
from fenics import *

mesh = Mesh()
with XDMFFile("/tmp/fenics/mesh.xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("/tmp/fenics/mf.xdmf") as infile:
    infile.read(mvc, "bc")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
domains = mf

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
cell_to_facet = mesh.topology()(mesh.topology().dim(), mesh.topology().dim()-1)
for cell in cells(mesh):
    for face in cell_to_facet(cell.index()):
        facet = Facet(mesh, face)
        if facet.exterior():
            if domains.array()[cell.index()]==1:
                boundaries.set_value(facet.index(), 1)

XDMFFile("bds.xdmf").write(boundaries)

When I look on bds.xdmf file in paraview I see only 0 values.

The most confusing part for me is this:

        if facet.exterior():
            if domains.array()[cell.index()]==1:
                boundaries.set_value(facet.index(), 1)

Why do we check for exterior? As I understand in last two lines we check if value in domains mesh is equal to 1 and then we set the value to the boundaries mesh. Is it right?