* fichier topoptim_10_Non-linear.dgibi
************************************************************************
* Section : Mathematiques Fonctions
************************************************************************

************************************************************************
** Non-linear topology optimization of a simple 2D structure in contact
** with a rigid sphere with prescribed displacement.
**
** Author:
** Guenhael Le Quilliec (LaMe - Polytech Tours)
**
** Version:
** V2.0 2025/11/13 Updated to make it compatible with TOPOPTIM V4.0
** V1.0 2017/04/18 Original version compatible with TOPOPTIM V2.1
************************************************************************

* Plot results
graph0 = FAUX ;

* General options
OPTI 'DIME' 2 'MODE' 'AXIS' 'ELEM' 'QUA4' ;

* Number of elements
nelp0 = 30 ;
nelr0 = 30 ;
nelz0 = 30 ;

* Dimensions
r0 = 20.0e-3 ;
l1 = nelr0 * 1.0e-3 ;
l2 = nelz0 * 1.0e-3 ;

* Bulk mesh
p0 = 0.0 0.0 ;
p1 = 0.0 (-1.0 * l2) ;
p2 = l1 (-1.0 * l2) ;
p3 = l1 0.0 ;
ln0 = DROI nelz0 p0 p1 ;
ln1 = DROI nelr0 p1 p2 ;
ln2 = DROI nelz0 p2 p3 ;
ln3 = DROI nelr0 p3 p0 ;
msh0 = SURF (ln0 ET ln1 ET ln2 ET ln3) ;

* Indent mesh
p5 = 0.0 r0 ;
p6 = r0 r0 ;
arc0 = CERC nelp0 p0 p5 p6 ;

* Model and material
mod0 = MODE msh0 'MECANIQUE' 'ELASTIQUE' ;
mat0 = MATE mod0 'YOUN' 0.1e9 'NU' 0.4 ;

* Time steps
nbstp0 = 20 ;
time0 = 2.0 ;
lst_tps0 = PROG 0.0 'PAS' (time0 / (FLOT nbstp0)) time0 ;

* Boundary conditions
bc1 = (BLOQ 'UR' ln0) ET (BLOQ 'UZ' ln1) ET (BLOQ 'UR' arc0) ;
bc2 = BLOQ 'UZ' arc0 ;

* Loading
evol0 = EVOL 'MANU' 'Temps' lst_tps0 'Coef'
              (PROG 0.0 'PAS' (1.0 / (FLOT nbstp0)) 1.0) ;
load0 = CHAR 'DIMP' (DEPI bc2 -2.0e-3) evol0 ;

* Contact
modct0 = MODE arc0 'CONTACT' 'FROTTANT' 'COULOMB'  ln3;
matct0 = MATE modct0 'MU' 0.5 ;

* Factor evolutions
pfct0 = EVOL 'MANU' 'CYCLES' (PROG 1  6) 'P' (PROG 1.0  3.0) ;
qfct0 = EVOL 'MANU' 'CYCLES' (PROG 2 32) 'Q' (PROG 1.0 10.0) ;

* Finite element model table
mdl0 = TABL ;
mdl0.'MODELE'                 = mod0 ET modct0 ;
mdl0.'CARACTERISTIQUES'       = mat0 ET matct0 ;
mdl0.'BLOCAGES_MECANIQUES'    = bc1 ET bc2 ;
mdl0.'CHARGEMENT'             = load0 ;
mdl0.'TEMPS_CALCULES'         = lst_tps0 ;

* Optimization table
tab0 = TABL ;
tab0.'TRAC'                   = graph0 ;
tab0.'RESOLUTION_PASAPAS'     = mdl0 ;
tab0.'FRACTION_VOLUME_LIMITE' = 0.2 ;
tab0.'PENALISATION_P'         = pfct0 ;
tab0.'OC_Q'                   = qfct0 ;
tab0.'SEUIL_ELEMENTS_ACTIFS'  = 1.0e-2 ;
tab0.'MES_SAUVEGARDES'        = TABL ;
tab0.'MES_SAUVEGARDES'.'RESOLUTION' = VRAI ;
*tab0.'FILTRE'                 = MOT 'EDP' ;
*tab0.'FILTRE_EDP_RAYON'       = 2.0 * l1 / (FLOT nelr0) ;
*tab0.'FILTRE'                 = MOT 'CHAPEAU' ;
*tab0.'FILTRE_CHAPEAU_RAYON'   = 2.5 * l1 / (FLOT nelr0) ;
*tab0.'FILTRER'                = MOT 'DENSITE' ;
*tab0.'OPTIMISEUR'             = MOT 'MMA' ;

* Optimization
TOPOPTIM tab0 ;

* Output data
topo0 = tab0.'TOPOLOGIE'.(tab0.'CYCLE') ;
topomsh0 = tab0.'MAILLAGE'.(tab0.'CYCLE') ;
res0 = tab0.'RESOLUTION'.(tab0.'CYCLE') ;

* Deform the mesh
DEPL (msh0 ET arc0) 'PLUS' res0.'DEPLACEMENTS'.(1) ;

* Plot final deformed mesh
* and final topology (physical density)
SI graph0 ;
    TRAC (topomsh0 ET (arc0 COUL 'BLEU')) ;
    TRAC (REDU topo0 topomsh0) (REDU mod0 topomsh0)
         (PROG 0.0 'PAS' (1.0 / 56.0) 1.0)
         'TITR' 'Topologie finale' ;
FINS ;

* Plot output evolutions
SI graph0 ;
    REPE loop0 (DIME tab0.'EVOLUTIONS_SAUVEES') ;
        evoname0 = EXTR tab0.'EVOLUTIONS_SAUVEES' &loop0 ;
        DESS tab0.'EVOLUTIONS'.evoname0 'POSY' 'EXCE'
             'TITR' (CHAI 'Evolution de' ' ' evoname0
                          ' au cours des cycles d''optimisation') ;
    FIN loop0 ;
FINS ;

FIN ;