Coefficient defined by two functions, one for each subdomain (with a piecewise function to identify domains)

Hello people,
I try to define a variable D so that D=f(x,C) with x the spatial location and C a field variable
So far, i have defined two different functions:
D = Da(C)
D = Db(C)
Also i have a piece wise function I(x) which is equal to 1 for cells within mesh A, and equal to 0 for cells within mesh B
then i do:
result = I(x)Da(C) + [1-I(x)]Db(C)

It works, except for nodes that belong to both mesh A and mesh B.

Here is a minimum code of my problem (Fenics 2016.1)

# Impose / is the decimal division (must be at the start of the code)
from __future__ import division

# Dolfin library
from dolfin import *
# Numpy library
import numpy as np

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

# Mesh
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(10, 10, 10), 20, 20, 20)

# Create subdomains object
subdomains_mesh = MeshFunction('size_t', mesh, 3)

# Defines subdomains
class OmegaA(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] <= 5 else False 
class OmegaB(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[0] >= 5 else False

# Mark subdomains
id_meshA = 10
id_meshB = 20
OmegaA().mark(subdomains_mesh, id_meshA)
OmegaB().mark(subdomains_mesh, id_meshB)

# Field that indicates if the cell are within the mesh A or within the mesh B
# Create function space
V_identify = FunctionSpace(mesh, 'DG', 0)
# Then, we define a function on this finite element function space
Identify_=Function(V_identify)
# Values are set (within the mesh A =1, within the mesh B =0)
Localtion_ = np.zeros(1000)
Localtion_[id_meshA]=1
Localtion_[id_meshB]=0
# Attribute value locally on each cell
# We perform a loop on each cells of the all the subdomain
for cell_number in range(len(subdomains_mesh.array())): # Select cell
    subdomain_number = subdomains_mesh.array()[cell_number] # Get the id of the subdomain that contains the cell
    Identify_.vector()[cell_number]=Localtion_[subdomain_number] # Attribute function value on each cell (piecewise constant function)
plot(Identify_, title='Within meshA (=1) or meshB (=0)', interactive=True)  

# Define functions: value of D we want for the mesh A and for the mesh B
def D_A(C):
    D=C
    return(D) 
def D_B(C):
    D=C/2
    return(D)

# Define function: value of D that change depending on if the cell are within mesh A or mesh B
def DA_DB (Identify, C):
    # Calculate the two variable field
    field_DA = D_A(C)
    field_DB = D_B(C)
    # Attribute spatially C
    # Identify = 1 if within mesh A, =0 if within mesh B
    # So that, DD = D_A(C) for cell within mesh A
    # So that, DD = D_B(C) for cell within mesh B
    # AND FOR NODES THAT BELONGS BOTH TO MESH A AND B ? HERE IS MY ISSUE.
    # I WOULD LIKE DD = (D_A(C)+D_B(C))/2 for nodes at the interface between mesh A and B
    DD = Identify*field_DA + (1-Identify)*field_DB
    # Calculate the flux
    return DD       

# Create variable field C
C_scalar=10

# Set function space for the variable field C
V = FunctionSpace(mesh, 'Lagrange', 1)

# Project C on the function space
C_field = project(C_scalar,V)
plot(C_field, title='C',interactive=True) # Value display is correct

# Calculate
field_DA = D_A(C_field)
field_DB = D_B(C_field)
plot(field_DA, title='D with expression we want only for mesh A',interactive=True) # Value display is correct
plot(field_DB, title='D with expression we want only for mesh B',interactive=True) # Value display is correct

# Calculate
DD_ = DA_DB(Identify_,C_field) # WRONG
plot(DD_, title='D with expression both for mesh A and B',interactive=True) # Value display is incorrect at the interface

Any help would be greatly appreciate.

Okay, i have found another way to solve my problem

Instead of using a piece wise function that FEniCS (and me) doesnt hand well, i have done this:

# Get spatial coordinnates
x_mod = SpatialCoordinate(mesh)

# Create x0 function
x0_coordinate_function = project(x_mod[0],V)
plot(x0_coordinate_function, title='x0_coordinate_function',interactive=True)
# The function is converted in an array
x0_array = x0_coordinate_function.vector().array()
# Apply a mask
np.place(x0_array, x0_array<=4.9, 1) # half the length minus a tolerance (to avoid to be on a node)
np.place(x0_array, x0_array>1, 0)
# Insert the new array in the fuction
x0_coordinate_function.vector()[:] = x0_array   

# Define function: value of D that change depening on if the cell are within mesh A or mesh B
def DA_DB_x (C, x):
    # Calculate the two variable field
    field_DA = D_A(C)
    field_DB = D_B(C)
    # Attribute spatially C
    DD = x*field_DA + (1-x)*field_DB
    # Calculate the flux
    return DD       

# Calculate
DD_ = DA_DB_x(C_field,x0_coordinate_function)
plot(DD_, title='D with expression both for mesh A and B',interactive=True) # Value display is incorrect at the interface

Although, if you know to do it with the piece wise function...