Setting free surface BC

H Marcos,

I have a little bit of trouble understanding your question. You have a complex geometry and want to use a mesh function to describe the boundaries. Right? And your code below is not working because…? Have you tried meshbuilder? I think it's for situations like this, where you need to mark a boundary of a complex geometry, but I haven't used it myself.

Mikael

Den Mar 20, 2012 kl. 10:30 PM skrev Marcos Samudio:

> New question #191248 on CBC.PDESys:
> https://answers.launchpad.net/cbcpdesys/+question/191248
>
> 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.
>
> 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).
>
> The code I append is a little long, but the key is really only in the on_free_surface function, so you may skip the rest.
> I am trying to do it this way:
>
> #First define class to gather surface nodes data
> class FS_coordinates:
> def __init__(self,x,y,z):
> self.x = x
> self.y = y
> self.z = z
> FS_coor = FS_coordinates(nodes_x, nodes_y, nodes_z)
>
>
> #TEST WITH SIMPLE MESH. Gather points on surface x = 1
> mesh = UnitCube(3,3,3)
> coor = mesh.coordinates()
> 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)
>
> #Define inside_function
> def on_free_surface(x, on_boundary):
> fsc = FS_coor
> tol = 1e-6
> is_on_FS = False
> print x
> 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
> dist = sqrt(dist)
> if dist <tol:
> is_on_FS = True
> break
> return on_boundary and is_on_FS
>
> #Define additional BC
> def walls_bcs(x, on_boundary):
> return on_boundary and(near(x[0], 0.) or near(x[0], 1.) or near(x[1], 0.)\
> or near(x[1],1.) or near(x[2],0.))
>
> #Create mesh function
> mf = FaceFunction('uint', mesh)
> mf.set_all(0)
>
> #Mark boundaries
> inlet = FlowSubDomain(on_free_surface, bc_type='Wall', func = {'u': Constant((1., 0., 0.))},mf = mf)
> walls = FlowSubDomain(walls_bcs, bc_type = 'Wall', mf = mf)
>
> #Save file for further visualisation
> file = File("free_surface.pvd")
> file << mf
>
> Thanks,
>
> Marcos
>
> --
> You received this question notification because you are an answer
> contact for CBC.PDESys.

Hello Mikael,

I have a complex geometry, but I manage to get all nodes on the free surface boundary using Paraview.

What I want to do next is something similar to what you do with simple geometries, which is checking whether x is close to some value and return true if it is. In this case, the values to check against are the coordinates of all the nodes on the surface.

I calculate the distance between a given x and the surface nodes, and if it is close enough, then it is on the free surface. This would be my inside function for the FlowSubDomain class.

def on_free_surface(x, on_boundary):
    fsc = FS_coor #list of free surface coordinates
    tol = 1e-6
    is_on_FS = 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

I think this should work, but for some reason it does not!

This question was expired because it remained in the 'Open' state without activity for the last 15 days.