$$$$ TOPORSTR
* TOPORSTR  PROCEDUR  FD218221  25/12/18    21:15:13     12429          

************************************************************************
** Procedure called by TOPOOBJE to apply geometric restrictions to an
** input field.
**
** Remark: the input field must be an intensive quantity
**         and defined on the design zone (D).
**
** Author:
** Guenhael Le Quilliec (LaMe - Polytech Tours)
**
** Version:
** 2.0 2025/11/13 Made compatible with TOPOPTIM V4.0
** 1.1 2021/03/19 Fixed a minor issue by adding 'SCALAIRE' type for the
**                sensitivity
** 1.0 2017/04/18 Initial version
************************************************************************

DEBP TOPORSTR tab*'TABLE' fldD*'MCHAML' ;

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

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

* If there is nothing to do
* *************************

SI (NON Ltab.'RESTRICTIONS') ;
*   Return input field
    QUIT TOPORSTR ;
FINS ;

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

* Criterion to fuse the nodes
crt = tab.'ELIM_CRITERE' ;
* Node mesh of the centroids
ctrmshD = Wtab.'MAILLAGE_GRAVITE_D' ;

* Preprocessing
* *************

* Switch from element field to nodal field
fldD = PROI ctrmshD fldD ;

* Apply restrictions
* ******************

* Loop over each restriction
REPE itr (DIME tab.'RESTRICTIONS') ;
*   Find the next restriction
    SI (EXIS tab.'RESTRICTIONS' &itr) ;
        rst = tab.'RESTRICTIONS'.&itr ;
    SINO ;
        QUIT itr ;
    FINS ;
*   Get node mesh of centroids of the restricted zone
*   and its initial sensitivity
    Lrstzon = EXIS rst 'ZONE' ;
    SI Lrstzon ;
        ctrmsh1 = (REDU Wtab.'UNIT_D' rst.'ZONE') POIN 'SUPE' -1.0 ;
        ELIM ctrmshD ctrmsh1 crt ;
        fld1 = REDU fldD ctrmsh1 ;
    SINO ;
        ctrmsh1 = ctrmshD ;
        fld1 = fldD ;
    FINS ;
*   Make copies with prescribed geometric transformation
*   Fuse nodes of their supports
*   Compute mean sensitivity fields
    SI (EGA rst.'TYPE' 'SYME_POIN') ;
        msh2 fld2 = ctrmsh1 fld1 SYME 'POIN' rst.'POIN1' ;
        ELIM ctrmsh1 msh2 crt ;
        fld1 = (fld1 + fld2) / 2.0 ;
    FINS ;
    SI (EGA rst.'TYPE' 'SYME_DROI') ;
        msh2 fld2 = ctrmsh1 fld1 SYME 'DROI' rst.'POIN1' rst.'POIN2' ;
        ELIM ctrmsh1 msh2 crt ;
        fld1 = (fld1 + fld2) / 2.0 ;
    FINS ;
    SI (EGA rst.'TYPE' 'SYME_PLAN') ;
        msh2 fld2 = ctrmsh1 fld1 SYME 'PLAN' rst.'POIN1' rst.'POIN2' rst.'POIN3' ;
        ELIM ctrmsh1 msh2 crt ;
        fld1 = (fld1 + fld2) / 2.0 ;
    FINS ;
    SI (EGA rst.'TYPE' 'PERIODICITE_CIRCULAIRE') ;
        SI (rst.'PERIODES' EGA 2) ;
            SI ((VALE 'DIME') EGA 3) ;
                msh2 fld2 = ctrmsh1 fld1 TOUR 180.0 rst.'POIN1' rst.'POIN2' ;
            SINO ;
                msh2 fld2 = ctrmsh1 fld1 TOUR 180.0 rst.'POIN1' ;
            FINS ;
            ELIM ctrmsh1 msh2 crt ;
            fld1 = (fld1 + fld2) / 2.0 ;
        SINO ;
            mshtab = TABL ;
            mshtab.(0) = ctrmsh1 ;
            fldtab = TABL ;
            fldtab.(0) = fld1 ;
            p0 = rst.'PERIODES' ;
            a0 = 360.0 / p0 ;
            SI ((VALE 'DIME') EGA 3) ;
                p1 = rst.'POIN1' ;
                p2 = rst.'POIN2' ;
                REPE itr1 (p0 - 1) ;
                    old = &itr1 - 1 ;
                    mshtab.&itr1 fldtab.&itr1 = mshtab.old fldtab.old TOUR a0 p1 p2 ;
                FIN itr1 ;
            SINO ;
                p1 = rst.'POIN1' ;
                REPE itr1 (p0 - 1) ;
                    old = &itr1 - 1 ;
                    mshtab.&itr1 fldtab.&itr1 = mshtab.old fldtab.old TOUR a0 p1 ;
                FIN itr1 ;
            FINS ;
            REPE itr1 (p0 - 1) ;
                ELIM ctrmsh1 mshtab.&itr1 crt ;
                fld1 = fld1 + fldtab.&itr1 ;
            FIN itr1 ;
            fld1 = fld1 / p0 ;
        FINS ;
    FINS ;
    SI (EGA rst.'TYPE' 'PERIODICITE_AXIALE') ;
        mshtab = TABL ;
        mshtab.(0) = ctrmsh1 ;
        fldtab = TABL ;
        fldtab.(0) = fld1 ;
        p0 = rst.'PERIODES' ;
        v0 = rst.'POIN2' MOIN rst.'POIN1' ;
        REPE itr1 (p0 - 1) ;
            old = &itr1 - 1 ;
            mshtab.&itr1 fldtab.&itr1 = mshtab.old fldtab.old PLUS v0 ;
            itr2 = -1 * &itr1 ;
            old = itr2 + 1 ;
            mshtab.itr2 fldtab.itr2 = mshtab.old fldtab.old MOIN v0 ;
        FIN itr1 ;
        REPE itr1 (p0 - 1) ;
            ELIM ctrmsh1 mshtab.&itr1 crt ;
            fld1 = fld1 + fldtab.&itr1 ;
            itr2 = -1 * &itr1 ;
            ELIM ctrmsh1 mshtab.itr2 crt ;
            fld1 = fld1 + fldtab.itr2 ;
        FIN itr1 ;
        fld1 = fld1 / p0 ;
    FINS ;
*   Compute final sensitivity
    SI Lrstzon ;
        fldD = fldD - (REDU fldD ctrmsh1) + fld1 ;
    SINO ;
        fldD = fld1 ;
    FINS ;
FIN itr ;

* Postprocessing
* **************

* Switch sensitivity back from nodal field to element field
fldD = MANU 'CHML' Wtab.'MODELE_D' 'REPA' 'SCAL' (EXTR fldD 'VALE')
           'TYPE' 'SCALAIRE' 'GRAVITE' ;
* In thermal mode, 'GRAVITE' option is not respected
* This issue is resolved by using CHAN 'GRAVITE'
SI (NON Ltab.'MECANIQUE') ;
    fldD = CHAN 'GRAVITE' fldD Wtab.'MODELE_D' ;
FINS ;

FINP fldD ;

$$$$