Example: 2_load_cases_3_dim
===========================

.. contents::

This is a 3 dimensional example to show :py:mod:`ttrr`.

This example is also available as :download:`batch_2_load_cases_3_dim.py`.

First start python and import necessary modules::

  >>> import ttrr,numpy,scipy

design space
------------
To define the design space::

  >>> dim = 3 # dimension of the design space
  >>> nelx = 3 # number of nodal points in x-direction
  >>> nely = nelx # number of nodal points in y-direction
  >>> nelz = nelx # number of nodal points in z-direction

To let you easily change the design space, we use here a
calculation to define the grid. The measurement unit is meter::

  >>> n = max(nelx,nely,nelz)
  >>> gx0 = 0.0
  >>> gx1 = (nelx-1)*1.0/(n-1)
  >>> gy0 = 0.0
  >>> gy1 = (nely-1)*1.0/(n-1)
  >>> gz0 = 0.0
  >>> gz1 = (nelz-1)*1.0/(n-1)

bearings
--------
We want to fix all nodes on the ground. We can describe these
area (All nodes will become fixed bearings)::

  >>> mountingarea = numpy.array([
  ...     [ gx0 , gx1 , gy0 , gy1 , gz0 , gz0 ]
  ...     ])
  >>> mountingarea[:,[0,2,4]] = mountingarea[:,[0,2,4]] - 1e-6
  >>> mountingarea[:,[1,3,5]] = mountingarea[:,[1,3,5]] + 1e-6

get the parameter
-----------------
Now we can get the parameter for the calculation from barsnodesongrid3d.
Here we can also use the optional parameter maxstablaenge to decide that
bars longer than maxstablaenge should be ignored::

  >>> maxstablaenge = 3 # maximum allowed length of a bar
  >>> [nodes,bars,mounting] = ttrr.barsnodesongrid3d(gx0,gx1,gy0,gy1,gz0,gz1,nelx,nely,nelz,mountingarea=mountingarea,maxstablaenge=maxstablaenge)
  mounting wird aus 1 Zeilen von mountingarea erstellt
  9 Knoten sind gelagert; 27 Lager
  302 potenzielle Staebe gefunden; beruecksichtige 27 Lager
  es bleiben 274 Staebe, die nicht nur gelagert sind
  274 potenzielle Staebe; beruecksichtige maximale Stablaenge 3.000000
  es bleiben 274 Staebe, die kuerzer sind als 3.000000

We can look at the Parameter::

  >>> nodes
  array([[ 0. ,  0. ,  0. ],
         [ 0. ,  0.5,  0. ],
         [ 0. ,  1. ,  0. ],
         [ 0.5,  0. ,  0. ],
         [ 0.5,  0.5,  0. ],
         [ 0.5,  1. ,  0. ],
         [ 1. ,  0. ,  0. ],
         [ 1. ,  0.5,  0. ],
         [ 1. ,  1. ,  0. ],
         [ 0. ,  0. ,  0.5],
         [ 0. ,  0.5,  0.5],
         [ 0. ,  1. ,  0.5],
         [ 0.5,  0. ,  0.5],
         [ 0.5,  0.5,  0.5],
         [ 0.5,  1. ,  0.5],
         [ 1. ,  0. ,  0.5],
         [ 1. ,  0.5,  0.5],
         [ 1. ,  1. ,  0.5],
         [ 0. ,  0. ,  1. ],
         [ 0. ,  0.5,  1. ],
         [ 0. ,  1. ,  1. ],
         [ 0.5,  0. ,  1. ],
         [ 0.5,  0.5,  1. ],
         [ 0.5,  1. ,  1. ],
         [ 1. ,  0. ,  1. ],
         [ 1. ,  0.5,  1. ],
         [ 1. ,  1. ,  1. ]])
  >>> bars
  array([[ 0,  9],
         [ 0, 12],
         [ 0, 21],
         [ 0, 15],
         [ 0, 10],
         [ 0, 19],
         [ 0, 11],
         [ 0, 13],
         [ 0, 22],
         [ 0, 14],
         [ 0, 23],
         [ 0, 16],
         [ 0, 25],
         [ 0, 17],
         [ 9, 18],
         [ 9, 10],
         ...
         [17, 22]], dtype=uint16)
  >>> mounting
  array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
         17, 18, 19, 20, 21, 22, 23, 24, 25, 26])
  >>> (kn,kr) = nodes.shape
  >>> (sn,sr) = bars.shape
  >>> (ln,) = mounting.shape

