Procédures | Cast3M

* TOPOCONT  PROCEDUR  FD218221  25/12/18    21:15:04     12429          
 
************************************************************************
** Procedure called by TOPOPTIM to calculate:
**
** 1 - The constraints
**     Each constraint must be in the form G(X) = f(X) - f_lim <= 0,
**     with 1 < f_lim < 100, and G differentiable with respect to the
**     design variables X.
**
** 2 - Their volume-normalized sensitivities
**     Defined as dG/dX / vol1_e.
**
** Finally, any requested filtering is applied.
**
** Author:
** Guenhael Le Quilliec (LaMe - Polytech Tours)
**
** Version:
** 1.0 2025/11/10 Initial version
************************************************************************
 
DEBP TOPOCONT tab*'TABLE' ;
 
* General input data
* ******************
 
Wtab = tab.'WTABLE' ;
Ltab = Wtab.'LOGIQUE' ;
 
* Initialization
* **************
 
nbrG      = 0 ;
Glst      = PROG ;
NdGdXDtab = TABL ;
NdGdYDtab = TABL ;
 
**********************************************************************
*               Volume constraint and its sensitivity                *
**********************************************************************
 
SI Ltab.'CONTRAINTE_VOLUME' ;
 
*   Input data
*   **********
 
    uniD     = Wtab.'UNIT_D' ;
    modD     = Wtab.'MODELE_D' ;
    matD     = Wtab.'CARACTERISTIQUES_D' ;
    VtotDlim = Wtab.'VOLUME_TOTAL_D' * tab.'FRACTION_VOLUME_LIMITE' ;
 
    resB = Wtab.'RESOLUTION_B' ;
 
*   Logical for physical density dependent loading
    LdFdY = EXIS resB 'SENSIBILITE_NORMALISEE_CHARGEMENT' ;
 
*   TODO faire en sorte que ça soit compatible.
    SI LdFdY ;
        ERRE 'Actuellement la contrainte de volume n''est pas compatible avec un chargement dependant de la topologie.' ;
    FINS ;
 
*   Volume constraint (G)
*   *********************
 
*   Lets note y_e be the element value of Y.
*   The volume constraint function is:
*                    G =     vol_tot(Y)      / vol_tot_lim  - 1  <= 0
*                      =  Σ  vol_e(y_e)      / vol_tot_lim  - 1
 
*   Note: not calculated for OC optimizer, as it is not used.
    SI (NEG tab.'OPTIMISEUR' 'OC') ;
        G = ((INTG Wtab.'Y_D' modD matD) / VtotDlim) - 1.0 ;
        Glst = Glst ET (PROG G) ;
    FINS ;
 
*   Normalized sensitivity with respect to Y
*   ****************************************
 
*   Its sensitivity w.r.t. y_e is:
*              dG/dy_e =  d_vol_e(y_e)/dy_e  / vol_tot_lim
*                      = d(y_e*vol1_e)/dy_e  / vol_tot_lim
*                      =       vol1_e        / vol_tot_lim
 
*   Its normalized sensitivity w.r.t. y_e is:
*                       dG/dy_e / vol1_e = 1 / vol_tot_lim
 
    SI (TOPOMODI tab (CHAI 'UNIT_D_/_VTOTMAX') uniD VtotDlim) ;
        Wtab.'UNIT_D_/_VTOTMAX' = uniD / VtotDlim ;
    FINS ;
    NdGdYD = Wtab.'UNIT_D_/_VTOTMAX' ;
 
*   Output normalized constraint sensitivity
    nbrG = nbrG + 1 ;
    NdGdYDtab.(nbrG) = NdGdYD ;
*   TODO This line is optional but allows retrieving the name of the N-th constraint in the sensitivity table.
*    NdGdXDtab.(CHAI (DIME cst)) = CHAI 'CONTRAINTE_NORMALISEE_VOLUME' ;
 
FINS ;
 
 
**********************************************************************
*          Von Mises stress constraint and its sensitivity           *
**********************************************************************
 
*SI Ltab.'CONTRAINTE_SIGMA_VM' ;
SI (Ltab.'VERBART' OU Ltab.'CONIGLIO') ;
 
