mark internal facets

Question

Hi guys,

I have a mesh constituted of 3d cubes. I want to mark all facets on the plane x[0]=a that are located between x1min and x1max along the direction x[1] and between x2min and x2max along the direction x[2].

The code below,

Facets_domain = FacetFunction('size_t', mesh)

class Bndry_(SubDomain):
    def inside(self, x, on_boundary):
        return (abs(x[0]-a) < tolerance_boundary_location
                and x[1] > x1min-tolerance_boundary_location
                and x[1] < x1max+tolerance_boundary_location
                and x[2] > x2min-tolerance_boundary_location
                and x[2] < x2max+tolerance_boundary_location)
Bndry_().mark(Facets_domain,1)

marks correctly the facets as i want, except that only the facets at the exterior are marked. Actually, i want to mark the facets within the volume to specify later a flux as an internal boundary condition.

Someone can tell me what's wrong on the code?
Thanks.

Comments

I am actually wandering how to code an internal Neumann BC on precise locations.
Even if i can mark correctly the internal facets (with a value 1 in this example), i am wondering if i can modify the variational formulation as simple as:

F=... - desired_flux_value*v*ds(1)

Especially since i want to impose a zero flux condition on these facets (desired_flux_value=0), then... i actually add 0 to the variational formulation. I guess that's not the correct way to do this.
Does someone know how to correctly introduce internal Neumann BC on internal facets?
Any help would be appreciate.

Answer 1

Regarding the boundary condition, you have to use dS instead of ds for integrals over interior facets. Implementing the boundary condition is more involved than just adding another integral on the right hand side though. You also have to allow for discontinuity across the internal boundary, at least for the test functions. Take a look at the answer to 10927.

Comments

Ok, then, here is a simplified version of my code:
I impose a Dirichlet BC at the left and at the right boundary of a rectangle, and i want one Neumann BC at the middle of the domain (stationary diffusion)

from dolfin import *

print "The current version of FEniCS is",dolfin.dolfin_version()
# The current version of FEniCS is 2016.1.0

### Check machine precision and choose as a consequence the tolerance
print 'Machine precision={value}'.format(value=DOLFIN_EPS)
tolerance_boundary_location = DOLFIN_EPS*1000;
print 'Boundary detection tolerance={value}'.format(value=tolerance_boundary_location)

### Create simple geometry
length_0=2; length_1=1; length_2=1
n_0 = 6; n_1=3; n_2=3
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length_0, length_1, length_2), n_0, n_1,n_2)
plot(mesh, title='Mesh of the total domain', interactive=True)

### Set finite elements
V = FunctionSpace(mesh, 'Lagrange', 1)

### Set Dirichlet boundary location at left and rigth
# Location
def LeftBndry(x):
    return ( abs(x[0]-0) < tolerance_boundary_location )
def RightBndry(x):
    return ( abs(x[0]-length_0) < tolerance_boundary_location )
# Values of the concentrations field on this location
c_left = 50 
c_right = 100
BC_Dirichlet_left_value = Constant(c_left)
BC_Dirichlet_right_value = Constant(c_right)
# Apply over the finite element function space, the specified concentrations values, on the specified boundary
BC_Dirichlet_left = DirichletBC(V, BC_Dirichlet_left_value, LeftBndry)
BC_Dirichlet_right = DirichletBC(V, BC_Dirichlet_right_value, RightBndry)
# Total dirichlet boundaries
bc_dirichlet_tot = [BC_Dirichlet_left, BC_Dirichlet_right]

### Set Internal Neumann Boundary condition
# Define all the facets
all_facet_domains = FacetFunction("size_t", mesh)
# Set location of the bc 
class Neumann_internal_boundary(SubDomain):
        def inside(self, x, on_boundary):
            return abs(x[0]-1) < tolerance_boundary_location
# Mark the facets that satisfy the location condition           
Neumann_internal_boundary().mark(all_facet_domains, 1)
# Specify internal facets
dS_internal_neumann = dS(subdomain_data = all_facet_domains)
# Set flux value
flux_value = 1


### Write variational formulation

# Define trial function u on the finite element function space
u_trial = TrialFunction(V)
# Define test functions v on the finite element function space
v_test = TestFunction(V)

# Define function f
f = Constant(0)

# Integrands of a(u,v)
a = inner (nabla_grad(u_trial),nabla_grad(v_test))*dx

# Integrands of L(v)
L = f*v_test*dx

# Define the unknown function to be computed on the finite element function space
C = Function(V)

### Solve the PDE variational problem 
solve(a == L, C, bc_dirichlet_tot)

# Plot the solution
plot(C, title='Solution C', interactive=True)

Then, from this point, how can i add the internal Neumann Boundary condition ?
I believe i have correctly identified the internal factes with the keyword dS.
Then, i guess i have to change the variational formulation, but i don't know how to do it correctly. I have read in the Fenics book about some operator (jump) that can be used only with internal facets, but i don't understand how to use them, and even if i have to use them in my case.
I am trying to understand the example you notice in your answer, if i can use it in my problem, but if you know how to write corretly the internal, it would help me a lot.
Thanks.

PS: it seems in the example you sent me (Neumann BC depending on solution), that the jump operator is used only because he wants the flux at the internal equal to the local difference from the two domains at the interface, which is not my case : (i want to add on Γ: a∇u⋅n=k, and not =k[[u]], with his notation)

---

I have found in the forum, that this expression should be enough:

L = f*v*dx + g("+")*v("+")*dSs(0)

https://fenicsproject.org/qa/2744/how-to-make-interior-neumann-bc-dependent-on-fe-function

https://fenicsproject.org/qa/2532/how-to-set-up-interior-neumann-boundaries

I am trying it.