This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Solid Mechanics, Topology Optimization - Local minimum or error?

+3 votes

Hello!

I am still beginning with FEniCS, so I am sorry if my questions are too simple, but I can't figure out what is wrong.

I am working with Optimization of Solid Mechanics (using IPOPT). I wrote a simple code to minimize compliance of a cantilever beam subject to distributed load. It seems to be working, but the results doesn't seem to be "complete". I believe 2 things can be happening:

1) I am reaching a local minimum
2) My code is wrong

Can anybody help?

See the code below:

from mshr import *
from dolfin import *
from dolfin_adjoint import *
import pyipopt

## Geometry and elasticity
t, h, L = 2., 1., 5.                          # Thickness, height and length
Rectangle = Rectangle(Point(0,0),Point(L,h))  # Geometry
E, nu = 210*10**9, 0.3                        # Young Modulus
G = E/(2.0*(1.0 + nu))                        # Shear Modulus
lmbda = E*nu/((1.0 + nu)*(1.0 -2.0*nu))       # Lambda

## SIMP
V = Constant(0.4) # Volume constraint
p = 3             # Exponent
eps = Constant(1.0e-3)

## Mesh, Control and Solution Spaces
mesh = generate_mesh(Rectangle, 100) # Mesh
def plot_mesh():                     # Mesh visualization
    plot(mesh, axes = True)
    return

A = VectorFunctionSpace(mesh, "Lagrange", 1) # Displacements
C = FunctionSpace(mesh, "Lagrange", 1)       # Control

## Volumetric Load
q = -1000000.0/t
b = Constant((0.0, q))

## Dirichlet BC em x[0] = 0
def Left_boundary(x, on_boundary):
    return on_boundary and abs(x[0]) < DOLFIN_EPS
u_L = Constant((0.0, 0.0))
bc = DirichletBC(A, u_L, Left_boundary)

## Forward problem solution
def forward(rho_opt):
    u = TrialFunction(A)  ## Trial and test functions
    w = TestFunction(A)
    sigma = lmbda*tr(grad(u))*Identity(2) + 2*G*sym(grad(u)) ## Stress
    F = (eps+(1-eps)*rho_opt**p)*inner(sigma, grad(w))*dx - dot(b, w)*dx
    a, L = lhs(F), rhs(F)
    u = Function(A)
    problem = LinearVariationalProblem(a, L, u, bcs=bc)
    solver = LinearVariationalSolver(problem)
    solver.parameters["symmetric"] = True
    solver.solve()
    return u

## Plot displacements
def plot_u():
    plot(u, mode = "displacement")
    interactive()
    return

## MAIN
if __name__ == "__main__":
    rho_opt = interpolate(V, C)    # Initial value interpolated in V
    u = forward(rho_opt)           # Forward problem

    J = Functional(dot(b, u)*dx)   # Functional Min Compliance
    m = Control(rho_opt)           # Control
    Jhat = ReducedFunctional(J, m) # Reduced Functional

    lb = 0.0  # Inferior
    ub = 1.0  # Superior

    class VolumeConstraint(InequalityConstraint):
        """A class that enforces the volume constraint g(a) = V - a*dx >= 0."""
        def __init__(self, V):
            self.V  = float(V)
            self.smass  = assemble(TestFunction(C) * Constant(1) * dx)
            self.tmpvec = Function(C)

        def function(self, m):
            self.tmpvec.vector()[:] = m
            integral = self.smass.inner(self.tmpvec.vector())
            if MPI.rank(mpi_comm_world()) == 0:
                print "Current control integral: ", integral
            return [self.V - integral]

        def jacobian(self, m):
            return [-self.smass]

        def output_workspace(self):
            return [0.0]

        def length(self):
            """Return the number of components in the constraint vector (here, one)."""
            return 1

    problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V))
    parameters = {"acceptable_tol": 1.0e-200, "maximum_iterations": 100}

    solver = IPOPTSolver(problem, parameters=parameters)
    u_opt = solver.solve()

    File("1_dist_load/control_solution_1.xdmf") << u_opt
asked Jan 6, 2016 by C. Okubo FEniCS Novice (470 points)

1 Answer

+2 votes
 
Best answer

You need some kind of regularisation. Either length scale control by means of a filter or a penalty term related to the boundary. The latter is the easiest option, as you can just do