*   Input data
*   **********
 
    modMD  = Wtab.'MECANIQUE'.'MODELE_D' ;
    resA   = Wtab.'RESOLUTION_A' ;
    YD     = Wtab.'Y_D' ;
    vnmlim = tab.'SIGMA_VM_LIMITE' ;
 
*   Mesh, model, material, and densities belonging to I
*   ***************************************************
 
    mshI  = Wtab.'MAILLAGE_I' ;
    modMI = Wtab.'MECANIQUE'.'MODELE_I' ;
    matMI = Wtab.'MECANIQUE'.'CARACTERISTIQUES_I' ;
    YI    = REDU YD mshI ;
    VMI   = REDU Wtab.'VOLUME_D' mshI ;
 
*   Elastic strains on I (at centroids)
*   ***********************************
 
    SI Ltab.'PASAPAS' ;
        tmp = REDU resA.'CONTRAINTES'.(resA.'MAX_ID') mshI ;
        epsI = ELAS modMI tmp matMI ;
    SINO ;
        epsI = REDU (resA.'DEFORMATIONS') mshI ;
    FINS ;
    epsI = CHAN 'GRAVITE' epsI modMI ;
 
*   Normalized microscopic von Mises stresses on D
*   **********************************************
 
*   Microscopic stresses using full-density stiffness E(X=1) on I
    micsigI = ELAS modMI epsI matMI ;
*    micsigI = CHAN 'GRAVITE' micsigI modMI ;
*   Microscopic von Mises stresses on I
    micvnmI = CHAN 'TYPE' (VMIS modMI micsigI matMI) 'SCALAIRE' ;
*   TODO Ne conserver que le max par élément? Ça va compliquer le calcule de sa sensibilité peut être
    micvnmI = CHAN 'GRAVITE' micvnmI modMI ;
*   Normalized microscopic von Mises stresses on I
    NmicvnmI = micvnmI / vnmlim ;
*   Normalized microscopic von Mises stresses on D
*   Remark: values of any element of I missing in D are set to 0.0
    NmicvnmD = REDU NmicvnmI modMD ;
 
*   Local constraints (g) and their derivatives with respect to Y on D
*   ******************************************************************
 
*   Local constraints on the microscopic von Mises stress ratio
*   Note: use Nmicvnm instead of micvnm to avoid infinite output
    gmicvnmD = NmicvnmD - 1.0 ;
*   Local constraints to be aggregated, as proposed by Verbart et al. 2017
    gD = YD * gmicvnmD ;
*   Its maximum in D
    maxg = MAXI gD ;
 
*   If the Coniglio constraints are requested instead of Verbart
    SI Ltab.'CONIGLIO' ;
*       Local constraints to be aggregated, as proposed by Coniglio et al. 2018
        gD = gD - maxg ;
    FINS ;
 
*   Their derivatives with respect to Y
    dgdYD = gmicvnmD ;
 
*   Local constraint aggregation and its derivative with respect to g on D
*   ***********************************************************************
 
    agg daggdgD = TOPOAGRE tab gD modMD ;
 
*   Maximum von Mises stress
*   ************************
 
*   Compute maximum von Mises stress for on-the-fly display
    maxvnm = (maxg + 1.0) * vnmlim ;
 
*   **************************************************************
*   **************************************************************
 
*   Global constraint (G) on D
*   **************************
 
    SI Ltab.'VERBART' ;
*       Approximation of max(g), as proposed by Verbart et al. 2017
        G = agg ;
    SINO ;
*       Global constraint as proposed by Coniglio et al. 2018
        G = agg + maxg ;
    FINS ;
 
    Glst = Glst ET (PROG G) ;
 
*   Sensitivity / Gradient / Total derivative of G on D
*   ***************************************************
 
*   Computing the total derivative of G with respect to X
*   (using adjoint vector method)
 
*   1) Partial derivative of G with respect to the physical densities Y on D
*   ------------------------------------------------------------------------
 
*   ∂G/∂Y = ∂agg/∂g * ∂g/∂Y
    dGdYD = daggdgD * dgdYD ;
 
*   2) Partial derivative of G with respect to the displacements U on I
*   -------------------------------------------------------------------
 
*   ∂micvnm/∂U =    ∂micvnm/∂σ_mic            * ∂σ_mic/∂ε   *   ∂ε/∂U
*           = ∂((σ_micᵀ*V*σ_mic)**0.5)/∂σ_mic * ∂(D_1 ε)/∂ε * ∂(B*U)/∂U
*           =      (V*σ_mic)ᵀ / micvnm        *    D_1      *     B
 
