Three-dimensional limit analysis problem using Semi-Definite Programming

In this demo, we consider a classical limit analysis problem, namely a slope stability problem for a cohesive-frictional material described by a Mohr-Coulomb criterion. The geometry being three-dimensional in this example, the corresponding problem will be a Semi-Definite Programming (SDP) problem. We show how to formulate such constraints using the fenics_optim package.

In the following, we will denote by \(\mathbb{S}_n\) the set of symmetric \(n\times n\) matrices.

Problem formulation

We consider a soil domain \(\Omega = [0;L]\times [0;W] \times [0;H]\) with homogeneous Dirichlet boundary conditions \(\boldsymbol{u}=0\) on the right \(x=L\) and bottom \(z=0\). The remaining boundaries have homogeneous Neumann boundary conditions. The loading consists of a graviatational body force \(\boldsymbol{f}=(0,0,-\gamma)\) with \(\gamma\) being the soil self-weight. The soil obeys a Mohr-Coulomb criterion of cohesion \(c\) and internal friction angle \(\phi\), i.e. the stress state \(\boldsymbol{\sigma}\in \mathbb{S}_3\) must statisfy \(\boldsymbol{\sigma}\in G\) where:

\begin{equation} G=\left\{\boldsymbol{\sigma}\in \mathbb{S}_3 \text{ s.t. } \sigma_M - a\sigma_m \leq \dfrac{2c\cos\phi}{1+\sin\phi}\right\} \end{equation}

where \(a=\dfrac{1-\sin\phi}{1+\sin\phi}\), \(\sigma_M = \max_{I} \{\sigma_I\}\) and \(\sigma_m = \min_I \{\sigma_I\}\) with \(\sigma_I\) being the eigenvalues of \(\boldsymbol{\sigma}\).

The limit analysis problem amounts to finding the slope stability factor given by \(SF=\lambda^+\dfrac{\gamma H}{c}\) where \(\lambda^+\) is obtained from solving:

\begin{equation} \begin{array}{rl} \displaystyle{\lambda^+ = \inf_{\boldsymbol{u}}} & \displaystyle{\int_\Omega \pi(\nabla^s \boldsymbol{u}) \text{ dx}}\\ \text{s.t.} & \displaystyle{\int_{\Omega} \boldsymbol{f}\cdot\boldsymbol{u} \text{ dx} = 1} \end{array} \end{equation}

in which \(\nabla^s \boldsymbol{u} = \frac{1}{2}(\nabla \boldsymbol{u} + \nabla \boldsymbol{u}^T)\) is the symmetric gradient and \(\pi\) is the support function of the convex set \(G\):

\begin{align*} \pi(\boldsymbol{d}) &= \sup_{\boldsymbol{\sigma}\in G} \{\boldsymbol{\sigma}:\boldsymbol{d}\}\\ &= \begin{cases} c \cot\phi \operatorname{tr}(\boldsymbol{d}) & \text{if } \displaystyle{\operatorname{tr}(\boldsymbol{d}) \geq \sin\phi\left(\sum_I |d_I|\right)} \\ +\infty & \text{otherwise} \end{cases} \end{align*}

Conic reformulation

Following [MAR08], the above support function can be expressed equivalently in a conic-representable fashion using two auxilary SDP variables \(\boldsymbol{Y}_1,\boldsymbol{Y}_2\) as follows:

\begin{equation} \begin{array}{rl} \displaystyle{\pi(\boldsymbol{d}) = \inf_{\boldsymbol{Y}_1,\boldsymbol{Y}_2\in \mathbb{S}_3}} & \dfrac{2c\cos\phi}{1+\sin\phi}\operatorname{tr}(\boldsymbol{Y}_1)\\ \text{s.t.} & \boldsymbol{d} = \boldsymbol{Y}_1 - \boldsymbol{Y}_2\\ & a\operatorname{tr}(\boldsymbol{Y}_1)=\operatorname{tr}(\boldsymbol{Y}_2)\\ & \boldsymbol{Y}_1 \succeq 0, \boldsymbol{Y}_2\succeq 0 \end{array} \end{equation}

Implementation

We first define the conic representation of the Mohr-Coulomb support function:

[1]:

