Exemples de jeux de données Cast3M

************************************************************************
* Example of topological optimization                                  *
* Density method + penalization (SIMP)                                 *
*                                                                      *
* This test case is based on the following article and reproduces      *
* the case in fig. 5 :                                                 *
* Wu, Aage, Westermann, Sigmund (2018)                                 *
* "Infill Optimization for Additive Manufacturing - Approaching        *
*  Bone-like Porous Structures"                                        *
* IEEE Transactions on Visualization and Computer Graphics,            *
* vol. 24, no. 2, pp. 1127-1140                                        *
* https://doi.org/10.1109/TVCG.2017.2655523                            *
*                                                                      *
* Application to a beam under bending in 2D plane stresses             *
* with an additional infill constraint                                 *
*                                                                      *
* Optimization algorithm: Method of Moving Asymptotes                  *
*                                                                      *
* Problem: minimize the compliance with:                               *
*          a constraint on the global volume fraction                  *
*          a constraint on the local  volume fraction                  *
*                                                                      *
*          min   C(x)  = u^T.F                                         *
*         s.t.   G1(x) = vf(x)/volfrac - 1    < 0                      *
*                G2(x) = p-norm(vf)/alpha - 1 < 0                      *
*                                                                      *
*          /+----------------------------------------------+           *
*          /|                                              |           *
*          /|                                              |           *
*          /|                                              |           *
*          /|                                              |           *
*          /|                                           p2 + |         *
*          /|                                              | |         *
*          /|                                              | |         *
*          /|                                              | v Force   *
*          /|                                              |           *
*          /+----------------------------------------------+           *
*                                                                      *
************************************************************************
 
 
** Filtering procedure
DEBP HFILT cham1*'MCHAML' mo*'MMODEL' mf*'RIGIDITE' mpt*'MAILLAGE' ;
  chp1     = MANU 'CHPO' mpt 'SCAL' (EXTR cham1 'VALE' 'SCAL') ;
  chp2     = mf * chp1 ;
  cham2    = MANU 'CHML' mo 'REPA' 'SCAL' (EXTR chp2 'VALE' mpt) 'TYPE' 'SCALAIRE' 'GRAVITE' ;
FINP cham2 ;
 
** Global parameters
itrac    = FAUX ;
pal1     = PALE 'GRIS' ;
OPTI 'DIME' 2 'MODE' 'PLAN' 'CONT' 'ELEM' 'QUA4' 'PALE' pal1 'ECHO' 0 ;
 
** Geometrical parameters (width and height) / choose 400 and 200 to match the case in the article
l        = 120. ;
h        = 60. ;
e        = 1. ;
 
** Mesh parameters / choose 400 and 200 to match the case in the article
nelx     = 120 ;
nely     = 60 ;
 
** Material parameters (isotropic elasticity)
e0       = 1. ;
emin     = e0 / 1.E9 ;
nu       = 0.3 ;
 
** Topology optimization parameters
*  penal     : SIMP penalization coefficient
*  volfrac   : minimal volume fraction
*  ft        : = 1 apply filter on the compliance sensitivity field
*              = 2 apply spatial filter on the density field
*              = 3 apply spatial filter on the density field + thresholding by heaviside function
*  rmin      : filter radius (in m)
*  beta      : initial value of the slope of the heaviside function (only for ft = 3)
*              will be doubled every 50 iterations
*  move      : limit of the maximal increment of density
*  changmax  : optimization stop criterion
*  xmin/xmax : min and max density bounds
*  rinfill   : radius of the local volume fraction constraint
*  ainfill   : upper bound for the local volume fraction
*  pinfill   : aggregation coefficient of the local volume fraction constrainsts
penal    = 3. ;
volfrac  = 0.5 ;
ft       = 2 ;
rmin     = 1.5 * l / nelx ;
beta     = 1 ;
move     = 0.2 ;
changmax = 0.01 ;
xmin     = 0. ;
xmax     = 1. ;
rinfill  = 6. * l / nelx ;
ainfill  = 0.6 ;
pinfill  = 16. ;
 
** Mesh
p0       = 0. 0. ;
p1       = 0. h ;
ll       = DROI nely p0 p1 ;
mesh     = ll TRAN nelx (l 0.) ;
con      = CONT mesh ;
p2       = con POIN 'PROC' (l (h / 2.)) ;
 
** Mechanical model
mod      = MODE mesh 'MECANIQUE' ;
 
