The solution of topology optimization problem in 2D formulation

There is a constrained optimization problem based on finite element method. We should find the optimal distribution of density inside the our region for minimazing the energy of deformation. The energy of deformation or compliance can be described by two ways:


Where C - global stiffness matrix, U - displacement vector;

Or by the next way:


Where F - global force vector;

So, we should minimize this value. Our variables are the densities of finite elements.

There is a solution of this problem realized in Mathematica in 3D formulation:


Launching the FEM package


Defining the 3D stress operator:

op = (λ + μ)*
    Grad[Div[{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}], {x, y, 
      z}] + μ*
    Laplacian[{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}]; 
{μ, λ} = {Y/(2*(1 + ν)), ν*(Y/((1 + ν)*(1 - 
           2*ν)))} /. {ν -> 1/3, Y -> 5}; 

Defining the geometrical region

a = 1.; 
reg1 = Parallelepiped[{0, 0, 0}, {{0, 0, 1}, {0, 1, 0}, {1, 0, 0}}]; 
Ω = reg1; 
mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 0.001, 
   "MeshOrder" -> 2]; 

Defining the boundary conditions:

dX = DirichletCondition[u[x, y, z] == 0, z == 0]; 
dY = DirichletCondition[v[x, y, z] == 0, z == 0]; 
dZ = DirichletCondition[w[x, y, z] == 0, z == 0]; 
nZ = NeumannValue[0.01, z == 1 && (x - 0.5)^2 + (y - 0.5)^2 <= 0.5]; 

Discretizing our PDE:

{state} = 
  NDSolve`ProcessEquations[{op == {0, 0, nZ}, dX, dY, dZ}, {u, v, w}, 
   Element[{x, y, z}, mesh], 
   Method -> {"FiniteElement", "PrecomputeGeometryData" -> False}]; 
femdata = state["FiniteElementData"]; 
pded = femdata["PDECoefficientData"]; 
bcd = femdata["BoundaryConditionData"]; 
md = femdata["FEMMethodData"]; 
sd = state["SolutionData"][[1]]; 
vd = state["VariableData"]; 
discretePDE = DiscretizePDE[pded, md, sd]; 
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"]; 
discreteBCs = DiscretizeBoundaryConditions[bcd, md, sd]; 
ele = DiscretizePDE[pded, md, sd, "SaveFiniteElements" -> True]; 

Modification of local stiffness matrices. We attend to each finite element a dependence on its density value by power law:

{stiffnessElements, nonzeroPos} = ele["StiffnessElements"]; 
density = Table[ρ[i], {i, 1, Length[stiffnessElements]}]; 
modifiedElements = density^3*stiffnessElements; 
dof = md["DegreesOfFreedom"]; 
inci = Join @@ md["Incidents"]; 
modifiedstiffness = 
  AssembleMatrix[{inci, inci}, modifiedElements, {dof, dof}, 
   "Sparse" -> nonzeroPos]; 
displacement = Table[δ[i], {i, 1, Length[stiffness]}];

Extracting the data for further computations:

coords = mesh["Coordinates"]; 
elems = mesh["MeshElements"]; 
volume = Join @@ mesh["MeshElementMeasure"];

Deploying the boundary conditions:

DeployBoundaryConditions[{load, stiffness}, discreteBCs]; 
DeployBoundaryConditions[{load, modifiedstiffness}, discreteBCs]; 
loadF = Flatten[load];  

Defining the objective function, constraints and initial point of method. This equations are defined by SIMP method of topology optimization. Initial point was chosen by empirical observations:

goalFunction = displacement . modifiedstiffness . displacement; 
massEquation = 
  Sum[volume[[k]]*density[[k]], {k, 1, Length[volume]}] <= 
   Sum[volume[[k]], {k, 1, Length[volume]}]*0.85; 
densityEq = Thread[(0. <= #1 <= 1 & )[density]]; 
equation = 
  Table[modifiedstiffness[[i]] . displacement == loadF[[i]], {i, 1, 
sysOfEq = Flatten[{goalFunction, massEquation, equation, densityEq}]; 
variables = Flatten[{density, displacement}]; 
initDis = Thread[({#1, 0} & )[displacement]]; 
initDen = Thread[({#1, 1.} & )[density]]; 
initCond = Flatten[{initDen, initDis}, 1]; 

Using of FindMinimum:

AbsoluteTiming[resultFMin = Quiet[FindMinimum[sysOfEq, initCond]]; ]
minimizedDen = density /. resultFMin[[2]]; 

Extracting the values of densities:

Evaluate[Table[ρ[i], {i, 1, Length[stiffnessElements]}]] = 

It is the code of topology optimization. And I have some questions about using FindMinimum.

1) Is it a justified decision to use FindMinimum? I tried to use NMinimize and it didn't give me any results. But I can't control anything inside the interior point which is used for solving constrained optimization problems.

2) I have no differences in processor usage when solving optimization problem with different value of MaxCellMeasure. This value near the 15-20%. But usage of operative memory is very high. How I can set up the upper edge of operative memory usage? And Is there any way to speedup solving of this problem? Maybe there is some low level functions for constrained optimization can be used?

3) Why using the second form of objective function leads to divergence of optimization problem?

If you need a full code of this project, you can look it by this link.

| improve this question | | | | |
  • $\begingroup$ The influence of the density on the PDE is not clear to me. The density should enter the equation also right between Grad and Div in the Grad-Div and the Div-Grad-term (the latter is the Laplacian). Moreover, it is most likely a good idea to write your one (projected?) gradient descent algorithm or Newton algorithm. In any case, it is a good idea to implement the differential of your objective function. $\endgroup$ – Henrik Schumacher Jun 5 '18 at 16:20
  • $\begingroup$ @HenrikSchumacher The influence of density is introduced in stiffness. $E=E_{min}+(E_{max}-E_{min})*\rho$ Thanks for idea of the gradient algorithm. $\endgroup$ – Андрей Кротких Jun 6 '18 at 12:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.