Télécharger topoboot.procedur

Retour à la liste

* TOPOBOOT  PROCEDUR  FD218221  25/12/18    21:15:02     12429          
 
************************************************************************
** Procedure called by TOPOPTIM for preprocessing input data before
** starting the topology optimization loop.
**
** Author:
** Guenhael Le Quilliec (LaMe - Polytech Tours)
**
** Version:
** 3.0 2025/12/01 Various updates and improvements for TOPOPTIM V4.0.
** 2.0 2021/04/22 Made compatible with Cast3M-2020 for linear
**                thermo-mechanical resolution
** 1.0 2021/04/15 Initial version
************************************************************************
*
*   TOPOBOOT sections:
*   ------------------
*
*   - Initialization
*   - TOPOPTIM version < 4.0 input data verification
*   - General options
*   - Base (B) resolution model (mod) and material characteristics (mat)
*   - Mechanism (mcn)
*   - Base (B), frozen (F), design (D) and contact (cnt) meshes (msh)
*   - Initialization of active (A) and intersection (I) meshes
*   - Design (D) and frozen (F) model and material characteristics
*   - Constraints (g & G)
*   - Objective (O)
*   - Optimizer
*   - Design variables (X)
*   - Restrictions
*   - Filtering of some sensitivities or the intermediate density (S)
*   - Binarization by projection of the physical density (Y)
*   - Physical density penalization (Z)
*   - Convergence
*   - Parameters evolving over optimization cycles
*   - On-the-fly monitoring
*   - Output data
*   - Unitary (uni) and volume (V) element fields
*   - Default model, characteristics, unitary and volume element fields
*   - Design zone (D) various data: zero field, centroids & total volume
*   - Restart (R)
*   - Initial design variables (X) and physical densities (Y)
*
**********************************************************************
 
DEBP TOPOBOOT tab*'TABLE' ;
 
**********************************************************************
*                           Initialization                           *
**********************************************************************
 
* Start timer
TEMP 'ZERO' ;
 
* Working (W) table (tab)
tab.'WTABLE' = TABL ;
Wtab = tab.'WTABLE' ;
 
* Boolean/Logiocal (L) table (tab)
Wtab.'LOGIQUE' = TABL ;
Ltab = Wtab.'LOGIQUE' ;
 
* Define "Not a Number"
Wtab.'NAN' = LOG -1.0 ;
 
* Define empty variable
Wtab.'VIDE' = VIDE 'MMODEL' ;
 
* Table of pointers used by TOPOMODI
Wtab.'POINTEURS' = TABL ;
 
* Number of the last completed cycle
SI (NON (EXIS tab 'CYCLE')) ;
    tab.'CYCLE' = 0 ;
FINS ;
cyc = (ENTI (tab.'CYCLE')) + 1 ;
Wtab.'CYCLE' = cyc ;
 
* Number of the current iteration and its timer
Wtab.'ITERATION' = 0 ;
Wtab.'HORLOGE' = 0.0 ;
 
**********************************************************************
*           TOPOPTIM version < 4.0 input data verification           *
**********************************************************************
 
* Check old input names for compatibility
 
oldtab = TABL ;
oldtab.'CRITERE'            = CHAI 'CONVERGENCE_CRITERE' ;
oldtab.'CYCLES_SAUVEGARDES' = CHAI 'CYCLES_SAUVES' ;
oldtab.'FRACTION_VOLUME'    = CHAI 'FRACTION_VOLUME_LIMITE' ;
oldtab.'FACTEUR_D'          = CHAI 'OC_D' ;
oldtab.'FACTEUR_P'          = CHAI 'PENALISATION_P' ;
oldtab.'FACTEUR_Q'          = CHAI 'OC_Q' ;
oldtab.'SEUIL'              = CHAI 'SEUIL_ELEMENTS_ACTIFS' ;
oldtab.'TOPOLOGIE_MAX_INC'  = CHAI 'X_CHANGE_SEUIL' ;
oldtab.'TOPOLOGIE_MIN'      = CHAI 'X_MIN_SEUIL' ;
 
tmp = INDE oldtab ;
REPE itr (DIME tmp) ;
    SI (EXIS tab tmp.&itr) ;
        ERRE (CHAI 'Attention, depuis TOPOPTIM V4.0, le parametre' ' '
                   tmp.&itr ' a ete renomme' ' ' oldtab.(tmp.&itr) '.') ;
    FINS ;
FIN itr ;
 
* Check old name of hat filter
SI (EXIS tab 'FILTRE') ;
    SI (EGA tab.'FILTRE' 'MATRICE') ;
        ERRE (CHAI 'Attention, depuis TOPOPTIM V4.0, le filtre MATRICE'
                   ' a ete renomme CHAPEAU.') ;
    FINS ;
FINS ;
 
* Check old index of initial topology
SI (EXIS tab 'TOPOLOGIE') ;
    SI (EXIS tab.'TOPOLOGIE' 0) ;
        ERRE (CHAI 'Attention, depuis TOPOPTIM V4.0, la topologie'
                   ' initiale est donnee a l''indice 1.') ;
    FINS ;
FINS ;
 
**********************************************************************
*                          General options                           *
**********************************************************************
 
* Numerical precision
SI (NON (EXIS tab 'PRECISION')) ;
    tab.'PRECISION' = 1.0e-9 ;
FINS ;
 
* Division by zero protection
SI (NON (EXIS tab 'ZERO_DIVISION')) ;
    tab.'ZERO_DIVISION' = 0.001 ;
FINS ;
 
* Maximum distance for node fusion
SI (NON (EXIS tab 'ELIM_CRITERE')) ;
    tab.'ELIM_CRITERE' = 1.0e-9 ;
FINS ;
 
* User-defined procedure called just before the resolution, which
* can be used to update the loading and any other data
* if they (partly or totally) depend on topology
SI (NON (EXIS tab 'PROCEDURE_TOPOUPDT')) ;
    tab.'PROCEDURE_TOPOUPDT' = FAUX ;
FINS ;
 
* Personal procedure that can be used, among others, to
* calculate/modify the sensitivity field and objective value
SI (NON (EXIS tab 'PROCEDURE_TOPOPERS')) ;
    tab.'PROCEDURE_TOPOPERS' = FAUX ;
FINS ;
 
 
**********************************************************************
* Base (B) resolution model (mod) and material characteristics (mat) *
**********************************************************************
 
* Table of the resolution (res) data of the base (B) finite element
* model
Ltab.'PASAPAS' = EXIS tab 'RESOLUTION_PASAPAS' ;
SI Ltab.'PASAPAS' ;
    resB = tab.'RESOLUTION_PASAPAS' ;
SINO ;
    resB = tab.'RESOLUTION_LINEAIRE' ;
FINS ;
Wtab.'RESOLUTION_B' = resB ;
 
* Base (B) mechanical (M), thermal (T), and contact (cnt) models (mod)
modB = resB.'MODELE' ;
LM = EXIS modB 'FORM' 'MECANIQUE' ;
Ltab.'MECANIQUE' = LM ;
SI LM ;
    Wtab.'MECANIQUE' = TABL ;
    modMB = EXTR modB 'FORM' 'MECANIQUE' ;
    Wtab.'MECANIQUE'.'MODELE_B' = modMB ;
FINS ;
LT = EXIS modB 'FORM' 'THERMIQUE' ;
Ltab.'THERMIQUE' = LT ;
SI LT ;
    Wtab.'THERMIQUE' = TABL ;
    modTB = EXTR modB 'FORM' 'THERMIQUE' ;
    Wtab.'THERMIQUE'.'MODELE_B' = modTB ;
FINS ;
Lcnt = EXIS modB 'FORM' 'CONTACT' ;
Ltab.'CONTACT' = Lcnt ;
SI Lcnt ;
    Wtab.'CONTACT' = TABL ;
    modcntB = EXTR modB 'FORM' 'CONTACT' ;
    Wtab.'CONTACT'.'MODELE_B' = modcntB ;
