Discontinuous type Argument must be restricted

Errors
1

Hello, I’m working with an internal boundary integration as the right side of a variational system:

V = FunctionSpace(mesh, 'CG', 1)
n = FacetNormal(mesh)
S = 2.0*mu*inner(e(u),e(u)) + lmbda*div(u)**2

theta,xi = [TrialFunction(V), TestFunction(V)]     

left_side = assemble((inner(grad(theta),grad(xi)) + inner(theta,xi))*dX)
solverav = LUSolver(left_side)
right_side = assemble(inner(jump(S)*n('+'),n('+'))*xi('+')*dS_int + Constant(0)*xi*dx(domain=mesh, subdomain_data=domains))

th = Function(V)                  
solverav.solve(th.vector(), right_side)

Here u is the solution of a linear elasticity problem and e(u) is the strain. I’m aware that e(u) is discontinuous, since involves \nabla u, and also that its values are associated with each cell. Therefore it makes sense to me that when integrating over facets I need to restrict this expression to either the positive or negative side of the cell. However, when running the code I still get the error “Discontinuous type Argument must be restricted”, so I understand that I have also to restrict my test function. But I don’t get why. By restricting it to (’+’) or (’-’) won’t It affect the solution of the problem?

Could anyone give me some help on this? I’ve seen a lot of posts here that helped me understanding the restrictions, but I don’t get this issue with the test function.
Thank you!

2

If you make this a MWE, we can help debug it.

3

Hi Nate! Here it is. Thank you!

from dolfin import *
import numpy as np

lx,ly = [2.0,1.0]; Nx, Ny = [62, 31]  

def phi(x):
    return (x[0] - 1.0)**2 + (x[1]-0.5)**2 - 0.2**2

class Omega(SubDomain):
    def inside(self, x, on_boundary):
        return phi(x) > 0

mesh  = RectangleMesh(Point(0.0,0.0),Point(lx,ly),Nx,Ny,'crossed') 
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dX = Measure('dx')
omega = Omega() 
omega.mark(domains, 1)
dx = Measure("dx", domain=mesh, subdomain_data=domains)
V = FunctionSpace(mesh, 'CG', 1)
Vvec = VectorFunctionSpace(mesh, 'CG', 1)

# Boundary Conditions
class DirBd(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],.0)     
dirBd = DirBd()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)    
dirBd.mark(boundaries, 1) 
ds = Measure("ds")(subdomain_data=boundaries) 
bcd  = [DirichletBC(Vvec, (0.0,0.0), boundaries, 1)] 
Load = Point(lx, 0.5)

eps_er, E, nu = [0.001, 1.0, 0.3] 
mu,lmbda = Constant(E/(2*(1 + nu))),Constant(E*nu/((1+nu)*(1-2*nu)))
lmbda = Constant(2*mu*lmbda/(lmbda+2*mu)) #plane stress       

# Solving Elasticity
u,v = [TrialFunction(Vvec), TestFunction(Vvec)]
S1 = 2.0*mu*inner(sym(grad(u)),sym(grad(v))) + lmbda*div(u)*div(v)
A = assemble( S1*eps_er*dx(0) + S1*dx(1) )   
b = assemble( inner(Expression(('0.0', '0.0'),degree=2) ,v) * ds(2))    
U = Function(Vvec)
delta = PointSource(Vvec.sub(1), Load, -1.0)
delta.apply(b) 
for bc in bcd: bc.apply(A,b)    
solver = LUSolver(A)
solver.solve(U.vector(), b);  

# Marking Internal Boundary
int_boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
tdim = mesh.topology().dim()
mesh.init(tdim-1, tdim)
facet_to_cell = mesh.topology()(tdim-1, tdim)
num_facets = facet_to_cell.size()
domain_values = domains.array()
facet_values = int_boundary.array()
for facet in range(num_facets):
    # Check if interior
    cells = facet_to_cell(facet) 
    if len(cells == 2):
        # Check if facet is on the boundary between two cells with different markers
        values = np.unique(domain_values[cells])
        if len(values) == 2:
            facet_values[facet] = 1        
    else:
        continue

dS_int = Measure("dS", domain=mesh, subdomain_data=int_boundary) 
File("facets.pvd") << int_boundary

# Variational problem of interest (where the error appears)
n = FacetNormal(mesh)
S = 2.0*mu*inner(sym(grad(U)),sym(grad(U))) + lmbda*div(U)**2
theta,xi = [TrialFunction(V), TestFunction(V)]   
left_side = assemble((inner(grad(theta),grad(xi)) + inner(theta,xi))*dX)
solverav = LUSolver(left_side)
right_side = assemble(inner(jump(S)*n('+'),n('+'))*xi*dS_int + Constant(0)*xi*dx(domain=mesh, subdomain_data=domains))
th = Function(V)                  
solverav.solve(th.vector(), right_side)
4

You need to choose on which side of the facet to evaluate xi, e.g. xi("+").

5

Thank you! This works, but I’m still a bit confused on what this restriction means, since I’m interested in the solution th defined on the whole domain. If I choose (’+’), will I somehow be concentrating the solution on “this side”?

Also, could you maybe give me some insight on the element dS? Is it an oriented element? I mean, does it have a normal vector “included” already?