$$$$ TOPOOBJE
* TOPOOBJE  PROCEDUR  FD218221  25/12/18    21:15:09     12429          

************************************************************************
** Procedure called by TOPOPTIM to calculate:
**
** 1 - The objective function, which can be:
**     - Mechanism output displacement maximization (based on
**       compliance);
**     - Compliance minimization (or more precisely, "2 * the total
**       elastic strain energy" in mechanics and "total thermal
**       compliance" in thermal problems);
**     - Total volume fraction.
**
** 2 - The volume-normalized objective sensitivity (w.r.t X)
**     Defined as ∂obj/∂X / vol1_e
**
**     For compliance:
**      - ∂cpl/∂X / vol1_e = - dZ/dX * q1_e / vol1_e
**     Only the explicit part is included (ignoring U(X) dependence).
**     The minus sign compensates for the missing implicit contribution,
**     ensuring the optimizer moves densities in the correct direction.
**
** Finally, any requested restrictions and/or filtering are applied.
**
** Author:
** Guenhael Le Quilliec (LaMe - Polytech Tours)
**
** Version:
** 2.0 2025/11/17 Renamed TOPOOBJE and made compatible with
**                TOPOPTIM V4.0. Added total volume fraction objective.
**                Call TOPORSTR and filter the objective sensitivity at
**                the end of TOPOOBJE.
** 1.0 2021/04/15 Initial version of TOPOSENS
************************************************************************

DEBP TOPOOBJE tab*'TABLE' ;

* General input data
* ******************

Wtab = tab.'WTABLE' ;
Ltab = Wtab.'LOGIQUE' ;


**********************************************************************
*                   Compliance and its sensitivity                   *
**********************************************************************

SI Ltab.'OBJECTIF_COMPLIANCE' ;

*   Input data
*   **********

    wgtmcn = tab.'POIDS_MECANISME' ;
    wgtM   = tab.'POIDS_ENERGIE_DEFO' ;
    wgtT   = tab.'POIDS_TEMPERATURE' ;
    ctrmshD = Wtab.'MAILLAGE_GRAVITE_D' ;
    resA    = Wtab.'RESOLUTION_A' ;
    SI Ltab.'CAS_MULTIPLES' ;
        nbrcas = Wtab.'NB_CAS' ;
    FINS ;
    SI Ltab.'MECANIQUE' ;
        modMD = Wtab.'MECANIQUE'.'MODELE_D' ;
    FINS ;
    SI Ltab.'THERMIQUE' ;
        modTD = Wtab.'THERMIQUE'.'MODELE_D' ;
    FINS ;

*   Mesh, model, material, and "physical density derivative factor" (dZ/dY / Z) belonging to I
*   ******************************************************************************************

    mshI = Wtab.'MAILLAGE_I' ;
    SI Ltab.'MECANIQUE' ;
        modMI    = Wtab.'MECANIQUE'.'MODELE_I' ;
        ZmatMI   = Wtab.'MECANIQUE'.'CARACTERISTIQUES_Z_I' ;
        dZdY_ZMI = Wtab.'MECANIQUE'.'SENSIBILITE_Z_I' /
                   Wtab.'MECANIQUE'.'Z_I' ;
    FINS ;
    SI Ltab.'THERMIQUE' ;
        modTI    = Wtab.'THERMIQUE'.'MODELE_I' ;
        ZmatTI   = Wtab.'THERMIQUE'.'CARACTERISTIQUES_Z_I' ;
        dZdY_ZTI = Wtab.'THERMIQUE'.'SENSIBILITE_Z_I' /
                   Wtab.'THERMIQUE'.'Z_I' ;
    FINS ;

*   Initialization
*   **************

    O = 0.0 ;
    dOdYI = Wtab.'VIDE' ;

*   Compliance and its sensitivity for mechanisms
*   *********************************************

    SI Ltab.'MECANISME' ;
        SI Ltab.'CAS_MULTIPLES' ;
            REPE itr nbrcas ;
*               Handle and response stress and strain
                SI Ltab.'PASAPAS' ;
                    tmp = resA.&itr.'MECANISME_ENTREE'.'MAX_ID' ;
                    siginpI = REDU (resA.&itr.'MECANISME_ENTREE'.'CONTRAINTES'.tmp) mshI ;
                    epsinpI = ELAS modMI siginpI ZmatMI ;
                    tmp = resA.&itr.'MECANISME_SORTIE'.'MAX_ID' ;
                    sigoutI = REDU (resA.&itr.'MECANISME_SORTIE'.'CONTRAINTES'.tmp) mshI ;
*                    epsoutI = ELAS modMI sigoutI ZmatMI ;
                SINO ;
*                    siginpI = REDU (resA.&itr.'MECANISME_ENTREE'.'CONTRAINTES') mshI ;
                    epsinpI = REDU (resA.&itr.'MECANISME_ENTREE'.'DEFORMATIONS') mshI ;
                    sigoutI = REDU (resA.&itr.'MECANISME_SORTIE'.'CONTRAINTES') mshI ;
*                    epsoutI = REDU (resA.&itr.'MECANISME_SORTIE'.'DEFORMATIONS') mshI ;
                FINS ;
                cpl dcpldYI = TOPOCMPL modMI ZmatMI dZdY_ZMI sigoutI epsinpI ;
                O = O + ((wgtmcn / nbrcas) * cpl) ;
                SI (EGA dOdYI Wtab.'VIDE') ;
                    dOdYI = ((wgtmcn / nbrcas) * dcpldYI) ;
                SINO ;
                    dOdYI = dOdYI + ((wgtmcn / nbrcas) * dcpldYI) ;
                FINS ;
            FIN itr ;
        SINO ;
*           Handle and response stress and strain
            SI Ltab.'PASAPAS' ;
                tmp = resA.'MECANISME_ENTREE'.'MAX_ID' ;
                siginpI = REDU (resA.'MECANISME_ENTREE'.'CONTRAINTES'.tmp) mshI ;
                epsinpI = ELAS modMI siginpI ZmatMI ;
                tmp = resA.'MECANISME_SORTIE'.'MAX_ID' ;
                sigoutI = REDU (resA.'MECANISME_SORTIE'.'CONTRAINTES'.tmp) mshI ;
*                epsoutI = ELAS modMI sigoutI ZmatMI ;
            SINO ;
*                siginpI = REDU (resA.'MECANISME_ENTREE'.'CONTRAINTES') mshI ;
                epsinpI = REDU (resA.'MECANISME_ENTREE'.'DEFORMATIONS') mshI ;
                sigoutI = REDU (resA.'MECANISME_SORTIE'.'CONTRAINTES') mshI ;
*                epsoutI = REDU (resA.'MECANISME_SORTIE'.'DEFORMATIONS') mshI ;
            FINS ;
            cpl dcpldYI = TOPOCMPL modMI ZmatMI dZdY_ZMI sigoutI epsinpI ;
            O = O + (wgtmcn * cpl) ;
            SI (EGA dOdYI Wtab.'VIDE') ;
                dOdYI = (wgtmcn * dcpldYI) ;
            SINO ;
                dOdYI = dOdYI + (wgtmcn * dcpldYI) ;
            FINS ;
        FINS ;
    FINS ;

*   Compliance and its sensitivity
*   ******************************

*    SI (Ltab.'MECANIQUE' ET ((ABS wgtM) > tab.'PRECISION')) ;
    SI (Ltab.'MECANIQUE' ET (NEG wgtM 0.0)) ;
        SI Ltab.'CAS_MULTIPLES' ;
            REPE itr nbrcas ;
                SI Ltab.'PASAPAS' ;
                    tmp = resA.&itr.'MAX_ID' ;
                    sigI = REDU (resA.&itr.'CONTRAINTES'.tmp) mshI ;
                    epsI = ELAS modMI sigI ZmatMI ;
                SINO ;
                    sigI = REDU (resA.&itr.'CONTRAINTES') mshI ;
                    epsI = REDU (resA.&itr.'DEFORMATIONS') mshI ;
                FINS ;
                cpl dcpldYI = TOPOCMPL modMI ZmatMI dZdY_ZMI sigI epsI ;
                O = O + ((wgtM / nbrcas) * cpl) ;
                SI (EGA dOdYI Wtab.'VIDE') ;
                    dOdYI = ((-1.0 * wgtM / nbrcas) * dcpldYI) ;
                SINO ;
                    dOdYI = dOdYI - ((wgtM / nbrcas) * dcpldYI) ;
                FINS ;
            FIN itr ;
        SINO ;
            SI Ltab.'PASAPAS' ;
                tmp = resA.'MAX_ID' ;
                sigI = REDU (resA.'CONTRAINTES'.tmp) mshI ;
                epsI = ELAS modMI sigI ZmatMI ;
            SINO ;
                sigI = REDU (resA.'CONTRAINTES') mshI ;
                epsI = REDU (resA.'DEFORMATIONS') mshI ;
            FINS ;
            cpl dcpldYI = TOPOCMPL modMI ZmatMI dZdY_ZMI sigI epsI ;
            O = O + (wgtM * cpl) ;
            SI (EGA dOdYI Wtab.'VIDE') ;
                dOdYI = ((-1.0 * wgtM) * dcpldYI) ;
            SINO ;
                dOdYI = dOdYI - (wgtM * dcpldYI) ;
            FINS ;
        FINS ;
    FINS ;

*   Thermal compliance and its sensitivity
*   **************************************

*    SI (Ltab.'THERMIQUE' ET ((ABS wgtT) > tab.'PRECISION')) ;
    SI (Ltab.'THERMIQUE' ET (NEG wgtT 0.0)) ;
        SI Ltab.'CAS_MULTIPLES' ;
            REPE itr nbrcas ;
*               Temperature nodal field
                SI Ltab.'PASAPAS' ;
                    tmp = resA.&itr.'MAX_ID' ;
                    TA = REDU (resA.&itr.'TEMPERATURES'.tmp) mshI ;
                SINO ;
                    TA = REDU (resA.&itr.'TEMPERATURES') mshI ;
                FINS ;
                cpl dcpldYI = TOPOCMPL modTI ZmatTI dZdY_ZTI TA ;
*               Switch the model of the element field eneK from Therm to Mecha
                SI Ltab.'MECANIQUE' ;
*                   Remark: values of any element of I missing in D will be set to 0.0
                    dcpldYD = TOPOCHAN (REDU dcpldYI modTD) modMD ctrmshD ;
                    dcpldYI = REDU dcpldYD modMI ;
                FINS ;
                O = O + ((wgtT / nbrcas) * cpl) ;
                SI (EGA dOdYI Wtab.'VIDE') ;
                    dOdYI = ((-1.0 * wgtT / nbrcas) * dcpldYI) ;
                SINO ;
                    dOdYI = dOdYI - ((wgtT / nbrcas) * dcpldYI) ;
                FINS ;
            FIN itr ;
        SINO ;
*           Temperature nodal field
            SI Ltab.'PASAPAS' ;
                tmp = resA.'MAX_ID' ;
                TA = REDU (resA.'TEMPERATURES'.tmp) mshI ;
            SINO ;
                TA = REDU (resA.'TEMPERATURES') mshI ;
            FINS ;
            cpl dcpldYI = TOPOCMPL modTI ZmatTI dZdY_ZTI TA ;
*           Switch the model of the element field eneK from Therm to Mecha
            SI Ltab.'MECANIQUE' ;
*               Remark: values of any element of I missing in D will be set to 0.0
                dcpldYD = TOPOCHAN (REDU dcpldYI modTD) modMD ctrmshD ;
                dcpldYI = REDU dcpldYD modMI ;
            FINS ;
            O = O + (wgtT * cpl) ;
            SI (EGA dOdYI Wtab.'VIDE') ;
                dOdYI = ((-1.0 * wgtT) * dcpldYI) ;
            SINO ;
                dOdYI = dOdYI - (wgtT * dcpldYI) ;
            FINS ;
        FINS ;
    FINS ;

*   Get the objective sensitivity with respect to physical densities Y on D
*   ***********************************************************************

*   If it is empty
    SI (EGA dOdYI Wtab.'VIDE') ;
*       Return a zero field
        dOdYD = Wtab.'ZERO_D' ;
    SINO ;
*       Use REDU to get O on D from I.
*       Remark: values of any element of D missing in I will be set to 0.0
        dOdYD = REDU dOdYI Wtab.'MODELE_D' ;
    FINS ;

*   Normalization
*   *************

*   Volume-normalized objective sensitivity
    NdOdYD = dOdYD / Wtab.'VOLUME_D' ;

FINS ;


**********************************************************************
*                  Volume ratio and its sensitivity                  *
**********************************************************************

SI Ltab.'OBJECTIF_VOLUME' ;

*   Input data
*   **********

    uniD  = Wtab.'UNIT_D' ;
    modD  = Wtab.'MODELE_D' ;
    matD  = Wtab.'CARACTERISTIQUES_D' ;
    VtotD = Wtab.'VOLUME_TOTAL_D' ;

*   Objective to be minimized
*   *************************

*   Total material volume normalized by the total element volume
    SI (TOPOMODI tab (CHAI 'FRACTION_VOLUME') Wtab.'Y_D' modD matD) ;
        Wtab.'FRACTION_VOLUME' = (INTG Wtab.'Y_D' modD matD) / VtotD ;
    FINS ;
    O = Wtab.'FRACTION_VOLUME' ;

*   Volume-normalized objective sensitivity with respect to Y
    SI (TOPOMODI tab (CHAI 'UNIT_D_/_VTOT_D') uniD VtotD) ;
        Wtab.'UNIT_D_/_VTOT_D' = uniD / VtotD ;
    FINS ;
    NdOdYD = Wtab.'UNIT_D_/_VTOT_D' ;

FINS ;


**********************************************************************
*        Postprocessing of the objective and its sensitivity         *
**********************************************************************

* Correct the sensitivity
* ***********************

* Modify sensitivity to get derivative with respect to X instead of Y
NdOdXD = TOPODYDX tab NdOdYD ;

* Add any requested restrictions
* ******************************

NdOdXD = TOPORSTR tab NdOdXD ;

* Objective sensitivity filtering
* *******************************

SI Ltab.'FILTRE_SENSIBILITE_OBJECTIF' ;
    tmp = NdOdXD * Wtab.'Y_D' ;
    tmp = TOPOFILT tab tmp ;
    NdOdXD = tmp / (BORN Wtab.'Y_D' 'SCAL' 'MINI' tab.'ZERO_DIVISION') ;
FINS ;

* Save output data
* ****************

* Save the compliance value and its sensitivity field before applying restrictions
Wtab.'OBJECTIF' = O ;
Wtab.'SENSIBILITE_NORMALISEE_OBJECTIF_D' = NdOdXD ;

FINP ;

$$$$