FINS ;
 
* Base (B) mechanical (M), thermal (T), and contact (cnt) material
* characteristics (mat)
matB = resB.'CARACTERISTIQUES' ;
SI LM ;
    matMB = REDU matB modMB ;
    Wtab.'MECANIQUE'.'CARACTERISTIQUES_B' = matMB ;
FINS ;
SI LT ;
    matTB = REDU matB modTB ;
    Wtab.'THERMIQUE'.'CARACTERISTIQUES_B' = matTB ;
FINS ;
SI Lcnt ;
    matcntB = REDU matB modcntB ;
    Wtab.'CONTACT'.'CARACTERISTIQUES_B' = matcntB ;
FINS ;
 
 
**********************************************************************
*                          Mechanism (mcn)                           *
**********************************************************************
 
* Compliant mechanism synthesis topology optimization
Lmcn = EXIS tab 'MECANISME' ;
Ltab.'MECANISME' = Lmcn ;
 
* Solve the problem also without springs
SI (NON (EXIS tab 'MECANISME_ZERO_SPRING')) ;
    tab.'MECANISME_ZERO_SPRING' = FAUX ;
FINS ;
 
* Verification
* !!!!!!!!!!!!
 
* Check that a mechanical model exists
SI (Lmcn ET (NON LM)) ;
    MESS 'Pour effectuer une optimisation de mecanisme deformable,' ;
    MESS 'il faut imperativement preciser un modele mecanique.' ;
    ERRE 'Veuillez preciser le modele mecanique a utiliser.' ;
FINS ;
 
 
**********************************************************************
*  Base (B), frozen (F), design (D) and contact (cnt) meshes (msh)   *
**********************************************************************
 
* Base (B) mesh
* *************
* Remark: without contact mesh, only mechanical or thermal
SI LM ;
    mshB = EXTR modMB 'MAIL' ;
SINO ;
    mshB = EXTR modTB 'MAIL' ;
FINS ;
Wtab.'MAILLAGE_B' = mshB ;
 
* Frozen area (F) mesh
* ********************
SI (NON (EXIS tab 'ZONE_FIGEE')) ;
    tab.'ZONE_FIGEE' = VIDE 'MAILLAGE' ;
FINS ;
LF = (NBEL tab.'ZONE_FIGEE') > 0 ;
Ltab.'ZONE_FIGEE' = LF ;
SI LF ;
    mshF = tab.'ZONE_FIGEE' ;
    mshD = DIFF mshB mshF ;
SINO ;
    mshD = mshB ;
FINS ;
 
* Design (D) mesh
* ***************
SI LF ;
    mshD = DIFF mshB mshF ;
SINO ;
    mshD = mshB ;
FINS ;
 
* Contact (cnt) mesh
* ******************
SI Lcnt ;
    Wtab.'CONTACT'.'MAILLAGE_B' = EXTR modcntB 'MAIL' ;
FINS ;
 
* Verification
* !!!!!!!!!!!!
 
* Check that the number of elements of thermal and mechanical meshes
* equals that of the base mesh, otherwise the procedure can't work
SI (LM ET LT) ;
    SI (NEG (NBEL (EXTR modMB 'MAIL')) (NBEL (EXTR modTB 'MAIL'))) ;
        MESS 'Les modeles thermique et mecanique doivent s''appuyer sur un meme maillage.' ;
        ERRE 'Veuillez adapter vos modeles en consequence.' ;
    FINS ;
FINS ;
 
 
**********************************************************************
*      Initialization of active (A) and intersection (I) meshes      *
**********************************************************************
 
* Initialize active (A) & intersection (I) zones with empty meshes
* Remark: this is needed by TOPOACTI at first iteration
Wtab.'MAILLAGE_A' = DIFF (mshD ELEM 1) (mshD ELEM 1) ;
Wtab.'MAILLAGE_I' = Wtab.'MAILLAGE_A' ;
 
 
**********************************************************************
*    Design (D) and frozen (F) model and material characteristics    *
**********************************************************************
 
SI LF ;
    SI LM ;
        modMD = REDU modMB mshD ;
        modMF = REDU modMB mshF ;
        Wtab.'MECANIQUE'.'MODELE_F' = modMF ;
        matMD = REDU matMB mshD ;
        matMF = REDU matMB mshF ;
        Wtab.'MECANIQUE'.'CARACTERISTIQUES_F' = matMF ;
    FINS ;
    SI LT ;
        modTD = REDU modTB mshD ;
        modTF = REDU modTB mshF ;
        Wtab.'THERMIQUE'.'MODELE_F' = modTF ;
        matTD = REDU matTB mshD ;
        matTF = REDU matTB mshF ;
        Wtab.'THERMIQUE'.'CARACTERISTIQUES_F' = matTF ;
    FINS ;
SINO ;
    SI LM ;
        modMD = modMB ;
        matMD = matMB ;
    FINS ;
    SI LT ;
        modTD = modTB ;
        matTD = matTB ;
    FINS ;
FINS ;
Wtab.'MAILLAGE_D' = mshD ;
SI LM ;
    Wtab.'MECANIQUE'.'MODELE_D' = modMD ;
    Wtab.'MECANIQUE'.'CARACTERISTIQUES_D' = matMD ;
FINS ;
SI LT ;
    Wtab.'THERMIQUE'.'MODELE_D' = modTD ;
    Wtab.'THERMIQUE'.'CARACTERISTIQUES_D' = matTD ;
FINS ;
 
 
**********************************************************************
*                        Constraints (g & G)                         *
**********************************************************************
 
* List of constraints
Wtab.'LISTE_CONTRAINTES' = MOTS ;
 
* Maximum volume fraction constraint
* Remark: No longer optional since TOPOPTIM V4, as this data now
*         indicates not only the maximum volume but also whether a
*         maximum volume constraint is desired.
Ltab.'CONTRAINTE_VOLUME' = EXIS tab 'FRACTION_VOLUME_LIMITE' ;
SI Ltab.'CONTRAINTE_VOLUME' ;
    Wtab.'LISTE_CONTRAINTES' = Wtab.'LISTE_CONTRAINTES' ET (MOTS 'VOLUME') ;
FINS ;
 
* Maximum von Mises stress constraint
* TODO could be either a scalar or a field?
Ltab.'VERBART'  = FAUX ;
Ltab.'CONIGLIO' = FAUX ;
Ltab.'CONTRAINTE_SIGMA_VM' = EXIS tab 'SIGMA_VM_LIMITATION' ;
SI Ltab.'CONTRAINTE_SIGMA_VM' ;
    SI (NON (EXIS tab 'SIGMA_VM_LIMITE')) ;
        ERRE 'La contrainte de von Mises limite SIGMA_VM_LIMITE doit etre renseignee.' ;
    FINS ;
    SI (NON (EXIS (MOTS 'VERBART' 'CONIGLIO') tab.'SIGMA_VM_LIMITATION')) ;
        ERRE 'Seuls les approches VERBART et CONIGLIO sont disponibles.' ;
    FINS ;
    SI (EGA tab.'SIGMA_VM_LIMITATION' 'VERBART') ;
        Ltab.'VERBART'  = VRAI ;
    SINO ;
        Ltab.'CONIGLIO' = VRAI ;
    FINS ;
SINO ;
    SI (EXIS tab 'SIGMA_VM_LIMITE') ;
        Ltab.'CONTRAINTE_SIGMA_VM' = VRAI ;
        Ltab.'VERBART'             = VRAI ;
    FINS ;
FINS ;
SI Ltab.'CONTRAINTE_SIGMA_VM' ;
    Wtab.'LISTE_CONTRAINTES' = Wtab.'LISTE_CONTRAINTES' ET (MOTS 'SIGMA_VM_LIMITE') ;
FINS ;
 
* Aggregation
SI (EXIS tab 'AGREGATION') ;
    Wtab.'AGREGATION' = tab.'AGREGATION' ;
