Integration on interior boundary

Hi,
I am solving a solid mechanics problem - deflection of cantilever beam with a membrane embedded. The problem is as shown in figure below,

Problem definition"


The grey colored block is the solid, and the membrane is along the green plane. I impose displacement = 0 bc at the left face, and a neuman traction bc on the top face of block (load).

Governing equation:
The standard solid mechanics equation for this problem would look like, (at steady state)
image + load term
Now, I will add the membrane contribution as follows, (adding red term)
image + load term
where FS = deformation gradient*2nd Piola Kirchoff stress = First Piola Kirchoff stress. In my case I have used Neo Hookean material model to derive stresses.

What I tried:
Initially I solved the simple solid problem without membrane and found no issue. Now, when i add the membrane term (red) I am facing issues. For debugging purpose, I have made the membrane term to be ‘zero’ as following,

I = Identity(dim)
K_m = Constant(0.0)*I
res_solid = inner(F*S,grad(du))*dx - inner(f,du)*ds(24) + inner(K_m,grad(du)("+"))*dS

In the above line, since I made K_m = zero, i expect the results to be same as standard problem, as if there is no membrane. But, I get different results (I am looking at the LHS matrix of assembled form for comparison). I do not understand what is happening here.
Can anyone help me with this?

Thanks.

MWE:

from dolfin import *
import numpy as np
from scipy.sparse import csr_matrix
from mshr import *

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(16,16,500), 2, 4, 10)

# Define subdomains
class Domain1(SubDomain):
    def inside(self, x, on_boundary):
        # Assuming the interface is at the halfway point in the x-direction
        return x[1] <= 8

class Domain2(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 8

domain1 = Domain1()
domain2 = Domain2()

# Mark subdomains with different markers
domains = MeshFunction('size_t', mesh, mesh.topology().dim())
domains.set_all(0)
domain1.mark(domains, 1)
domain2.mark(domains, 2)

# Define and mark boundaries
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 0)

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 16)

class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 8)

left_boundary = LeftBoundary()
top_boundary = TopBoundary()
interface = Interface()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left_boundary.mark(boundaries, 23)  # Mark the left boundary as 23
top_boundary.mark(boundaries, 24)   # Mark the top boundary as 24
interface.mark(boundaries, 25)      # Mark the interface as 25

# Define measures for integrating over subdomains and boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=boundaries)  # for interfaces

nsd = mesh.geometry().dim()
I = Identity(nsd)

# define function space
V = VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)
du = TestFunction(V)

# material properties
nu = Constant(0.3)      # poisson ratio
mu = Constant(3.84615)      # shear modulus
E = 2*mu*(1+nu)     # youngs modulus
kappa = E*nu/((1+nu)*(1-2*nu)) + (2/3)*E/(2*(1 + nu))
rho = Constant(pow(10,-6))

# solid model
F = I + grad(u)
J = det(F)
C = F.T*F
E = 0.5*(C-I)
f = Constant((0,0.00005,0)) # traction force on top face
S = mu*(J**(-2/3))*(I-(tr(C) + (3-nsd))/3*inv(C)) + kappa/2*((J**2)-1)*inv(C) # neo-hookean

# Boundary conditions
bcs = [
    DirichletBC(V, Constant((0,0,0)), boundaries, 23),
]

# define the variational problem
res_solid = inner(F*S,grad(du))*dx - inner(f,du)*ds(24)
K_m = Constant(0.0)*I
res_solid = res_solid + inner(K_m,grad(du)("+"))*dS(25)    # line i am talking about
Dres_solid = derivative(res_solid,u)  # Jacobian

# problem 
problem_m = NonlinearVariationalProblem(res_solid,u,bcs,Dres_solid)
solver_m = NonlinearVariationalSolver(problem_m)
solver_m.parameters['newton_solver']\
                   ['maximum_iterations'] = 15
solver_m.parameters['newton_solver']\
                   ['relative_tolerance'] = 1e-6
solver_m.parameters['newton_solver']['linear_solver'] = 'mumps'

#  following lines are for debugging
A, b = assemble_system(Dres_solid, -res_solid, bcs)
A_csr = csr_matrix(A.array())
np.savetxt("AK.csv", A_csr.toarray(), delimiter=",")
b_np = b.get_local()
np.savetxt("bK.csv", b_np, delimiter=",")

# solve
solver_m.solve()  

# output u
ufile = File("u.pvd")
ufile << u

Note: Paraview files might visually look same for both cases, but when i look into my matrices, I can see significant differences.

It is unclear to me how you added that term to your variational formulation.
This term should contribute to cells on either side of the membrane, not just a single side.

As a side note the residual that you are assembling should not touch the lhs at all, as it is not a function of u. Should start by looking at the solution vector, as it has significant differences between using the term and not using it.

Thanks for the reply @dokken
I found something interesting. When i said the results are different, I have also looked at the difference. It seems like only the DOFs corresponding to Dirichlet BC nodes are affected. This is the reason why my final results are unchanged.
I am not sure how fenics is assembling my matrix, and where this difference is coming from. And should I bother about it or not?

What do you meant by,

This term should contribute to cells on either side of the membrane, not just a single side.

Is it because i do not have the dummy term (Constant(0.0)*dx) that guides fenics with “+” and “-” sides? The issue remains same even with the dummy term.

In physical sense, I would like to stiffen the solid only at the interface. Hence I add the membrane term (red).
I derive a new stress term (red term) for membrane that has higher stiffness. The derivation of this new term involves 2D parametrization of the interface and hence usage of FacetJacobians. So I must restrict side of this Jacobian.
All of this did not work! Hence I am trying to debug it by making the additional term be zero.

Just adding a term, not stemming from the variational formulation might not be a good idea.

That means that the rhs and lhs is scaled consistently in both cases I guess. This might be due to how legacy dolfin sets the diagonal.