Télécharger deduadap.procedur

Retour à la liste

* DEDUADAP  PROCEDUR  GOUNAND   24/09/18    21:15:02     12010          
************************************************************************
* NOM         : DEDUADAP
*
*    Opérateur DEDU option ADAP
*    --------------------------
*    CHPO2 = 'DEDU' 'ADAP' MAIL (CHAM1) (RIG1 (CHPO1)) ;
*
*    Objet :
*    _______
*
*    Génère un champ de déplacement permettant de régulariser un
*    maillage ou de l'adapter suivant une métrique sans changer
*    sa topologie.
*
*    Commentaire :
*    _____________
*
*    MAIL   : maillage à régulariser ou adapter
*
*    CHAM1  : champ par élément aux noeuds donnant une métrique :
*             tenseur symétrique de composantes G11,G22,G12,...
*             (par défaut, le tenseur unité)
*
*    RIG1   : Conditions sur les déplacements
*    CHPO1    (par défaut, on bloque les noeuds frontières de MAIL)
*
*    CHPO2  : champ de déplacement.
*
*    Notes :
*    _______
*
*    L'option 'ADAP' est censée fonctionner sans conditions sur le
*    maillage.
*
*    Référence principale :
*    ______________________
*
*@Article{huang4,
*  author =       {Weizhang Huang},
*  title =        {Variational Mesh Adaptation:
*                  Isotropy and Equidistribution},
*  journal =      {JCP},
*  year =         {2001},
*  volume =       {174},
*  pages =     {903-924},
*  endroit = {Classeur Mesh Movement (VIIIb)}
*}
*
*  Je dois faire également un rapport technique...
*
*
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 05/04/2006, version initiale
* HISTORIQUE : v1, 05/04/2006, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
*
*
'DEBPROC' DEDUADAP ;
*
* Lecture des arguments
*
'ARGUMENT' mail*'MAILLAGE' ;
*
'ARGUMENT' rblo/'RIGIDITE' ;
lrb = 'EXISTE' rblo ;
'SI' lrb ;
   'ARGUMENT' cblo/'CHPOINT' ;
   lcb = 'EXISTE' cblo ;
'SINON' ;
   lcb = FAUX ;
'FINSI' ;
*
lmotcle = 'MOTS' 'METR' 'DISG' 'THET' 'GAMM' 'NITM' 'ACVG' 'METG'
   'IDIR' 'TINV' 'DENS' ;
lchad = FAUX ; lchp = FAUX ;
ldisg   = FAUX ;
theta   = 0.2D0 ;
gamma   = 2.D0 ;
*theta   = 0.2D0 ;
*gamma   = 2.D0 ;
itmax   = 40 ;
acccvg  = VRAI ;
lmg = FAUX ;
lidir = FAUX ;
ltinv = FAUX ;
*
'REPETER' imotcle ;
   'ARGUMENT' motcle/'MOT' ;
   'SI' ('NON' ('EXISTE' motcle)) ;
      'QUITTER' imotcle ;
   'FINSI' ;
   lmc = 'EXISTE' lmotcle motcle ;
   'SI' ('NON' lmc) ;
*1052 2
*Mot-cle incorrect "%M1:4". Voici la liste des valeurs admises : %M5:40
      'ERRE' 1052 'AVEC' ('CHAI' motcle
'METR DISG THET GAMM NITM ACVG METG IDIR TINV') ;
   'FINSI' ;