SINO ;
    Wtab.'AGREGATION' = MOT 'KSL' ;
*    SI (Ltab.'VERBART' OU Ltab.'CONIGLIO') ;
*        Wtab.'AGREGATION' = MOT 'KSL' ;
*    FINS ;
FINS ;
* And its parameter
SI (EXIS tab 'AGREGATION_P') ;
    Wtab.'AGREGATION_P' = tab.'AGREGATION_P' ;
SINO ;
    SI (EXIS Wtab 'AGREGATION') ;
        SI (EGA Wtab.'AGREGATION' 'KSL') ;
            Wtab.'AGREGATION_P' = 8.0 ;
        FINS ;
    FINS ;
FINS ;
 
* Verification
* !!!!!!!!!!!!
 
* Check that a constraint is set
SI (EGA (DIME Wtab.'LISTE_CONTRAINTES') 0);
    MESS 'Il manque une contrainte d''optimisation.'
    ERRE 'Par exemple, indiquez une valeur pour FRACTION_VOLUME_LIMITE.' ;
FINS ;
 
 
**********************************************************************
*                           Objective (O)                            *
**********************************************************************
 
* Type of objective function
* **************************
 
* Remark: mechanism synthesis uses compliance objective function type
SI (NON (EXIS tab 'OBJECTIF')) ;
    SI Ltab.'CONTRAINTE_VOLUME' ;
        tab.'OBJECTIF' = MOT 'COMPLIANCE' ;
    SINO ;
        tab.'OBJECTIF' = MOT 'FRACTION_VOLUME' ;
    FINS ;
FINS ;
 
Ltab.'OBJECTIF_COMPLIANCE' = EGA tab.'OBJECTIF' 'COMPLIANCE' ;
Ltab.'OBJECTIF_VOLUME'     = EGA tab.'OBJECTIF' 'FRACTION_VOLUME' ;
 
* Weights for multi-objective applications
* ****************************************
 
SI (NON (EXIS tab 'POIDS_MECANISME')) ;
    tab.'POIDS_MECANISME' = 1.0 ;
FINS ;
 
SI (NON (EXIS tab 'POIDS_ENERGIE_DEFO')) ;
    SI Lmcn ;
        tab.'POIDS_ENERGIE_DEFO' = 0.0 ;
    SINO ;
        tab.'POIDS_ENERGIE_DEFO' = 1.0 ;
    FINS ;
FINS ;
 
SI (NON (EXIS tab 'POIDS_TEMPERATURE')) ;
    SI Lmcn ;
        tab.'POIDS_TEMPERATURE' = 0.0 ;
    SINO ;
        tab.'POIDS_TEMPERATURE' = 1.0 ;
    FINS ;
FINS ;
 
* Verification
* !!!!!!!!!!!!
 
SI (NON (EXIS (MOTS 'COMPLIANCE' 'FRACTION_VOLUME') tab.'OBJECTIF')) ;
    ERRE (CHAI 'Il n''est pas possible de minimiser l''OBJECTIF :' ' ' tab.'OBJECTIF') ;
FINS ;
 
 
**********************************************************************
*                             Optimizer                              *
**********************************************************************
 
* Default optimization algorithm
* ******************************
 
SI (NON (EXIS tab 'OPTIMISEUR')) ;
    tab.'OPTIMISEUR' = MOT 'OC' ;
FINS ;
 
* Optimality criterion (OC)
* *************************
 
* Damping factor (D)
SI (NON (EXIS tab 'OC_D')) ;
*   Default value of D
    tab.'OC_D' = 0.5 ;
FINS ;
 
* Grey Scale Factor (Q)
SI (NON (EXIS tab 'OC_Q')) ;
*   Default value of Q
    tab.'OC_Q' = 1.0 ;
FINS ;
 
* Initial lower bound for the Lagrangian multiplier
SI (NON (EXIS tab 'OC_L1')) ;
    tab.'OC_L1' = 0.0 ;
FINS ;
 
* Initial upper bound for the Lagrangian multiplier
*SI (NON (EXIS tab 'OC_L2')) ;
*    tab.'OC_L2' = 1.0e5 ;
*FINS ;
 
* Minimal value for the upper bound
SI (NON (EXIS tab 'OC_L2_MIN')) ;
    tab.'OC_L2_MIN' = 1.0e-40 ;
FINS ;
 
* Convergence criterion for the bi-sectioning algorithm
SI (NON (EXIS tab 'OC_CRITERE')) ;
    tab.'OC_CRITERE' = 1.0e-8 ;
FINS ;
 
* Maximum number of iterations for the bi-sectioning algorithm
SI (NON (EXIS tab 'OC_MAX_IT')) ;
    tab.'OC_MAX_IT' = 100000 ;
FINS ;
 
* Lower bound for B (where B**D is the ratio between the new and the
* previous topology fields)
SI (NON (EXIS tab 'OC_B_MIN')) ;
    tab.'OC_B_MIN' = 1.0e-10 ;
FINS ;
 
* Method of Moving Asymptotes (MMA)
* *********************************
 
* Coefficient(s) C of MMA
* TODO à adapter le jour où on proposerra d'appliquer plus d'une contrainte
* TODO cette valeur par défaut est-elle judicieuse en cas de contrainte de
*      sur le volume?
SI (NON (EXIS tab 'MMA_C')) ;
    tab.'MMA_C' = PROG 1.0e4 ;
SINO ;
    SI (EGA (TYPE tab.'MMA_C') 'FLOTTANT') ;
        tab.'MMA_C' = PROG tab.'MMA_C' ;
    FINS ;
FINS ;
 
* Coefficient to be applied on the objective function so that it remains
* between 1 and 100
SI (NON (EXIS tab 'MMA_COEF_OBJECTIF')) ;
    tab.'MMA_COEF_OBJECTIF' = 1.0 ;
    SI (EGA tab.'OBJECTIF' 'COMPLIANCE') ;
*        tab.'MMA_COEF_OBJECTIF' = Wtab.'NAN' ;
        tab.'MMA_COEF_OBJECTIF' = 0.0 ;
    FINS ;
    SI (EGA tab.'OBJECTIF' 'FRACTION_VOLUME') ;
        tab.'MMA_COEF_OBJECTIF' = 100.0 ;
    FINS ;
FINS ;
 
* Verification
* !!!!!!!!!!!!
 
* With OC, only maximum volume fraction must be set
SI ((EGA tab.'OPTIMISEUR' 'OC') ET
     (((DIME Wtab.'LISTE_CONTRAINTES') > 1) OU
      (NON Ltab.'CONTRAINTE_VOLUME')));
    ERRE 'Seule une contrainte sur le volume est autorisee avec OC. Tentez d''utiliser l''optimiseur MMA.' ;
FINS ;
 
* Stress constraint is compatible only with MMA and a mechanical model should be defined
SI Ltab.'CONTRAINTE_SIGMA_VM' ;
    SI (NEG tab.'OPTIMISEUR' 'MMA') ;
        ERRE 'Seul l''optimiseur MMA est compatible avec une contrainte de von Mises maximale.' ;
    FINS ;
    SI (NON Ltab.'MECANIQUE') ;
        ERRE 'Un modele mecanique est indispensable pour imposer une contrainte de von Mises maximale.' ;
    FINS ;
FINS ;
 
 
**********************************************************************
*                        Design variables (X)                        *
**********************************************************************
 
* Move (absolute value of the increment) threshold of the design
* variable between two cycles
SI (NON (EXIS tab 'X_CHANGE_SEUIL')) ;
    tab.'X_CHANGE_SEUIL' = 0.2 ;
FINS ;
 
* Design variable minimum threshold
SI (NON (EXIS tab 'X_MIN_SEUIL')) ;
    SI Lmcn ;
        tab.'X_MIN_SEUIL' = 0.001 ;
    SINO ;
        tab.'X_MIN_SEUIL' = 0.0 ;
    FINS ;