** Material properties field with a unit Young modulus
maun     = MATE mod 'YOUN' 1. 'NU' nu ;
 
** Boundary conditions
blo      = BLOQ 'UX' 'UY' ll ;
 
** Load (local force)
f        = FORC (0. -1.) p2 ;
 
** Volume of each element (ve) and total volume (vtot)
un       = MANU 'CHML' mod 'SCAL' 1. 'GRAVITE' ;
ve       = INTG 'ELEM' mod un maun ;
vtot     = INTG mod un maun ;
 
** Initial density field
xini     = volfrac ;
x        = MANU 'CHML' mod 'SCAL' xini 'GRAVITE' ;
change   = 1. ;
 
** Physical density field
SI ((ft EGA 1) OU (ft EGA 2)) ;
  xphys    = x ;
FINSI ;
SI (ft EGA 3) ;
  xtilde   = x ;
  xphys    = 1. - (EXP (-1. * beta * xtilde)) + (xtilde * (EXP (-1. * beta))) ;
FINSI ;
 
** Initial volume fraction
vfx      = (INTG mod xphys maun) / vtot ;
 
** Gravity centers and filtering matrix
*  the weight depends on the volume of each element
ptg      = un POIN 'SUPERIEUR' 0. ;
wg       = MANU 'CHPO' ptg 'SCAL' (EXTR ve 'VALE' 'SCAL') ;
kfil     = MFIL wg rmin 1. 0. ;
ung      = MANU 'CHPO' ptg 1 'SCAL' 1. ;
ks       = kfil * ung ;
kfil     = NFIL kfil ks ;
 
** Filtering matrix for the infill constraint
kin      = MFIL wg rinfill 0. 0. ;
ks       = kin * ung ;
kin      = NFIL kin ks ;
 
** Initialization of the table for the MMA
nx       = NBEL mesh ;
tmma     = TABL ;
* initial values of the design variables x
lx       = EXTR x 'VALE' 'SCAL' ;
tmma . 'X' = lx ;
* bounds for x
tmma . 'XMIN' = xmin ;
tmma . 'XMAX' = xmax ;
* other parameters
tmma . 'MOVE' = move ;
 
** Lets's start a timer
TEMP 'ZERO' ;
 
** Topology optimization loop
liso     = PROG 0. 'PAS' 0.05 1. ;
loopbeta = 0 ;
liter    = PROG ;
lobj     = PROG ;
lvf      = PROG ;
lchange  = PROG ;
REPE b1 150 ;
  loopbeta = loopbeta + 1 ;
* penalization of the stiffness matrix (modified SIMP)
  ep       = emin + ((xphys ** penal) * (e0 - emin)) ;
  map      = MATE mod 'YOUN' ep 'NU' nu 'DIM3' e ;
  k        = RIGI mod map ;
* resolution of the FE problem
  kbc      = k ET blo ;
  u        = RESO kbc f ;
* objective function: compliance = u^T.F
* value
  c        = MAXI (RESU (PSCA u f (MOTS 'UX' 'UY') (MOTS 'FX' 'FY'))) ;
* sensitivity (gradient with respect to the physical variable)
  eps      = EPSI 'LINE' mod u ;
  sigun    = ELAS mod maun eps ;
  eneun    = ENER mod eps sigun ;
  eneun    = INTG 'ELEM' mod eneun ;
  dc       = -1. * penal * (xphys ** (penal - 1.)) * (e0 - emin) * eneun ;
* constraint function: g1 = vf(x)/volfrac - 1
* value
  g1       = (vfx / volfrac) - 1. ;
* sensitivity (gradient with respect to the physical variable)
  dg1      = ve / vtot / volfrac ;
* constraint function: g2 = p-norm(vf(x))/alpha - 1
* value
  xbar     = HFILT xphys mod kin ptg ;
  lxbar    = EXTR xbar 'VALE' 'SCAL' ;
  gpmoy    = AGRE 'PMOY' lxbar pinfill ;
  g2       = (gpmoy / ainfill) - 1. ;
* sensitivity (gradient with respect to the physical variable)
  tmp      = gpmoy / ((AGRE 'PMOM' lxbar pinfill) / nx) ;
  dg2      = (1. / (nx * ainfill)) * tmp * (xbar ** (pinfill - 1.)) ;
