Exemples de jeux de données Cast3M

************************************************************************
* Exemple de methode d'optimisation topologique                        *
* Methode a densite + penalisation SIMP                                *
*                                                                      *
* Ce cas test est une reproduction en Gibiane du code de reference :   *
* Andreassen, Clausen, Schevenels, Lazarov, & Sigmund (2011)           *
* "Efficient topology optimization in MATLAB using 88 lines of code"   *
* Structural and Multidisciplinary Optimization, 43(1), 1-16           *
* https://doi.org/10.1007/s00158-010-0594-7                            *
*                                                                      *
* Application a une poutre en flexion en 2D contraintes planes         *
*                                                                      *
* Algorithme d'optimisation : Optimality Criteria                      *
*                                                                      *
* Probleme : minimisation de la compliance sous contrainte de volume   *
*                                                                      *
*          min   C(x) = u^T.F                                          *
*      tel que   V(x) < V0      (avec V0 = volfrac * Vplein)           *
*                                                                      *
*       |                                                              *
*       | Force                                                        *
*       v                                                              *
*     O>+------------------------------------------------------+       *
*       | p1                                                   |       *
*       |                                                      |       *
*     O>|                                                      |       *
*       |                                                      |       *
*       |                                                      |       *
*     O>+------------------------------------------------------+ p2    *
*                                                              Ʌ       *
*                                                              O       *
*                                                                      *
************************************************************************
 
 
** Procedure de filtrage
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 ;
 
** Parametres globaux
itrac    = FAUX ;
OPTI 'DIME' 2 'MODE' 'PLAN' 'CONT' 'ELEM' 'QUA4' 'ECHO' 0 ;
 
** Parametres geometriques (longueur et hauteur)
l        = 60. ;
h        = 20. ;
e        = 1. ;
 
** Parametres du maillage (nombre d'elements sur x et y)
nelx     = 60 ;
nely     = 20 ;
 
** Parametres materiau (elasticite, isotrope)
e0       = 1. ;
emin     = e0 / 1.E9 ;
nu       = 0.3 ;
 
** Parametres de l'optimisation topologique
*  penal     : coefficient de penalisation SIMP
*  volfrac   : fraction volumique pour la contrainte sur le volume
*  ft        : = 1 filtrage de la sensibilite de la compliance
*              = 2 filtrage de la densite
*              = 3 filtrage de la densite + seuillage par heaviside
*  rmin      : rayon de filtrage (en m)
*  beta      : valeur initiale du coefficient beta pour le seuillage heaviside
*              celui-ci sera double toutes les 50 iterations d'optimisation
*              utile seulement pour ft = 3
*  move      : limite d'increment pour la densite
*  changmax  : critere d'arret de l'optimisation
*  xmin/xmax : densites min et max
penal    = 3. ;
volfrac  = 0.5 ;
ft       = 2 ;
rmin     = 0.04 * l ;
beta     = 1 ;
move     = 0.2 ;
changmax = 0.01 ;
xmin     = 0. ;
xmax     = 1. ;
 
** Maillage
p0       = 0. 0. ;
p1       = 0. h ;
lg       = DROI nely p0 p1 ;
mail     = lg TRAN nelx (l 0.) ;
con      = CONT mail ;
p2       = con POIN 'PROC' (l 0.) ;
nx       = NBEL mail ;
 
** Modele mecanique
mod      = MODE mail 'MECANIQUE' ;
 
** Champ materiau pour un module d'Young unitaire
maun     = MATE mod 'YOUN' 1. 'NU' nu 'DIM3' e ;
 
** Blocages mecaniques
blo      = (BLOQ 'UX' lg) ET (BLOQ 'UY' p2) ;
 
** Chargement de force ponctuelle
f        = FORC (0. -1.) p1 ;
 
** Champ unitaire
un       = MANU 'CHML' mod 'SCAL' 1. 'GRAVITE' ;
 
** Volume total (vtot) et volume cible (v0)
vtot     = INTG mod un ;
v0       = vtot * volfrac ;
 
** Initialisation de la densite
xini     = volfrac ;
x        = MANU 'CHML' mod 'SCAL' xini 'GRAVITE' ;
change   = 1. ;
 
** Calcul de la densite physique
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 ;
 
** Volume et fraction volumique de la topologie
vx       = INTG mod xphys ;
fvx      = vx / vtot ;
 
** Champ de sensibilite de la fonction contrainte (volume)
dvun     = INTG mod un 'ELEM' ;
 
** Maillage des centres de gravite et matrice de filtrage
*  ici, les poids des cellules dependent de leur volume
ptg      = un POIN 'SUPERIEUR' 0. ;
wg       = MANU 'CHPO' ptg 'SCAL' (EXTR dvun 'VALE' 'SCAL') ;
kfil     = MFIL wg rmin 1. 0. ;
 
** On demarre un chrono
TEMP 'ZERO' ;
 
** Boucle d'optimisation topologique
liso     = PROG 0. 'PAS' 0.05 1. ;
loopbeta = 0 ;
liter    = PROG ;
lobj     = PROG ;
lfv      = PROG ;
lchange  = PROG ;
REPE b1 500 ;
  loopbeta = loopbeta + 1 ;
* penalisation de la rigidite (SIMP modifiee)
  ep       = emin + ((xphys ** penal) * (e0 - emin)) ;
  map      = MATE mod 'YOUN' ep 'NU' nu 'DIM3' e ;
  k        = RIGI mod map ;
* resolution du probleme mecanique
  kbc      = k ET blo ;
  u        = RESO kbc f ;
* fonction objectif : compliance = u^T.K.u
  c        = MAXI (RESU (PSCA u f (MOTS 'UX' 'UY') (MOTS 'FX' 'FY'))) ;
* sensibilite de c
  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 ;
* sensibilite de v
  dv       = dvun ;
* filtrage de la sensibilite de la compliance
  SI (ft EGA 1) ;
    dc       = (HFILT (x * dc) mod kfil ptg) / (BORN x 'MINIMUM' 1.E-3) ;
  FINSI ;
* changement des sensibilites due au filtrage de la densite
  SI (ft EGA 2) ;
    dc       = HFILT dc mod kfil ptg ;
    dv       = HFILT dv 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 ;
    dvtilde  = dv * dx ;
    dv       = HFILT dvtilde mod kfil ptg ;
  FINSI ;
* infos sur la topologie courante
  info     = CHAI 'It:' (&b1 - 1) / 5 'Obj:' / 10 c > 1 'Fvol:' > 4 fvx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
  MESS info ;
  SI itrac ;
    def1     = DEFO mail u 0.05 ;
    TRAC xphys mod con def1 liso 'TITR' info 'NCLK' ;
  FINSI ;
* optimisation d'une nouvelle topologie (methode du critere d'optimalite)
  l1       = 0. ;
  l2       = 1.E9 ;
  REPE b2 200 ;
    SI (((l2 - l1) / (l1 + l2)) < 0.001) ;
      QUIT b2 ;
    FINSI ;
    lmid     = 0.5 * (l1 + l2) ;
    b        = -1. * dc / (lmid * dv) ;
    xinf     = BORN (x - move) 'MINIMUM' xmin ;
    xsup     = BORN (x + move) 'MAXIMUM' xmax ;
    xnew     = x * (b ** 0.5) ;
    minf     =  (xnew - xinf) MASQ 'INFERIEUR' 0. ;
    mmil     = ((xnew - xinf) MASQ 'SUPERIEUR' 0.) * ((xnew - xsup) MASQ 'INFERIEUR' 0.) ;
    msup     =  (xnew - xsup) MASQ 'SUPERIEUR' 0. ;
    xnew     = (xinf * minf) + (xnew * mmil) + (xsup * msup) ;
