Conditional function in the variational formula through stress tensor

variational formulation
1

Hello!

As the title says, I get some trouble with the splitting function lhs() and rhs() when I try to split the followed variational formulation :

f_u = inner( sigma(u,g , G, K0,ndim) , epsilon(v) ) * dx + dot( tn, v ) * ds
a_u = lhs(f_u)
L_u = rhs(f_u)

with function define here :

# degradation function
def gk( alpha,k = 1e-6):
    return (1 - alpha )**2 + k

#shape function
def w(alpha):
    return alpha**2

def epsilon(u):
    return sym(nabla_grad(u))

def ElasticEnergyDensityPlus(u,  K0, G):
    return 0.5* K0*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2 + G *inner(dev(epsilon(u)),dev(epsilon(u)))

def sigma(u ,g , G, K0,ndim):
    """
   Ko = initial stiffness
   G = shear modulus
   g = degradation function  = gk()
   u = displacement
    """
    Kd = conditional(gt(tr(epsilon(u)),0),g*K0, K0)
    return (2. * Kd  * tr(epsilon(u))*Identity(ndim) + G *g * dev(epsilon(u)))

Through the stress expression, I try to define the compression/traction distinction of the material behavior. In function of the sign of tr(epsilon), the stiffness coefficient Kd is weigted by a degradation function g

When I run, I get this error : Found Argument in , this is an invalid expression.
Traceback (most recent call last):
File “…/Main_SinglNotc_E2.py”, line 216, in
a_u = lhs(f_u)

This is the beginning of my code :

from fenics import *
import numpy as np
from Function import *
import matplotlib.pyplot as plt

# ------------------
# Parameters
# ------------------
ERROR = 40
set_log_level(ERROR) # log level
parameters.parse()   # read paramaters from command line
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs" # quadrature is deprecated

"""
=========================================================
                        Time Discretization
=========================================================
"""

T = 1.0           # final time
num_steps = 1000    # number of time steps
dt = T / num_steps # time step size
time_step = np.linspace(0, 1.0, 300)

"""
=========================================================
                        Create mesh
=========================================================
"""

nx = 50
ny = 50
mesh = UnitSquareMesh(nx, ny,'crossed')

"""
=========================================================
                        Space discretization
=========================================================
"""
degreeU   = 1
degreeAlpha = 1
degreeH   = 0
V         = VectorFunctionSpace(mesh,'Lagrange', degreeU)
V_alpha       = FunctionSpace(mesh,'Lagrange', degreeAlpha)
V_H       = FunctionSpace(mesh,'DG', degreeH)
Tens      = TensorFunctionSpace(mesh, 'Lagrange', 1)


"""
=========================================================
    Define Trial and Test Function
    Define functions for solutions at current times steps
=========================================================
"""
u,v,u_            = TrialFunction(V), TestFunction(V), Function(V)
alpha,beta,alpha_ = TrialFunction(V_alpha), TestFunction(V_alpha), Function(V_alpha)


"""
===================================================================
                        Physical parameters
Parameters
___________________________________________________________________
E= Yound modulus            nu : poisson coef
mu=shear modulus            K0 : stiffness
mu, lmbda = lamé coeff      Gc : critical energy released rate J/mm²
lw= regularization parameter (internal length) mm

cw=normalized parameter (depends on the shape function w)
for w = phi1² cw = 2
    w = phi1  cw = 3/8
===================================================================
"""

E, nu     = Constant(1000), Constant(0.3) #[MPa]
G    = Constant(E / (2.0*(1.0 + nu))) # [MPa]
lmbda = Constant (E * nu / ((1.0 + nu)* (1.0 - 2.0*nu)))
K0, Gc    = Constant(lmbda+2./3.*G), Constant(0.3)
cw, lw    = Constant(2), Constant(0.1) # [mm]
Y_seuil = 0#Constant(Gc/(2*lw))

"""
=========================================================
                        Boundary condition
=========================================================
"""
#loading : linear displacement (ux,uy) = (0,t/4) 
Disp=Expression( ('0.0', 't/4.'),t = 0.0, degree = 1 )

# boolean boundary function
boundary_bottom  = 'near(x[1], 0)'
boundary_up      = 'near(x[1], 1)'
#Crach_boundary      = 'near(x[1], 0.5)' && 'between(x[0], (0., 0.3)) ' 
class Precrack(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (0., 0.3))  and near(x[1], 0.5) ) # and : &&

# Initialize mesh function for exterior domains / corresponding tags to the boundaries

#Up, bottom and crack domain
crack = Precrack()
bottom = CompiledSubDomain(boundary_bottom)
up = CompiledSubDomain(boundary_up)
boundaries = MeshFunction("size_t", mesh,1)
boundaries.set_all(0)
bottom.mark(boundaries, 1) # mark left as 1
up.mark(boundaries, 2) # mark right as 2
crack.mark(boundaries, 3) # mark right as 2
ds = Measure("ds",subdomain_data=boundaries) # left: ds(1), right: ds(2), crack : ds(3), ds(0) all the edges

