Télécharger topoadjo.procedur

Retour à la liste

* TOPOADJO  PROCEDUR  FD218221  25/12/18    21:15:02     12429          
 
************************************************************************
** Procedure called by TOPOCONT to compute the total derivative of the
** total constraint (G) with respect to the design variables, using the
** partial derivatives with respect to the physical densities and the
** adjoint method.
**
** If the loading is physical density dependent, its volume-normalized
** sensitivity with respect to the design variables should also be
** provided on the design zone (D).
**
** Author:
** Guenhael Le Quilliec (LaMe - Polytech Tours)
**
** Version:
** 1.0 2025/11/21 Initial version
************************************************************************
 
DEBP TOPOADJO tab*'TABLE' dGdYD*'MCHAML' dGdUI*'CHPOINT' ;
 
* Input data
* **********
 
Wtab   = tab.'WTABLE' ;
*Ltab   = Wtab.'LOGIQUE' ;
 
VD     = Wtab.'VOLUME_D' ;
 
mshI   = Wtab.'MAILLAGE_I' ;
modMD  = Wtab.'MECANIQUE'.'MODELE_D' ;
modMA  = Wtab.'MECANIQUE'.'MODELE_A' ;
modMI  = Wtab.'MECANIQUE'.'MODELE_I' ;
ZmatMA = Wtab.'MECANIQUE'.'CARACTERISTIQUES_Z_A' ;
bdcMB  = Wtab.'MECANIQUE'.'BLOCAGES' ;
 
resB = Wtab.'RESOLUTION_B' ;
 
UI     = REDU Wtab.'RESOLUTION_A'.'DEPLACEMENTS' mshI ;
 
* Initialization
* **************
 
* Derivative of the material characteristics with respect to Y
dZmatMI = TOPOMATE modMI Wtab.'MECANIQUE'.'CARACTERISTIQUES_I'
                         Wtab.'MECANIQUE'.'SENSIBILITE_Z_I'
                         tab.'COMPOSANTES' ;
 
* Logical for physical density dependent loading
LdFdY = EXIS resB 'SENSIBILITE_NORMALISEE_CHARGEMENT' ;
* If loading is physical density dependent
SI LdFdY ;
    NdFdYI = REDU resB.'SENSIBILITE_NORMALISEE_CHARGEMENT' mshI ;
FINS ;
 
* Vector of adjoint variables λ on I
* **********************************
 
* λ = - K_Z⁻¹ * ∂G/∂U
* Remark1: -λ is calculated on mshA, then restricted to I
* Remark2: minus sign is compensated later
SI (EGA tab.'SEUIL_ELEMENTS_ACTIFS' 0.0) ;
* Penalized rigidity matrix K on A
ZKA = (RIGI modMA ZmatMA) ET bdcMB ;
lamI = REDU (RESO ZKA dGdUI) mshI ;
OUBL ZKA ;
SINO ;
ERRE 'Actuellement non compatible avec SEUIL_ELEMENTS_ACTIFS <> 0' ;
* TODO à rendre compatible en cas de zone active non connexe.
*      En attendant on fait la résolution sur B
*      Appliquer un seuil sur l'aggrégation des contraintes locales? (c.f. Verbart Fig.14)
tmp1 = Wtab.'MECANIQUE'.'MODELE_B' ;
tmp2 = Wtab.'MECANIQUE'.'CARACTERISTIQUES_B' ;
tmp2 = TOPOMATE tmp1 tmp2 Wtab.'MECANIQUE'.'Z_B' tab.'COMPOSANTES' ;
ZKB = (RIGI tmp1 tmp2) ET bdcMB ;
lamI = REDU (RESO ZKB dGdUI) mshI ;
OUBL ZKB ;
FINS ;
 
* Rename components and switch to mass integration points
tmp = EXCO (EXTR lamI 'COMP') lamI (EXTR dGdUI 'COMP') ;
lamI1 = CHAN 'CHAM' tmp modMI 'MASSE' ;
lamI0 = CHAN 'ATTRIBUT' lamI 'NATURE' 'DISCRET' ;
 
* Sensitivities of the constraints term C with respect to physical densities Y on I
* *********************************************************************************
 
* As the constraints are KU - F = 0, the constraints term is C = λᵀ * (KU - F).
* Their sensitivities are ∂C/∂Y = λᵀ * (∂K/∂Y U - ∂F/∂Y)
 
