Exemples de jeux de données Cast3M

************************************************************************
* Exemple de méthode d'optimisation topologique                        *
* Méthode à densité + pénalisation SIMP                                *
* Algorithme d'optimisation : Optimality Criteria - Bendsoe            *
* Application à une poutre en flexion                                  *
************************************************************************
 
 
** Paramètres globaux
itrac  = FAUX ;
OPTI 'DIME' 2 'MODE' 'PLAN' 'DEFO' 'ELEM' 'QUA4' 'ECHO' 0 ;
l      = 1. ;
h      = 0.2 ;
yo     = 200.E9 ;
nu     = 0.3 ;
 
** Maillage (controlé par sa taille de maille)
den1   = l / 120. ;
OPTI 'DENS' den1 ;
p1     = 0. 0. ;
p2     = l  0. ;
lb     = DROI p1 p2 ;
mail   = lb TRAN (0. h) ;
p3     = mail POIN 'PROC' ((l / 2.) h) ;
con    = CONT mail ;
 
** Modèle
mod    = MODE mail 'MECANIQUE' ;
ma0    = MATE mod 'YOUN' yo 'NU' nu ;
 
** Blocages
blo    = (BLOQ 'UY' p1) ET (BLOQ 'UX' 'UY' p2) ;
 
** Chargement
f      = FORC (0. -1.E5) p3 ;
 
** Paramètres d'optimisation
p      = 3. ;
fv     = 0.4 ;
eta    = 0.5 ;
rmin   = 1.5 * den1 ;
m      = 0.1 ;
xmin   = 0.001 ;
xmax   = 1. ;
 
** Initialisation de la topologie (avec la fraction volumique cible)
x      = MANU 'CHML' mod 'SCAL' fv 'GRAVITE' ;
 
** Matrice de filtrage
un     = MANU 'CHML' mod 'SCAL' 1. 'GRAVITE' ;
vole   = INTG mod un 'ELEM' ;
mcg    = un POIN 'SUPERIEUR' 0. ;
wg     = PROI mcg vole ;
kfil   = MFIL wg rmin 1. 0. ;
 
** Volume plein, initial et cible
v0     = INTG mod un ;
vx     = INTG mod x ;
fvx    = vx / v0 ;
vcib   = v0 * fv ;
chgx   = 0. ;
 
** Boucle d'optimisation topologique
liso   = PROG 0. 'PAS' 0.05 1. ;
REPE b1 100 ;
* pénalisation de la rigidité
  yop    = (x ** p) * yo ;
  map    = MATE mod 'YOUN' yop 'NU' nu ;
* résolution du problème mécanique
  rip    = RIGI mod map ;
  u      = RESO (rip ET blo) f ;
* fonction objectif : compliance = uT.K.u = Int(sig:eps)dV
  eps    = EPSI 'LINE' mod u ;
  sig    = ELAS mod map eps ;
  psi    = INTG mod (ENER mod sig eps) ;
* sensibilité
  sig0   = ELAS mod ma0 eps ;
  ene0   = ENER mod eps sig0 ;
  ene    = CHAN ene0 mod 'GRAVITE' ;
  dpsi   = (-1. * p * (x ** (p - 1.)) * ene) ;
* filtrage de la sensibilité
  xdpsi  = x * dpsi ;
  xdpsi0 = PROI mcg xdpsi ;
  xdpsi1 = kfil * xdpsi0 ;
  xdpsi  = MANU 'CHML' mod 'REPA' 'SCAL' (EXTR xdpsi1 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
  dpsi   = xdpsi / x ;
* infos sur la topologie courante
  info   = CHAI 'It:' (&b1 - 1) / 5 'Obj:' / 10  psi > 1 'Fvol:' > 4 fvx > 1 'Change:' > 4 chgx > 1 ;
  MESS info ;
  SI itrac ;
    def1   = DEFO mail u 1000. ;
    TRAC x mod con def1 liso 'TITR' info 'NCLK' ;
  FINSI ;
* optimisation d'une nouvelle topologie (méthode du critère d'optimalité)
  l1     = 0. ;
  l2     = 100000000. ;
  REPE b2 100 ;
    SI ((l2 - l1) < 0.0001) ;
      QUIT b2 ;
    FINSI ;
    lmid   = (l1 + l2) / 2. ;
    b      = -1. * dpsi / (lmid * vole) ;
    xinf   = BORN (x - m) 'MINIMUM' xmin ;
    xsup   = BORN (x + m) 'MAXIMUM' xmax ;
    xnew   = x * (b ** eta) ;
    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) ;
    vxnew  = INTG mod xnew ;
    SI (vxnew > vcib) ;
      l1     = lmid ;
    SINON ;
      l2     = lmid ;
    FINSI ;
  FIN b2 ;
* bilan de l'itération
  fvx    = vxnew / v0 ;
  chgx   = MAXI 'ABS' (x - xnew) ;
  SI (chgx < 0.01) ;
    info   = CHAI 'It:' &b1 / 5 'Obj:' / 10  psi > 1 'Fvol:' > 4 fvx > 1 'Change:' > 4 chgx > 1 ;
    MESS info ;
    QUIT b1 ;
  FINSI ;
* préparation de la nouvelle itération
  x      = xnew ;
FIN  b1 ;
 
* Tracé de la dernière itération
SI itrac ;
  def1   = DEFO mail u 1000. ;
  TRAC xnew mod con def1 liso 'TITR' info ;
FINSI ;
 
FIN ;