# Determine boundaries by cell-by-cell inspection

Asked by Marcos Samudio on 2012-03-20

I am trying to perform a simulation of a 3D free surface flow.

I already have the mesh, and now I have to set the BCs for it.
For the no-slip surfaces I would just use on_boundary.
I can easily identify the inflow and outflow faces, as those are in planes, so I would do something like:

def inside(x, on_boundary):
on_boundary and (x in the inflow/outflow plane)

For the free surface BC I do have some additional work to do. I extracted the coordinates of all points lying on the free surface with Paraview, stored them in a Python array(FS_coor), performed an iteration through all cells in the mesh and stored in a list all cells(FS_cells) that have at least three of their nodes among the surface nodes (I am using CG degree 1 elements).

Is it possible to create a MeshFunction that can get the cells a face belongs to, and mark that face if it is both on_boundary and the cell it belongs to (the only cell it belongs to) is among FS_cells?

I hope the explanation was understandable, and if any of you can think of an easier way to do it, please let me know!

Thanks,

Marcos

## Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Mikael Mortensen
Solved:
Last query:
 Marcos Samudio (marcossamudio) said on 2012-03-20: #1

I think there might be an easier way to achieve this, but I am still missing something.

I append a piece of code that pretends marking faces of x = 1 on an unitary cube. The input would be the coordinates of the points on that boundary. This could be generalised to any boundary later.

class FS_coordinates:
def __init__(self,x,y,z):
self.x = x
self.y = y
self.z = z
from dolfin import *

#TEST WITH SIMPLE MESH
mesh = UnitCube(3,3,3)
coor = mesh.coordinates()
#Gather points on surface x = 1
tol = 1e-5
nodes_x = []
nodes_y = []
nodes_z = []
for i in coor:
if abs(i) > 1-tol:
nodes_x.append(i)
nodes_y.append(i)
nodes_z.append(i)
FS_coor = FS_coordinates(nodes_x, nodes_y, nodes_z)

class free_surf_boundary(SubDomain):
def __init__(self, free_surface_coordinates, *args):
SubDomain.__init__(self, *args)
self.free_surf_coor = free_surface_coordinates
def inside(self, x, on_boundary):
fsc = self.free_surf_coor
tol = 1e-4
is_on_FS = False
for y in range(len(fsc.x)):
p = [fsc.x[y], fsc.y[y], fsc.z[y]]
dist = 0.
for m in range(3): dist += (p[m] - x[m])**2
if dist <tol:
is_on_FS = True
break
return on_boundary and is_on_FS

mf = MeshFunction("uint",mesh, mesh.topology().dim()-1)
mf.set_all(0)
surface_boundary = free_surf_boundary(FS_coor)
surface_boundary.mark(mf, 1)
plot(mf, interactive=True)

I get a black box! What am I doing wrong?

Thanks again!

Marcos

 Anders Logg (logg) said on 2012-03-20: #2

Is it not possible to mark the free surface as the nodes which do
*not* belong to noslip, inflow our outflow?

--
Anders

On Tue, Mar 20, 2012 at 03:05:43PM -0000, Marcos Samudio wrote:
> New question #191226 on DOLFIN:
>
> I am trying to perform a simulation of a 3D free surface flow.
>
> I already have the mesh, and now I have to set the BCs for it.
> For the no-slip surfaces I would just use on_boundary.
> I can easily identify the inflow and outflow faces, as those are in planes, so I would do something like:
>
> def inside(x, on_boundary):
> on_boundary and (x in the inflow/outflow plane)
>
> For the free surface BC I do have some additional work to do. I extracted the coordinates of all points lying on the free surface with Paraview, stored them in a Python array(FS_coor), performed an iteration through all cells in the mesh and stored in a list all cells(FS_cells) that have at least three of their nodes among the surface nodes (I am using CG degree 1 elements).
>
> Is it possible to create a MeshFunction that can get the cells a face belongs to, and mark that face if it is both on_boundary and the cell it belongs to (the only cell it belongs to) is among FS_cells?
>
> I hope the explanation was understandable, and if any of you can think of an easier way to do it, please let me know!
>
> Thanks,
>
> Marcos
>

 Anders Logg (logg) said on 2012-03-20: #3