* ∂K/∂Y on I
dKdYI = RIGI modMI dZmatMI ;
 
* ∂K/∂Y * U on I
**tmp = dKdYI * UI ;
* TODO version approchée en attendant d'avoir une version Fortran ?
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* -λ * ∂K/∂Y * U on I
SI FAUX ;
tmp = PROG ;
REPE itr (NBEL mshI) ;
    tmp1 = mshI ELEM &itr ;
    tmp2 = (REDU dKdYI tmp1) * UI ;
    tmp3 = REDU lamI0 tmp1 ;
    tmp4 = PSCA tmp2 tmp3 (EXTR tmp2 'COMP') (EXTR tmp3 'COMP') ;
    tmp4 = RESU tmp4 ;
    tmp = tmp ET (EXTR tmp4 'VALE' 'SCAL') ;
FIN itr ;
tmp = MANU 'CHML' modMI 'REPA' 'SCAL' tmp 'TYPE' 'SCALAIRE' 'GRAVITE' ;
SINO ;
* Création de maillages d'éléments dissociés à partir d'un maillage
* initial dont tous les éléments reposent sur une grille régulière
SI (NON (EXIS Wtab 'TOPOADJO_TEMPORAIRE')) ;
    MESS 'Pareparation du maillage temporaire pour TOPOADJO...' ;
    crt = tab.'ELIM_CRITERE' ;
    msh = Wtab.'MAILLAGE_B' ;
    SI (EGA (VALE 'MODE') 'TRID') ;
        env = ENVE msh ;
        edg = (ARET env) COUL 'ROUG' ;
        c1 = edg COOR 1 ;
        c2 = edg COOR 2 ;
        c3 = edg COOR 3 ;
        min1 = MINI c1 ;
        max1 = MAXI c1 ;
        min2 = MINI c2 ;
        max2 = MAXI c2 ;
        min3 = MINI c3 ;
        max3 = MAXI c3 ;
        el1 = msh ELEM 1 ;
        cel1 = el1 COOR 1 ;
        cel2 = el1 COOR 2 ;
        cel3 = el1 COOR 3 ;
        minel1 = MINI cel1 ;
        maxel1 = MAXI cel1 ;
        minel2 = MINI cel2 ;
        maxel2 = MAXI cel2 ;
        minel3 = MINI cel3 ;
        maxel3 = MAXI cel3 ;
        el1 = el1 PLUS ((min1 min2 min3) MOIN (minel1 minel2 minel3)) ;
        s1 = maxel1 - minel1 ;
        s2 = maxel2 - minel2 ;
        s3 = maxel3 - minel3 ;
        nb1 = ENTI 'PROC' ((max1 - min1) / s1) ;
        nb2 = ENTI 'PROC' ((max2 - min2) / s2) ;
        nb3 = ENTI 'PROC' ((max3 - min3) / s3) ;
        raw0 = el1 ;
        REPE itr (nb1 - 1) ;
            raw0 = raw0 ET (el1 PLUS ((&itr * s1) 0.0 0.0)) ;
        FIN itr ;
        ELIM raw0 crt ;
        raw1 = el1 ;
        REPE itr ((nb1 - 1) / 2) ;
            raw1 = raw1 ET (el1 PLUS ((s1 * (&itr * 2)) 0.0 0.0)) ;
        FIN itr ;
        ELIM raw0 raw1 crt ;
        raw2 = DIFF raw0 raw1 ;
        layer1 = raw1 ;
        layer2 = raw2 ;
        REPE itr ((nb2 - 1) / 2) ;
            layer1 = layer1 ET (raw1 PLUS (0.0 (s2 * (&itr * 2)) 0.0)) ;
            layer2 = layer2 ET (raw2 PLUS (0.0 (s2 * (&itr * 2)) 0.0)) ;
        FIN itr ;
        raw3 = raw1 PLUS (0.0 s2 0.0) ;
        raw4 = raw2 PLUS (0.0 s2 0.0) ;
        layer3 = raw3 ;
        layer4 = raw4 ;
        REPE itr ((nb2 - 2) / 2) ;
            layer3 = layer3 ET (raw3 PLUS (0.0 (s2 * (&itr * 2)) 0.0)) ;
            layer4 = layer4 ET (raw4 PLUS (0.0 (s2 * (&itr * 2)) 0.0)) ;
        FIN itr ;
        msh1 = layer1 ;
        msh2 = layer2 ;
        msh3 = layer3 ;
        msh4 = layer4 ;
        REPE itr ((nb3 - 1) / 2) ;
            msh1 = msh1 ET (layer1 PLUS (0.0 0.0 (s3 * (&itr * 2)))) ;
            msh2 = msh2 ET (layer2 PLUS (0.0 0.0 (s3 * (&itr * 2)))) ;
            msh3 = msh3 ET (layer3 PLUS (0.0 0.0 (s3 * (&itr * 2)))) ;
            msh4 = msh4 ET (layer4 PLUS (0.0 0.0 (s3 * (&itr * 2)))) ;
        FIN itr ;
        layer5 = layer1 PLUS (0.0 0.0 s3) ;
        layer6 = layer2 PLUS (0.0 0.0 s3) ;
        layer7 = layer3 PLUS (0.0 0.0 s3) ;
        layer8 = layer4 PLUS (0.0 0.0 s3) ;
        msh5 = layer5 ;
        msh6 = layer6 ;
        msh7 = layer7 ;
        msh8 = layer8 ;
        REPE itr ((nb3 - 2) / 2) ;
            msh5 = msh5 ET (layer5 PLUS (0.0 0.0 (s3 * (&itr * 2)))) ;
            msh6 = msh6 ET (layer6 PLUS (0.0 0.0 (s3 * (&itr * 2)))) ;
            msh7 = msh7 ET (layer7 PLUS (0.0 0.0 (s3 * (&itr * 2)))) ;
            msh8 = msh8 ET (layer8 PLUS (0.0 0.0 (s3 * (&itr * 2)))) ;
        FIN itr ;
        msh1to8 = msh1 ET msh2 ET msh3 ET msh4 ET msh5 ET msh6 ET msh7 et msh8 ;
        ELIM msh msh1to8 crt ;
        mshD = Wtab.'MAILLAGE_D' ;
        mshtab = TABL 'ESCLAVE' ;
        mshtab.(1) = mshD INTE msh1 ;
        mshtab.(2) = mshD INTE msh2 ;
        mshtab.(3) = mshD INTE msh3 ;
        mshtab.(4) = mshD INTE msh4 ;
        mshtab.(5) = mshD INTE msh5 ;
        mshtab.(6) = mshD INTE msh6 ;
        mshtab.(7) = mshD INTE msh7 ;
        mshtab.(8) = mshD INTE msh8 ;
    MESS 'Maillage temporaire genere pour TOPOADJO.' ;
    SINO ;
        cnt = CONT msh ;