*
   'SI' ('EGA' motcle 'METR') ;
      'ARGUMENT' chad/'MCHAML' ;
      lchad = 'EXISTE' chad ;
      'SI' ('NON' lchad) ;
         'ARGUMENT' chp/'CHPOINT' ;
         lchp = 'EXISTE' chp ;
         'SI' lchp ;
            'ARGUMENT' metdisc/'MOT' ;
         'FINS' ;
         lmetdisc = 'EXIS' metdisc ;
      'SINON' ;
         argu kmodl1*'MMODEL' ;
         lchp = FAUX ;
      'FINSI' ;
   'SINON' ;
     'SI' ('EGA' motcle 'DENS') ;
       'ARGUMENT' chp*'CHPOINT' ;
       chp    = 'REDU' chp mail ;
       kmodl1 = 'MODE' mail 'THERMIQUE' ;
       chad   = 'CHAN' 'CHAM' kmodl1 chp ;
       chad   = (chad ** -2) 'NOMC' 'G11' ;
       'SI' (('VALE' DIME) '>' 1) ;
         chad0  = 0.*chad ;
         chad   = chad 'ET' (chad 'NOMC' 'G22') 'ET' (chad0 'NOMC' 'G21') ;
       'FINS' ;
       'SI' (('VALE' DIME) '>' 2) ;
         chad   = chad 'ET' (chad 'NOMC' 'G33') 'ET' (chad0 'NOMC' 'G31') 'ET' (chad0 'NOMC' 'G32') ;
       'FINS' ;
       lchad  = VRAI ;
       lchp   = FAUX ;
     'FINSI' ;
   'FINSI' ;
*
   'SI' ('EGA' motcle 'DISG') ;
      'ARGUMENT' gdisc*'MOT' ;
      ldisg = VRAI ;
   'FINSI' ;
*
   vprec = ('VALE' 'PREC') '**' 0.75 ;
   mvprec = '*' vprec -1. ;
   'SI' ('EGA' motcle 'THET') ;
      'ARGUMENT' theta*'FLOTTANT' ;
      lok = 'ET' ('>' theta mvprec) ('<' theta ('+' 1. vprec)) ;
      'SI' ('NON' lok) ;
*  42 2   %m1:8 = %r1 non compris entre %r2 et %r3
         'ERRE' 42 'AVEC' 'FLOT1' ('PROG' theta 0. 1.) ;
      'FINSI' ;
   'FINSI' ;
*
   'SI' ('EGA' motcle 'GAMM') ;
      'ARGUMENT' gamma*'FLOTTANT' ;
      lok = ('>' gamma ('-' 1.D0 vprec)) ;
      'SI' ('NON' lok) ;
*  41 2 %m1:8 = %r1 inferieur a %r2
         'ERRE' 41 'AVEC' 'FLOT2' ('PROG' gamma 1.) ;
      'FINSI' ;
   'FINSI' ;
*
   'SI' ('EGA' motcle 'NITM') ;
      'ARGUMENT' itmax*'ENTIER' ;
      lok = ('>' itmax 0) ;
      'SI' ('NON' lok) ;
*  41 2 %m1:8 = %r1 inferieur a %r2
         'ERRE' 41 'AVEC' 'ENTI1' ('PROG' itmax 1.) ;
      'FINSI' ;
   'FINSI' ;
*
   'SI' ('EGA' motcle 'ACVG') ;
      'ARGUMENT' acccvg*'LOGIQUE' ;
   'FINSI' ;
*
   'SI' ('EGA' motcle 'METG') ;
      'ARGUMENT' methgau*'MOT' ;
      lmg = VRAI ;
   'FINSI' ;
*
   'SI' ('EGA' motcle 'IDIR') ;
      'ARGUMENT' idir*'ENTIER' ;
      lidir = VRAI ;
   'FINSI' ;
*
   'SI' ('EGA' motcle 'TINV') ;
      'ARGUMENT' tinv*'TABLE' ;
      ltinv = VRAI ;
   'FINSI' ;
*
'FIN' imotcle ;
ladap = 'OU' lchad lchp ;
*
'ARGUMENT' debug/'LOGIQUE' ;
'SI' ('NON' ('EXISTE' debug)) ;
   debug = FAUX ;
'FINSI' ;
*
* Initialisations
*
idim = 'VALEUR' 'DIME' ;
imod = 'VALEUR' 'MODE' ;
vdim = DEADUTIL 'DIMM' mail ;
vtyp = DEADUTIL 'TYPM' mail ;
laxi = DEADUTIL 'AXI?' ;
lsph = DEADUTIL 'SPH?' ;
*
'SI' ('OU' ('<' idim 1) ('>' idim 3)) ;
* 709 2 Fonction indisponible en dimension %i1.
   'ERREUR' 709 'AVEC' idim ;