This means: We have kn=27 nodes and sn=274 bars. ln=27 node-directions are fixed. sr=2 is the number of nodes for a bar and kr=3 is the dimension of the design space.

load cases
----------
The loads are stored in a sparse matrix. For every node we have dim directions. So we have dim*kn rows and 2 cols for the 2 load cases::

  >>> loads = scipy.sparse.lil_matrix((dim*kn,2),dtype=numpy.float64) # 2 load cases

The number k1=19 of a node can be calculated from its coordinates. The direction r1=57 is in z-direction. The magnitude -100000.0 descibes that the direction is the negative z-direction::

  >>> k1 = ttrr.nodenumberfromcoordinates(nodes,gx0,(gy0+gy1)/2,gz1) # position of the first force
  >>> r1 = dim*k1 # direction of the first force
  >>> loads[r1,0] = -100000.0 # Newton = ca. 1000 kg
  >>> k1,r1
  (19, 57)

The same for the other load case::

  >>> k2 = ttrr.nodenumberfromcoordinates(nodes,gx1,(gy0+gy1)/2,gz1) # position of the second force
  >>> r2 = dim*k2 # direction of the second force
  >>> loads[r2,1] = 100000.0 # Newton = ca. 1000 kg
  >>> k2,r2
  (25, 75)

We can look at the result::

  >>> print loads
    (57, 0)       -100000.0
    (75, 1)       100000.0

calculation
-----------
We need a few parameters for the calculation::

  >>> smax = 1 # meter^2
  >>> sigma = 1e8 # Pa = 100 MPa = 100 N / mm^2 (yield strength for aluminium)
  >>> sigmamin = - sigma # compressive strength
  >>> sigmamax = sigma # tensile strength
  >>> density = 2.7e3 # kg / m^3 (density for aluminium)
  >>> name="2_load_cases_3_dim"

Now we can do the classical topology optimization::

  >>> [s,w,lengths] = ttrr.ttrr(dim,nodes,bars,mounting,loads,smax,name=name,sigmamin=sigmamin,sigmamax=sigmamax,density=density,opt=0)

To calculate a redundant solution we need (This takes about 47 seconds on a computer with 3.6 GHz)::

  >>> redundancy = 3
  >>> name="2_load_cases_3_dim_redundant"
  >>> [s,w,lengths] = ttrr.ttrr_redundant(dim,nodes,bars,mounting,loads,smax,name=name,sigmamin=sigmamin,sigmamax=sigmamax,density=density,redundancy=redundancy,opt=3,sort=1)

visualization
-------------
To get pictures, we can use :doc:`ttrr_tools` on the command line. To get help::

  $ ttrr_tools.py -h

We get quick and dirty pictures with gnuplot::

  $ ttrr_tools.py --png_gnuplot -f 2_load_cases_3_dim*.tar.bz2

Nice pictures are rendered by povray::

  $ ttrr_tools.py --png_pov -f 2_load_cases_3_dim*.tar.bz2

Because we are here viewing a 3 dimensional structure the 2 dimensional picture
is insufficient. Therefore it is possible to get a movie animation::

  $ ttrr_tools.py --ogv -f 2_load_cases_3_dim*.tar.bz2

If you only want a small animation you can do one of the following::

  $ ttrr_tools.py --apng -xm 320 -ym 240 -f 2_load_cases_3_dim*.tar.bz2
  $ ttrr_tools.py --mng -xm 320 -ym 240 -f 2_load_cases_3_dim*.tar.bz2
  $ ttrr_tools.py --gif -xm 320 -ym 240 -f 2_load_cases_3_dim*.tar.bz2

You can also combine these options::

  $ ttrr_tools.py -xm 320 -ym 240 -f 2_load_cases_3_dim*.tar.bz2 --png_pov --apng --mng --gif