That's more code than I have time to read...

Why does it need to be that complex? It should be enough with
something like

class Boundary(SubDomain):
def inside(x, on_boundary):
return near(x, 0.0)

f = FacetFunction(mesh)
b = Boundary()
b.set_all(0)
b.mark(f, 1)

--
Anders

On Tue, Mar 20, 2012 at 05:30:53PM -0000, Marcos Samudio wrote:
> Question #191226 on DOLFIN changed:
>
> I think there might be an easier way to achieve this, but I am still
> missing something.
>
> I append a piece of code that pretends marking faces of x = 1 on an
> unitary cube. The input would be the coordinates of the points on that
> boundary. This could be generalised to any boundary later.
>
> class FS_coordinates:
> def __init__(self,x,y,z):
> self.x = x
> self.y = y
> self.z = z
> from dolfin import *
>
> #TEST WITH SIMPLE MESH
> mesh = UnitCube(3,3,3)
> coor = mesh.coordinates()
> #Gather points on surface x = 1
> tol = 1e-5
> nodes_x = []
> nodes_y = []
> nodes_z = []
> for i in coor:
> if abs(i) > 1-tol:
> nodes_x.append(i)
> nodes_y.append(i)
> nodes_z.append(i)
> FS_coor = FS_coordinates(nodes_x, nodes_y, nodes_z)
>
> class free_surf_boundary(SubDomain):
> def __init__(self, free_surface_coordinates, *args):
> SubDomain.__init__(self, *args)
> self.free_surf_coor = free_surface_coordinates
> def inside(self, x, on_boundary):
> fsc = self.free_surf_coor
> tol = 1e-4
> is_on_FS = False
> for y in range(len(fsc.x)):
> p = [fsc.x[y], fsc.y[y], fsc.z[y]]
> dist = 0.
> for m in range(3): dist += (p[m] - x[m])**2
> if dist <tol:
> is_on_FS = True
> break
> return on_boundary and is_on_FS
>
> mf = MeshFunction("uint",mesh, mesh.topology().dim()-1)
> mf.set_all(0)
> surface_boundary = free_surf_boundary(FS_coor)
> surface_boundary.mark(mf, 1)
> plot(mf, interactive=True)
>
> I get a black box! What am I doing wrong?
>
> Thanks again!
>
> Marcos
>

 Marcos Samudio (marcossamudio) said on 2012-03-20: #4

Thanks for the fast response Anders.
I did not get what you said in #3. I am not trying to solve for a cube, but for a general domain. The cube is just a simpler example.
About #2, how can I get those faces (the ones that do not belong either to the noslip, inflow or outflow) marked in a MeshFunction?

Thanks,

Marcos

 Anders Logg (logg) said on 2012-03-20: #5

On Tue, Mar 20, 2012 at 07:15:46PM -0000, Marcos Samudio wrote:
> Question #191226 on DOLFIN changed:
>
>
> Marcos Samudio is still having a problem:
> Thanks for the fast response Anders.
> I did not get what you said in #3. I am not trying to solve for a cube, but for a general domain. The cube is just a simpler example.
> About #2, how can I get those faces (the ones that do not belong either to the noslip, inflow or outflow) marked in a MeshFunction?

They will be the ones you have not marked.

Mark everything as say 0, then mark noslip as 1, inflow as 2, outflow
as 3. The free surface should then be the ones marked 0.

--
Anders

 Marcos Samudio (marcossamudio) said on 2012-03-20: #6

The mesh geometry is really complex, so I mark all exterior faces as no_slip first, then I mark the inflow and outflow, and the free surface remains to be marked (it was previously marked as no_slip).

A shorter code would be:

#fsc holds the coordinates of the surface nodes
def on_free_surface(x, on_boundary):
tol = 1e-4
is_on_surface = False
for i in range(len(fsc.x)):
p = [fsc.x[i], fsc.y[i], fsc.z[i]]
dist = 0.
for m in range(3): dist += (p[m] - x[m])**2
dist = sqrt(dist)
if dist <tol:
is_on_FS = True
break
return on_boundary and is_on_FS