'FINSI' ;
'SI' (('EGA' imod 'AXIS') 'OU' ('EGA' imod 'UNIDAXIS') 'OU'
   ('EGA' imod 'FOUR') 'OU' ('EGA' imod 'SPHE')) ;
*-105 0 Mode de calcul actuel %m1:32
   'ERRE'  -105 'AVEC' imod ;
* 710 2 Fonction indisponible pour ce mode de calcul
   'ERRE' 710 ;
'FINSI' ;
*
vquaf = ('EGA' vtyp 'QUAF') ;
'SI' ('ET' ldisg ('NON' vquaf)) ;
   'MESS' 'DISG option :' ;
*  66 2 L'objet %m1:8 doit etre de type %m9:16
   'ERRE' 66 'AVEC' 'MAIL    QUAF' ;
'FINSI' ;
*'SI' ('ET' lchp ('NON' vquaf)) ;
*   'MESS' 'METR+CHPO option :' ;
**  66 2 L'objet %m1:8 doit etre de type %m9:16
*   'ERRE' 66 'AVEC' 'MAIL    QUAF' ;
*'FINSI' ;
*
* Paramètres du solveur non-linéaire
*
* Evaluation de la matrice tangente
*   iktan = 1 : Matrice tangente exact (chère)
*   iktan = 2 :  " approchée 1 (on néglige les termes extradiagonaux
*                                    (UX,FY)...)
*                             (par défaut)
*   iktan = 3 :  " approchée 2 (on néglige en plus les dérivées
*                                 extradiagonales (ddxi, ddeta))
*                             (ne fonctionne pas)
*
* Accélération de convergence
*   acccvg = VRAI : accélération à la PV
*
* Période de recalcul de la matrice tangente
*   rktan = i (toutes les i itérations non-linéaires)
*
*   itmax : Nombre maxi d'itérations
*   rfonc : critère d'arrêt sur la variation de la fonctionnelle
*   fback : facteur divisif pour la relaxation lors d'un backtracking
*   nback : nombre maxi de backtracking
*   fvdet : critère pour le backtracking, variation maxi de la valeur
*           du jacobien
 
iktan = 2 ;
*!!!iktan = 1 ;
rktan  = 1 ;
'SI' debug ;
   rfonc = 5.D-8 ;
'SINON' ;
   rfonc = 5.D-2 ;
'FINSI' ;
fback = 2.D0 ; fvdet = 4.D0 ;
nback = 10 ;   damp  = 1.0D0  ;
*
* Maillage
*
'SI' vquaf ;
   _mt = mail ;
'SINON' ;
   _mt = 'CHANGER' mail 'QUAF' ;
'FINSI' ;
*
* Inconnus et discrétisation
*
'SI' ('NON' lmg) ;
   'SI' ('EGA' vtyp 'LINE') ;
      methgau = 'GAR1' ;
   'SINON' ;
      methgau = 'GAR2' ;
   'FINSI' ;
'FINSI' ;
'SI' ('NON' ldisg) ;
   gdisc   = vtyp ;
'FINSI' ;
*
'SI' ('OU' laxi lsph) ;
   lu = 'MOTS' 'UR' 'UZ' ;
   lf = 'MOTS' 'FR' 'FZ' ;
'SINON' ;
   lu = 'MOTS' 'UX' 'UY' 'UZ' ;
   lf = 'MOTS' 'FX' 'FY' 'FZ' ;
'FINSI' ;
*
lextr = 'LECT' 1 'PAS' 1 idim ;
*
lpr = 'EXTRAIRE' lu lextr ; ldu = 'EXTRAIRE' lf lextr ;
*
* Métrique
*
'SI' ladap;
   'SI' lchad ;