The rendering by povray is time-consuming. Therefore it's done in parallel.

The above animation as apng:

  +---------------------------------------------------------------------+-------------------------------------------------------------------------------+
  | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_apng.png | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_redundant_apng.png |
  |    :width: 7cm                                                      |    :width: 7cm                                                                |
  +---------------------------------------------------------------------+-------------------------------------------------------------------------------+

.. only:: html

  The above animation as gif:

    +----------------------------------------------------------------+--------------------------------------------------------------------------+
    | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim.gif | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_redundant.gif |
    |    :width: 7cm                                                 |    :width: 7cm                                                           |
    +----------------------------------------------------------------+--------------------------------------------------------------------------+

  The above animation as mng:

    +----------------------------------------------------------------+--------------------------------------------------------------------------+
    | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim.mng | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_redundant.mng |
    |    :width: 7cm                                                 |    :width: 7cm                                                           |
    +----------------------------------------------------------------+--------------------------------------------------------------------------+

conclusion
----------
Many people are wondering because the asymmetrical result. At last it's a
symmetric problem. But we have a linear problem. Therefore there possibly
exists more than one edge of the polytope as an optimal solution --- the
optimal value is the same. If we use an interior point algorithm we get a
result between these edges.

The interior point algorithm of glpk has a big problem with scaling.
Therefore we do a small scaling per hand (Note: This is only a scaling
factor. The topology is the same.)::

  >>> name="2_load_cases_3_dim_interior"
  >>> [s,w,lengths] = ttrr.ttrr(dim,nodes,bars,mounting,loads,smax,name=name,sigmamin=sigmamin,sigmamax=sigmamax,density=density,opt=1,scale_sigma=1e-8,scale_forces=1e-8)

Now we can look at the result:

  +---------------------------------------------------------------------+------------------------------------------------------------------------------+
  | simplex                                                             | interior point                                                               |
  +---------------------------------------------------------------------+------------------------------------------------------------------------------+
  | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_apng.png | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_interior_apng.png |
  |    :width: 7cm                                                      |    :width: 7cm                                                               |
  +---------------------------------------------------------------------+------------------------------------------------------------------------------+

I hope you get the picture of the 2 edges of the polytope in this problem.

Another interesting point is the fact that we see a 2 dimensional structure in a 3 dimension space. Because the loads are only acting in one plane this is possible. Sure a small perturbation would destroy this behavior. Therefore it is not a practical solution. If you want a practical solution think of robust::

  >>> name="2_load_cases_3_dim_robust"
  >>> [s,w,lengths] = ttrr.ttrr_robust(dim,nodes,bars,mounting,loads,smax,name=name,sigmamin=sigmamin,sigmamax=sigmamax,density=density,p=0.3,art=0)

And an animation::

  $ ttrr_tools.py -xm 320 -ym 240 -f 2_load_cases_3_dim_robust.tar.bz2 --png_pov --apng --mng --gif

The redundant solution is for a redundancy of 1/3. Every color represent a part
of the structure. Every part is allow to drop out. The other parts are able to
support the given loads. The interesting point is that the red or blue structure
are not connected.

+-------------------------------------------------------------------------------+----------------------------------------------------------------------------+
| redundant                                                                     | robust                                                                     |
+-------------------------------------------------------------------------------+----------------------------------------------------------------------------+
| .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_redundant_apng.png | .. image:: ../ttrr_src/examples_results/2_load_cases_3_dim_robust_apng.png |
|    :width: 7 cm                                                               |    :width: 7 cm                                                            |
+-------------------------------------------------------------------------------+----------------------------------------------------------------------------+

We can also compare the optimal values (This is the volume of the structures.):

+------------------+--------------------+------------------+-----------------+
| optimal solution | redundant solution | optimal solution | robust solution |
+------------------+--------------------+------------------+-----------------+
| simplex          | branch and bound   | interior point   | simplex         |
+------------------+--------------------+------------------+-----------------+
| 0.000875         | 0.001563           | 0.000875         | 0.001154        |
+------------------+--------------------+------------------+-----------------+