urgent help needed: did I do wrong to set up boundaries with Dirichlet and Neumann boundaries with mesh from an xml file

Question

Hi, all experts,

I am trying to solve a flow in porous media problem. After simplification, it becomes a transient diffusion problem for pressure. The 3D domain has 7 surfaces, and the 3D mesh was generated using NETGEN and exported as diffpack format, then it was converted using "dolfin-convert". However, when I run my project, at the first time step, it gives "Segmentation fault (core dumped)".

I tested using a 2-D rectangle domain and it works fine and the result is good.

So I suspect that the way to set up the boundaries for my 3D case is not right. The piece of code is attached below. Any help is highly appreciated.

The boundary surfaces are numbered from 1 to 7 in NETGEN. Surface #7 is the Dirichlet boundary with Pressure = constant. Surface #5 is the Neumann boundary, dP/dn = beta = constant (non-zero). For all other surfaces, dP/dn=0.

As you can see, beta is a large value. Thus in both 2D case and 3D case, similar high mesh resolution is used along the normal direction at the Neumann boundary with dP/dn = beta.

Here is part of my code:

mesh = Mesh("xld_r5d05_v2.xml")
V = FunctionSpace(mesh, "Lagrange",1)

alpha = 5.48359766E-02
beta = 9.22675065E+07
Pf = Constant(0.0)

boundaries = MeshFunction('size_t',mesh)

bcs = DirichletBC(V, Pf, boundaries, 7)

ds = Measure('ds')[boundaries]

t_stop = 5 # second
dt = 1e-3
theta = 1

p_1 = interpolate(Pf, V)

p = TrialFunction(V)
v = TestFunction(V)

a = pvdx + thetadtalphainner(nabla_grad(v), nabla_grad(p))dx

L = (p_1v - (1-theta)dtalphainner(nabla_grad(v), nabla_grad(p)))dx - dtalphabetav*ds(5)

Comments

boundaries = MeshFunction('size_t',mesh)

you need specify the type of mesh function here, cell, facet or vertices?

---

I intended to obtain the correct boundary surfaces in order to set up the boundary conditions.

As I mentioned, my mesh is generated using NETGEN then converted to xml using dolfin-convert. For example, my mesh file is "xld_r5d05_v2.xml", and dolfin-convert also generated several xml file such as "xld_r5d05_v2_marker_1.xml" and so on. My guess is that these files represent the surfaces, but I don't know if there is a way to directly get the right boundary from these files.

I need to set Dirichlet boundary on one surface, and Neumann bc on another one. All other surfaces with Neumann condition of zero-flux.

I searched around and can't ready find any clue.

For your comment, I guess that the type should be facet.

Your help is highly appreciated.

Answer 1

try:
boundaries = MeshFunction('size_t',mesh,2)