from dolfin import (
    BoxMesh,
    Point,
    near,
    dx,
    Constant,
    DirichletBC,
    FunctionSpace,
    VectorFunctionSpace,
    XDMFFile,
)
from ufl import sym, dot, grad, pi, sin, cos, tr
from fenics_optim import MosekProblem, ConvexFunction, to_mat, to_vect, SDP

c = Constant(1.0)
phi = Constant(pi / 6.0)


class MohrCoulomb(ConvexFunction):
    """SDP implementation of Mohr-Coulomb criterion."""
    def conic_repr(self, X):
        Y1 = self.add_var(6, cone=SDP(3))
        Y2 = self.add_var(6, cone=SDP(3))
        a = (1-sin(phi))/(1+sin(phi))
        self.add_eq_constraint(X - Y1 + Y2)
        self.add_eq_constraint(tr(to_mat(Y2))-a*tr(to_mat(Y1)))
        self.set_linear_term(2*c*cos(phi)/(1+sin(phi))*tr(to_mat(Y1)))

In the above, symmetric \(n\times n\) matrix variables are represent using a vector of \(n(n+1)/2\) values of the upper-triangular part. The SDP constraint is enforced through the cone SDP(3) in which it is check that the dimension \(n=3\) of the cone matches the dimension of the vector \(n(n+1)/2=6\). The to_mat function enables to give the matrix representation from its vector representation, enabling to use predefined UFL operators such as the trace tr.

We now define the mesh, loading, function spaces and boundary conditions:

[2]:

L, W, H = (1.2, 2., 1.)
Nx, Ny, Nz = (20, 1, 20)
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, H), Nx, Ny, Nz)

gamma = 10.
f = Constant((0, 0, -gamma))

def border(x, on_boundary):
    return near(x[0], L) or near(x[2], 0)


V = VectorFunctionSpace(mesh, "CG", 2)
R = FunctionSpace(mesh, "R", 0)

bc = DirichletBC(V, Constant((0.,0.,0.)), border)

Note that we used on purpose only 1 element in the \(y\)-direction with quite a large width in order to reproduce a 2D plane-strain situation for which we have a good approximation of the exact solution.

We now initiate the MosekProblem object and first add the linear equality constraint:

[3]:

prob = MosekProblem("3D limit analysis")
u = prob.add_var(V, bc=bc, name="Collapse mechanism")

def Pext(lamb):
    return lamb*dot(f,u)*dx
prob.add_eq_constraint(R, A=Pext, b=1)

We now add the convex term corresponding to the support function. Note that we use the to_vect function to obtain the vector representation of the symmetric gradient tensor. We select a vertex integration scheme ensuring an upper bound status of the computed stability factor (see [MAK07]).

[4]:

crit = MohrCoulomb(to_vect(sym(grad(u))), quadrature_scheme="vertex")
prob.add_convex_term(crit)

The problem can then be solved and results are exported to Paraview.

[5]:

prob.optimize()

print("-------------------------------------")
print("2D plane strain factor [Chen] (for phi=30°):", 6.69)
print("Computed factor:", prob.pobj*float(gamma*H/c))

with XDMFFile("slope_3D.xdmf") as ffile:
    ffile.write(u)

Matrix shape: (67930, 15129)
Number of cones: 0
Problem
  Name                   :
  Objective sense        : min
  Type                   : CONIC (conic optimization problem)
  Constraints            : 67930
  Cones                  : 0
  Scalar variables       : 15129
  Matrix variables       : 19200
  Integer variables      : 0

Optimizer started.
Presolve started.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00
Lin. dep.  - tries                  : 0                 time                   : 0.00
Lin. dep.  - number                 : 0
Presolve terminated. Time: 0.06
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
  Name                   :
  Objective sense        : min
  Type                   : CONIC (conic optimization problem)
  Constraints            : 67930
  Cones                  : 0
  Scalar variables       : 15129
  Matrix variables       : 19200
  Integer variables      : 0

