New link to C++SCPIP and more stuff

I consider to publish more C++ code like C++SCPIP. I hope for
a C++ILUPACK version, maybe a C++MMA and then there is a FEM code in C++ I wrote 2 years ago together with Michael Tesch, which could serve educational purpose (C++FEM).

Therefore I renamed the Sourforge project to C++Math. As a matter of fact most academic state of the art code is written in Fortran, runs
excellent but often lacks ease of use. On the other side is this my contribution to bring more object oriented techniques into the academic community.

The next post will have some pictures again! :)

Battle of the tools

As I worked for 5 years as a software developer after my
diploma degree before I did my master in computational engineering, I got used to work with tools. I shall add that I do a lot of object oriented stuff - always in the seek of a balance between framework generalization and pragmatism.

What I miss is Together, one of the best (an for non-academic used most expensive) UML tools for visual supported designing,
refactoring or simply navigation within the code. I was not able yet to make it parse the huge C++ finite element simulation code from my institute.

What I really find a big boost ist Eclipse with the C/C++ plugin. It is state of the arte and has a lot of features lightyears away from emacs and vi. For users I suggest the blog of Doug, the CDT project lead. Contact me first if you want to try it with Linux on a 64bit system (Code2, Opteron, ...)!

It integraties smooth with our version control sytem subversion. And if you ever worked with one you will never want to miss it. It is open source, fast to set-up and supported my many tools on all platforms.

For finding memory leaks we use valgrind, with the profiling tool chachegrind (which simulates the cache!). To be used with valkyrie (very important!) and kcachegrind.

For bugtracking, wiki and web-browsing our repository we use trac, a feature rich web-tool.

So it is at least my personal experience, that object oriented software and tools is good for this kind of academic, high-performance oriented software done at the universities. I simply don't belive, that FORTRAN is significantly faster. Hey, FORTRAN77 is 3 years younger than the Intel 4004! But it is on the other side true, that the learning curve for OO takes years, but this is also true for math actually for every profession.

C++SCPIP published

I finally managed to publish my open source C++ interface to the SQP/MMA optimizer SCPIP from Chrisitan Zillober. The current entry point is cppmath.wiki.sourceforge.net.

I want to thank Christian Zillober for providing me the code of SCPIP for my thesis and Andreas Wächter from whom I took a lot of ideas and code from his IPOPT.

I hope C++SCPIP is for use of some optimizers out there. If so, please let me know :)

DCAMM Advanced school

I participated at the DCAMM Advanced school from June 20 - 26, 2007. This was offered by the DTU in Denmark.

This international course for Ph.D. students provided a good introduction to the SIMP
method with exercises based on the 99 lines code. DTU is full of SIMP experts (to start
with Martin Bendsoe and Ole Sigmund - who are really nice!!)

It was really inspiring and it helped me a lot to talk to the people there. Especially
Jakob Søndergaard Jensen :)

I had a good time there, met nice people and I enjoyed it to talk to so many people working
on similar problems as I do. There is an picture of the final presentation of the group results.



So, if you have an option then I advice you to attend this course the next time. One learns a lot!

Maximizing the compliance

I was wondering what would happen if I maximize the compliance - in other words to minimize stiffness. Note that the compliance is defined as the scalar product of force (externaly applied in the center of this bridge) and the displacement. <f, u>. When minimizing compliance I mininimize the mechanical energy which is <u, Ku> what is equivalent. Here the first formulation is more interesting, as our force is applied in the example only at a single node. The simulation is done with IPOPT as optimizer, the animation shows, that in a first attempt IPOPT weakens the material at the fix points - but removing the material at the load point allows clearly the highest displacement at this single load point.

Comparing IPOPT with the Optimality Condition

Here I compare the optimality condition against IPOPT via the bridge well known example. Now I fix (horizontal and vertical) the lower right end but let the left end move horizontal (only the vertical displacement is fxed). I think this is a nice example about the impact of boundary conditions (compare with the 2D bridge example!).

The video using IPOPT is here IPOPT video, the video using optimality condition is optimality condition video. The final image for the the IPOPT and optimality condition optimization is:



Note, that the volume constraint (0.5) is always fulfilled! The "flickering" in the IPOPT animation can be better understood having a look at the iteration plot . IPOPT quickly runs into a local minima and it takes some time to leave it. Then one can see brave tries to escape the next minima (which could be close to the global minimum) several times. Note that I do not let IPOPT quit the optimization process but rather limit the maximal number of iterations.