*      $mt = 'MODE' _mt 'NAVIER_STOKES' 'QUAF' ;
*      chpmet = 'KCHA' $mt chad 'CHPO' 'QUAF' ;
*      metdisc = 'CSTE' ;
*      $mt = 'MODELISER' mail 'THERMIQUE' ;
      chpmet  = 'CHANGER' 'CHPO' kmodl1 chad ;
      metdisc = gdisc ;
   'FINSI' ;
   'SI' lchp ;
      chpmet = chp ;
      'SI' ('NON' lmetdisc) ;
         metdisc = gdisc ;
      'FINS' ;
   'FINSI' ;
'FINSI' ;
*
* Conditions aux limites sur les déplacements
*
'SI' ('NON' lrb) ;
   'SI' ('EGA' idim 1) ;
* Faute de mieux : contour ne fonctionne pas en 1D
      bord = 'ET' ('POIN' mail 'INITIAL')
                  ('POIN' mail 'FINAL') ;
   'FINSI' ;
   'SI' ('EGA' idim 2) ;
      bord = 'CONTOUR' mail ;
*      rblo = 'BLOQUE' 'UX' 'UY' bord ;
   'FINSI' ;
   'SI' ('EGA' idim 3) ;
      bord = 'ENVELOPPE' mail ;
*      rblo = 'BLOQUE' 'UX' 'UY' 'UZ' bord ;
   'FINSI' ;
   rblo = 'BLOQUE' 'DEPL' bord ;
'FINSI' ;
*
* Initilisations diverses
*
iniform = 'FORME' ;
tres = 'TABLE' ;
dx   = 'MANUEL' 'CHPO' mail lpr ('PROG') ;
det0 = DEADJACO _mt gdisc methgau ;
tyde = 'TYPE' det0 ;
'SI' ('EGA' tyde 'ENTIER') ;
   chinfo = 'CHAINE' 'DEDUADAP : the mesh is already skew !' ;
   'MESSAGE' chinfo ;
* 426 2 Maillage incorrect
   'ERRE' 426 ;
'FINSI' ;
'SI' ladap ;
   fonc0 = DEADFONC _mt gdisc methgau theta gamma
                    chpmet metdisc 'ELEM' ;
   'SI' lidir ;
      res   = DEADRESI _mt gdisc methgau theta gamma ldu
         chpmet metdisc idir ;
   'SINON' ;
      res   = DEADRESI _mt gdisc methgau theta gamma ldu
         chpmet metdisc ;
   'FINSI' ;
'SINON' ;
   fonc0 = DEADFONC _mt gdisc methgau theta gamma 'ELEM' ;
   'SI' lidir ;
      res   = DEADRESI _mt gdisc methgau theta gamma ldu idir ;
   'SINON' ;
      res   = DEADRESI _mt gdisc methgau theta gamma ldu ;
   'FINSI' ;
'FINSI' ;
dampi = damp ;
fonci = fonc0 ;
deti  = det0 ;
'SI' acccvg ;
   znacce = 3 ; itdep = 3 ;
   acfp1  =  'COPIER' ('*' res 0.) ;  acfp2  = acfp1 ;
   acfp3  = acfp1 ;  acfep1 = acfp1;
   acfep2 = acfp1 ;  correc=0.; freap = 0. ;
'FINSI' ;
*
* Algorithme non linéaire
*
cvgok = FAUX ;
'REPETER' it itmax ;
   dampi = 'MINIMUM' ('PROG' damp ('*' dampi fback)) ;
   'SI' debug ;
      ch = 'CHAINE' '  Itération : ' &it ;
      'MESSAGE' ch ;
      ch = 'CHAINE' '    dampi=' dampi ;
      'MESSAGE' ch ;
   'FINSI' ;
*
*  Calcul du résidu
*
*    partie liée à la fonctionnelle
   'SI' ladap ;
      'SI' lidir ;
         resf = DEADRESI _mt gdisc methgau theta gamma ldu
            chpmet metdisc idir ;
      'SINON' ;
         resf = DEADRESI _mt gdisc methgau theta gamma ldu
            chpmet metdisc ;
      'FINSI' ;
   'SINON' ;
      'SI' lidir ;
         resf = DEADRESI _mt gdisc methgau theta gamma ldu idir ;
      'SINON' ;
         resf = DEADRESI _mt gdisc methgau theta gamma ldu ;
      'FINSI' ;
   'FINSI' ;