bcu_bottom = DirichletBC( V , Constant((0.,0.)), boundary_bottom)
bcu_up     = DirichletBC( V , Disp, boundary_up )
bcu        = [bcu_bottom, bcu_up]

bcp_bottom = DirichletBC(V_alpha, Constant(0.), boundary_bottom)
bcp_up     = DirichletBC(V_alpha, Constant(0.), boundary_up)
bcp_crack     = DirichletBC(V_alpha, Constant(1.), boundaries,3)
bcp        = [bcp_bottom, bcp_up, bcp_crack]


"""
=========================================================
                        Define Variational problem
=========================================================
"""


# INITIAL CONDITION
u_n            = interpolate(Constant((0.,0.)),V)
alpha_n        = interpolate(Constant(0.), V_alpha)
H_n            = interpolate(Expression("0.0",degree=0), V_H)
H_n_prev       = interpolate(Expression("0.0",degree=0), V_H)

#---------------- Mechanical equation

g = gk(alpha_n,1.0e-6) # type : ufl.algebra.Sum
tn = Constant((0., 0.)) # external stress

def sigma(u ,g , G, K0,ndim):
    Kd = conditional(gt(tr(epsilon(u)),0),g*K0, K0)
    return  (2. * Kd * tr(epsilon(u))*Identity(ndim) + G *g * dev(epsilon(u)))

#--------------- Phase Field equation
f_u = inner( sigma(u,g , G, K0,ndim) , epsilon(v) ) * dx + dot( tn, v ) * ds
a_u = lhs(f_u)
L_u = rhs(f_u)

I know that the error comes from the conditional expression, but I don’t know what can I do to replace it.

Thanks for your help !

2

The problem here is that lhs and rhs are intended to be applied to residuals of linear problems, which are sums of forms of arity 1 (i.e., linear forms in a single TestFunction argument, which go on the right-hand side (RHS)) and forms of arity 2 (i.e., bilinear forms in a TrialFunction and a TestFunction, which go on the left-hand side (LHS)). Your formulation is nonlinear in the TrialFunction, so it is an invalid argument to lhs and rhs.

The way to formulate a nonlinear problem is to define a single residual form of arity 1 that is linear in a TestFunction, in which the unknown solution is a Function (not a TrialFunction). The (Gateaux) derivative of that form with respect to the solution Function, in the direction of an arbitrary TrialFunction, will then be a bilinear form of arity 2, which serves as the LHS in each step of Newton’s method. See the following demo for an example of solving a nonlinear problem:

https://fenicsproject.org/olddocs/dolfin/latest/python/demos/nonlinear-poisson/demo_nonlinear-poisson.py.html

1 Like
3

Thanks very much for the explication, it works well now !
But I wonder if I am dependent of nonlinear solver in this case or should I change the variational formulation to use a linear solver ?

4

If the physical phenomenon you’re modeling is inherently nonlinear, then there’s no way to capture that nonlinearity with a linear formulation. However, solution strategies for nonlinear formulations typically do solve a sequence of linear problems. The default behavior of solve(F==0,...) is to apply Newton’s method, where the linear problem to be solved in each step can be derived automatically using computer algebra in UFL (viz., the derivative function). However, in specific applications, you may instead want to implement some problem-specific solution algorithm yourself. For instance, it looks like you’re solving a phase field fracture problem, where some researchers alternate between solving elasticity with a fixed phase field and the phase field equation with a fixed displacement field (to reduce the size of the system of algebraic equations and/or to re-use existing solvers for elasticity and reaction–diffusion).

5

Indeed, I am agreed. Before getting into trouble with my code, I used the followed expression of the stress and the results were pretty good with a linear solver.

def sigma(u ,g, G, lmbda,ndim):
     return g * ( 2 * G * epsilon(u) + lmbda*tr(epsilon(u))*Identity(ndim))

To get linear problems, I used, as you said, the alternate minimization scheme (ie fixing one of the two variables for each problem). But I wanted to explore other splitting formulations of the elastic energy density that affects, in consequence, the stress expression too. And it seems that the first stress expression that generates an error makes the mechanical problem nonlinear due to the tr(epsilon(u)) in the conditional function if I am not wrong. If I assess the trace according to the previous time step, I could solve the problem with a linear solver.

def sigma(u_old, u ,g , G, K0,ndim):
    Kd = conditional(gt(tr(epsilon(u_old)),0),g*K0, K0)
    return  (2. * Kd * tr(epsilon(u))*Identity(ndim) + G *g * dev(epsilon(u)))
6

If I assess the trace according to the previous time step, I could solve the problem with a linear solver.

You can also do something similar as a fixed point iteration within each time step, which doesn’t introduce any splitting error, assuming the fixed point iteration converges for each step.

7

Yes, you right. Thanks for the advice !