p2 = Constant(0.01) #experiment with this
J += p2*dot(grad(rho_opt),grad(rho_opt))
answered Jan 6, 2016 by KristianE FEniCS Expert (12,900 points)
selected Jan 10, 2016 by C. Okubo

Hi KristianE! Thanks for your help!

My code was lacking what you comment, but even with that, it is not working the way it should yet.

I tried to solve the example in the following link with my code, but the results is still weird.

https://bitbucket.org/dolfin-adjoint/da-applications/src/217862f7d1f8bf16f2ccd9ecd751cb793d7210fc/structural_optimisation/?at=default

I am working on the code, but if you see some more errors or have any more comments, I'll be glad to see!

Here is an example without IPopt that should work

### PLANE STRESS BY KRISTIAN EJLEBJERG JENSEN
from dolfin import *
from dolfin_adjoint import *
from numpy import ones, array
from numpy import dot as pydot
from numpy import abs as pyabs
from numpy import sqrt as pysqrt

set_log_level(INFO+1)
def XXline(Lmin_in=5e-2, Tvol=0.5,Emin_in=1e-3, nu_in=0.3, mvlimit=0.1, Lx=2., Ly=1.0, Lz=1.0, \
meshsz=1, rhomin=1e-2, simpP_in=1.0, maxN = None, do3D=False, itertotal=100):
    def boundary(x):
        return x[0] < DOLFIN_EPS
    if do3D:
        def boundary2a(x):
            return x[1]-0.5*Ly < DOLFIN_EPS
        def boundary2b(x):
            return x[2]-0.5*Lz < DOLFIN_EPS
        class Loadbnd(SubDomain):
            def inside(self, x, on_boundary):
               return Lx-DOLFIN_EPS < x[0]  and 0.5*Ly-0.05-DOLFIN_EPS < x[1] and x[1] < 0.5*Ly+0.05+DOLFIN_EPS and \
                                                0.5*Lz-0.05-DOLFIN_EPS < x[2] and x[2] < 0.5*Lz+0.05+DOLFIN_EPS
    else:
        def boundary2(x):
            return x[1]-0.5 < DOLFIN_EPS
        class Loadbnd(SubDomain):
            def inside(self, x, on_boundary):
               return Lx-DOLFIN_EPS < x[0]  and 0.5*Ly-0.05-DOLFIN_EPS < x[1] and x[1] < 0.5*Ly+0.05+DOLFIN_EPS
    if do3D:
        if dolfin.__version__ > '1.5.0':
            mesh = BoxMesh(Point(0.,Ly/2.,Lz/2.),Point(Lx,Ly,Lz),40*meshsz,10*meshsz,10*meshsz)
        else:
            mesh = BoxMesh(0.,Ly/2.,Lz/2.,Lx,Ly,Lz,40*meshsz,10*meshsz,10*meshsz)
    else:
        if dolfin.__version__ > '1.5.0':
            mesh = RectangleMesh(Point(0.,Ly/2.),Point(Lx,Ly),40*meshsz,10*meshsz)
        else:
            mesh = RectangleMesh(0.,Ly/2.,Lx,Ly,40*meshsz,10*meshsz)
    Vd = FunctionSpace(mesh,'DG',0) #CG, 1
    rho = interpolate(Constant(Tvol),Vd)
    Vrho = (assemble(interpolate(Expression("1"),Vd)*TestFunction(Vd)*dx)).array()
    ### CONSTANTS
    simpP = Constant(simpP_in)
    Emin = Constant(Emin_in)
    Emax = Constant(1.0)
    nu = Constant(nu_in)
    Lmin = Constant(Lmin_in)
    if do3D:
        u0 = Constant((0.0, 0.0, 0.0))
        g = Constant((0.0, -1.0, 0.0))
    else:
        u0 = Constant((0.0, 0.0))
        g = Constant((0.0,-1.0))
    def forward_model(designANDp, annotate=True):
     rho = designANDp[0]
     simpP = designANDp[1]
     mesh = rho.function_space().mesh()

     ### INITIALIZE VARIABLES
     V = VectorFunctionSpace(mesh, "CG" ,1)
     u_trl = TrialFunction(V)
     u_tst = TestFunction(V)
     u = Function(V)

     boundaries = FacetFunction("size_t",mesh)
     loadbnd = Loadbnd()
     boundaries.set_all(0)
     loadbnd.mark(boundaries, 1)
     bc1 = DirichletBC(V, u0, boundary)
     ds = Measure("ds")[boundaries]
     if do3D:
         [Vx,Vy,Vz] = V.split()
         bc2a = DirichletBC(Vx, Constant(0.), boundary2a)
         bc2b = DirichletBC(Vz, Constant(0.), boundary2b)
         bc2c = DirichletBC(Vz, Constant(0.), boundary2a)
         bcs = [bc1,bc2a,bc2b,bc2c]   
     else:
         [Vx,Vy] = V.split()
         bc2 = DirichletBC(Vx, Constant(0.), boundary2)
         bcs = [bc1,bc2]

     E = Emin+(Emax-Emin)*pow(rho,simpP)
     ## FIND RESPONSE
     if do3D:
       TrGradU = tr(grad(u_trl)) 
     else:
       TrGradU = tr(grad(u_trl))-nu*div(u_trl)/(1.0-nu)
     sigma = E*(1.0/(1.0+nu)*sym(grad(u_trl)) + nu/(1.0+nu)/(1.0-2.0*nu)*TrGradU*Identity(u_trl.cell().d))
     a = inner(sigma, sym(grad(u_tst)))*dx
     L = inner(g,u_tst)*ds(1)
     solve(a == L, u, bcs, annotate=annotate)
     return [u,ds]

    ### ITERATE TO FIND DESIGN NEW DESIGN   
    fid  = File("rho.pvd")
    pmax = 3.0+do3D; deltap = 2.0/itertotal*1.5;
    for iii in range(itertotal):
     ## RUN FORWARD MODEL