*   mise a jour de la densite physique (selon le filtrage)
    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 ;
*   mise a jour du volume et de la fraction volumique
    vx       = INTG mod xphys ;
    fvx      = vx / vtot ;
    SI (vx > v0) ;
      l1       = lmid ;
    SINON ;
      l2       = lmid ;
    FINSI ;
  FIN b2 ;
* bilan de l'iteration
  change   = MAXI 'ABS' (xnew - x) ;
  liter    = liter ET &b1 ;
  lobj     = lobj ET c ;
  lfv      = lfv ET fvx ;
  lchange  = lchange ET change ;
* preparation de la nouvelle iteration
  x        = xnew ;
* mise a jour du parametre de seuillage beta (si ft = 3)
  SI (ft EGA 3) ;
    SI ((beta < 512) ET ((loopbeta >EG 50) OU (change &lt;EG changmax))) ;
      beta     = 2 * beta ;
      loopbeta = 0 ;
      change   = 1. ;
    FINSI ;
  FINSI ;
* critere d'arret de l'optimisation
  SI (change < changmax) ;
    info     = CHAI 'It:' &b1 / 5 'Obj:' / 10  c > 1 'Fvol:' > 4 fvx > 1 'Change:' > 4 change > 1 'Beta:' > 4 beta > 2 ;
    MESS info ;
    QUIT b1 ;
  FINSI ;
FIN  b1 ;
 
** Mesure du temps ecoule
***TEMP 'IMPR' 'SOMM' 'HORL' ;
 
** Trace de la topologie finale
evobj    = EVOL 'ROUG' 'MANU' 'Iter' liter 'Compliance' lobj ;
evfv     = EVOL 'ORAN' 'MANU' 'Iter' liter 'Frac. vol.' lfv ;
evchange = EVOL 'VERT' 'MANU' 'Iter' liter 'Max. change' lchange ;
SI itrac ;
  info     = CHAI '[Topologie finale]' ' ' info ;
  TRAC xphys mod con liso 'TITR' info ;
  DESS evobj    'TITR' 'Fonction objectif' ;
  DESS evfv     'TITR' 'Fraction volumique' ;
  DESS evchange 'TITR' 'Changement max.' ;
FINSI ;
 
** Maillage de la topologie finale
tab1     = TABL ;
tab1 . 'EPAISSEUR' = 0 ;
tab1 . 'MODELE'    = mod ;
tab1 . 'TOPOLOGIE' = xphys ;
mailf    = TOPOSURF tab1 ;
SI itrac ;
  TRAC mailf 'TITR' 'Maillage de la topologie finale' ;
FINSI ;
 
FIN ;