* ft = 1  --> filtering the complicance sensitivity
  SI (ft EGA 1) ;
    dc       = (HFILT (x * dc) mod kfil ptg) / (BORN x 'MINIMUM' 1.E-3) ;
  FINSI ;
* ft = 2,3 --> updating the sensitivies (to become gradients with respect to the design variable)
  SI (ft EGA 2) ;
    dc       = HFILT dc  mod kfil ptg ;
    dg1      = HFILT dg1 mod kfil ptg ;
    dg2      = HFILT dg2 mod kfil ptg ;
  FINSI ;
  SI (ft EGA 3) ;
    dx       = (beta * (EXP (-1. * beta * xtilde))) + (EXP (-1. * beta)) ;
    dctilde  = dc * dx ;
    dc       = HFILT dctilde mod kfil ptg ;
    dg1tilde = dg1 * dx ;
    dg1      = HFILT dg1tilde mod kfil ptg ;
    dg2tilde = dg2 * dx ;
    dg2      = HFILT dg2tilde mod kfil ptg ;
  FINSI ;
* information about the current topology
  info     = CHAI 'It:' (&b1 - 1) / 5 'Obj:' / 10 c > 1 'Vol. frac:' > 4 vfx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
  MESS info ;
  SI itrac ;
    TRAC xphys mod con liso 'TITR' info 'NCLK' ;
  FINSI ;
* update the topology by optimization
  tmma . 'F0VAL' = c ;
  tmma . 'DF0DX' = EXTR dc 'VALE' 'SCAL' ;
  tmma . 'FVAL'  = PROG g1 g2 ;
  tmma . 'DFDX'  = ENUM (EXTR dg1 'VALE' 'SCAL') (EXTR dg2 'VALE' 'SCAL') ;
  MMA tmma ;
  lxnew    = tmma . 'X' ;
  xnew     = MANU 'CHML' mod 'REPA' 'SCAL' lxnew 'TYPE' 'SCALAIRE' 'GRAVITE' ;
* update the physical density (by filtering and thresholding)
  SI (ft EGA 1) ;
    xphys    = xnew ;
  FINSI ;
  SI (ft EGA 2) ;
    xphys    = HFILT xnew mod kfil ptg ;
  FINSI ;
  SI (ft EGA 3) ;
    xtilde   = HFILT xnew mod kfil ptg ;
    xphys    = 1. - (EXP (-1. * beta * xtilde)) + (xtilde * (EXP (-1. * beta))) ;
  FINSI ;
* updating the volume fraction
  vfx      = (INTG mod xphys maun) / vtot ;
* summary of the current iteration
  change   = MINI (MAXI 'ABS' (lxnew - (tmma . 'XOLD1')))
                  (MAXI 'ABS' (lxnew - (tmma . 'XOLD2'))) ;
  liter    = liter ET &b1 ;
  lobj     = lobj ET c ;
  lvf      = lvf ET vfx ;
  lchange  = lchange ET change ;
* preparing the next iteration
  x        = xnew ;
  lx       = lxnew ;
* ft = 3 --> updatind the slope of the thresholding function
  SI (ft EGA 3) ;
    SI ((beta < 128) ET ((loopbeta >EG 50) OU (change &lt;EG changmax))) ;
      beta     = 2 * beta ;
      loopbeta = 0 ;
      change   = 1. ;
    FINSI ;
  FINSI ;
* stop criterion
  SI (change < changmax) ;
    info     = CHAI 'It:' &b1 / 5 'Obj:' / 10  c > 1 'Vol. frac:' > 4 vfx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
    MESS info ;
    QUIT b1 ;
  FINSI ;
FIN  b1 ;
 
** Elapsed time
TEMP 'IMPR' 'SOMM' 'HORL' ;
 
** Plotting the final topology
evobj    = EVOL 'ROUG' 'MANU' 'Iterations' liter 'Compliance' lobj ;
evvf     = EVOL 'ORAN' 'MANU' 'Iterations' liter 'Frac. vol.' lvf ;
evchange = EVOL 'VERT' 'MANU' 'Iterations' liter 'Max. change' lchange ;
SI itrac ;
  info     = CHAI '[Final topology]' ' ' info ;
  TRAC xphys mod con liso 'TITR' info ;
  DESS evobj    'TITR' 'Objective function' ;
  DESS evvf     'TITR' 'Volume fraction' ;
  DESS evchange 'TITR' 'Max. change' ;
FINSI ;
 
 
FIN ;