*       cnt = cnt COUL 'ROUG' ;
        c1 = cnt COOR 1 ;
        c2 = cnt COOR 2 ;
        min1 = MINI c1 ;
        max1 = MAXI c1 ;
        min2 = MINI c2 ;
        max2 = MAXI c2 ;
        el1 = msh ELEM 1 ;
        cel1 = el1 COOR 1 ;
        cel2 = el1 COOR 2 ;
        minel1 = MINI cel1 ;
        maxel1 = MAXI cel1 ;
        minel2 = MINI cel2 ;
        maxel2 = MAXI cel2 ;
        el1 = el1 PLUS ((min1 min2) MOIN (minel1 minel2)) ;
        s1 = maxel1 - minel1 ;
        s2 = maxel2 - minel2 ;
        nb1 = ENTI 'PROC' ((max1 - min1) / s1) ;
        nb2 = ENTI 'PROC' ((max2 - min2) / s2) ;
        raw0 = el1 ;
        REPE itr (nb1 - 1) ;
            raw0 = raw0 ET (el1 PLUS ((&itr * s1) 0.0)) ;
        FIN itr ;
        ELIM raw0 crt ;
        raw1 = el1 ;
        REPE itr ((nb1 - 1) / 2) ;
            raw1 = raw1 ET (el1 PLUS ((s1 * (&itr * 2)) 0.0)) ;
        FIN itr ;
        ELIM raw0 raw1 crt ;
        raw2 = DIFF raw0 raw1 ;
        msh1 = raw1 ;
        msh2 = raw2 ;
        REPE itr ((nb2 - 1) / 2) ;
            msh1 = msh1 ET (raw1 PLUS (0.0 (s2 * (&itr * 2)))) ;
            msh2 = msh2 ET (raw2 PLUS (0.0 (s2 * (&itr * 2)))) ;
        FIN itr ;
        raw3 = raw1 PLUS (0.0 s2) ;
        msh3 = raw3 ;
        REPE itr ((nb2 - 2) / 2) ;
            msh3 = msh3 ET (raw3 PLUS (0.0 (s2 * (&itr * 2)))) ;
        FIN itr ;
        msh1234 = msh1 ET msh2 ET msh3 ;
        ELIM msh msh1234 crt ;
        msh4 = DIFF msh msh1234 ;
        msh1234 = msh1234 ET msh4 ;
        mshD = Wtab.'MAILLAGE_D' ;
        mshtab = TABL 'ESCLAVE' ;
        mshtab.(1) = mshD INTE msh1 ;
        mshtab.(2) = mshD INTE msh2 ;
        mshtab.(3) = mshD INTE msh3 ;
        mshtab.(4) = mshD INTE msh4 ;
    FINS ;
    Wtab.'TOPOADJO_TEMPORAIRE' = mshtab ;
