How to solve PDE in subdomain with mesh import from Gmsh?


332
views
1
5 months ago by
Hello all,
I have a 2 phase problem to solve, it is like the figure below:

Now I want to solve a navier-stoke equation in the liquid part, which is labelled as 2.
And its boundary is labelled as 3.
I generate my mesh from Gmsh.
The code is like:
mesh = Mesh("filename.xml")

cd  = MeshFunction('size_t', mesh, "filename_physical_region.xml")
cdl = MeshFunction('size_t', mesh, "filename_facet_region.xml")

subdomain = CellFunction("size_t", mesh)
subdomain.set_all(0)
subdomain.array()[cd.array()==1] = 1
subdomain.array()[cd.array()==2] = 2

boundary = FacetFunction("size_t", mesh)
boundary.set_all(0)
boundary.array()[cdl.array()==3] = 3

# Now I guess in order to solve the navier stokes equation only in subdomain 2,
# I need to create a submesh and a FunctionSpace in the subdomain like:

submesh = SubMesh(mesh, subdomain, 2)
W = FunctionSpace(submesh, 'P', 1)

# but how do I set the Dirichlet boundary condition on the curve 3?
# I tried:

bcw = Dirichlet(W, Constant(0.0), cdl, 3)

# of course I got the error's reason: User MeshFunction and FunctionSpace meshes are different​

Thank you!!!

Community: FEniCS Project
The code like below does not work neither, there is no error report but the result makes no sense:
cdl2 = MeshFunction("size_t", submesh, "filename_facet_region.xml")

bcw = DirichletBC(W, Constant(0.0), cdl2, 3)
​
written 5 months ago by SICHENG SUN  

1 Answer


2
5 months ago by
Hi SICHENG SUN,
if I understood correctly and you succeeded in creating the submesh for subdomain 2, a straightforward solution could be:
def nolisp(x, on_boundary):
    if on_boundary:
        if near(x[0], Right, tol) or near(x[0], Left, tol) or near(x[1], Bottom, tol):
            return False
        else:
            return True
    else:
        return False​

bcw = DirichletBC(V, Constant(0.0), nolisp)

Thank you!
Well, my no slip boundary is only on the red line, which is labelled as 3.
Other boundaries are just Neuman boundary with flux = 0.(So we can leave them alone.)
----------------------------
But I actually have solved it in another way.
I solve the velocity in the whole domain rather than the sub domain. And set the velocity in the solid equal to 0.
written 5 months ago by SICHENG SUN  
2
The noslip function returns False for the bottom, left and right boundaries.
It should (probably) return True only for the red boundary.
written 5 months ago by Igor Baratta  
Oh, I see it now. That's such a clever way to define to boundary. I have never thought like this.
Thx a lot!
Although I have already solved this problem, you provided me another different and intelligent method. I am really appreciate it.
written 5 months ago by SICHENG SUN  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »