Conditional of jump has no effect?

variational formulation
1

I’m trying to solve a problem with DG-method this includes a conditional depending on the vectorfunction u:

def tcohesive(u):
dn = jump(u,n)
tn = conditional(lt(dn,0.002) , 1.0 , 0.0)
return as_vector([0.0,tn]

n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, ‘DG’, 1)
v = TestFunction(V)
u_ = Function(V)

F = inner(sigma(u_),epsilon(v))*dx + dot(jump(v),tcohesive(u_))*dS
J = derivative(F, u_)
problem = NonlinearVariationalProblem(F, u_, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()

The problem is that the conditional never changes to its False condition. Being new to FeniCS, I appreciate any help.

2

Could you take some time to carefully format the code in your question correctly, along with completing it such that it provides a minimal working example so that we may debug?

1 Like
3

I’m really sorry. I hope this serves the purpose:

#!/usr/bin/python2
from dolfin import *
tol     = DOLFIN_EPS*1E3
youngs  = 10E3
nu	    = 0.3
##################################
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return ( (abs(x[1])  < tol or abs(x[1]-1)  < tol) and on_boundary)   
##################################
class InternalEdges(SubDomain):
    def inside(self, x, on_boundary):
        return (not on_boundary)
##################################
class InterfaceEdges(SubDomain):
    def inside(self, x, on_boundary):
         return (abs(x[1]-0.5)<tol)
##################################
def sigma(u):
     return (youngs/(1.0-nu*nu))*((1.0-nu)*sym(grad(u))+nu*Identity(2)*tr(sym(grad(u))))
##################################
def tcohesive(x):
    dn = jump(x,n)
    tn = conditional(lt(dn,0.002),1.0,0.0)*dn
    return as_vector([0.0,tn])
#################################

mesh		= UnitSquareMesh(10,10)

boundaries  = MeshFunction('size_t', mesh,mesh.topology().dim()-1)
DirichletBoundary().mark(boundaries, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

internals = MeshFunction('size_t', mesh,mesh.topology().dim()-1)
InternalEdges().mark(internals,2)
InterfaceEdges().mark(internals,3)
dS = Measure('dS', domain=mesh, subdomain_data=internals)

# Some geometric entities
n     = FacetNormal(mesh)
h     = CellSize(mesh)
h_avg = 0.5*(h('+')+h('-'))
alpha = 10.0*youngs/(4.0*h_avg*(1.0+nu)*(1.0-2.0*nu))

# Function spaces and functions
U           = VectorElement('DG',mesh.ufl_cell(),1)
V           = VectorFunctionSpace(mesh, 'DG', 1)
v   		= TestFunction(V)
u_  		= Function(V)

bcs         = [DirichletBC(V,Expression(('0.0','x[1]*0.004'), degree=2),DirichletBoundary(),"geometric")]

#bilinear form. the problem is in tcohesive(u_)
aVol        	= inner(sigma(u_),sym(grad(v)))
aBulk       	= -inner(outer(jump(v),n('+')),avg(sigma(u_)))-inner(outer(jump(u_),n('+')),avg(sigma(v)))+ alpha*dot(jump(u_),jump(v))
aCohesive   	= dot(jump(v),tcohesive(u_))

F           = aVol*dx+aBulk*dS(2)+aCohesive*dS(3)
J  		    = derivative(F, u_) 
problem 	= NonlinearVariationalProblem(F, u_, bcs, J)
solver  	= NonlinearVariationalSolver(problem)
solver.solve()
T   = TensorFunctionSpace(mesh,'DG',1)
s   = project(sigma(u_),T)
File('./Stress.pvd')<<s

The example is rather long, because its a pretty special problem. To narrow it down: The problem clearly lies in the conditional in the function tcohesive. When I replace it with a simple mathematical expression (for example sin(dn)) the result is as expected. But when I use the conditional the result has nothing to do with the values I plugin or expect. The rest of the code works fine with no error messages.

1 Like
4

One thing I’ll mention is that, based on the hashbang line specifying python2 (which has not been supported in the last few FEniCS releases), you appear to be using an old version of the software. I tried to reproduce the issue with a more minimal example, and it actually fails with a compilation error when I force it to use a deprecated form compiler (commented line near top) while working as expected with the default UFLACS:

from dolfin import *

# Compilation fails with an error using the old "quadrature"
# representation in FFC, which is deprecated.

#parameters["form_compiler"]["representation"] = "quadrature"

# Two element mesh:
mesh = UnitSquareMesh(1,1)
n = FacetNormal(mesh)

# Get a vector function f whose 0-th component is
# 1 on one element and zero on the other.
V = VectorFunctionSpace(mesh,"DG",1)
x = SpatialCoordinate(mesh)
f0 = conditional(lt(x[1],x[0]),1,0)
f = project(as_vector([f0,0]),V)

# One should be 0, the other should be sqrt(2); which is
# which depends on the arbitrary normal orientation.
print(assemble(conditional(lt(jump(f,n),0),1,0)*dS))
print(assemble(conditional(gt(jump(f,n),0),1,0)*dS))

("tsfc" representation works correctly as well, but that is not part of standard FEniCS distributions.)

3 Likes
5

Thank you very much for your efforts and answer,
I have to use python2 on the pc at work but I use python 3 and fenics 2019.2.0.dev0 on my private mashine resulting in the same problem. Even adding

 parameters["form_compiler"]["representation"] = "uflacs"

had no effect. Your example works fine. Could the problem lie in the combination of conditional and
NonlinearVariationalSolver or derivative?

6

Not a solution to the original problem, but an addendum regarding the old FFC (non quadrature / non UFLACS) representation. There’s a module level variable that provides a workaround. I think UFLACS became the default in the 2017.0.1 version so it’s long been unnecessary, but I’ll add this here for anyone who absolutely must use very old versions of FEniCS.

import ufl
ufl.algorithms.apply_derivatives.CONDITIONAL_WORKAROUND = True
1 Like
7

Looking at your formulation, it looks like you’re trying to model some kind of contact/fracture problem? If I play around with the parameters and correct for the lack of coercivity on the boundary dS(3) I get reasonable results (in the eyeball norm). What is it that you’re expecting which you’re not seeing? Perhaps you could link us to the model equations you’re trying to discretise?

1 Like
8

You are absolutely right. It’s a test case for a cohesive zone model. The conditional I used here is pretty extreme because I wanted to see if the conditional works. What I expect is a change in the stress in vertical direction (s22)

  T = TensorFunctionSpace(mesh , 'DG' , 1)
  s22 = project(sigma(u_),T)(0.5,0.25)[3]

when the jump in u passes the threshold. But no matter how large I make the bc

bcs         = [DirichletBC(V,Expression(('0.0','x[1]*0.004'), degree=2),DirichletBoundary(),"geometric")]

it never happens.

I expect the red line but I get the blue.

path12472
path124725492×4458 267 KB

9

I think I’ve got it! The definition of sides (’+’ and ‘-’) and therefore the definition of jump is different then it is usually in cohesive zone modelling. In fact it is turned around. So if I put

  dn = -jump(x,n)

instead of

 dn = jump(x,n)

it works as expected. Sorry for keeping you busy with such a trivial mistake. Thank you very much for your efforts and kind help.

1 Like