Procédures | Cast3M

* TOPOLOGY  PROCEDUR  FD218221  25/12/18    21:15:08     12429          
 
************************************************************************
** Procedure called by TOPOPTIM for updating the topology using:
**   - Optimality Criteria method (OC);
**   - Method of Moving Asymptotes (MMA).
**
** This procedure also calculates the maximum increment of the design
** variables.
**
** Author:
** Guenhael Le Quilliec (LaMe - Polytech Tours)
**
** Version:
** 3.0 2025/11/13 Improved OC, added MMA, update MAX_X_CHANGE
** 2.2 2021/03/19 Modified default value of 'OC_L2'
** 2.1 2018/01/29 Fixed issue with quadratic elements
** 2.0 2017/11/11 Initial version
************************************************************************
 
DEBP TOPOLOGY tab*'TABLE' ;
 
* General input data
* ******************
 
Wtab = tab.'WTABLE' ;
Ltab = Wtab.'LOGIQUE' ;
 
modD      = Wtab.'MODELE_D' ;
oldXD     = Wtab.'X_D' ;
oldYD     = Wtab.'Y_D' ;
NdOdXD    = Wtab.'SENSIBILITE_NORMALISEE_OBJECTIF_D' ;
NdGdXDtab = Wtab.'SENSIBILITES_NORMALISEES_CONTRAINTES_D' ;
 
* Compliance minimization under total volume constraint
* using the Optimality Criteria (Bendsoe 1995)
* *****************************************************
 
SI (EGA tab.'OPTIMISEUR' 'OC') ;
 
*   Input fields
*   ************
    matD     = Wtab.'CARACTERISTIQUES_D' ;
    uniD     = Wtab.'UNIT_D' ;
    XmaxD    = BORN (oldXD + tab.'X_CHANGE_SEUIL') 'SCAL' 'MAXI' 1.0 ;
    XminD    = BORN (oldXD - tab.'X_CHANGE_SEUIL') 'SCAL' 'MINI'
                    tab.'X_MIN_SEUIL' ;
    VtotDlim = Wtab.'VOLUME_TOTAL_D' * tab.'FRACTION_VOLUME_LIMITE' ;
*   OC allows only one single volume constraint (verified in TOPOBOOT)
    NdGdXD   = NdGdXDtab.(1) ;
 
*   Compute "specific objective sensitivity"
*       dO/dX / dG/dX = dcpl/dX / dV/dX
*   Remark: constraint derivative is expressed as dV/dX / V_total,
*           Lagrange multiplier λ adapts automatically
    dO_dGD = NdOdXD / NdGdXD ;
 
*   OC specific input data
*   **********************
    eta  = Wtab.'OC_D' ;
    gsf  = Wtab.'OC_Q' ;
    Lgsf = gsf > 1 ;
****    Lmcn = Ltab.'MECANISME' ;
*   Bounds of Lagrange multiplier
    lam1 = tab.'OC_L1' ;
    SI (EXIS tab 'OC_L2') ;
        lam2 = tab.'OC_L2' ;
    SINO ;
        lam2 = MAXI 'ABS' dO_dGD ;
    FINS ;
    lam2min = tab.'OC_L2_MIN' ;
    crt  = tab.'OC_CRITERE' ;
****     SI Lmcn ;
        Bmin = tab.'OC_B_MIN' ;
****     FINS ;
 
*   OC optimization (Bisection algorithm)
*   *************************************
 
    REPE itr tab.'OC_MAX_IT' ;
 
*       Update value of the Lagrange multiplier
        lam = 0.5 * (lam2 + lam1) ;
 
*       Compute scaled sensitivity B
        B = dO_dGD * (-1.0 / lam) ;
*       Negative values may occur for compliant mechanism synthesis
*       or elastic strain energy minimization with quadratic elements
****        SI Lmcn ;
            B = BORN B 'SCAL' 'MINI' Bmin ;
****        FINS ;
 
