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
