Procédures | Cast3M

* TOPOFILT  PROCEDUR  FD218221  25/12/18    21:15:06     12429          
 
************************************************************************
** Procedure called by verious TOPOPTIM subprocedures for filtering
** an input field.
**
** Author:
** Guenhael Le Quilliec (LaMe - Polytech Tours)
**
** Version:
** 4.0 2025/12/16 Field to be filtered is given as input, and modified
**                sensitivity filtering is now available. The field can
**                also be normalized if it is an extensive quantity.
**                Added a third filtering function: Helmholtz (PDE)
** 3.0 2021/04/19 Added a second filtering function using a filtering
**                matrix generated with the MFIL operator
** 2.0 2017/11/11 Initial version: default Gibiane filter
************************************************************************
 
DEBP TOPOFILT tab*'TABLE' fldD*'MCHAML' LdX/'LOGIQUE' LN/'LOGIQUE' ;
 
* Input data
* **********
 
Wtab = tab.'WTABLE' ;
Ltab = Wtab.'LOGIQUE' ;
 
* Default value of modified filter logical (for X derivative)
SI (NON (EXIS LdX)) ;
    LdX = FAUX ;
FINS ;
 
* Default value of normalization logical (for non-extensive input fields)
* Remark: Never used by TOPOPTIM, but may be useful elsewhere.
SI (NON (EXIS LN)) ;
    LN = FAUX ;
FINS ;
 
modD = Wtab.'MODELE_D' ;
 
* Pre-processing
* **************
 
* If the field is extensive and needs normalization
SI LN ;
    VD = Wtab.'VOLUME_D' ;
*   Transform it into an intensive field
    fldD = fldD / VD ;
FINS ;
 
* Helmholtz (PDE) filter type (Lazarov et Sigmund 2010)
* *****************************************************
 
SI (EGA tab.'FILTRE' 'EDP') ;
 
*   Build the Helmholtz (PDE) filter matrix
*   ***************************************
 
    oldtab = Wtab.'FILTRE_EDP_ID' ;
 
*   If the mesh (or more precisely, the model and the frozen area)
*   or other filtering parameters have changed, then we update H and Hm.
*   Note: for 'MODEL' and 'ZONE_FIGEE', EGA compare their pointers.
    SI  ((NEG oldtab.'MODEL'      Wtab.'RESOLUTION_B'.'MODELE')
      OU (NEG oldtab.'ZONE_FIGEE' tab.'ZONE_FIGEE')
      OU (NEG oldtab.'RAYON'      tab.'FILTRE_EDP_RAYON')) ;
*       Verification only due to the division by the number of nodes
        SI ((DIME (Wtab.'MAILLAGE_D' ELEM 'TYPE')) > 1) ;
            ERRE (CHAI 'Actuellement le filtre EDP n''est compatible'
                       ' qu''avec des maillages constitues d''un seul'
                       ' et unique type d''element') ;
        FINS ;
*       Update identification model, frozen area, and parameters
        oldtab.'MODEL'      = Wtab.'RESOLUTION_B'.'MODELE' ;
        oldtab.'ZONE_FIGEE' = tab.'ZONE_FIGEE' ;
        oldtab.'RAYON'      = tab.'FILTRE_EDP_RAYON' ;
*       Helmholtz filter equivalent radius
        Hval = tab.'FILTRE_EDP_RAYON' / 2.0 / (3.0**0.5) ;
*       Helmholtz filter model
        HmodD = MODE Wtab.'MAILLAGE_D' 'DIFFUSION' 'FICK'
                    'INCO' 'CO' 'SCAL' ;
*       Helmholtz filter caracteristics
        HmatD = MATE HmodD 'KD' (Hval ** 2) 'CDIF' 1.0 ;
*       Helmholtz filtering matrix
        Wtab.'FILTRE_EDP' = (COND HmodD HmatD) ET (CAPA HmodD HmatD) ;
*       Some cleaning
        OUBL HmodD ;
        OUBL HmatD ;
    FINS ;
 
*   Filter the input field
*   **********************
 
*   Get the node field to be filtered
*   TODO find a way to remove that limitation...
    tmp = NBNO (ELEM Wtab.'MAILLAGE_D' 1) ;
    fldD = CHAN 'CHPO' modD (fldD / tmp) 'SOMM' ;
*   Resolution and rename the component
    fldD = EXCO 'CO' (RESO Wtab.'FILTRE_EDP' fldD) 'SCAL' ;
*   Get the element field from node field
    fldD = CHAN 'CHAM' fldD modD 'GRAVITE' 'SCALAIRE' ;
*   In thermal mode, 'GRAVITE' option is not respected
*   This issue is resolved by using CHAN 'GRAVITE'
    SI (NON Ltab.'MECANIQUE') ;
        fldD = CHAN 'GRAVITE' fldD modD ;
    FINS ;
 
FINS ;
 
* Hat filter matrix type (Bruns et Tortorelli 2001)
* *************************************************
 