*       Update design variables X on D
        XD = oldXD * (B**eta) ;
 
*       Grey scale factor Q
        SI Lgsf ;
            XD = XD**gsf ;
        FINS ;
 
*       Bound X: void < X < 1, and |X - Xold| < Xmax_increment
        msk = XD MASQ 'INFE' XmaxD ;
        XD  = (XD * msk) - (XmaxD * (msk - uniD)) ;
        msk = XD MASQ 'SUPE' XminD ;
        XD  = (XD * msk) - (XminD * (msk - uniD)) ;
        Wtab.'X_D' = XD ;
 
*       Update physical densities (Y) from design variables (X) on D
        TOPOPHYS tab ;
 
*       Update new bounds of Lagrange multiplier for total volume constraint
        SI ((INTG Wtab.'Y_D' modD matD) > VtotDlim) ;
            lam1 = lam ;
        SINO ;
            lam2 = lam ;
        FINS ;
 
*       Convergence test
        SI (NON ((((lam2 - lam1) / (lam1 + lam2)) > crt) ET (lam2 > lam2min))) ;
            QUIT itr ;
        FINS ;
 
    FIN itr ;
 
FINS ;
 
* Minimization under constraints using
* Method of Moving Asymptotes (Svanberg 1987)
* *******************************************
 
SI (EGA tab.'OPTIMISEUR' 'MMA') ;
 
*   MMA specific input data
*   ***********************
 
    VD   = Wtab.'VOLUME_D' ;
    O    = Wtab.'OBJECTIF' ;
    Glst = Wtab.'CONTRAINTES' ;
 
*   MMA input table
*   ***************
 
*   Extracting list of current design variables
    Xlst = EXTR oldXD 'VALE' 'SCAL' ;
*   Number of design variables
    nbrX = DIME Xlst ;
 
*   Existing MMA table
    SI (EXIS Wtab 'MMA') ;
        mmatab = Wtab.'MMA' ;
*       If the mesh (or more precisely, the model and the frozen area)
*       have changed, then we reset the MMA table
        SI  ((NEG mmatab.'ID'.'MODEL'      Wtab.'RESOLUTION_B'.'MODELE')
          OU (NEG mmatab.'ID'.'ZONE_FIGEE' tab.'ZONE_FIGEE')) ;
            OTER Wtab 'MMA' ;
        FINS ;
    FINS ;
 
*   Initialize table for MMA operator
    SI (NON (EXIS Wtab 'MMA')) ;
        mmatab        = TABL ;
        mmatab.'XMAX' = 1.0 ;
        mmatab.'LOW'  = PROG nbrX*1.0 ;
        mmatab.'UPP'  = PROG nbrX*1.0 ;
        mmatab.'A0'   = 1.0 ;
        mmatab.'A'    = PROG 0.0 ;
        mmatab.'D'    = PROG 0.0 ;
        Wtab.'MMA' = mmatab ;
*       TODO renvoie FAUX alors que les deux valent NaN
*        SI (EGA tab.'MMA_COEF_OBJECTIF' Wtab.'NAN') ;
        SI (EGA tab.'MMA_COEF_OBJECTIF' 0.0) ;
            tab.'MMA_COEF_OBJECTIF' = 100.0 / O ;
        FINS ;
*       Save identification data
        mmatab.'ID' = TABL ;
        mmatab.'ID'.'MODEL'      = Wtab.'RESOLUTION_B'.'MODELE' ;
        mmatab.'ID'.'ZONE_FIGEE' = tab.'ZONE_FIGEE' ;
    FINS ;
 
*   If the number of design variables has changed (e.g. remeshing)
    SI (NEG nbrX (DIME mmatab.'LOW')) ;
        mmatab.'LOW' = PROG nbrX*1.0 ;
        mmatab.'UPP' = PROG nbrX*1.0 ;
    FINS ;
 