*    REPE itr 4 ;
*        TRAC (mshtab.(&itr) ET cnt) ;
*    FIN itr ;
    SI (NEG (NBEL (DIFF mshD (ET mshtab))) 0) ;
        ERRE (CHAI 'Le maillage ne repose pas sur une grille reguliere.'
                   ' Pour le moment, seuls les maillages reguliers sont'
                   ' compatibles avec une contrainte limite de von Mises.') ;
    FINS ;
FINS ;
mshtab = Wtab.'TOPOADJO_TEMPORAIRE' ;
matMI = Wtab.'MECANIQUE'.'CARACTERISTIQUES_I' ;
VI = REDU VD mshI ;
nbr = 4 ;
SI (EGA (VALE 'MODE') 'TRID') ;
    nbr = 8 ;
FINS ;
REPE itr nbr ;
    tmp1 = mshI INTE mshtab.(&itr) ;
    tmp2 = (REDU dKdYI tmp1) * UI ;
    tmp3 = REDU lamI0 tmp1 ;
    tmp4 = PSCA tmp2 tmp3 (EXTR tmp2 'COMP') (EXTR tmp3 'COMP') ;
    tmp5 = REDU modMI tmp1 ;
    tmp6 = REDU matMI tmp5 ;
    tmp4 = CHAN 'CHAM' tmp4 tmp5 'NOEUD' 'SCALAIRE' ;
*   Je ne peux pas faire la résultante par élément, à la place
*   j'intègre, je divise par le volume puis je multiplie par le Nb de
*   noeuds. Ça diffère en axi car le volume que "représente" chaque noeud
*   est différent mais cette différence s'estompe loin de l'axe Z
*   et si les éléments sont radialement étroits.
*   Ça va aussi différer si les éléments sont quadratiques.
*   Bref c'est une rustine temporaire et qui aide à avoir des performances
*   raisonables.
    tmp4 = INTG 'ELEM' tmp4 tmp5 tmp6 ;
    tmp4 = tmp4 / (REDU VI tmp5) ;
    SI (EGA &itr 1) ;
       tmp = tmp4 ;
    SINO ;
       tmp = tmp Et tmp4 ;
    FINS ;
FIN itr ;
tmp = tmp * (NBNO (mshtab.(1) ELEM 1)) ;
FINS ;
*TRAC tmp modMI ;
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
* dC/dY on I
* Remark: minus sign on λ is compensated here
SI LdFdY ;
*   If the loading is physical density dependent
*   ∂C/∂Y = λᵀ * (∂K/∂Y * U - ∂F/∂Y)
    tmp1 = PSCA lamI1 NdFdYI (EXTR NdFdYI 'COMP') (EXTR NdFdYI 'COMP') ;
    dCdYI = (INTG 'ELEM' tmp1 modMI) - tmp ;
SINO ;
*   Otherwise
*   ∂C/∂Y = λᵀ * (∂K/∂Y * U)
    dCdYI = -1.0 * tmp ;
FINS ;
 
* Total derivative of G with respect to physical densities Y on D
* ***************************************************************
 
* Total derivative dG/dY
totdGdYD = dGdYD + (REDU dCdYI modMD) ;
 
* Volume-normalization
* ********************
 
NtotdGdYD = totdGdYD / VD ;
 
FINP NtotdGdYD ;
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales