Determine boundaries by cell-by-cell inspection

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[0] = 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[0]) > 1-tol:
        nodes_x.append(i[0])
        nodes_y.append(i[1])
        nodes_z.append(i[2])
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

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:
> https://answers.launchpad.net/dolfin/+question/191226
>
> 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
>

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.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:
> https://answers.launchpad.net/dolfin/+question/191226
>
> Marcos Samudio gave more information on the question:
> 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[0] = 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[0]) > 1-tol:
> nodes_x.append(i[0])
> nodes_y.append(i[1])
> nodes_z.append(i[2])
> 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
>

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

On Tue, Mar 20, 2012 at 07:15:46PM -0000, Marcos Samudio wrote:
> Question #191226 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/191226
>
> Status: Answered => Open
>
> 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

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.

On Tue, Mar 20, 2012 at 09:25:42PM -0000, Marcos Samudio wrote:
> Question #191226 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/191226
>
> Status: Answered => Open
>
> 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.
>

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.

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

Thanks Mikael Mortensen, that solved my question.