*   Set parameters
*   (or update if tab has changed between different calls to TOPOPTIM)
    mmatab.'XMIN' = tab.'X_MIN_SEUIL' ;
    mmatab.'MOVE' = tab.'X_CHANGE_SEUIL' ;
    mmatab.'C'    = tab.'MMA_C' ;
 
*   List of current design variables
    mmatab.'X' = Xlst ;
 
*   Objective function value and its sensitivity
*   (MMA objective value should not exceed 100)
    mmatab.'F0VAL' = O * tab.'MMA_COEF_OBJECTIF' ;
    SI ((mmatab.'F0VAL' > 100) OU (mmatab.'F0VAL' < 1));
        SAUT 1 'LIGN' ;
        SI (mmatab.'F0VAL' > 100) ;
            MESS (CHAI ' /!\ F0 = MMA_COEF_OBJECTIF * OBJECTIF =' ' ' mmatab.'F0VAL' ' > 100') ;
        SINO ;
            MESS (CHAI ' /!\ F0 = MMA_COEF_OBJECTIF * OBJECTIF =' ' ' mmatab.'F0VAL' ' < 1') ;
        FINS ;
        MESS (CHAI '     avec MMA_COEF_OBJECTIF =' ' ' tab.'MMA_COEF_OBJECTIF') ;
        MESS (CHAI '     or il est preferable que 1 <= F0 <= 100.') ;
*       ERRE 'Veuillez adapter la valeur du coefficient MMA_COEF_OBJECTIF.' ;
        MESS (CHAI '     Veuillez adapter la valeur du coefficient MMA_COEF_OBJECTIF.') ;
        SAUT 1 'LIGN' ;
    FINS ;
*   mmatab.'DF0DX' = (EXTR (NdOdXD * VD) 'VALE' 'SCAL') * tab.'MMA_COEF_OBJECTIF' ;
    SI (TOPOMODI tab (CHAI 'VOLUME_D_x_MMA_COEF_OBJECTIF')
                      VD tab.'MMA_COEF_OBJECTIF') ;
        Wtab.'VOLUME_D_x_MMA_COEF_OBJECTIF' = VD * tab.'MMA_COEF_OBJECTIF' ;
    FINS ;
    mmatab.'DF0DX' = (EXTR (NdOdXD * Wtab.'VOLUME_D_x_MMA_COEF_OBJECTIF') 'VALE' 'SCAL') ;
*   Constraints (fi < 0) and their sensitivity
    mmatab.'FVAL' = Glst ;
    mmatab.'DFDX' = TABL ;
    REPE itr (DIME Glst) ;
        mmatab.'DFDX'.&itr = EXTR (NdGdXDtab.&itr * VD) 'VALE' 'SCAL' ;
    FIN itr ;
 
*   Call MMA operator to update design variables
*   ********************************************
 
    MMA mmatab ;
 
*   Postprocess MMA output data
*   ***************************
 
*   MCHAML of new design variables
    XD = MANU 'CHML' modD 'REPA' 'SCAL' mmatab.'X' 'TYPE' 'SCALAIRE' 'GRAVITE' ;
 
*   In thermal mode, 'GRAVITE' option is not respected
*   This issue is resolved by using CHAN 'GRAVITE'
    SI (NON Ltab.'MECANIQUE') ;
        XD = CHAN 'GRAVITE' XD modD ;
    FINS ;
    Wtab.'X_D' = XD ;
 
*   Update physical densities from design variables
    TOPOPHYS tab ;
 
SINO ;
 
*   If MMA is not used but a table exists
    SI (EXIS Wtab 'MMA') ;
*       Delete it
        OTER Wtab 'MMA' ;
    FINS ;
 
FINS ;
 
* Other output
* ************
 
* Update maximum increment of design variables
Wtab.'MAX_X_CHANGE' = MAXI (ABS (XD - oldXD)) ;
 
* Save previous design variables and previous physical densities
Wtab.'X_D_PRECED' = oldXD ;
Wtab.'Y_D_PRECED' = oldYD ;
 
FINP ;