FINS ;
 
 
**********************************************************************
*                            Restrictions                            *
**********************************************************************
 
* Procedure to add a geometric restrictions
Ltab.'RESTRICTIONS' = EXIS tab 'RESTRICTIONS' ;
 
 
**********************************************************************
*  Filtering of some sensitivities or the intermediate density (S)   *
**********************************************************************
 
* Type of filter
SI (NON (EXIS tab 'FILTRE')) ;
    tab.'FILTRE' = MOT 'GIBIANE' ;
FINS ;
 
* Parameters for Helmholtz (PDE) filter
SI (EGA tab.'FILTRE' 'EDP');
*   Filter radius
    SI (NON (EXIS tab 'FILTRE_EDP_RAYON')) ;
        ERRE 'Veuillez preciser un rayon de filtrage.' ;
    FINS ;
FINS ;
* Default Helmholtz (PDE) filter identification data
Wtab.'FILTRE_EDP_ID' = TABL ;
Wtab.'FILTRE_EDP_ID'.'MODEL'      = VIDE 'MMODEL' ;
Wtab.'FILTRE_EDP_ID'.'ZONE_FIGEE' = VIDE 'MAILLAGE' ;
Wtab.'FILTRE_EDP_ID'.'RAYON'      = 0.0 ;
 
* Parameters for hat filter
SI (EGA tab.'FILTRE' 'CHAPEAU');
*   Filter radius
    SI (NON (EXIS tab 'FILTRE_CHAPEAU_RAYON')) ;
        ERRE 'Veuillez preciser un rayon de filtrage.' ;
    FINS ;
*   Filter exponent
    SI (NON (EXIS tab 'FILTRE_CHAPEAU_EXPOSANT')) ;
        tab.'FILTRE_CHAPEAU_EXPOSANT' = 1.0 ;
    FINS ;
*   Filter weight minimum threshold
    SI (NON (EXIS tab 'FILTRE_CHAPEAU_CRITERE')) ;
        tab.'FILTRE_CHAPEAU_CRITERE' = 0.0 ;
    FINS ;
FINS ;
* Default hat filter identification data (compatible with SAUV/REST)
Wtab.'FILTRE_CHAPEAU_ID' = TABL ;
Wtab.'FILTRE_CHAPEAU_ID'.'MODEL'      = VIDE 'MMODEL' ;
Wtab.'FILTRE_CHAPEAU_ID'.'ZONE_FIGEE' = VIDE 'MAILLAGE' ;
Wtab.'FILTRE_CHAPEAU_ID'.'RAYON'      = 0.0 ;
Wtab.'FILTRE_CHAPEAU_ID'.'EXPOSANT'   = 0.0 ;
Wtab.'FILTRE_CHAPEAU_ID'.'CRITERE'    = 0.0 ;
 
* Parameters for Gibiane filter type
SI (EGA tab.'FILTRE' 'GIBIANE');
*   Filter rate
    SI (NON (EXIS tab 'FILTRE_GIBIANE_TAUX')) ;
        tab.'FILTRE_GIBIANE_TAUX' = 1 ;
    FINS ;
FINS ;
 
* Field to be filtered
Ltab.'FILTRE_SENSIBILITE_OBJECTIF' = FAUX ;
Ltab.'FILTRE_SENSIBILITES_CONTRAINTES' = FAUX ;
Ltab.'FILTRE_DENSITE' = FAUX ;
SI (EXIS tab 'FILTRER') ;
    SI (EXIS (MOTS 'DENSITE'
                   'SENSIBILITE_OBJECTIF'
                   'SENSIBILITES_CONTRAINTES') tab.'FILTRER') ;
        SI (EGA tab.'FILTRER' 'DENSITE') ;
            Ltab.'FILTRE_DENSITE' = VRAI ;
        FINS ;
        SI (EGA tab.'FILTRER' 'SENSIBILITE_OBJECTIF') ;
            SI (NON Ltab.'CONTRAINTE_VOLUME') ;
                ERRE (CHAI 'L''option FILTRER ne peut valoir SENSIBILITE_OBJECTIF'
                     ' qu''en cas de contrainte sur le volume.') ;
            FINS ;
            Ltab.'FILTRE_SENSIBILITE_OBJECTIF' = VRAI ;
        FINS ;
        SI (EGA tab.'FILTRER' 'SENSIBILITES_CONTRAINTES') ;
            SI (Ltab.'CONTRAINTE_VOLUME') ;
                ERRE (CHAI 'L''option FILTRER ne peut valoir SENSIBILITES_CONTRAINTES'
                     ' qu''en cas de contrainte sur von Mises.') ;
            FINS ;
            Ltab.'FILTRE_SENSIBILITES_CONTRAINTES' = VRAI ;
        FINS ;
    SINO ;
        ERRE 'Le champ a filtrer n''est pas reconnu.' ;
    FINS ;
SINO ;
    SI Ltab.'CONTRAINTE_VOLUME' ;
        tab.'FILTRER' = MOT 'SENSIBILITE_OBJECTIF' ;
        Ltab.'FILTRE_SENSIBILITE_OBJECTIF' = VRAI ;
    SINO ;
        tab.'FILTRER' = MOT 'SENSIBILITES_CONTRAINTES' ;
        Ltab.'FILTRE_SENSIBILITES_CONTRAINTES' = VRAI ;
    FINS ;
FINS ;
 
* Verification
* !!!!!!!!!!!!
 
*SI (Ltab.'FILTRE_DENSITE' ET (NEG tab.'FILTRE' 'CHAPEAU')) ;
**   If the density is filtered then the "modified filter" is needed
**   to get the X derivative. This is possible only with the hat
**   filter.
*    ERRE (CHAI 'Seul le filtre de type CHAPEAU permet de filtrer la'
*               ' DENSITE. Veuillez adapter les parametres d''entree.') ;
*FINS ;
 
 
**********************************************************************
*       Binarization by projection of the physical density (Y)       *
**********************************************************************
 
* Physical density projection with a hyperbolic sigmoid function
* **************************************************************
 
* TODO Créer une erreure ou une alerte si OC_Q est actif en même temps que la projection (car ils ont le même but)
*      Idem s'il est actif avec MMA?...
 
Ltab.'PROJECTION' = EXIS tab 'PROJECTION_RAIDEUR' ;
Ltab.'PROJECTION_VOLUME_CONSTANT' = FAUX ;
 
SI (NON (EXIS tab 'PROJECTION')) ;
    tab.'PROJECTION' = MOT 'GUEST_2004' ;
FINS ;
 
SI (EGA tab.'PROJECTION' 'GUEST_2004')  ;
   Wtab.'PROJECTION_SEUIL' = 0.0 ;
FINS ;
 
* Constant density threshold projection or constant volume projection
SI ((EGA tab.'PROJECTION' 'SIGMOID_TANH') ET
    (NON (EXIS tab 'PROJECTION_SEUIL'))) ;
    Wtab.'PROJECTION_SEUIL' = 0.5 ;
    Ltab.'PROJECTION_VOLUME_CONSTANT' = VRAI ;
FINS ;
SI (EXIS tab 'PROJECTION_SEUIL') ;
    Wtab.'PROJECTION_SEUIL' = tab.'PROJECTION_SEUIL' ;
FINS ;
* Maximum iterations of the volume-preserving bisection algorithm
SI (NON (EXIS tab 'PROJECTION_MAX_IT')) ;
    tab.'PROJECTION_MAX_IT' = 100 ;
FINS ;
* Bisection algorithm criterion on the volume-preserving threshold value
SI (NON (EXIS tab 'PROJECTION_SEUIL_CRITERE')) ;
    tab.'PROJECTION_SEUIL_CRITERE' = 1.0e-7 ;