*    partie liée aux blocages
   resb = '*' rblo dx ;
   'SI' lcb ;
      resb = resb '-' cblo ;
   'FINSI' ;
   res = resf '+' resb ;
*  acceleration de convergence
   'SI' acccvg ;
      correcp = correc; correc = 0;
      acfp0 = 'ENLEVER' (res - freap) 'FLX' ;
      acfep0 = acfp0;
      acfep0 = acfep0 - correcp ;
      'SI' ('MULT'  &it znacce)   ;
         'SI' ('>' &it itdep) ;
            'SI' ('>' dampi ('*' damp 0.9)) ;
               'SI' debug ;
                  'MESSAGE' 'Convergence acceleration' ;
               'FINSI' ;
                correc = act3  acfep2  acfep1 acfep0
                              acfp3 acfp2 acfp1  acfp0 ;
               res = '-' res correc;
            'FINSI' ;
         'FINSI' ;
      'FINSI' ;
      'SI' ('>' &it 3) ;
         'DETRUIT' acfp3 ;
         'DETRUIT' acfep2 ;
      'FINSI' ;
      acfp3  = acfp2  ; acfp2  = acfp1  ; acfp1  = acfp0 ;
      acfep2 = acfep1 ; acfep1 = acfep0 ;
      'MENAGE' ;
   'FINSI' ;
*
*  Calcul de la matrice tangente
*
   'SI' ('MULT' ('-' &it 1) rktan) ;
      'SI' ladap ;
         'SI' lidir ;
            jac  = DEADKTAN _mt gdisc methgau theta gamma
                            lpr ldu chpmet metdisc iktan idir ;
         'SINON' ;
            jac  = DEADKTAN _mt gdisc methgau theta gamma
                            lpr ldu iktan chpmet metdisc ;
         'FINSI' ;
      'SINON' ;
         'SI' lidir ;
            jac  = DEADKTAN _mt gdisc methgau theta gamma
                            lpr ldu iktan idir ;
         'SINON' ;
            jac  = DEADKTAN _mt gdisc methgau theta gamma
                            lpr ldu iktan ;
         'FINSI' ;
      'FINSI' ;
      jact = 'ET' jac rblo ;
   'FINSI' ;
*
* Résolution du problème linéarisé
*
   'SI' ltinv ;
      mat smb = jact ('*' -1.D0 res) ;
      'SI' ('EGA' (tinv . 'TYPINV') 0) ;
         'OPTI' impi 0 ;
         sol = 'RESOUD' mat smb 'NOID' ;
         'OPTI' impi 0 ;
      'SINON' ;
         'SI' ('EXISTE' tinv 'LTIME') ;
            ltime = tinv . 'LTIME' ;
         'SINON' ;
            ltime = FAUX ;
         'FINSI' ;
*
         'SI' ('EGA' ltime vrai) ;
            sol tt  = 'KRES' mat smb 'TYPI' tinv ;
            'LISTE' tt ;
         'SINON' ;
            sol  = 'KRES' mat smb 'TYPI' tinv ;
         'FINSI' ;
      'FINSI' ;
      ddx   = sol ;
   'SINON' ;
      ddx   = 'RESO' jact ('*' -1.D0 res) ;
   'FINSI' ;
*   'LISTE' ddx ;
*
*
*
   'SI' acccvg ;
      freap = reac rblo ddx ;
*   resx  = 'ENLEVER' ('-' res (reac rblo ddx)) 'FLX' ;
*   resx2  = 'ENLEVER' (res) 'FLX' ;
   'FINSI' ;
*
* Backtracking
*
   backok = FAUX ;
   'REPETER' iback nback ;
     'SI' ('>' &iback 1) ;
        dampi = '/' dampi fback ;
        'SI' debug ;
           ch = 'CHAINE' '    dampi=' dampi ;
           'MESSAGE' ch ;
         'FINSI' ;
      'FINSI' ;
      ddxi = '*' ddx dampi ;
      depx = 'EXCO' lpr ddxi lpr 'NOID' ;