*   ∂g/∂U   = Y * ∂gmicvnm/∂U
*           = Y/σvm_lim * ∂micvnm/∂U
*           = Y/(σvm_lim*micvnm) * (V * σ_mic)ᵀ * D_1 * B
*           = Y/(σvm_lim*micvnm) * Bᵀ * D_1 * V * σ_mic
 
*   ∂G/∂U   = ∂agg/∂g * ∂g/∂U
*           = ∂agg/∂g * Y/(σvm_lim*micvnm) * Bᵀ * D_1 * V * σ_mic
 
*   With:
*       - B     : the strain-displacement matrix
*       - D_1   : the full-density stiffness matrix (Y=1)
*       - V     : the von Mises matrix such that vnm² = σᵀ V σ
*       - σ_mic : the microscopic stresses
 
*   Remarks:
*         1) Using BSIG instead of Bᵀ produces integrated output.
*            Divide by vol1_e for expected values.
*         2) Calculated on I; zero for void elements of D.
*         3) ∂G/∂U is a nodal field. Non-zero values belong to I,
*            but can be applied to any mesh containing I.
 
    tmp = TOPODVSA micsigI modMI matMI (((REDU daggdgD mshI) * YI) /
                                        (vnmlim * micvnmI * VMI)) ;
    dGdUI = BSIG modMI tmp ;
 
*   3) Normalized total derivative of G with respect to the physical densities Y on D
*   ---------------------------------------------------------------------------------
 
*   Compute the volume-normalized total derivative of G w.r.t. Y,
*   via the adjoint method
*   dG/dY / vol1_e = (∂G/∂Y + ∂G/∂U * ∂U/∂Y) / vol1_e
    NtotdGdYD = TOPOADJO tab dGdYD dGdUI ;
 
*   Output normalized constraint sensitivity
    nbrG = nbrG + 1 ;
    NdGdYDtab.(nbrG) = NtotdGdYD ;
*   TODO This line is optional but allows retrieving the name of the N-th constraint in the sensitivity table.
*    NdGdXDtab.(CHAI (DIME cst)) = CHAI 'CONTRAINTE_NORMALISEE_SIGMA_VM' ;
 
FINS ;
*FINS ;
 
 
**********************************************************************
 
* Constraint sensitivities with respect to design variables X on D
* ****************************************************************
 
*   As Y=X when X is not projected nor filtered, dG/dx_e = dG/dy_e.
*   Otherwise Y≠X and we need to calculate dG/dx_e from dG/dy_e:
*              dG/dx_e = dG/dy_e * dy_e/dx_e
*                      = dG/dy_e * dy_e/ds_e * ds_e/dx_e
 
*   dG/dy_e / vol1_e = dG/dy_e * dy_e/dx_e / vol1_e
REPE itr nbrG ;
    NdGdXDtab.&itr = TOPODYDX tab NdGdYDtab.&itr ;
FIN itr ;
 
* Constraint sensitivities filtering
* **********************************
 
SI Ltab.'FILTRE_SENSIBILITES_CONTRAINTES' ;
    REPE itr nbrG ;
        tmp = NdGdXDtab.&itr * Wtab.'Y_D' ;
        tmp = TOPOFILT tab tmp ;
        NdGdXDtab.&itr = tmp /
               (BORN Wtab.'Y_D' 'SCAL' 'MINI' tab.'ZERO_DIVISION') ;
    FIN itr ;
FINS ;
 
* Output data
* ***********
 
Wtab.'CONTRAINTES' = Glst ;
Wtab.'SENSIBILITES_NORMALISEES_CONTRAINTES_D' = NdGdXDtab ;
 
* In the specific case of stress constraint
SI Ltab.'CONTRAINTE_SIGMA_VM' ;
*   Save the maximum von Mises stress
    Wtab.'SIGMA_VM_MAX' = maxvnm ;
*   Save Microscopic von Mises stresses on I for on-the-fly field
*   visualization if TRAC is activated
    SI tab.'TRAC' ;
        Wtab.'MICROSCOPIC_VM_I' = micvnmI ;
    FINS ;
FINS ;
 
FINP ;