Télécharger la plaquette Cast3M


fr en

Télécharger dynamic.procedur

Retour à la liste

Numérotation des lignes : 

  1. * DYNAMIC   PROCEDUR  FANDEUR   22/03/01    21:15:03     11301          
  2. * ======================================================================
  3. * CALCUL DE LA REPONSE DYNAMIQUE PAR INTEGRATION TEMPORELLE
  4. * (SCHEMA DE NEWMARK OU HHT OU ALPHA-GENERALISE)
  5. * ======================================================================
  6. * Creation : ???, 1988 ???
  7. * Refonte  : Benoit Prabel, 23/07/2015
  8. * Modifications : JCARDO 15/02/2017 => post-traitement tableur + VTK
  9. * ======================================================================
  10. DEBP DYNAMIC ETAB*'TABLE' ;
  11.  
  12.    IECHO = VALE 'ECHO' ;
  13.  
  14. *  +-----------------------------------------------------------------+
  15. *  |                                                                 |
  16. *  |           L E C T U R E   D E S   P A R A M E T R E S           |
  17. *  |                                                                 |
  18. *  +-----------------------------------------------------------------+
  19.  
  20. *  ================
  21. *  OPTIONS INTERNES
  22. *  ================
  23.  
  24.    KMENAG = VRAI ;
  25.  
  26.  
  27. *  ===================================================
  28. *  PARAMETRES DE L'ALGORITHME D'INTEGRATION TEMPORELLE
  29. *  ===================================================
  30.  
  31.    SI ((EXIS ETAB 'ALPHA_F') ET (EXIS ETAB 'RHO_INF')) ;
  32.       MESS 'Veuillez fournir 1 seul parametre parmi :' ;
  33.       MESS '  ALPHA_F => algorithme HHT ' ;
  34.       MESS '  RHO_INF => algorithme alpha-generalise' ;
  35.       ERRE 21 ;
  36.    FINS ;
  37.  
  38. *  => Newmark acceleration moyenne (IALGO = 0)
  39.    IALGO  = 0  ;
  40.    ALPHAM = 0. ;
  41.    ALPHAF = 0. ;
  42.  
  43. *  => HHT (IALGO = 1)
  44.    SI (EXIS ETAB 'ALPHA_F') ;
  45.       ALPHAF = ETAB.'ALPHA_F' ;
  46.       SI (IECHO >EG 0) ;
  47.          MESS 'ALGO = HHT (ALPHA_F=' ALPHAF ')' ;
  48.       FINS ;
  49.       IALGO = 1 ;
  50.    FINS ;
  51.  
  52. *  => alpha-generalise (IALGO = 2) 
  53.    SI (EXIS ETAB 'RHO_INF') ;
  54.       RHOINF = ETAB . 'RHO_INF' ;
  55.       ALPHAM = ((2.*RHOINF) - 1.) / (1.+RHOINF) ;
  56.       ALPHAF = RHOINF / (1.+RHOINF) ;
  57.       SI (IECHO >EG 0) ;
  58.          MESS 'ALGO = alpha-generalise (ALPHA_M=' ALPHAM 
  59.                                     ' & ALPHA_F=' ALPHAF ')' ;
  60.       FINS ;
  61.       IALGO = 2 ;
  62.    FINS ;
  63.  
  64.  
  65. *  ====================
  66. *  CONDITIONS INITIALES
  67. *  ====================
  68.  
  69.    U0 = ETAB.'DEPL' ;
  70.    V0 = ETAB.'VITE' ;
  71.  
  72.  
  73. *  =======================
  74. *  MATRICES ET CHARGEMENTS
  75. *  =======================
  76.  
  77. *  Matrices et chargements constants
  78.    K = ETAB.'RIGI' ;
  79.    M = ETAB.'MASS' ;
  80.    KAMOR = EXIS ETAB 'AMOR' ;
  81.    SI KAMOR ;
  82.       C = EXTR (ETAB.'AMOR') 'RIGI' 'NOMU' ;
  83.       SI (IECHO >EG 0) ;
  84.          MESS 'Presence d une matrice d amortissement' ;
  85.       FINS ;
  86.    FINS ;
  87.    CHAR1 = ETAB.'CHAR' ;
  88.  
  89. *  Matrices et chargements "utilisateur" via la procedure CHARMECA
  90.    KCHARM = FAUX ;
  91.    SI (EXIS ETAB 'CHARMECA') ;
  92.       KCHARM = ETAB.'CHARMECA' ;
  93.    FINS ;
  94.  
  95.  
  96. *  ====================================
  97. *  APPEL A VITEUNIL POUR LES CONTACTS ?
  98. *  ====================================
  99.  
  100.    SI (EXIS ETAB 'VITEUNIL') ;
  101.       KVITUN = ETAB.'VITEUNIL' ;
  102.    SINON ;
  103.       KVITUN = (EGA IALGO 0) ;
  104.    FINS ;
  105.  
  106.    SI (KVITUN ET (NEG IALGO 0)) ;
  107.       MESS 'VITEUNIL compatible uniquement avec le schema de Newmark ' 
  108.            'acceleration moyenne' ;
  109.       ERRE 21 ;
  110.    FINS ;
  111.  
  112.  
  113. *  ====================
  114. *  PARAMETRES TEMPORELS
  115. *  ====================
  116.  
  117. *  Nouvelle syntaxe :
  118.    SI (EXIS ETAB 'TEMPS_CALCULES') ;
  119.       TPROG = ETAB.'TEMPS_CALCULES' ;
  120.  
  121.       SI ((MINI ((ENLE TPROG 1) - (ENLE TPROG (DIME TPROG)))) <EG 0.) ;
  122.          MESS 'La liste des TEMPS_CALCULES doit etre strictement ' 
  123.               'croissante !' ;
  124.          ERRE 21 ;
  125.       FINS ;
  126.  
  127. *  Ancienne syntaxe : on avertit et on traduit
  128.    SINON ; SI (EXIS ETAB  'FREQ') ;
  129.       MESS '!!!  SYNTAXE BIENTOT OSBOLETE : CONSULTEZ LA NOTICE  !!!' ;
  130.       MESS '!!!  (LE CALCUL SE POURSUIT MALGRE TOUT...)          !!!' ;
  131.       DTFREQ = 0.25 / ETAB.'FREQ' ;
  132.       SI (EXIS ETAB 'DEBU') ;
  133.          T0 = ETAB.'DEBU' ;
  134.       SINON ;
  135.          T0 = 0. ;
  136.       FINS ;
  137.  
  138.       TSORT = ORDO ETAB.'INST' ;
  139.       ETAB.'TEMPS_SAUVES' = TSORT ;
  140.  
  141.       TFIN = DTFREQ * (ENTI 'SUPERIEUR' (MAXI TSORT  / DTFREQ)) ;
  142.       TPROG = PROG T0 PAS DTFREQ TFIN ;
  143.       ETAB.'TEMPS_CALCULES' = TPROG ;
  144.  
  145.    SINON ;
  146.       MESS 'Liste des TEMPS_CALCULES absente de la table !' ;
  147.       ERRE 21 ;
  148.    FINS ; FINS ;
  149.  
  150.    NPAS = DIME TPROG ;
  151.  
  152.  
  153. *  ========================
  154. *  PARAMETRES DE SAUVEGARDE
  155. *  ========================
  156.  
  157.    IGIB = 0 ;
  158.    ICSV = 0 ;
  159.    IVTK = 0 ;
  160.  
  161.  
  162. *  SAUVEGARDE VTK
  163. *  --------------
  164.  
  165. *  Instants de sauvegarde
  166.    SI (EXIS ETAB 'PAS_SAUVES_VTK') ;
  167.       SI (EGA (TYPE ETAB.'PAS_SAUVES_VTK') 'ENTIER') ;
  168.          IPASVTK = ETAB.'PAS_SAUVES_VTK' ;
  169.       SINON ;
  170.          SI (EGA (ETAB.'PAS_SAUVES_VTK') 'TOUS')  ;
  171.             IPASVTK = 1    ;
  172.          SINON ; SI (EGA (ETAB.'PAS_SAUVES_VTK') 'FINAL') ;
  173.             IPASVTK = NPAS ;
  174.          SINON ;
  175.             MESS 'PAS_SAUVES_VTK est de type incorrect !' ;
  176.             ERRE 21 ;
  177.          FINS ; FINS ;
  178.       FINS ;
  179.       IVTK = 1 ;
  180.    SINON ; SI (EXIS ETAB 'TEMPS_SAUVES_VTK') ;
  181.       SI (NEG (TYPE ETAB.'TEMPS_SAUVES_VTK') 'LISTREEL') ;
  182.          MESS 'TEMPS_SAUVES_VTK doit etre de type LISTREEL !' ;
  183.          ERRE 21 ;
  184.       FINS ;
  185.       IVTK = 2 ;
  186.    FINS ; FINS ;
  187.  
  188. *  Emplacement de sauvegarde
  189.    SI (EXIS ETAB 'FICHIER_VTK') ;
  190.       FICVTK = ETAB.'FICHIER_VTK' ;
  191.    SINON ;
  192.       FICVTK = 'DYNAMIC' ;
  193.    FINS ;
  194.  
  195. *  Liste des maillages a sauvegarder
  196.    SI (IVTK > 0) ;
  197.       SI (EXIS ETAB 'MAILLAGE_VTK') ;
  198.          MAIVTK = ETAB.'MAILLAGE_VTK' ;
  199.       SINON ;
  200.          MESS 'L indice MAILLAGE_VTK est obligatoire !' ;
  201.          ERRE 21 ;
  202.       FINS ;
  203.    FINS ;
  204.  
  205.  
  206. *  SAUVEGARDE CSV (TABLEUR)
  207. *  ------------------------
  208.  
  209. *  Instants de sauvegarde
  210.    SI (EXIS ETAB 'PAS_SAUVES_CSV') ;
  211.       SI (EGA (TYPE ETAB.'PAS_SAUVES_CSV') 'ENTIER') ;
  212.          IPASCSV = ETAB.'PAS_SAUVES_CSV' ;
  213.       SINON ;
  214.          SI (EGA (ETAB.'PAS_SAUVES_CSV') 'TOUS')  ;
  215.             IPASCSV = 1    ;
  216.          SINON ; SI (EGA (ETAB.'PAS_SAUVES_CSV') 'FINAL') ;
  217.             IPASCSV = NPAS ;
  218.          SINON ;
  219.             MESS 'PAS_SAUVES_CSV est de type incorrect !' ;
  220.             ERRE 21 ;
  221.          FINS ; FINS ;
  222.       FINS ;
  223.       ICSV = 1 ;
  224.    SINON ; SI (EXIS ETAB 'TEMPS_SAUVES_CSV') ;
  225.       SI (NEG (TYPE ETAB.'TEMPS_SAUVES_CSV') 'LISTREEL') ;
  226.          MESS 'TEMPS_SAUVES_CSV doit etre de type LISTREEL !' ;
  227.          ERRE 21 ;
  228.       FINS ;
  229.       ICSV = 2 ;
  230.    FINS ; FINS ;
  231.  
  232. *  Emplacement de sauvegarde (eventuellement a la fin du calcul)
  233.    SI (EXIS ETAB 'FICHIER_CSV') ;
  234.       KFICCSV = VRAI ;
  235.       FICCSV = ETAB.'FICHIER_CSV' ;
  236.    SINON ;      
  237.       KFICCSV = FAUX ;
  238.    FINS ;
  239.  
  240. *  Liste des noeuds et composantes a sauvegarder
  241.    SI (ICSV > 0) ;
  242.       SI (EXIS ETAB 'MAILLAGE_CSV') ;
  243.          MAICSV = CHAN 'POI1' ETAB.'MAILLAGE_CSV' ;
  244.          NELCSV = NBEL MAICSV ;
  245.          SI (NELCSV EGA 0) ;
  246.             MESS 'Le maillage MAILLAGE_CSV est vide !' ;
  247.             ERRE 21 ;
  248.          FINS ;
  249.       SINON ;
  250.          MESS 'L indice MAILLAGE_CSV est obligatoire !' ;
  251.          ERRE 21 ;
  252.       FINS ;
  253.  
  254.       KLISCO = FAUX ;
  255.       SI (EXIS ETAB 'COMPOSANTES_CSV') ;
  256.          LISCO = ETAB.'COMPOSANTES_CSV' ;
  257.          KLISCO = VRAI ;
  258.          SI ((DIME LISCO) EGA 0) ;
  259.             MESS 'La liste COMPOSANTES_CSV est vide !' ;
  260.             ERRE 21 ;
  261.          FINS ;
  262.       FINS ;
  263.    FINS ;
  264.  
  265.  
  266. *  SAUVEGARDE GIBIANE
  267. *  ------------------
  268.  
  269. *  Instants de sauvegarde
  270.    SI (EXIS ETAB 'PAS_SAUVES') ;
  271.       SI (EGA (TYPE ETAB.'PAS_SAUVES') 'ENTIER') ;
  272.          IPASGIB = ETAB.'PAS_SAUVES' ;
  273.       SINON ;
  274.          SI (EGA (ETAB.'PAS_SAUVES') 'TOUS')  ;
  275.             IPASGIB = 1    ;
  276.          SINON ; SI (EGA (ETAB.'PAS_SAUVES') 'FINAL') ;
  277.             IPASGIB = NPAS ;
  278.          SINON ;
  279.             MESS 'PAS_SAUVES est de type incorrect !' ;
  280.             ERRE 21 ;
  281.          FINS ; FINS ;
  282.       FINS ;
  283.       IGIB = 1 ;
  284.    SINON ; SI (EXIS ETAB 'TEMPS_SAUVES') ;
  285.       SI (NEG (TYPE ETAB.'TEMPS_SAUVES') 'LISTREEL') ;
  286.          MESS 'TEMPS_SAUVES doit etre de type LISTREEL !' ;
  287.          ERRE 21 ;
  288.       FINS ;
  289.       IGIB = 2 ;
  290.    SINON ;
  291.       IPASGIB = 4 ;
  292.       IGIB = 1 ;
  293.    FINS ; FINS ;
  294.  
  295. *  Sauvegarde incrementale ?
  296. *  (appel a "SAUV MUET ETAB" => ecriture sur disque)
  297.    KINCR = FAUX ;
  298.    SI (EXIS ETAB 'SAUV') ;
  299.       KINCR = ETAB.'SAUV' ;
  300.       SI (KINCR ET (IECHO >EG 0)) ;
  301.          MESS 'SAUVegarde incrementale activee' ;
  302.       FINS ;
  303.    FINS ;
  304.  
  305. *  Fantomisation des resultats au fur et a mesure du calcul
  306.    KECON = FAUX ;
  307.    SI (EXIS ETAB 'ECON') ;
  308.       KECON = ETAB.'ECON' ;
  309.    FINS ;
  310.  
  311. *  Restriction de la geometrie a sauvegarder
  312.    KMAIGIB = EXIS ETAB 'MAILLAGE_SAUVE' ;
  313.    SI KMAIGIB ;
  314.       MAIGIB = ETAB.'MAILLAGE_SAUVE' ;
  315.    FINS ;
  316.  
  317.  
  318.  
  319. *  +-----------------------------------------------------------------+
  320. *  |                                                                 |
  321. *  |                  I N I T I A L I S A T I O N S                  |
  322. *  |                                                                 |
  323. *  +-----------------------------------------------------------------+
  324.  
  325. *  ===================================================
  326. *  PARAMETRES DE L'ALGORITHME D'INTEGRATION TEMPORELLE
  327. *  ===================================================
  328.  
  329.    ALPHA = ALPHAF - ALPHAM ;
  330.    GAMMA = 0.5 + ALPHA ;
  331.    BETA  = 0.25 * ((1.+ALPHA)**2) ;
  332.    1MAF  = 1. - ALPHAF ;
  333.    1MAM  = 1. - ALPHAM ;
  334.    GAMBET = GAMMA / BETA ;
  335.    SI (IECHO >EG 0) ;
  336.       MESS (CHAI 'Calcul avec GAMMA=' GAMMA ' BETA=' BETA
  337.                  ' (GAMMA/BETA=' GAMBET ')') ;
  338.    FINS ;
  339.  
  340.  
  341. *  ====================
  342. *  PARAMETRES TEMPORELS
  343. *  ====================
  344.  
  345.    IPAS = 1 ;
  346.    T0   = EXTR TPROG IPAS ;
  347.    DT0  = 0. ;
  348.    CHACHA = CHAI 'PAS #' IPAS ' TEMPS' T0 ;
  349.  
  350.  
  351. *  ===================
  352. *  SYSTEME D'EQUATIONS
  353. *  ===================
  354.  
  355.    OPER = VIDE 'RIGIDITE'/'RIGIDITE' ;
  356.    F0 = TIRE CHAR1 T0 ;
  357.  
  358.  
  359. *  =====================================================
  360. *  ACCELERATION INITIALE A0 (NOUVEAUX SCHEMAS TEMPORELS)
  361. *  =====================================================
  362.  
  363.    SI (EXIS ETAB 'ACCE') ;
  364.       A0 = ETAB.'ACCE' ;
  365.    SINON ;
  366.       K0 = K ;
  367.       F2ND = (F0 ENLE 'FLX') - (K * U0) ;
  368.       SI KAMOR ;
  369.          F2ND = F2ND - (C * V0) ;
  370.       FINS ;
  371.       SI KCHARM ;
  372.          TCHAR = CHARMECA WTAB ;
  373.          KCHAF = EXIS TCHAR 'ADDI_SECOND' ; 
  374.          SI KCHAF ;
  375.             F2ND = F2ND + TCHAR.'ADDI_SECOND' ;
  376.          FINS ;
  377.          KCHAK = EXIS TCHAR 'ADDI_MATRICE' ;
  378.          SI KCHAK ;
  379.             K0 = K0 ET TCHAR.'ADDI_MATRICE' ;
  380.          FINS ;
  381.      FINS ;
  382.      FREAC0 = REAC K0 U0 ;
  383.      F2ND   = F2ND + FREAC0 ;
  384. *    (1.*M) POUR EVITER PB HORODATAGE SI REUTILISATION DE M PLUS TARD
  385.      A0 = RESO (1.*M) F2ND ;
  386.    FINS ;
  387.  
  388.  
  389. *  ===================
  390. *  TABLES DE RESULTATS
  391. *  ===================
  392.  
  393.    INDVTK = 0 ;
  394.    INDCSV = 0 ;
  395.    INDGIB = 0 ;
  396.    NFANTO = 0 ;
  397.    EPOCH0 = DATE 'EPOCH';
  398.  
  399. *  SAUVEGARDE VTK
  400. *  --------------
  401.    SI (IVTK > 0) ;
  402.       SI (IVTK EGA 1) ;
  403.          KVTK = VRAI ;
  404.       SINON ;
  405.          KVTK = DANS ETAB.'TEMPS_SAUVES_VTK' T0 ;
  406.       FINS ;
  407.  
  408.       SI KVTK ;
  409.          INDVTK = INDVTK + 1 ;
  410.          OPTI 'SORT' FICVTK ;
  411.          SORT 'VTK' MAIVTK U0 'DEPL' V0 'VITE' 'TEMP' T0 ;
  412.          CHACHA = CHAI CHACHA '  => VTK #' INDVTK ;
  413.       FINS ;
  414.    FINS ;
  415.  
  416. *  SAUVEGARDE CSV
  417. *  --------------
  418.    SI (ICSV > 0) ;
  419.  
  420. *     Initialisation de la table TABCSV
  421.       TABCSV = TABL ;
  422.       ETAB.'RESULTATS_CSV' = TABCSV ;
  423.       TABCSV.'TEMPS' = PROG ;
  424.       LCOCSV = MOTS ;
  425.       LPOCSV = LECT ;
  426.  
  427. *     Titres de colonnes et initialisation des LISTREEL
  428.       CHPO1 = REDU (CHPO 'UNIFORME' 0. K) MAICSV ;
  429.       REPE I1 NELCSV ;
  430.          POIN1 = MAICSV POIN &I1 ;
  431.          COMP1 = EXTR (REDU CHPO1 POIN1) 'COMP' ;
  432.          SI KLISCO ;
  433.             COMP1 = COMP1 SAUF (COMP1 SAUF LISCO) ;
  434.          FINS ;
  435.          NCO1 = DIME COMP1 ;
  436.          SI (NCO1 EGA 0) ;
  437.             ITER I1 ;
  438.          FINS ;
  439.          LPOCSV = LPOCSV ET (LECT NCO1*&I1) ;
  440.          REPE I2 NCO1 ;
  441.             CHA1 = EXTR COMP1 &I2 ;
  442.             REPE I3 4 ;
  443.                J3 = 5 - &I3 ;
  444.                SI (NEG (EXTR CHA1 J3 J3) ' ') ;
  445.                   CHA1 = EXTR CHA1 1 J3 ;
  446.                   QUIT I3 ;
  447.                FINS ;
  448.             FIN I3 ;
  449.             LCOCSV = LCOCSV ET CHA1 ;
  450.             NOM1 = CHAI 'DEPL_' CHA1 '_' &I1 ;
  451.             NOM2 = CHAI 'VITE_' CHA1 '_' &I1 ;
  452.             TABCSV.NOM1 = PROG ;
  453.             TABCSV.NOM2 = PROG ;
  454.          FIN I2 ;
  455.       FIN I1 ;
  456.       _TABCSV = INDE TABCSV ;
  457.       NBCSV = ((DIME TABCSV) - 1) / 2 ;
  458.  
  459. *     Garder l'instant initial ?
  460.       SI (ICSV EGA 1) ;
  461.          KCSV = VRAI ;
  462.       SINON ;
  463.          KCSV = DANS ETAB.'TEMPS_SAUVES_CSV' T0 ;
  464.       FINS ;
  465.  
  466.       SI KCSV ;
  467.          INDCSV = INDCSV + 1 ;
  468.          TABCSV.'TEMPS' = TABCSV.'TEMPS' ET T0 ;
  469.          REPE I1 NBCSV ;
  470.             COMP1 = EXTR LCOCSV &I1 ;
  471.             POIN1 = MAICSV POIN (EXTR LPOCSV &I1) ;
  472.             X1 = EXTR U0 'VALE' COMP1 POIN1 'NOID' ;
  473.             X2 = EXTR V0 'VALE' COMP1 POIN1 'NOID' ;
  474.             K1 = 2*&I1 ;         
  475.             K2 = K1 + 1 ;
  476.             MOT1 = MOT _TABCSV.K1 ;
  477.             MOT2 = MOT _TABCSV.K2 ;
  478.             TABCSV.MOT1 = TABCSV.MOT1 ET X1 ;
  479.             TABCSV.MOT2 = TABCSV.MOT2 ET X2 ;
  480.          FIN I1 ;
  481.          CHACHA = CHAI CHACHA '  => CSV #' INDCSV ;
  482.       FINS ;
  483.    FINS ;
  484.  
  485. *  SAUVEGARDE GIBIANE
  486. *  ------------------
  487.    SI (IGIB > 0) ;
  488.  
  489. *     Initialisation de la table STAB (indicee par INDGIB)
  490.       STAB = TABL ;
  491.       ETAB.'RESULTATS' = STAB ;
  492.  
  493. *     Garder l'instant initial ?
  494.       SI (EGA IGIB 1) ;
  495.          KGIB = VRAI ;
  496.       SINON ;
  497.          KGIB = DANS ETAB.'TEMPS_SAUVES' T0 ;
  498.       FINS ;
  499.  
  500.       SI (KGIB OU (EGA IPAS NPAS)) ;
  501. *        (on force la sauvegarde si dernier pas)
  502.          INDGIB = INDGIB + 1 ;
  503.          STN = TABL ;
  504.          STAB.INDGIB = STN ;
  505.          STN.'TEMP' = T0 ;
  506.          SI KMAIGIB ;
  507.             STN.'DEPL' = REDU U0 MAIGIB ;
  508.             STN.'VITE' = REDU V0 MAIGIB ;
  509.          SINON ;
  510.             STN.'DEPL' = U0 ;
  511.             STN.'VITE' = V0 ;
  512.          FINS ;
  513.          CHACHA = CHAI CHACHA '  => GIBI #' INDGIB ;
  514.       FINS ;
  515.    FINS ;
  516.  
  517.  
  518. *  =======================
  519. *  TABLE DE CALCUL INTERNE
  520. *  =======================
  521.  
  522.    WTAB = TABL ;
  523.  
  524. *  On transmet a WTAB la table de PARAmetres utilisateurs si elle existe  
  525.    SI (EXIS ETAB 'PARA') ;
  526.       WTAB.'PARA' = ETAB.'PARA' ;
  527.    FINS ;
  528.    WTAB.'TEMP' = T0 ;
  529.    WTAB.'DEPL' = U0 ;
  530.    WTAB.'VITE' = V0 ;
  531.  
  532.  
  533.  
  534.  
  535. *  +-----------------------------------------------------------------+
  536. *  |                                                                 |
  537. *  |      B O U C L E   S U R   L E S   P A S   D E   T E M P S      |
  538. *  |                                                                 |
  539. *  +-----------------------------------------------------------------+
  540.  
  541.  
  542. *  AFFICHAGE DU MESSAGE D'INFORMATION POUR L'ETAT INITIAL
  543.    MESS CHACHA ;
  544.  
  545.    REPE BPAS (NPAS - 1) ;
  546.       IPAS = IPAS + 1 ;
  547.  
  548. *     ===========================
  549. *     PREPARATION DU PAS DE TEMPS
  550. *     ===========================
  551.  
  552. *     PARAMETRES TEMPORELS  
  553. *     --------------------
  554.       T1  = EXTR TPROG IPAS ;
  555.       DT  = T1 - T0 ;
  556.       DT2 = DT ** 2 ;
  557.       CHACHA = CHAI 'PAS #' IPAS ' TEMPS' T1 ;
  558.  
  559. *     QUELQUES CONSTANTES
  560. *     -------------------
  561.       XM = 1. / (BETA*DT2) ;
  562.       XC = GAMMA / (BETA*DT) ;
  563.  
  564. *     WTAB AVEC U_N ET T_N+1
  565.       WTAB.'TEMP' = T1 ;
  566.       WTAB.'DEPL' = U0 ;
  567.       WTAB.'VITE' = V0 ;
  568. *       WTAB.'ACCE' = A0 ;
  569.  
  570. *     PROCEDURE "UTILISATEUR" (CHARMECA)
  571. *     ----------------------------------
  572.       KCHAF = FAUX ;
  573.       KCHAK = FAUX ;
  574.       KCHAFNL = FAUX ;
  575.       KCHAKNL = FAUX ;
  576.       KCHACNL = FAUX ;
  577.       SI KCHARM ;
  578.          TCHAR = CHARMECA WTAB ;
  579.  
  580. *        Indices destines a l'ajout de relations (mult. de Lagrange)
  581. *        entre deplacements
  582.          KCHAF = EXIS TCHAR 'ADDI_SECOND' ; 
  583.          SI KCHAF ;
  584.             FCHAR1 = TCHAR.'ADDI_SECOND' ;
  585.          FINS ;
  586.          KCHAK = EXIS TCHAR 'ADDI_MATRICE' ;
  587.          SI KCHAK ;
  588.             KCHAR1 = TCHAR.'ADDI_MATRICE' ;
  589.          FINS ;
  590.  
  591. *        Indices destines aux autres comportements non lineaires
  592.          KCHAFNL = EXIS TCHAR 'ADDI_FNL' ;
  593.          SI KCHAFNL ;
  594.             FNL0 = TCHAR.'ADDI_FNL' ;
  595.          FINS ;
  596.          KCHAKNL = EXIS TCHAR 'ADDI_KNL' ;
  597.          SI KCHAKNL ;
  598.             KNL1 = TCHAR.'ADDI_KNL' ;
  599.          FINS ;
  600.          KCHACNL = EXIS TCHAR 'ADDI_CNL' ;
  601.          SI KCHACNL ;
  602.             CNL1 = TCHAR.'ADDI_CNL' ;
  603.          FINS ;
  604.       FINS ;
  605.  
  606. *     CALCUL DU SECOND MEMBRE
  607. *     -----------------------
  608.       F1 = TIRE CHAR1 T1 ;
  609. *     (aucune difference entre ADDI_SECOND et ADDI_FNL si ce n'est
  610. *      l'affichage dans CHACHA)
  611.       SI KCHAF ; 
  612.         F1 = F1 + FCHAR1 ; 
  613.         CHACHA = CHAI CHACHA ' +FLX' ; 
  614.       FINS ;
  615.       SI KCHAFNL ;
  616.         F1 = F1 + FNL0 ; 
  617.         CHACHA = CHAI CHACHA ' +FNL' ; 
  618.       FINS ; 
  619.       FEXT1 = (1MAF*(F1 ENLE 'FLX')) + (ALPHAF*(F0 ENLE 'FLX'))
  620.             + (EXCO 'FLX' F1 'NOID' 'FLX') ;
  621.  
  622. *     ANCIEN SCHEMA => COMPATIBLE AVEC VITEUNIL
  623.       SI KVITUN ; 
  624. *     SI (EGA IALGO 0) ;
  625.          K1 = 'EXTRAI' K 'RIGI' 'NOMU' ;                                        
  626.          K2 = 'EXTRAI' K 'RIGI' 'MULT' ;
  627.          MAICO = EXTR K 'MAIL' 'UNIL' ;
  628.          A1PRED = COLI V0 (1MAM / (BETA*DT)) ; 
  629.          FA = M * A1PRED ;
  630.          FB = K1 * (2. * (U0 - (EXCO 'LX' U0 'NOID' 'LX'))) ; 
  631.          SI ((NBEL MAICO) > 0) ;
  632.             U0PRED = U0 - (REDU U0 MAICO) ;
  633.          SINON ;
  634.             U0PRED = U0 ;   
  635.          FINS ;
  636.          FG = K2 * U0PRED ;
  637.          F2ND = FEXT1 + (F0 ENLE 'FLX') - FB + FA - FG ;
  638.  
  639. *     SYNTAXE ADAPTEE AU HHT ET ALPHA-GEN (MAIS NON COMPATIBLE AVEC
  640. *     RELATIONS UNILATERALES ET VITEUNIL)
  641.       SINON ;
  642.          A1PRED = COLI V0 ( 1MAM / (BETA*DT) )  
  643.                        A0 ( (1MAM * 0.5 / BETA) - 1. ) ;
  644.          V1PRED = COLI V0 ( (1MAF * GAMBET) - 1. )  
  645.                        A0 ( 1MAF * DT * ((0.5*GAMBET) - 1.) ) ;
  646.          U0PRED = ENLE U0 'LX' ;
  647.          F2ND   = FEXT1 - (K * (U0PRED))  + (M * A1PRED) ;
  648.          SI KAMOR ;
  649.             F2ND = F2ND + (C *  V1PRED) ;
  650.          FINS ;
  651.          SI KCHACNL ;
  652.             F2ND = F2ND + (CNL1 * (V1PRED + V0)) ;
  653.          FINS ;
  654. *        2nd membre complet = Fext_n+1 - Fint_n + Famor_n + M*(4/dt*V_n+A_n) 
  655.  
  656.       FINS ;
  657.  
  658. *     ON RECALCULE L'OPERATEUR...
  659. *     ...SI LE PAS DE TEMPS CHANGE OU SI UNE MATRICE CHANGE
  660. *     -----------------------------------------------------
  661.       DTTOL = 1.E-6 * (MAXI DT DT0);
  662.       SI ((DT NEG DT0 DTTOL) OU KCHAK) ;
  663.          SI (KMENAG) ;
  664.             DETR OPER ;
  665.          FINS ;
  666.          OPER = (1MAF * K) ET ((1MAM * XM) * M) ;
  667.          SI KAMOR ;
  668.             OPER = OPER ET ((1MAF * XC) * C) ;
  669.          FINS ;
  670.          SI KCHAK ; 
  671.             OPER = OPER ET KCHAR1 ; 
  672.             CHACHA = CHAI CHACHA ' +[LX]' ; 
  673.          FINS ;
  674.          SI KCHAKNL ; 
  675.             OPER = OPER ET (1MAF * KNL1) ;
  676.             CHACHA = CHAI CHACHA ' +[KNL]' ; 
  677.          FINS ;
  678.          SI KCHACNL ; 
  679.             OPER = OPER ET ((1MAF * XC) * CNL1) ; 
  680.             CHACHA = CHAI CHACHA ' +[CNL]' ; 
  681.          FINS ;
  682.       FINS ;
  683.  
  684.  
  685. *     ======================
  686. *     CALCUL DU PAS DE TEMPS
  687. *     ======================
  688.  
  689. *     RESOLUTION DU SYSTEME LINEAIRE
  690.       DU = RESO OPER F2ND ;
  691.  
  692. *     MISE A JOUR DES VARIABLES AU TEMPS N+1
  693. *     U1 = U0 + DU ;
  694. *     A priori inutile d'enlever les LX, mais histoire d'etre sur...
  695.       U1 = (ENLE U0 'LX') + DU ;
  696.       DU1 = ENLE DU 'LX' ;
  697.       V1 = COLI DU1 XC V0 (1. - GAMBET)   A0 (DT*(1. - (0.5*GAMBET))) ;
  698.       A1 = COLI DU1 XM V0 (-1./(BETA*DT)) A0 (1. - (0.5/BETA)) ;
  699.  
  700. * *     Calcul du Residu (a n+1 ou avec les alphas pour NL?)
  701. * *     rem : inutile en lineaire
  702. * *           K*U1 contient deja les reactions
  703. *       Res1 = F1 - (K * U1) - (M * A1);
  704. *       MaxRes1 = MAXI Res1 'ABS' 'SANS' (mots 'FLX');
  705. *       chacha = chai chacha ' max(R)=' MaxRes1;
  706. *       
  707. * *     Calcul de l'energie (a n+1)
  708. *       Wmeca1 = 0.5 * (XTMX K (U1 enle LX) + (XTMX M V1)) ;
  709. *       chacha = chai chacha ' Wmeca=' FORMAT '(F8.1)' Wmeca1 ;
  710.  
  711. *     CORRECTION DES VITESSES (SEULEMENT POUR NEWMARK ACC. MOY.)
  712.       SI KVITUN ; 
  713.  
  714. *        Eventuelle mise a jour de l'indicateur de presence de relations
  715. *        unilaterales
  716.          SI (KCHAK OU (IPAS <EG 2)) ; 
  717.             MTA = EXTR OPER 'CONTACT' ;
  718.             KCONTA = ((DIME MTA) > 0) ;
  719.          FINS ;
  720.  
  721.          SI KCONTA ;
  722.             TBID = TABL ;                                                       
  723.             V1 = VITEUNIL OPER M V1 DU U0 DT F2ND TBID ; 
  724.             SI (EXIS TBID 'RATE_VITEUNIL') ;
  725.                MESS 'Erreur fatale dans VITEUNIL' ;
  726.                ERRE 5 ;
  727.             FINS ;
  728.          FINS ;
  729.       FINS ;
  730.  
  731.  
  732. *     ==========================
  733. *     SAUVEGARDE DU PAS DE TEMPS
  734. *     ==========================
  735.  
  736. *     SAUVEGARDE VTK
  737. *     --------------
  738.       SI (IVTK > 0) ;
  739.          SI (IVTK EGA 1) ;
  740.             KVTK = MULT (IPAS - 1) IPASVTK ;
  741.          SINON ;
  742.             KVTK = DANS ETAB.'TEMPS_SAUVES_VTK' T1 ;
  743.          FINS ;
  744.  
  745.          SI KVTK ;
  746.             INDVTK = INDVTK + 1 ;
  747.             OPTI 'SORT' FICVTK ;
  748.             SORT 'VTK' MAIVTK U1 'DEPL' V1 'VITE' 'TEMP' T1 ;
  749.             CHACHA = CHAI CHACHA '  => VTK #' INDVTK ;
  750.          FINS ;
  751.       FINS ;
  752.  
  753. *     SAUVEGARDE CSV
  754. *     --------------
  755.       SI (ICSV > 0) ;
  756.          SI (ICSV EGA 1) ;
  757.             KCSV = MULT (IPAS - 1) IPASCSV ;
  758.          SINON ;
  759.             KCSV = DANS ETAB.'TEMPS_SAUVES_CSV' T1 ;
  760.          FINS ;
  761.  
  762.          SI KCSV ;
  763.             INDCSV = INDCSV + 1 ;
  764.             TABCSV.'TEMPS' = TABCSV.'TEMPS' ET T1 ;
  765.             REPE I1 NBCSV ;
  766.                COMP1 = EXTR LCOCSV &I1 ;
  767.                POIN1 = MAICSV POIN (EXTR LPOCSV &I1) ;
  768.                X1 = EXTR U1 'VALE' COMP1 POIN1 'NOID' ;
  769.                X2 = EXTR V1 'VALE' COMP1 POIN1 'NOID' ;
  770.                K1 = 2*&I1 ;         
  771.                K2 = K1 + 1 ;
  772.                MOT1 = MOT _TABCSV.K1 ;
  773.                MOT2 = MOT _TABCSV.K2 ;
  774.                TABCSV.MOT1 = TABCSV.MOT1 ET X1 ;
  775.                TABCSV.MOT2 = TABCSV.MOT2 ET X2 ;
  776.             FIN I1 ;
  777.             CHACHA = CHAI CHACHA '  => CSV #' INDCSV ;
  778.          FINS ;
  779.  
  780.          SI ((EGA IPAS NPAS) ET KFICCSV) ;
  781.             OPTI 'SORT' FICCSV ;
  782.             SORT 'EXCE' TABCSV ;
  783.          FINS ;
  784.       FINS ;
  785.  
  786. *     SAUVEGARDE GIBIANE
  787. *     ------------------
  788.       SI (IGIB > 0) ;
  789.          SI (IGIB EGA 1) ;
  790.             KGIB = MULT (IPAS - 1) IPASGIB ;
  791.          SINON ;
  792.             KGIB = DANS ETAB.'TEMPS_SAUVES' T1 ;
  793.          FINS ;
  794.  
  795.          EPOCH1='DATE' 'EPOCH';
  796.          SI ((IPAS 'EGA' NPAS) 'OU' KGIB);
  797. *           (on force la sauvegarde si dernier pas)
  798.             INDGIB = INDGIB + 1 ;
  799.             STN = TABL ;
  800.             STAB.INDGIB = STN ;
  801.             STN.'TEMP' = T1 ;
  802.             SI KMAIGIB ;
  803.                STN.'DEPL' = REDU U1 MAIGIB ;
  804.                STN.'VITE' = REDU V1 MAIGIB ;
  805.             SINON ;
  806.                STN.'DEPL' = U1 ;
  807.                STN.'VITE' = V1 ;
  808.             FINS ;
  809.             CHACHA = CHAI CHACHA '  => GIBI #' INDGIB ;
  810.  
  811. *           Ecriture incrementale sur disque dur
  812.             SI (KINCR 'ET'
  813.                 (((EPOCH1 - EPOCH0) '>EG' 300.D0) 'OU'
  814.                  ((EPOCH0 EGA 0.D0) 'ET' (EPOCH1 EGA 0.D0)) 'OU'
  815.                  (IPAS 'EGA' NPAS))) ;
  816. *               (Il est necessaire de sauver aussi OPER : RIGIDITE)
  817.                 ETAB.'OPER' = OPER ;
  818.                 SAUV 'MUET' ETAB   ;
  819.                 EPOCH0='DATE' 'EPOCH';
  820.  
  821. *               Fantomisation des resultats sauves precedemment
  822. *               (attention : on a encore besoin de ceux en cours pour le
  823. *                            prochain pas)
  824.                 SI (KECON) ;
  825.                   NFANTOL=(INDGIB - 1) - NFANTO ;
  826.                   SI(NFANTOL > 0);
  827.                     REPE SURI NFANTOL;
  828.                       JJ=NFANTO + &SURI;
  829.                       SI(JJ <EG 0);
  830.                         ITER SURI;
  831.                       FINS;
  832.                       'FANT' (STAB.JJ) 'DEPL' ;
  833.                       'FANT' (STAB.JJ) 'VITE' ;
  834.                     FIN SURI;
  835.                     NFANTO=NFANTO+NFANTOL;
  836.                   FINS ;
  837.                 FINS ;
  838.             FINS ;
  839.          FINS ;
  840.       FINS ;
  841.  
  842.  
  843. *     ===============================
  844. *     PASSAGE AU PAS DE TEMPS SUIVANT
  845. *     ===============================
  846.  
  847. *     AFFICHAGE DU MESSAGE D'INFORMATION DU PAS DE TEMPS COURANT
  848.       MESS CHACHA ;
  849.  
  850.       T0 = T1 ;
  851.       DT0 = DT ;
  852. *       U0 = U1 ENLE 'LX' ;
  853.       U0 = U1 ;
  854.       V0 = V1 ;
  855.       DETR A0 ;
  856.       A0 = A1 ;
  857.       DETR F0 ;
  858.       F0 = F1 ;
  859.       SI (KMENAG) ;
  860.          DETR A1PRED ;
  861.          DETR DU ;
  862.       FINS ;
  863.  
  864.    FIN BPAS ;
  865.  
  866.  
  867. *  Sauvegarde de la derniere acceleration (pour reprise eventuelle)
  868.    STN.'ACCE' = A1 ; 
  869.  
  870. FINP STAB ;
  871.  
  872.  
  873.  
  874.  
  875.  

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