* Test si le déplacement calculé inverse un jacobien
* ou le change trop
      oldconf = 'FORME' ;
      'FORME' depx ;
      detip = DEADJACO _mt gdisc methgau ;
      'FORME' oldconf ;
      tyde = 'TYPE' detip ;
      'SI' ('EGA' tyde 'ENTIER') ;
         'SI' debug ;
            ch = 'CHAINE' '    Warning : inv. loc. jacobien !' ;
            'MESSAGE' ch ;
         'FINSI' ;
      'SINON' ;
         vardet = ('/' detip deti) ;
         mivd = 'MINIMUM' vardet ;
         mavd = 'MAXIMUM' vardet ;
         'SI' debug ;
            'MESSAGE' ('CHAINE' 'Mini var jaco = ' mivd ) ;
            'MESSAGE' ('CHAINE' 'Maxi var jaco = ' mavd ) ;
         'FINSI' ;
         bigvar = 'OU' ('>' mavd fvdet) ('<' mivd ('/' 1.D0 fvdet)) ;
         'SI' bigvar ;
            'SI' debug ;
               ch = 'CHAINE'
               '    Warn : trop grande variation du jaco !' ;
               'MESSAGE' ch ;
            'FINSI' ;
         'SINON' ;
            backok = VRAI ;
            'QUITTER' iback ;
         'FINSI' ;
      'FINSI' ;
   'FIN' iback ;
   'SI' ('NON' backok) ;
      chinfo1 = 'CHAINE'
            'DEDUADAP : Backtracking failed to converge !' ;
      'MESSAGE' chinfo1 ;
      chinfo2 = 'CHAINE' 'Please check the output displacement'
      'MESSAGE' chinfo2 ;
      dep = 'EXCO' lpr dx lpr ;
      'RESPRO' dep ;
      'QUITTER' DEDUADAP ;
   'FINSI' ;
   dx = '+' dx ddxi ;
   'FORME' depx ;
   'SI' ladap ;
      foncip = DEADFONC _mt gdisc methgau theta gamma
                        chpmet metdisc 'ELEM' ;
   'SINON' ;
      foncip = DEADFONC _mt gdisc methgau theta gamma 'ELEM' ;
   'FINSI' ;
*
* Critère de convergence
*
   vrrfon = '/' ('-' foncip fonci) fonci ;
   critconv = 'MAXIMUM' vrrfon 'ABS' ;
   tres . &it = critconv ;
   cvgok = ('<' critconv rfonc) ;
   'SI' debug ;
      chmess = 'CHAINE' '  critere    = ' critconv ;
      'MESSAGE' chmess ;
   'FINSI' ;
   deti = detip ;
   fonci = foncip ;
   'SI' debug ;
      tit   = 'CHAINE' 'Maillage ; it=' &it ' crit=' critconv ;
      TRAC 'CACH' mail 'TITR' tit 'NCLK' ;
   'FINSI' ;
   'SI' cvgok ;
      'QUITTER' it ;
   'FINSI' ;
'FIN' it ;
'FORME' iniform ;
dep = 'EXCO' lpr dx lpr ;
*
itfin = ('-' &it 1) ;
'SI' debug ;
   ch = 'CHAINE' '  Iteration finale : ' itfin ;
   'MESSAGE' ch ;
   lit = 'PROG' 1. 'PAS' 1. 'NPAS' ('-' itfin 1) ;
   lres = 'PROG' ;
   'REPETER' it itfin ;
      lres = 'ET' lres ('PROG'
         ('/' ('LOG' ('+' (tres . &it ) 1.D-16))
              ('LOG' 10.D0) )) ;
   'FIN' it ;
   titev = 'CHAINE' 'LOG10 Critere (iterations)' ;
   evres = 'EVOL' 'MANU' lit lres 'TITR' titev ;
   'DESSIN' evres ;
'FINSI' ;
*
'RESPRO' dep ;
*
* End of procedure file DEDUADAP
*
'FINPROC' ;
 
 

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