FINS ;
* Bisection algorithm criterion on the volume-preserving difference rate
SI (NON (EXIS tab 'PROJECTION_VOLUME_CRITERE')) ;
    tab.'PROJECTION_VOLUME_CRITERE' = 1.0e-3 ;
FINS ;
 
* Verification
* !!!!!!!!!!!!
 
SI ((NEG tab.'PROJECTION' 'SIGMOID_TANH') ET
    (EXIS tab 'PROJECTION_SEUIL')) ;
    ERRE (CHAI 'La seuil de projection n est pas compatible'
               ' avec la fonction de projection'
               ' ' tab.'PROJECTION' '.') ;
FINS ;
 
**********************************************************************
*                 Physical density penalization (Z)                  *
**********************************************************************
 
* Penalization function
SI (NON (EXIS tab 'PENALISATION')) ;
    tab.'PENALISATION' = MOT 'SIMP' ;
FINS ;
 
* Penalty factor (P)
SI (NON (EXIS tab 'PENALISATION_P')) ;
    tab.'PENALISATION_P' = 3.0 ;
FINS ;
 
* List of material components that are proportional to the penalized
* physical densities (Z)
SI (NON (EXIS tab 'COMPOSANTES')) ;
    tab.'COMPOSANTES' = MOTS 'YOUN' 'K' 'RHO' ;
FINS ;
 
* Ratio between the mechanical stiffness of the voids and the material
SI (NON (EXIS tab 'RAPPORT_RAIDEURS_MECANIQUES')) ;
    SI (Lmcn OU LT) ;
        tab.'RAPPORT_RAIDEURS_MECANIQUES' = 1.0e-8 ;
    SINO ;
        tab.'RAPPORT_RAIDEURS_MECANIQUES' = 0.0 ;
    FINS ;
FINS ;
 
* Ratio between the thermal stiffness of the voids and the material
SI (NON (EXIS tab 'RAPPORT_RAIDEURS_THERMIQUES')) ;
    tab.'RAPPORT_RAIDEURS_THERMIQUES' = 0.001 ;
FINS ;
 
* Minimal threshold value of the penalized physical densities (Z) to
* save mesh elements
SI (NON (EXIS tab 'SEUIL_ELEMENTS_ACTIFS')) ;
    tab.'SEUIL_ELEMENTS_ACTIFS' = 1.0e-9 ;
FINS ;
 
 
**********************************************************************
*                            Convergence                             *
**********************************************************************
 
* Maximum number of optimization cycles
SI (NON (EXIS tab 'MAX_CYCLES')) ;
    tab.'MAX_CYCLES' = 100 ;
FINS ;
 
* Convergence criterion
SI (NON (EXIS tab 'CONVERGENCE_CRITERE')) ;
    tab.'CONVERGENCE_CRITERE' = 0.01 ;
FINS ;
 
* Initialize or reset convergence logical
tab.'CONVERGENCE' = FAUX ;
 
* Stop optimization trigger
Ltab.'STOP' = FAUX ;
 
* Verification
* !!!!!!!!!!!!
 
* Make sure the optimization loop will start
SI ((Wtab.'CYCLE') > tab.'MAX_CYCLES') ;
    SAUT 1 'LIGN' ;
    MESS ' /!\ ALERTE | Reprise de TOPOPTIM impossible.                 |' ;
    MESS '     ====== | Le nombre maximum de cycles a deja ete atteint. |' ;
    SAUT 1 'LIGN' ;
FINS ;
 
 
**********************************************************************
*            Parameters evolving over optimization cycles            *
**********************************************************************
 
* Input parameters that may be evolutions over optimization cycles
* (otherwise they are scalar values)
Wtab.'EVOLUTIONS' = MOTS 'PENALISATION_P' 'OC_D' 'OC_Q' 'PROJECTION_RAIDEUR' ;
 
* Initialization call to TOPOEVAL
* Remark: this can be needed for the call to TOPOPHYS at the end of
*         TOPOBOOT.
TOPOEVOL tab ;
 
**********************************************************************
*                       On-the-fly monitoring                        *
**********************************************************************
 
* Plot results on-the-fly
SI (NON (EXIS tab 'TRAC')) ;
    tab.'TRAC' = VRAI ;
FINS ;
 
* Table of data for visualization
Wtab.'VISU' = TABL ;
 
* User list of variable names to be printed on the fly
SI (EXIS tab 'MESS') ;
    SI (NEG (TYPE tab.'MESS') 'LISTMOTS') ;
        tab.'MESS' = MOTS tab.'MESS' ;
    FINS ;
SINO ;
    tab.'MESS' = MOTS ;
FINS ;
 
* Full list of variable names to be printed on the fly
* Remark: order will be preserved when printed
tmp = MOTS ;
tmp = tmp ET (MOTS 'OBJECTIF') ;
SI Ltab.'CONTRAINTE_SIGMA_VM' ;
    tmp = tmp ET (MOTS 'SIGMA_VM_MAX') ;
FINS ;
SI (Ltab.'CONTRAINTE_VOLUME' ET (NEG tab.'OPTIMISEUR' 'OC')) ;
*   Not shown in the case of volume minimization as the objective
*   corresponds to volume fraction and is already displayed
*   OU Ltab.'OBJECTIF_VOLUME'
    tmp = tmp ET (MOTS 'FRACTION_VOLUME') ;
FINS ;
tmp = tmp ET (MOTS 'MAX_X_CHANGE') ;
*   Add any requested variables
tmp = tmp ET tab.'MESS' ;
*   Time / cycle
tmp = tmp ET (MOTS 'TEMPS_CYCLE') ;
*   Number of active elements
SI  (NON ((LM ET (tab.'RAPPORT_RAIDEURS_MECANIQUES' > 
                  tab.'SEUIL_ELEMENTS_ACTIFS')) OU
          (LT ET (tab.'RAPPORT_RAIDEURS_THERMIQUES' > 
                  tab.'SEUIL_ELEMENTS_ACTIFS')))) ;
    tmp = tmp ET (MOTS 'NB_ACTIVE_ELEM') ;
FINS ;
Wtab.'VISU'.'VARIABLES' = tmp ;
 
 
**********************************************************************
*                            Output data                             *
**********************************************************************
 
* Logical to save data on a file
SI (NON (EXIS tab 'SAUV')) ;
    tab.'SAUV' = FAUX ;
FINS ;
 
* Save data that might be needed in case of restart
SI (NON (EXIS tab 'SAUVEGARDE_DE_REPRISE')) ;
    tab.'SAUVEGARDE_DE_REPRISE' = VRAI ;
FINS ;
 
* Cycles to be saved
SI (NON (EXIS tab 'CYCLES_SAUVES')) ;
    tab.'CYCLES_SAUVES' = LECT ;
FINS ;
 
* Data to be saved
SI (NON (EXIS tab 'MES_SAUVEGARDES')) ;
    tab.'MES_SAUVEGARDES' = TABL ;
FINS ;
SI (NON (EXIS tab.'MES_SAUVEGARDES' 'X')) ;
    tab.'MES_SAUVEGARDES'.'X' = VRAI ;
FINS ;
SI (NON (EXIS tab.'MES_SAUVEGARDES' 'X_FILTRE')) ;
    tab.'MES_SAUVEGARDES'.'X_FILTRE' = VRAI ;
FINS ;
SI (NON (EXIS tab.'MES_SAUVEGARDES' 'TOPOLOGIE')) ;
    tab.'MES_SAUVEGARDES'.'TOPOLOGIE' = VRAI ;
FINS ;
SI (NON (EXIS tab.'MES_SAUVEGARDES' 'MAILLAGE')) ;
    tab.'MES_SAUVEGARDES'.'MAILLAGE' = VRAI ;
FINS ;
SI (NON (EXIS tab.'MES_SAUVEGARDES' 'RESOLUTION')) ;
    tab.'MES_SAUVEGARDES'.'RESOLUTION' = FAUX ;
FINS ;
 