SI (EGA tab.'FILTRE' 'CHAPEAU') ;
 
*   Build the hat filter matrix
*   ***************************
 
    oldtab = Wtab.'FILTRE_CHAPEAU_ID' ;
 
*   If the mesh (or more precisely, the model and the frozen area)
*   or other filtering parameters have changed, then we update H and Hm.
*   Note: for 'MODEL' and 'ZONE_FIGEE', EGA compare their pointers.
    SI  ((NEG oldtab.'MODEL'      Wtab.'RESOLUTION_B'.'MODELE')
      OU (NEG oldtab.'ZONE_FIGEE' tab.'ZONE_FIGEE')
      OU (NEG oldtab.'RAYON'      tab.'FILTRE_CHAPEAU_RAYON')
      OU (NEG oldtab.'EXPOSANT'   tab.'FILTRE_CHAPEAU_EXPOSANT')
      OU (NEG oldtab.'CRITERE'    tab.'FILTRE_CHAPEAU_CRITERE')) ;
*       Update identification model, frozen area, and parameters
        oldtab.'MODEL'      = Wtab.'RESOLUTION_B'.'MODELE' ;
        oldtab.'ZONE_FIGEE' = tab.'ZONE_FIGEE' ;
        oldtab.'RAYON'      = tab.'FILTRE_CHAPEAU_RAYON' ;
        oldtab.'EXPOSANT'   = tab.'FILTRE_CHAPEAU_EXPOSANT' ;
        oldtab.'CRITERE'    = tab.'FILTRE_CHAPEAU_CRITERE' ;
*       Prepare other data needed for generating the new filtering matrix
        VD = Wtab.'VOLUME_D' ;
        MESS 'Mise à jour de la matrice du filtre chapeau en cours...' ;
        ctrVD = MANU 'CHPO' Wtab.'MAILLAGE_GRAVITE_D' 'SCAL'
                     (EXTR VD 'VALE' 'SCAL') ;
*       Filtering matrix
        HD = MFIL ctrVD tab.'FILTRE_CHAPEAU_RAYON'
                        tab.'FILTRE_CHAPEAU_EXPOSANT'
                        tab.'FILTRE_CHAPEAU_CRITERE' ;
*       Unitary barycenter field
        uniD = MANU 'CHPO' Wtab.'MAILLAGE_GRAVITE_D' 'SCAL'
                    (EXTR Wtab.'UNIT_D' 'VALE' 'SCAL') ;
*       Weights (sum of each lines of the filtering matrix
        wgtD = uniD * HD ;
*       Normalized filtering matrix
        Wtab.'FILTRE_CHAPEAU_NORMALISE' = NFIL HD wgtD ;
*       Modified filtering matrix
        Wtab.'FILTRE_CHAPEAU_MODIFIE' = NFIL HD wgtD 'MODI' ;
*       Some cleaning as this matrix might be big
        OUBL HD ;
        MESS 'Fin de mise à jour de la matrice du filtre chapeau.' ;
    FINS ;
 
*   Filter the input field
*   **********************
 
*   Nodal field to be filtered
    fldD = MANU 'CHPO' Wtab.'MAILLAGE_GRAVITE_D' 'SCAL' (EXTR fldD 'VALE' 'SCAL') ;
*   Filter the nodal field
*   Remark: multiplication with the rigidity filtering affects the
*           order of the output values
    SI LdX ;
        fldD = fldD * Wtab.'FILTRE_CHAPEAU_MODIFIE' ;
    SINO ;
        fldD = fldD * Wtab.'FILTRE_CHAPEAU_NORMALISE' ;
    FINS ;
*   Filtered element field
*   Remark: Wtab.'MAILLAGE_GRAVITE_D' preserves the initial order to extract the values
    fldD = MANU 'CHML' modD
                  'REPA' 'SCAL' (EXTR fldD 'VALE' Wtab.'MAILLAGE_GRAVITE_D')
                  '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.'THERMIQUE'.'MODELE_D' ;
    FINS ;
 
FINS ;
 
* Default Gibiane filtering
* *************************
 
SI (EGA tab.'FILTRE' 'GIBIANE') ;
    SI (tab.'FILTRE_GIBIANE_TAUX' > 0) ;
*       Filtering
        REPE itr0 tab.'FILTRE_GIBIANE_TAUX' ;
            fldD = CHAN 'CHAM' (CHAN 'CHPO' modD fldD 'MOYE') modD 'GRAVITE' ;
        FIN itr0 ;
*       In thermal mode, 'GRAVITE' option is not respected
*       This issue is resolved by using CHAN 'GRAVITE'
        SI (NON Ltab.'MECANIQUE') ;
            fldD = CHAN 'GRAVITE' fldD modD ;
        FINS ;
    FINS ;
FINS ;
 
* Post-processing
* ***************
 
* Switch back to extensive quantity field
SI LN ;
    fldD = fldD * VD ;
FINS ;
 
FINP fldD ;