What I find noticable is the with the optimality condition, I really have almost no grayness (I will have to quantify this!) what actually makes sense (and is what we want). In the literature and with [1] find more grayness at the silhouette (as with IPOPT). I'm currently reimplementing the filtering of the gradient and now I also get more "aliasing". What do you get?

[Update 10.2010]
I found that this pretty old post is sometimes still read. The key point is, that sensitivity filtering has been used. It is well known that this disturbs the original gradient and general purpose optimizers have problems with that - SNOPT normally fails. IPOPT seems to be more robust here while MMA and OC are rock solid against it. So to compare IPOPT you need rigorous regularization like density filtering, slope constraints or feature size control (e.g. MOLE) to mention a relevant regularization techniques mentioned in THE BOOK.

Used Toolchain

Within this posting I will add more and more of the tools I use.

• As external optimizer I use IPOPT, to be replaced soon by a MMA implementation.
  A. Wächter and L. T. Biegler, On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming, Mathematical Programming 106(1), pp. 25-57, 2006

3D mechanic compliance minimization

The simulations in 3d look much better, but still need some hours to calculate. I hope to speed this up by a factor of 5-10, I will need it for real piezo optimization.

This is a bridge, fixed at the four lower corner points and three load points in the lower plane. I allowed half of the total volume. It is quite symmetric. The optimization is done with the optimiality criteria and the iterations for all this examples are between 15 and 40.


Here the bridge has only a single lower load point.


Here the single load point example but a volume fraction of 10%.


Here Tom (thanks again :)) suggested to optimize the eiffel tower. The top center force in positive z-direction, the lower four corners are fixed.


This is the most challenging example. For this bridge I apply a downwards force on every node on the lower plane. Currently it is still running, so this is a snapshot of iteration 15. The volume fraction is 20% and the penalization is 5 instead of 3. Still I don't think that there will ever be a nice black and white result. The green is for a pseudo density of 0.5. I'm curious if you have any idea to solve this optimization problem. I thought about adaptive penalization (increasingly higher depending on the greyness) - does anyone have experience in this field?


Finally the aqueduct which looks much better in 3d :) Here I use a penalization of 4. With the standard 3 I get nice big holes but I cannot get rid of the greyness. At leas here higher penalization works. The load applied on the upper long central line (a "pipe"). I fix the upper four corners and the front end back point of the feet. Note that there is a bead!

2D mechanic optimizations

From time to time, I calculate 2d and 3d mechanical minimizing compliance problems. This is some kind of leasure while I work on the piezo optimization - so don't take it serious, it's just for fun. On the other side, some optimizations don't really fit to what I expected from publications. So if you also do also thus optimizations, please let me know if you have the same or different results. I can give the parametrization details on request. I use quite fine discretizations and a filter regularization as given in [2] to make the images nicer .:)

The cantilever is in my opinion too thin on the right end which is loaded.


This bridge gets its load from the top. I miss a lower horizontal bar.


Here the bridge has a lower load in the center.


This is a aqueduct, it looks much cooler in 3d. Thanks to Tom for the suggestion for this optimization problem.


If you have any ideas what could interesting to optimize, just let me know.

Using IPOPT as Optimizer

I integrated IPOPT as an external optimizer, again still for the mechanical compliance minimization. I do not use the Hessian of the Lagrangian but let IPTOP to an L-BFGS approximation. The software is really nice - good documentation, open source with a good license (CPL) and a good C++ frontend with cool features.

The results are also good for 2D, the number of iterations is roughly comparable to the optimality criteria and there are about 50% more function evaluations for line search (determine the step length).

Bad news are, that in 3D the performance is much worse - but I still did not do examine what is really going on. Also I had to increase the mechanical load by a few magnitudes - otherwise I had nearly instant convergence with only slight modification of my inital guess for the desing parameters. I'll have to implement the break conditions by myself - what is supported by IPOPT

Scope of current work

Before doing structural optimization of a piezo based loudspeaker, I implement the standard approach in simple compliance minimization in mechanics.

More precisely I do SIMP, following the publications from Sigmund and Bensoe. This includes the simple heuristic optimality criteria update rule. I found the form given in the 99-lines code publication to suit better than the form in the book (see my publications reference list)