* Only for step-by-step resolution
SI Ltab.'PASAPAS' ;
*   Save time list
    SI (NON (EXIS resB 'TEMPS_SAUVES')) ;
        tmp = EXTR (LECT 1 (DIME resB.'TEMPS_CALCULES')) resB.'TEMPS_CALCULES' ;
    SINO ;
        tmp = resB.'TEMPS_SAUVES' ;
    FINS ;
    Wtab.'TEMPS_SAUVES' = tmp ;
FINS ;
 
* List of names of output evolutions
SI (NON (EXIS tab 'EVOLUTIONS_SAUVEES')) ;
*   Copy 'VISU'.'VARIABLES' without 'EVOLUTIONS'
*    tmp = MOTS ;
*    REPE itr (DIME Wtab.'VISU'.'VARIABLES') ;
*        tmp1 = EXTR Wtab.'VISU'.'VARIABLES' &itr ;
*        SI (NON (EXIS Wtab.'EVOLUTIONS' tmp1)) ;
*            tmp = tmp ET (MOTS tmp1) ;
*        FINS ;
*    FIN itr ;
    tmp = Wtab.'VISU'.'VARIABLES' ;
*   Save them by default
    tab.'EVOLUTIONS_SAUVEES' = tmp ;
SINO ;
    SI (EGA (TYPE tab.'EVOLUTIONS_SAUVEES') 'MOT') ;
        tab.'EVOLUTIONS_SAUVEES' = MOTS tab.'EVOLUTIONS_SAUVEES' ;
    FINS ;
FINS ;
 
* Table of output evolutions
SI (NON (EXIS tab 'EVOLUTIONS')) ;
    tab.'EVOLUTIONS' = TABL ;
FINS ;
 
* Table of evolution default colors
tmp = TABL ;
tmp.'TEMPS_CYCLE'         = MOT 'VERT' ;
tmp.'OBJECTIF'            = MOT 'ROUG' ;
tmp.'FRACTION_VOLUME'     = MOT 'JAUN' ;
tmp.'SIGMA_VM_MAX'        = MOT 'CORA' ;
tmp.'MAX_X_CHANGE'        = MOT 'ORAN' ;
tmp.'NB_ACTIVE_ELEM'      = MOT 'AZUR' ;
tmp.'PENALISATION_P'      = MOT 'ROSE' ;
tmp.'OC_Q'                = MOT 'GRIS' ;
tmp.'OC_D'                = MOT 'CYAN' ;
tmp.'PROJECTION_RAIDEUR'  = MOT 'BLAN' ;
Wtab.'EVOLUTIONS_COULEUR' = tmp ;
 
* Design variable output table
SI (NON (EXIS tab 'X')) ;
    tab.'X' = TABL ;
FINS ;
 
* Intermediate densities output table
SI (NON (EXIS tab 'X_FILTRE')) ;
    tab.'X_FILTRE' = TABL ;
FINS ;
 
* Topology output table
SI (NON (EXIS tab 'TOPOLOGIE')) ;
    tab.'TOPOLOGIE' = TABL ;
FINS ;
 
* Mesh output table
SI (NON (EXIS tab 'MAILLAGE')) ;
    tab.'MAILLAGE' = TABL ;
FINS ;
 
* Resolution output table
SI (NON (EXIS tab 'RESOLUTION')) ;
    tab.'RESOLUTION' = TABL ;
FINS ;
 
 
**********************************************************************
*            Unitary (uni) and volume (V) element fields             *
**********************************************************************
 
* For the mechanical (M) model
SI LM ;
*   Unitary field
    uniMD = MANU 'CHML' modMD 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
*   Volume of each element
    VMD = INTG 'ELEM' uniMD modMD matMB ;
*   Change the unitary field so that it is not considered as constant
*   to avoid problems with the function POIN 'SUPE'
    uniMD = VMD / VMD ;
    Wtab.'MECANIQUE'.'UNIT_D' = uniMD ;
    Wtab.'MECANIQUE'.'VOLUME_D' = VMD ;
    SI LF ;
*       Unitary field
        uniMF = MANU 'CHML' modMF 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
*       Volume of each element
        VMF = INTG 'ELEM' uniMF modMF matMB ;
        Wtab.'MECANIQUE'.'UNIT_F' = uniMF ;
        Wtab.'MECANIQUE'.'VOLUME_F' = VMF ;
    FINS ;
FINS ;
 
* For the thermal (T) model
SI LT ;
*   Unitary field
    uniTD = MANU 'CHML' modTD 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
*   In thermal mode, 'GRAVITE' option is not respected
*   This issue is resolved by using CHAN 'GRAVITE'
    uniTD = CHAN 'GRAVITE' uniTD modTD ;
*   Volume of each element
    VTD = INTG 'ELEM' uniTD modTD matTB ;
*   In thermal mode, INTG 'ELEM' is not given at 'GRAVITE' even if
*   input fields were defined at 'GRAVITE'
*   This problem is solved by using CHAN 'GRAVITE'
    VTD = CHAN 'GRAVITE' VTD modTD ;
*   Divide the element volume field by itself to get uniTD so that
*   it is not considered as constant and avoid problems with the
*   function POIN 'SUPE'
    uniTD = VTD / VTD ;
    Wtab.'THERMIQUE'.'UNIT_D' = uniTD ;
    Wtab.'THERMIQUE'.'VOLUME_D' = VTD ;
    SI LF ;
*       Unitary field
        uniTF = MANU 'CHML' modTF 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
*       In thermal mode, 'GRAVITE' option is not respected
*       This issue is resolved by using CHAN 'GRAVITE'
        uniTF = CHAN 'GRAVITE' uniTF modTF ;
*       Volume of each element
        VTF = INTG 'ELEM' uniTF modTF matTB ;
*       In thermal mode, INTG 'ELEM' is not given at 'GRAVITE' even if
*       input fields were defined at 'GRAVITE'
*       This problem is solved by using CHAN 'GRAVITE'
        VTF = CHAN 'GRAVITE' VTF modTF ;
*       Divide the element volume field by itself to get uniTF so that
*       it is not considered as constant and avoid problems with the
*       function POIN 'SUPE'
        uniTF = VTF / VTF ;
        Wtab.'THERMIQUE'.'UNIT_F' = uniTF ;
        Wtab.'THERMIQUE'.'VOLUME_F' = VTF ;
    FINS ;
FINS ;
 
**********************************************************************
* Default model, characteristics, unitary and volume element fields  *
**********************************************************************
 
SI LM ;
*   If a mechanical (M) model exists
    Wtab.'MODELE_B' = modMB ;
    Wtab.'MODELE_D' = modMD ;
    Wtab.'CARACTERISTIQUES_B' = matMB ;
    Wtab.'CARACTERISTIQUES_D' = matMD ;
    SI LF ;
        Wtab.'MODELE_F' = modMF ;
        Wtab.'CARACTERISTIQUES_F' = matMF ;
    FINS ;
    uniD = uniMD ;
    SI LF ;
        uniF = uniMF ;
    FINS ;
    Wtab.'VOLUME_D' = VMD ;
    Wtab.'VOLUME_F' = VMF ;
SINO ;
*   If only a thermal (T) model exists
    Wtab.'MODELE_B' = modTB ;
    Wtab.'MODELE_D' = modTD ;
    Wtab.'CARACTERISTIQUES_B' = matTB ;
    Wtab.'CARACTERISTIQUES_D' = matTD ;
    SI LF ;
        Wtab.'MODELE_F' = modTF ;
        Wtab.'CARACTERISTIQUES_F' = matTF ;
    FINS ;
    uniD = uniTD ;
    SI LF ;
        uniF = uniTF ;
    FINS ;
    Wtab.'VOLUME_D' = VTD ;
    Wtab.'VOLUME_F' = VTF ;
FINS ;
Wtab.'UNIT_D' = uniD ;
SI LF ;
    Wtab.'UNIT_F' = uniF ;