#     adj_reset()
     if iii == 0:
      [u,ds] = forward_model([rho,simpP])
      Phi  = Functional(dot(u,g)*ds(1))
      #control = [InitialConditionParameter(rho),ScalarParameter(simpP)]
      control = [Control(rho),ConstantControl(simpP)]
      rf  = ReducedFunctional(Phi,  control)
     df0dx = rf.derivative(forget=False)[0]
     f0val = rf([rho,simpP])
     xval = rho.vector().array()
     fval = pydot(xval,Vrho)/Vrho.sum()

     ## FILTER SENSITIVITY
     Vr = FunctionSpace(mesh, "CG", 1)
     f_trl = TrialFunction(Vr)
     f_tst = TestFunction(Vr)
     df0dx_ = Function(Vr)
     df0dx.vector().set_local(df0dx.vector().array()/Vrho*xval)
     a2 = (Lmin**2*dot(grad(f_trl),grad(f_tst))+f_trl*f_tst)*dx 
     L2 = df0dx*f_tst*dx
     solve(a2 == L2, df0dx_, [])
     df0dx = interpolate(df0dx_,FunctionSpace(mesh,'DG',0)).vector().array()*Vrho/xval

     ## UPDATE DESIGN
     if iii == 0:
      f0old = f0val
     xval = OptC(xval, df0dx, array([Vrho/Vrho.sum()]), Tvol, mvlimit, rhomin)     
     f0chg_pc  = 100*(f0val-f0old)/f0old
     f0old = f0val
     rho.vector().set_local(xval)

     # PRINT AND SAVE RESULTS 
     outstr = "%3.0f, C: %0.4f, Vol: %0.4f, p: %0.4f, Conv: %0.4f" % (iii, f0val, fval, float(simpP), f0chg_pc)
     log(INFO+1,outstr) #; fid.write(outstr+"\n")
     fid  << rho

     if float(simpP) <= pmax:
       simpP.assign(float(simpP)+deltap)
     #plot(rho,interactive=True)

def OptC(xval, df0dx, dfdx, Tvol, mvlimit, rhomin):
    i1=0; i2=1e6
    dfdx = dfdx[0,:]   
    N = len(df0dx)
    while i2 - i1 > 1e-4:
        lmid = (i1+i2)/2.
        xnew = array([array([array([array([xval+mvlimit,xval*pysqrt(pyabs(-df0dx)/lmid/dfdx)]).min(0),\
                     ones(N)]).min(0),xval-mvlimit]).max(0),ones(N)*rhomin]).max(0)
        if pydot(xnew,dfdx)/sum(dfdx) - Tvol > 0:
            i1 = lmid
        else:
            i2 = lmid
    return xnew