mf = FaceFunction('uint', mesh)
mf.set_all(0)
inlet = AutoSubDomain(on_free_surface)
inlet.mark(mf,1)

I tried this with the simple unitary cube, but the whole domain is being marked.

 Anders Logg (logg) said on 2012-03-20: #7

On Tue, Mar 20, 2012 at 09:25:42PM -0000, Marcos Samudio wrote:
> Question #191226 on DOLFIN changed:
>
>
> Marcos Samudio is still having a problem:
> The mesh geometry is really complex, so I mark all exterior faces as
> no_slip first, then I mark the inflow and outflow, and the free surface
> remains to be marked (it was previously marked as no_slip).

Do you have other boundaries in addition to noslip, inflow, outflow,
free surface?

--
Anders

> A shorter code would be:
>
> #fsc holds the coordinates of the surface nodes
> def on_free_surface(x, on_boundary):
> tol = 1e-4
> is_on_surface = False
> for i in range(len(fsc.x)):
> p = [fsc.x[i], fsc.y[i], fsc.z[i]]
> dist = 0.
> for m in range(3): dist += (p[m] - x[m])**2
> dist = sqrt(dist)
> if dist <tol:
> is_on_FS = True
> break
> return on_boundary and is_on_FS
>
> mf = FaceFunction('uint', mesh)
> mf.set_all(0)
> inlet = AutoSubDomain(on_free_surface)
> inlet.mark(mf,1)
>
> I tried this with the simple unitary cube, but the whole domain is being
> marked.
>

 Marcos Samudio (marcossamudio) said on 2012-03-21: #8

I have one inflow, two outflow boundaries (I am studying bifurcation behaviour), the noslip condition at the bottom and the free surface.

I am using the CBC.PDESys package, but I thought it would be useful to post the question here as well. The actual formulation I am using would be the one in https://answers.launchpad.net/cbcpdesys/+question/191248.

What I posted here is about the same, using only dolfin. Mikael Mortensen (mikael-mortensen) said on 2012-03-21: #9

Hi Marcos,

Your code is not working because the midpoint of the facet is not in you list of nodes. See towards the end of mesh/SubDomain.cpp. You can check out the section starting with // check midpoint and recompile dolfin. Or you can provide the midpoints in your list of points to check against as well.

Best regards

Mikael

Den Mar 21, 2012 kl. 2:10 PM skrev Marcos Samudio:

> Question #191226 on DOLFIN changed:
>
>
> Marcos Samudio is still having a problem:
> I have one inflow, two outflow boundaries (I am studying bifurcation
> behaviour), the noslip condition at the bottom and the free surface.
>
> I am using the CBC.PDESys package, but I thought it would be useful to
> post the question here as well. The actual formulation I am using would
>
> What I posted here is about the same, using only dolfin.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

 Marcos Samudio (marcossamudio) said on 2012-03-21: #10

Hi Mikael,

I was not aware of that. Anyway I think I found out a way to do it:

mf = FacetFunction("uint", mesh)
mf.set_all(0)
walls = FlowSubDomain(lambda x, on_boundary: on_boundary, bc_type = 'Wall', mf = mf)

mesh_coor = mesh.coordinates()
tol = 1e-4
fsc = vertices_on_free_surface
for facet in facets(mesh):
#Iterate only through facets on the boundary
if mf[facet] != 1:
continue
vert_on_bnd = 0
for vertex in vertices(facet):
for i in range(len(fsc.x)):
p = [fsc.x[i], fsc.y[i], fsc.z[i]]
dist = 0.
j = vertex.index()
for m in range(3):
dist += (mesh_coor[j][m] - p[m])**2
dist = sqrt(dist)
if dist <tol:
vert_on_bnd +=1
if vert_on_bnd == 3:
mf[facet] = 2

It worked out for the simple cube. The thing is I do not have the midpoints of facets, so this is faster in my case. If you notice anything wrong, please let me know!

Best regards,

Marcos

 Marcos Samudio (marcossamudio) said on 2012-03-21: #11

Thanks Mikael Mortensen, that solved my question.