FINS ;
 
 
**********************************************************************
* Design zone (D) various data: zero field, centroids & total volume *
**********************************************************************
 
* Zero element field
Wtab.'ZERO_D' = uniD * 0.0 ;
 
* Node mesh of the centroids
Wtab.'MAILLAGE_GRAVITE_D' = uniD POIN 'SUPE' -1.0 ;
 
SI LF ;
    SI LM ;
        VtotD = INTG uniD modMD matMD ;
*        VtotF = INTG uniF modMF matMF ;
    SINO ;
        VtotD = INTG uniD modTD matTD ;
*        VtotF = INTG uniF modTF matTF ;
    FINS ;
*    Wtab.'VOLUME_TOTAL_B' = VtotD + VtotF ;
SINO ;
    SI LM ;
        VtotD = INTG uniD modMD matMD ;
    SINO ;
        VtotD = INTG uniD modTD matTD ;
    FINS ;
*    Wtab.'VOLUME_TOTAL_B' = VtotD ;
FINS ;
Wtab.'VOLUME_TOTAL_D' = VtotD ;
 
 
**********************************************************************
*                            Restart (R)                             *
**********************************************************************
 
* Remark: The following inputs are not supposed to be provided by the
*         user. They are supposed to be generated by TOPOPTIM during a
*         previous call.
 
* List of restarts
* ****************
 
* Initialize or update the list of restarts
SI (EGA Wtab.'CYCLE' 1) ;
    tab.'REPRISES' = LECT ;
SINO ;
*   It should exist, but this prevents any error if it's missing
    SI (NON (EXIS tab 'REPRISES')) ;
        tab.'REPRISES' = LECT ;
    FINS ;
    tab.'REPRISES' = tab.'REPRISES' ET cyc ;
FINS ;
 
* Restart data table
* ******************
 
SI (EXIS tab 'RTABLE') ;
 
    Rtab = tab.'RTABLE' ;
 
*   Initialize TOPOMODI for loading and mechanism
    TOPOMODI tab (CHAI 'CHARGEMENT') Rtab.'CHARGEMENT' ;
    SI Ltab.'MECANISME' ;
        TOPOMODI tab (CHAI 'MECANISME') Rtab.'MECANISME' ;
    FINS ;
 
*   If any hat filter matrix exists
    SI (EXIS Rtab 'FILTRE_CHAPEAU_NORMALISE') ;
*       We will try to use it
        Wtab.'FILTRE_CHAPEAU_NORMALISE' = Rtab.'FILTRE_CHAPEAU_NORMALISE' ;
*       If a modified matrix exists too
        SI (EXIS Rtab 'FILTRE_CHAPEAU_MODIFIE') ;
*           We will try to use it
            Wtab.'FILTRE_CHAPEAU_MODIFIE' = Rtab.'FILTRE_CHAPEAU_MODIFIE' ;
        FINS ;
*       Get their corresponding identification data
        Wtab.'FILTRE_CHAPEAU_ID' = Rtab.'FILTRE_CHAPEAU_ID' ;
 
*       Finally, as Wtab.'MAILLAGE_GRAVITE_D' was recreated, the saved
*       hat filter matrix is not compatible. We need to fuse the
*       matrix mesh nodes with those of Wtab.'MAILLAGE_GRAVITE_D'
        tmp1 = Wtab.'MAILLAGE_GRAVITE_D' ;
        tmp2 = CHAN POI1 (EXTR Wtab.'FILTRE_CHAPEAU_NORMALISE' 'MAIL') ;
        ELIM tmp1 tmp2 tab.'ELIM_CRITERE' ;
        SI ((NBEL (DIFF tmp1 tmp2)) > 0) ;
*           This situation should never happen. Nevertheless, if it
*           does, then we will force the regeneration of the hat filter
*           matrix by resetting one of its identification data
            Wtab.'FILTRE_CHAPEAU_ID'.'MODEL' = VIDE 'MMODEL' ;
        FINS ;
    FINS ;
 
*   Reload the previous MMA table
    SI (EXIS Rtab 'MMA') ;
        Wtab.'MMA' = Rtab.'MMA' ;
    FINS ;
 
*   Use the last density threshold projection
*   if it has been saved in Rtab and if constant volume projection
    SI ((EXIS Rtab 'PROJECTION_SEUIL') ET Ltab.'PROJECTION_VOLUME_CONSTANT') ;
        Wtab.'PROJECTION_SEUIL' = Rtab.'PROJECTION_SEUIL' ;
    FINS ;
 
FINS ;
 
 
**********************************************************************
*          Initial design variables and physical densities           *
**********************************************************************
 
* If floats are given as initial design variables (X)
* or initial physical densities (Y)
* transform them into a homogeneous fields on B
* ------------------------------------------------------------------
 
SI (EXIS tab.'X' cyc) ;
    SI (EGA (TYPE tab.'X'.cyc) 'FLOTTANT') ;
        tab.'X'.cyc = uniD * tab.'X'.cyc ;
        SI LF ;
*           Combine the field on D with the unitary field on F
*           This creates a 2-zone field merged into a single zone
*           field, preserving the original element order using REDU
            tab.'X'.cyc = REDU (tab.'X'.cyc ET uniF) Wtab.'MODELE_B' ;
        FINS ;
    FINS ;
FINS ;
 
SI (EXIS tab.'X_FILTRE' cyc) ;
    SI (EGA (TYPE tab.'X_FILTRE'.cyc) 'FLOTTANT') ;
        tab.'X_FILTRE'.cyc = uniD * tab.'X_FILTRE'.cyc ;
        SI LF ;
*           Combine the field on D with the unitary field on F
*           This creates a 2-zone field merged into a single zone
*           field, preserving the original element order using REDU
            tab.'X_FILTRE'.cyc = REDU (tab.'X_FILTRE'.cyc ET uniF)
                                           Wtab.'MODELE_B' ;
        FINS ;
    FINS ;
FINS ;
 
SI (EXIS tab.'TOPOLOGIE' cyc) ;
    SI (EGA (TYPE tab.'TOPOLOGIE'.cyc) 'FLOTTANT') ;
        tab.'TOPOLOGIE'.cyc = uniD * tab.'TOPOLOGIE'.cyc ;
        SI LF ;
*           Combine the field on D with the unitary field on F
*           This creates a 2-zone field merged into a single zone
*           field, preserving the original element order using REDU
            tab.'TOPOLOGIE'.cyc = REDU (tab.'TOPOLOGIE'.cyc ET uniF)
                                             Wtab.'MODELE_B' ;
        FINS ;
    FINS ;
FINS ;
 
* If there are initial design variables (X)
* and/or initial intermediate densities (S)
* and/or initial physical densities (Y)
* and/or a restart table (that always contains X, S and Y)
* then we will use those data at best to initilize X, S and Y
* -----------------------------------------------------------
 
* Initialization
LX    = EXIS tab.'X' cyc ;
LS    = EXIS tab.'X_FILTRE' cyc ;
LY    = EXIS tab.'TOPOLOGIE' cyc ;
LRtab = EXIS tab 'RTABLE' ;
 
* If some initial values are provided
SI (LX OU LS OU LY OU LRtab) ;
 
*   Define at least X
*   -----------------
 
*   By default we will use X, S and Y from user defined input
*   otherwise from Rtab
    SI LX ;
        XB = tab.'X'.cyc ;
    FINS ;
    SI LS ;
        SB = tab.'X_FILTRE'.cyc ;
    FINS ;
    SI LY ;
        YB = tab.'TOPOLOGIE'.cyc ;
    FINS ;
    SI LRtab ;
        RXB = tab.'RTABLE'.'X_B' ;
        RSB = tab.'RTABLE'.'S_B' ;
        RYB = tab.'RTABLE'.'Y_B' ;
    FINS ;
 
*   If the only existing X, S and Y are from Rtab
    SI (LRtab ET (NON LX) ET (NON LS) ET (NON LY))  ;