Optimizer  - threads                : 2
Optimizer  - solved problem         : the primal
Optimizer  - Constraints            : 67201
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 14401             conic                  : 14401
Optimizer  - Semi-definite variables: 19200             scalarized             : 115200
Factor     - setup time             : 7.06              dense det. time        : 1.21
Factor     - ML order time          : 0.65              GP order time          : 3.98
Factor     - nonzeros before factor : 5.46e+06          after factor           : 6.16e+07
Factor     - dense dim.             : 2                 flops                  : 1.03e+11
Factor     - GP saved nzs           : 1.56e+07          GP saved flops         : 8.03e+10
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME
0   6.4e+01  1.0e+00  9.3e+00  0.00e+00   8.313843876e+00   0.000000000e+00   1.0e+00  7.31
1   6.0e+01  9.4e-01  4.2e+01  -2.01e+01  1.995899916e+01   2.161828566e+01   9.4e-01  12.63
2   5.5e+01  8.5e-01  6.5e+01  -3.86e+01  1.518431835e+00   6.617098782e+01   8.5e-01  17.53
3   4.2e+01  6.6e-01  6.3e+01  -2.09e+00  1.978801416e+00   1.064313946e+02   6.6e-01  22.32
4   3.0e+01  4.7e-01  5.5e+01  -1.73e+00  2.721072587e+00   1.556825262e+02   4.7e-01  27.62
5   1.4e+01  2.1e-01  2.9e+01  -1.17e+00  5.214853147e+00   2.115635460e+02   2.1e-01  32.48
6   3.1e+00  4.9e-02  4.0e+00  7.65e-02   1.282676009e+01   8.870108609e+01   4.9e-02  37.24
7   4.3e-01  6.6e-03  8.6e-02  9.97e-01   6.985128550e+00   8.879161141e+00   6.6e-03  43.68
8   1.9e-01  2.9e-03  4.2e-03  4.94e+00   1.667969811e+00   1.687221235e+00   2.9e-03  48.80
9   3.4e-02  5.3e-04  1.7e-04  3.11e+00   7.794060738e-01   7.802760766e-01   5.3e-04  54.46
10  7.3e-03  1.1e-04  8.5e-06  1.20e+00   7.225951896e-01   7.225976589e-01   1.1e-04  60.08
11  4.1e-03  6.4e-05  3.3e-06  1.05e+00   7.152058297e-01   7.152018064e-01   6.4e-05  64.49
12  2.3e-03  3.7e-05  1.3e-06  1.03e+00   7.110184087e-01   7.110128366e-01   3.7e-05  69.16
13  1.1e-03  1.7e-05  3.4e-07  1.02e+00   7.080158146e-01   7.080117406e-01   1.7e-05  73.64
14  1.6e-04  2.5e-06  2.5e-09  1.01e+00   7.058644281e-01   7.058630877e-01   2.5e-06  79.36
15  1.3e-04  2.0e-06  1.8e-09  1.00e+00   7.057887743e-01   7.057877034e-01   2.0e-06  84.07
16  5.2e-05  8.1e-07  4.3e-10  1.00e+00   7.056064811e-01   7.056060520e-01   8.1e-07  89.40
17  6.1e-06  1.2e-07  2.3e-11  1.00e+00   7.055029851e-01   7.055029206e-01   1.2e-07  94.69
Optimizer terminated. Time: 94.81


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 7.0550298511e-01    nrm: 1e+01    Viol.  con: 6e-09    var: 0e+00    barvar: 0e+00
  Dual.    obj: 7.0550292062e-01    nrm: 7e+02    Viol.  con: 0e+00    var: 3e-10    barvar: 7e-09
-------------------------------------
2D plane strain factor [Chen] (for phi=30°): 6.69
Computed factor: 7.055029851148466

[6]:

print("2D factor [Chen] (for phi=30°):", 6.69)
print("Computed factor:", prob.pobj*float(gamma*H/c))

2D factor [Chen] (for phi=30°): 6.69
Computed factor: 7.055029851148466

References

[MAK07] Makrodimopoulos, A., & Martin, C. M. (2007). Upper bound limit analysis using simplex strain elements and second‐order cone programming. International journal for numerical and analytical methods in geomechanics, 31(6), 835-865.

[MAR08] Martin, C. M., & Makrodimopoulos, A. (2008). Finite-element limit analysis of Mohr–Coulomb materials in 3D using semidefinite programming. Journal of engineering mechanics, 134(4), 339-347.

Cheeger problems Mixed 1D/2D modelling of a reinforced slope stability

fenics_optim 2.0.3 documentation

Three-dimensional limit analysis problem using Semi-Definite Programming

Problem formulation

Conic reformulation

Implementation

References