if __name__=="__main__":
 XXline(Lmin_in=1e-2, meshsz=3)

Hi Kristian!

I believe I found my mistake. Beyond what you said before about regularisation, the following line:

lmbdatr(sym(grad(u)))Identity(2) + 2Gsym(grad(u))

Was:

lmbdatr(grad(u))Identity(2) + 2Gsym(grad(u))

It was lacking "sym".

Thanks for your help!

Dear Okubo, could you please share your code? I've tried to fix sigma exactly like you did and also applying regularization, but no convergence still.

Hi gtt,

Sorry for taking so long to answer. It's been a while that I don't work on this. My research turned to another path. I put the complete code below, but I am not sure if it is working. It may need some adjusts to current versions. Hope it helps you.

from mshr import *
from dolfin import *
from dolfin_adjoint import *
import pyipopt

## Geometry and elasticity
t, h, L = 2., 1., 5.                          # Thickness, height and length
Rectangle = Rectangle(Point(0,0),Point(L,h))  # Geometry
E, nu = 210*10**3, 0.3                        # Young Modulus
G = E/(2.0*(1.0 + nu))                        # Shear Modulus
lmbda = E*nu/((1.0 + nu)*(1.0 -2.0*nu))       # Lambda

## SIMP
def simp(x):
    return eps+(1-eps)*x**p

V = Constant(0.4*L*h)  # Volume constraint
p = 4                  # Exponent
eps = Constant(1.0e-6) # Epsilon for SIMP

## Mesh, Control and Solution Spaces
nelx = 250
nely = 50
mesh = RectangleMesh(Point(0.0, 0.0), Point(L, h), nelx, nely)    
#mesh = generate_mesh(Rectangle, 100) # Mesh

A = VectorFunctionSpace(mesh, "Lagrange", 1) # Displacements
C = FunctionSpace(mesh, "Lagrange", 1)       # Control

## Volumetric Load
q = -10.0/t
b = Constant((0.0, q))

## Dirichlet BC em x[0] = 0
def Left_boundary(x, on_boundary):
    return on_boundary and abs(x[0]) < DOLFIN_EPS
u_L = Constant((0.0, 0.0))
bc = DirichletBC(A, u_L, Left_boundary)

## Forward problem solution
def forward(x):
    u = TrialFunction(A)  ## Trial and test functions
    w = TestFunction(A)
    sigma = lmbda*tr(sym(grad(u)))*Identity(2) + 2*G*sym(grad(u)) ## Stress
    F = simp(x)*inner(sigma, grad(w))*dx - dot(b, w)*dx
    a, L = lhs(F), rhs(F)
    u = Function(A)
    solve(a == L, u, bc)
    return u

## MAIN
if __name__ == "__main__":
    Vin = Constant(V/(L*h))
    x_in = interpolate(Vin, C)  # Initial value interpolated in V
    u = forward(x_in)           # Forward problem

    J = Functional(dot(b, u)*dx + Constant(1.0e-8)*dot(grad(x_in),grad(x_in))*dx)
    m = Control(x_in)               # Control
    Jhat = ReducedFunctional(J, m)  # Reduced Functional

    lb = 0.0  # Inferior
    ub = 1.0  # Superior

    class VolumeConstraint(InequalityConstraint):
        """A class that enforces the volume constraint g(a) = V - a*dx >= 0."""
        def __init__(self, V):
            self.V  = float(V)
            self.smass  = assemble(TestFunction(C) * Constant(1) * dx)
            self.tmpvec = Function(C)

        def function(self, m):
            self.tmpvec.vector()[:] = m
            integral = self.smass.inner(self.tmpvec.vector())
            if MPI.rank(mpi_comm_world()) == 0:
                print "Current control integral: ", integral
            return [self.V - integral]

        def jacobian(self, m):
            return [-self.smass]

        def output_workspace(self):
            return [0.0]

        def length(self):
            """Return the number of components in the constraint vector (here, one)."""
            return 1

    problem = MinimizationProblem(Jhat, bounds=(lb, ub), constraints=VolumeConstraint(V))
    parameters = {"acceptable_tol": 1.0e-16, "maximum_iterations": 100}

    solver = IPOPTSolver(problem, parameters=parameters)
    rho_opt = solver.solve()

    File("1_dist_load/control_solution_1.xdmf") << rho_opt
...