*       We will use them together
        XB = RXB ;
        LX = VRAI ;
        SB = RSB ;
        LS = VRAI ;
        YB = RYB ;
        LY = VRAI ;
    FINS ;
 
*   If X is given with Rtab
    SI (LRtab ET LX ET (NON LS) ET (NON LY))  ;
*       If it is the same as the one in Rtab
        SI (EGA XB RXB) ;
*           Then we can use S and Y from Rtab
            SB = RSB ;
            LS = VRAI ;
            YB = RYB ;
            LY = VRAI ;
        FINS ;
    FINS ;
 
*   Verification
*   !!!!!!!!!!!!
 
*   If X does not exist at this step
*   Remark: It means Rtab didn't exist either, otherwise X would have
*           been taken from Rtab.
    SI (NON LX) ;
        tmp = CHAI 'Vous avez fourni des valeurs imposees' ;
        SI LS ;
            tmp = CHAI tmp ' des densites intermediaires' ;
        FINS ;
        SI (LS ET LY) ;
            tmp = CHAI tmp ' et' ;
        FINS ;
        SI LY ;
            tmp = CHAI tmp ' de la topologie' ;
        FINS ;
        tmp = CHAI tmp ' pour le prochain cycle sans avoir fourni'
                 ' les variables de conceptions correspondantes' ;
        SI (NON LS) ;
            tmp = CHAI tmp ' ni les densites intermediaires,' ;
        FINS ;
        SI (NON LY) ;
            tmp = CHAI tmp ' ni la topologie,' ;
        FINS ;
        ERRE (CHAI tmp ' qui sont indispensables dans ce cas.') ;
    FINS ;
 
*   If Y is missing but not S
    SI ((LS ET (NON LY)) OU ((NON LS) ET LY));
        tmp = CHAI 'Vous avez fourni des valeurs imposees des'
                   ' variables de conceptions et' ;
        SI LS ;
            tmp = CHAI tmp ' des densites intermediaires' ;
        FINS ;
        SI LY ;
            tmp = CHAI tmp ' de la topologie' ;
        FINS ;
        tmp = CHAI tmp ' mais sans les valeurs imposees' ;
        SI (NON LS) ;
            tmp = CHAI tmp ' des densites intermediaires,' ;
        FINS ;
        SI (NON LY) ;
            tmp = CHAI tmp ' de la topologie,' ;
        FINS ;
        ERRE (CHAI tmp ' qui sont alors indispensables. Ou bien vous'
                       ' pouvez ne fournir que les variables de'
                       ' conceptions seules.') ;
    FINS ;
 
 
*   Projection of X, S and Y on D
*   -----------------------------
 
*   This projection is used for remeshing compatibility but
*   also in case the active zone was updated or the model etc...
*   but also because X, S and Y are given on B while we
*   need them only on D.
 
*   Logical for identical S and X
    LSegaX = FAUX ;
    SI (LS ET LX) ;
*       If they are identical
        LSegaX = EGA SB XB ;
    FINS ;
 
*   Logical for identical Y and S
    LYegaS = FAUX ;
    SI (LY ET LS) ;
*       If they are identical
        LYegaS = EGA YB SB ;
    FINS ;
 
*   Logical for identical Y and X
    LYegaX = FAUX ;
    SI (LY ET LX) ;
*       If they are identical
        LYegaX = EGA YB XB ;
    FINS ;
 
    SI LM ;
*       Project the input topology element field at the centroids of
*       the elements of the model modMD
        SI LX ;
            tmp = PROI Wtab.'MAILLAGE_GRAVITE_D' XB ;
            XD = MANU 'CHML' modMD 'REPA' 'SCAL' (EXTR tmp 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
        FINS ;
        SI (LS ET (NON LSegaX)) ;
            tmp = PROI Wtab.'MAILLAGE_GRAVITE_D' SB ;
            SD = MANU 'CHML' modMD 'REPA' 'SCAL' (EXTR tmp 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
        FINS ;
        SI (LY ET (NON LYegaS) ET (NON LYegaX)) ;
            tmp = PROI Wtab.'MAILLAGE_GRAVITE_D' YB ;
            YD = MANU 'CHML' modMD 'REPA' 'SCAL' (EXTR tmp 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
        FINS ;
    SINO ;
*       Project the input topology element field at the centroids of
*       the elements of the model modTD
        SI LX ;
            tmp = PROI Wtab.'MAILLAGE_GRAVITE_D' XB ;
            XD = MANU 'CHML' modTD 'REPA' 'SCAL' (EXTR tmp 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
*           In thermal mode, 'GRAVITE' option is not respected
*           This issue is resolved using CHAN 'GRAVITE'
            XD = CHAN 'GRAVITE' XD modTD ;
        FINS ;
        SI (LS ET (NON LSegaX)) ;
            tmp = PROI Wtab.'MAILLAGE_GRAVITE_D' SB ;
            SD = MANU 'CHML' modTD 'REPA' 'SCAL' (EXTR tmp 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
*           In thermal mode, 'GRAVITE' option is not respected
*           This issue is resolved using CHAN 'GRAVITE'
            SD = CHAN 'GRAVITE' SD modTD ;
        FINS ;
        SI (LY ET (NON LYegaS) ET (NON LYegaX)) ;
            tmp = PROI Wtab.'MAILLAGE_GRAVITE_D' YB ;
            YD = MANU 'CHML' modTD 'REPA' 'SCAL' (EXTR tmp 'VALE') 'TYPE' 'SCALAIRE' 'GRAVITE' ;
*           In thermal mode, 'GRAVITE' option is not respected
*           This issue is resolved using CHAN 'GRAVITE'
            YD = CHAN 'GRAVITE' YD modTD ;
        FINS ;
    FINS ;
 
    SI LSegaX ;
        SD = XD ;
    FINS ;
    SI LYegaS ;
        YD = SD ;
    FINS ;
    SI LYegaX ;
        YD = XD ;
    FINS ;
 
*   Add requested restrictions on X in case they are not already respected
*   ----------------------------------------------------------------------
 
*   It must be done on initial X as further calls to TOPORSTR during
*   optimisation loop will affect only the increments of X, not X itself.
*   Remark: Y does not respect the restrictions in the case of density
*   filtering and does not need to (but it usuallly close to respect them).
    XD = TOPORSTR tab XD ;
 
*   Save X
*   ------
 
    Wtab.'X_D' = XD ;
 
*   Save S and Y or generate them from X if they do not exist
*   ---------------------------------------------------------
 
    SI (LS ET LY) ;
        Wtab.'S_D' = SD ;
        Wtab.'Y_D' = YD ;
    SINO ;
*       Generate S and Y from X on D
*       Remark: this is done using the filtering and projection
*               parameters of the current cycle while S and Y shoud
*               have been generated at the end of the previous cycle
*               with the previous parameters.
        TOPOPHYS tab ;
    FINS ;
 
SINO ;
    SI (NEG cyc 1) ;
        MESS (CHAI 'Le cycle de depart est le' ' ' cyc '.') ;
        MESS (CHAI 'Or la topologie' ' ' cyc ' n''existe pas.') ;
        ERRE 'Veuillez preciser la topologie du nouveau cycle.' ;
    SINO ;
        SI Ltab.'CONTRAINTE_VOLUME' ;
            Wtab.'X_D' = uniD * tab.'FRACTION_VOLUME_LIMITE' ;
        SINO ;
            Wtab.'X_D' = uniD * 0.5 ;
        FINS ;
*       Get Y from X on D
        TOPOPHYS tab ;
    FINS ;
FINS ;
 
* Remark : After TOPOPHYS, we may not respect the volume fraction of
*          Y in the case of volume constraint at first iteration.
*          It should not be a problem either with OC nor MMA.
*          We could scale the fields to fit the volume fraction but
*          it should be done preserving X_MIN_SEUIL<X<1...
 
FINP ;
 
 
 

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