Télécharger la plaquette Cast3M


fr en

Télécharger xfem04.dgibi

Retour à la liste

Numérotation des lignes : 

  1. optio norm auto;
  2. * fichier : xfem04.dgibi
  3. ************************************************************************
  4. ************************************************************************
  5. ************************************************************************
  6. *
  7. * Mots-clés : Mécanique de la rupture, XFEM, contact, frottement, LATIN
  8. *
  9. * Etude de reponse d'une plaque en compression avec fissure inclinee
  10. * via une formulation mixte a trois champs. Le contact frottant et la
  11. * resolution par la LATIN sont pris en compte par des procedures perso. 
  12. * test des operateurs : RELA ACCRO FAIBLE, G_THETA
  13. * creation : 2013-02-13, Benoit Trolle (these SNCF-Lamcos)
  14. * integration dans Cast3m : 2013-12-17, Benoit Prabel (#)
  15. *
  16. ************************************************************************
  17.  
  18. OPTION DIME 2 MODE PLAN DEFO ELEM QUA4;
  19.  
  20.  
  21. * Masquer lignes excutées
  22. *  'OPTION' echo 0;
  23.  
  24. ************************************************************************
  25. *                           PROCEDURES INTERNES                        *
  26. ************************************************************************
  27.  
  28. *-----------------------------------------------------------------------
  29. * PLOTITER : Affichage des champs locaux au cours des itérations de la LATIN
  30. *----------------------------------------------------------------------
  31. 'DEBPROC' PLOTITER Tplus*CHPOINT Tmoins*CHPOINT Wplus*CHPOINT
  32.                      Wmoins*CHPOINT meshint*MAILLAGE;
  33.    fac2 = 20;
  34.    fac1 = 1E-12;
  35.    VTplus = vect Tplus forc rouge;
  36.    VTmoins = vect Tmoins forc  vert;
  37.    def1 = defo wplus meshint   VTplus rouge;                          
  38.    def2 = defo wmoins meshint  VTmoins vert;
  39.    def3 = def1 'ET' def2;
  40.    'TRACER' def3 'NCLK';
  41. 'FINPROC';
  42.  
  43.  
  44. *----------------------------------------------------------------------
  45. * LATIN  : contient les étapes du solveur non linéaires LATIN 
  46. *
  47. * ENTREE :
  48. *  - crack : maillage de l'interface
  49. *  - RHS0  : Chargement mécanique
  50. *  - K1tot : matrice totale sans Kuwt (avec C.L.)
  51. *  - Kuwt  : matrice issu de RELA ACCRO FAIBLE
  52. *  - Ktot1 : Kuwt et K1tot
  53. *  - mod1tot : modele de tout le massif
  54. *  - WPglo, WMglo, WPglo, WMglo : w+,w-,t+,t- solution du pas précédent
  55. *  - f_crack : Coefficient de frottement entre les lèvres de la fissure
  56. *  - ItMax : Nombre d'itération maximum
  57. *  - eps1 : Précision demandée
  58. *
  59. * SORTIE :  
  60. *  - (U, W, T) 
  61. *  - erreur à chaque itération
  62. *----------------------------------------------------------------------
  63. 'DEBPROC' LATIN Tdisplay*TABLE crack*MAILLAGE K1tot*RIGIDITE
  64.  mod1tot*MMODEL RHS0*CHPOINT IdPasT*ENTIER Ktot1*RIGIDITE Kuwt*RIGIDITE
  65.  cht1*CHPOINT chn*CHPOINT ItMax*ENTIER eps1*FLOTTANT k*FLOTTANT
  66.  f_crack*FLOTTANT mai1tot*MAILLAGE cht2/CHPOINT TPglo/CHPOINT
  67.  TMglo/CHPOINT WPglo/CHPOINT WMglo/CHPOINT;
  68.  
  69.  
  70. * initialisation avec second menbre incomplet uniquement pour le premier pas de de temps
  71. * pour les suivants on initialise avec les champs locaux du pas précédents
  72.  
  73. *'SI' ((IdPasT 'EGA' 1) 'OU' ouv1);
  74. 'SI'  (IdPasT 'EGA' 1);
  75.  
  76. *-----------------------------< Résolution initiale: étape globale
  77.    ma1 = 'EXTRAIRE' mod1 'MAIL';
  78.    cl2 =  ('BLOQUE' 'AX'ma1) 'ET' ('BLOQUE' 'B1X'ma1)  'ET' 
  79.    ('BLOQUE' 'C1X'ma1) 'ET' ('BLOQUE' 'D1X'ma1) 'ET' 
  80.    ('BLOQUE' 'E1X'ma1) 'ET' ('BLOQUE' 'AY'ma1)  'ET' 
  81.    ('BLOQUE' 'B1Y'ma1) 'ET' ('BLOQUE' 'C1Y'ma1) 'ET'  
  82.    ('BLOQUE' 'D1Y'ma1) 'ET' ('BLOQUE' 'E1Y'ma1) 'ET' 
  83.    ('BLOQUE' 'AZ'ma1)  'ET' ('BLOQUE' 'B1Z'ma1) 'ET' 
  84.    ('BLOQUE' 'C1Z'ma1) 'ET' ('BLOQUE' 'D1Z'ma1) 'ET' 
  85.    ('BLOQUE' 'E1Z'ma1) ;
  86.    ktot2 = ktot1 'ET' cl2;
  87. *   u1 = RESO Ktot2 RHS0;
  88.  
  89.    u1 = RESO Ktot1 RHS0;
  90.  
  91. *-----------------------------< extrait w et t fissure de u1
  92. * (dans leur formalisme global = moyenne + saut)
  93.    wsauv fsauv = EXTCRACK u1 crack dim1; 
  94.  
  95. *-----------------------------< recuperation champ + et -
  96.    Wmoins Wplus Tmoins Tplus = SPLIT crack fsauv wsauv;
  97.    'SI' (Tdisplay.iteration);
  98.          'TITRE' 'apres etape GLOBALE INITIALE';
  99.       PLOTITER Tplus Tmoins  Wplus Wmoins crack;
  100.    'FINSI';
  101. *-----------------------------< Projection dans la base locale xxx_b = xxx_before (local step)
  102.    'SI' (('VALEUR' 'DIME')  'EGA' 3);
  103.    TlocP_b TlocM_b WlocP_b WlocM_b = GLO2LOC cht1 chn Tplus Tmoins Wplus
  104.    Wmoins cht2;
  105.    'SINON' ;
  106.    TlocP_b TlocM_b WlocP_b WlocM_b = GLO2LOC cht1 chn Tplus Tmoins Wplus
  107.    Wmoins ;
  108.    'FINSI';
  109.  
  110.  
  111. *-----------------------------< Etape locale pour la première itération 
  112.    TlocP TlocM WlocP WlocM  = LOCSTEPI k f_crack crack TlocP_b TlocM_b
  113.                                           WlocP_b WlocM_b;
  114.  
  115. *-----------------------------< Projection des w+- et  t+- (i +1/2) dans le repère global
  116.    'SI' (('VALEUR' 'DIME')  'EGA' 3);
  117.       WMglo WPglo TMglo TPglo = LOC2GLO cht1 chn TlocP TlocM WlocP WlocM
  118.                                                      crack cht2;    
  119.     'SINON';
  120.       WMglo WPglo TMglo TPglo = LOC2GLO cht1 chn TlocP TlocM WlocP WlocM
  121.                                                      crack;
  122.     'FINSI';
  123.  
  124.  'FINSI';
  125. *---> initialisation itérations
  126.    it1 = 1 ;
  127.  
  128. * utlise solutions locales au pas de temps précédent pour initialiser
  129.  'SI' (IdPasT > 1);
  130.  
  131. * 'SI' (idpasT  < 6); itmax = 10; 'FINSI';
  132.  
  133. *-----------------------------< Projection dans la base locale xxx_b = xxx_before (local step) 
  134.    'SI' (('VALEUR' 'DIME')  'EGA' 3);
  135.       TlocP_b TlocM_b WlocP_b WlocM_b = GLO2LOC cht1 chn TPglo TMglo
  136.                                                    WPglo WMglo cht2;
  137.    'SINON' ;
  138.       TlocP_b TlocM_b WlocP_b WlocM_b = GLO2LOC cht1 chn TPglo TMglo
  139.                                                    WPglo WMglo ;
  140.    'FINSI'; 
  141.  
  142.    TlocP = TlocP_b;
  143.    TlocM = TlocM_b;
  144.    WlocP = WlocP_b;
  145.    WlocM = WlocM_b;
  146.  'FINSI';   
  147. * sauvegarde champ locaux après étape locale
  148.  ttlocp = table;  ttlocm = table;  twlocp = table;  twlocm = table;
  149.  ttlocp.it1 = tlocp_b; ttlocm.it1 = tlocm_b;
  150.  twlocp.it1 = wlocp_b; twlocm.it1 = wlocm_b;
  151.    'SI' (Tdisplay.iteration);
  152.       'TITRE' 'apres etape LOCALE';
  153.       PLOTITER TPglo TMglo  WPglo WMglo crack;
  154.    'FINSI'; 
  155. *-----------------------------< Recombine les champs
  156.  Wglo Tglo = LIPSRECO WMglo WPglo TMglo TPglo;
  157.  
  158. *-----------------------------< Mise à jour du second membre
  159.  rhs1 = UPDATRHS Kuwt Tglo Wglo crack it1 rhs0;
  160.  
  161. *-----------------------------< Etape globale
  162.  u1 = 'RESOUD' ktot1 rhs1;
  163.  
  164. *-----------------------------< extrait w et t fissure de u1
  165. * (dans leur formalisme global = moyenne + saut)
  166.  wfiss fcon = EXTCRACK u1 crack dim1;
  167.  
  168. *-----------------------------< recuperation champ + et -
  169.  Wmoins Wplus Tmoins Tplus = SPLIT Crack fcon wfiss;
  170.  'SI' (Tdisplay.iteration);
  171.        'TITRE' 'apres etape GLOBALE';
  172.    PLOTITER Tplus Tmoins  Wplus Wmoins crack;
  173.  'FINSI';
  174.  
  175.  
  176. *---< Projection dans la base locale
  177.  'SI' (('VALEUR' 'DIME')  'EGA' 3);
  178.    TlocP_n TlocM_n WlocP_n WlocM_n = GLO2LOC cht1 chn Tplus Tmoins
  179.    Wplus Wmoins cht2;
  180.  'SINON' ;
  181.    TlocP_n TlocM_n WlocP_n WlocM_n = GLO2LOC cht1 chn Tplus Tmoins
  182.    Wplus Wmoins ;
  183.  'FINSI';
  184.  
  185.  
  186.  'REPETER' blo7 (Itmax);
  187.    it1 = it1 '+' 1;
  188.  
  189.    TlocP_b = TlocP_n;
  190.    TlocM_b = TlocM_n; 
  191.    WlocP_b = WlocP_n;
  192.    WlocM_b = WlocM_n;
  193.  
  194. *-----------------------------< Etape locale (_a = after ; _b = before)
  195.    TlocP_a TlocM_a WlocP_a WlocM_a ouv1 = LOCSTEP1
  196.    k f_crack crack TlocP_b TlocM_b WlocP_b WlocM_b it1
  197.    (twlocp.(it1-1)) (twlocm.(it1-1));
  198.  
  199.    'SI' (ouv1 'ET' (it1 > 15));
  200.       'QUITTER' blo7;
  201.    'FINSI';   
  202.  
  203.  
  204. *-----------------------------< Projection des des w+- et  t+- (i +1/2) dans le repère global
  205.    'SI' (('VALEUR' 'DIME')  'EGA' 3);
  206.       WMglo WPglo TMglo TPglo = LOC2GLO cht1 chn TlocP_a TlocM_a WlocP_a
  207.                                        WlocM_a crack cht2;    
  208.    'SINON';
  209.       WMglo WPglo TMglo TPglo = LOC2GLO cht1 chn TlocP_a TlocM_a WlocP_a
  210.                                        WlocM_a crack;
  211.    'FINSI';
  212.  
  213.    'SI' (Tdisplay.iteration);
  214.          'TITRE' 'apres etape LOCALE';
  215.       PLOTITER TPglo TMglo  WPglo WMglo crack;
  216.    'FINSI';
  217. *-----------------------------< Recombine les champs
  218.    Wglo Tglo = LIPSRECO WMglo WPglo TMglo TPglo;
  219.  
  220. *-----------------------------< Mise à jour du second membre
  221.    rhs1 = UPDATRHS Kuwt Tglo Wglo crack it1 rhs0;
  222.  
  223. *-----------------------------< Etape globale
  224.    u1 = 'RESOUD' ktot1 rhs1;
  225.  
  226. *-----------------------------< extrait w et t fissure de u1
  227. * (dans leur formalisme global = moyenne + saut)
  228.    wfiss fcon = EXTCRACK u1 crack dim1;
  229.  
  230. *-----------------------------< recuperation champ + et -
  231.    Wmoins Wplus Tmoins Tplus = SPLIT crack fcon wfiss;
  232.    'SI' (Tdisplay.iteration);
  233.          'TITRE' 'apres etape GLOABLE';
  234.       PLOTITER Tplus Tmoins Wplus Wmoins crack;
  235.    'FINSI';
  236.  
  237. *-----------------------------< Projection dans la base locale
  238.    'SI' (('VALEUR' 'DIME')  'EGA' 3);
  239.       TlocP_n TlocM_n WlocP_n WlocM_n = GLO2LOC cht1 chn Tplus Tmoins
  240.                Wplus Wmoins cht2;
  241.    'SINON' ;
  242.       TlocP_n TlocM_n WlocP_n WlocM_n = GLO2LOC cht1 chn Tplus Tmoins
  243.                Wplus Wmoins ;
  244.    'FINSI';
  245.  
  246. * sauvegarge de l'histoire au cours des itérations 
  247. *   ttlocp.it1 = tlocp_n;
  248. *   ttlocm.it1 = tlocm_n;
  249. *   twlocp.it1 = wlocp_n;
  250. *   twlocm.it1 = wlocm_n;
  251.  
  252.    ttlocp.it1 = tlocp_a;
  253.    ttlocm.it1 = tlocm_a;
  254.    twlocp.it1 = wlocp_a;
  255.    twlocm.it1 = wlocm_a;
  256.  
  257.    'SI' (it1 'EGA' 2);
  258. *-----------------------------< Calcul de l'indicateur d'erreur INITIAL
  259.       TotalErr denN denT etaT etaN = ERRORINI TlocP_n TlocM_n WlocP_n
  260.       WlocM_n TlocP TlocM WlocP WlocM k;
  261.       lerror = 'PROG' TotalErr;letaT = 'PROG'etaT;letaN = 'PROG' etaN;
  262.    'SINON';
  263. *-----------------------------< Calcul de l'indicateur d'erreur (dénominateur fixé)
  264.       TotalErr etaT etaN = GETERROR TlocP_n TlocM_n WlocP_n WlocM_n
  265.       TlocP_a TlocM_a WlocP_a WlocM_a denN denT k;
  266.       lerror = lerror'ET'TotalErr;letaT=letaT'ET'etaT;
  267.       letaN=letaN'ET'etaN;
  268.    'FINSI';
  269. *   'MESSAGE' '';
  270. *   'MESSAGE' 'denN =' denN;
  271. *   'MESSAGE' 'denT =' denT;
  272. *   'MESSAGE' '';
  273.  
  274.    SI (TotalErr < eps1) ;
  275.       MESS 'La precision demandee a ete atteinte';
  276.       QUITTER blo7 ;
  277.    FINSI ;
  278.  
  279.    SI (TotalErr > 1000000) ;
  280.       'MESSAGE' ' ';
  281.       'MESSAGE' ' ';
  282.       'MESSAGE' ' ';
  283.       MESS '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!';
  284.       MESS '!!!                                                 !!!';
  285.       MESS '!!!     LE CALCUL DIVERGE !! -> ARRET DU CALCUL     !!!';
  286.       MESS '!!!                                                 !!!';
  287.       MESS '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!';
  288.       'MESSAGE' ' ';
  289.       'MESSAGE' ' ';
  290.       'MESSAGE' ' ';
  291.       QUITTER blo7 ;
  292.    FINSI ;
  293.  
  294.    MESS ' ';
  295.    MESS ' ';
  296.    MESS '-----------> iteration :' it1 ;
  297.    MESS 'Erreur :' totalerr;
  298.  
  299.  'FIN' blo7;
  300. * fin de la boucle sur les itération du solveur non-linéaire
  301.  
  302.  'MESSAGE' 'fin du cas fermeture';
  303. * 'FINSI';
  304. * fin du cas fermeture
  305.  
  306. * cas fissure complètement ouverte -> résoluion directe -> w = u sur la fissure et t = 0
  307.  'SI' ouv1;
  308.    'MESSAGE'  ' ' ;
  309.    'MESSAGE' 'Fisssure COMPLETEMENT ouverte -> résolution DIRECTE';
  310.    'MESSAGE'  ' ' ;
  311.    u2 = 'RESOUD' k1tot RHS0;
  312.    WMglo WPglo = XFEM 'FISS' crack u2 mod1tot;
  313.    TPglo = 0.*Tplus;   TMglo = 0.*Tplus;lerror = 'PROG' 0.;   
  314.  'SINON';
  315.     u2 = 'REDU' u1 mai1tot;
  316.  'FINSI';
  317.  
  318.  
  319. 'FINPROC' u2 TPglo TMglo WPglo WMglo lerror;
  320.  
  321.  
  322. *----------------------------------------------------------------------
  323. *
  324. * SAVFIELD : sauvegarde les champs dans les tables
  325. *
  326. *----------------------------------------------------------------------
  327.  
  328. 'DEBPROC' SAVFIELD u1*CHPOINT tplus*CHPOINT tmoins*CHPOINT wplus*CHPOINT
  329.             wmoins*CHPOINT it1*ENTIER;
  330.  
  331. 'SI' (it1 'EGA' 0);
  332.   ttplus = 'TABLE';
  333.   ttmoins = 'TABLE';
  334.   twplus = 'TABLE';
  335.   twmoins = 'TABLE';
  336.   tu = 'TABLE';
  337. 'FINSI';
  338.   ttplus.it1 = tplus;
  339.   ttmoins.it1 = tmoins;
  340.   twplus.it1 = wplus;
  341.   twmoins.it1 = wmoins;
  342.   tu.it1 = u1;
  343. 'FINPROC' ttplus ttmoins twplus twmoins tu;
  344.  
  345.  
  346. *----------------------------------------------------------------------
  347. *
  348. * EXTCRACK : extrait les inconnus relatives à la fissure de u1
  349. *
  350. *----------------------------------------------------------------------
  351.  
  352. 'DEBPROC' EXTCRACK u1*CHPOINT meshint*MAILLAGE dim1*ENTIER;
  353.  
  354. 'SI' (dim1 'EGA' 2); 
  355.  lmot1 = 'MOTS' 'UX  ' 'UY  ' 'AX  ' 'AY  ';
  356.  lmot2 = 'MOTS' 'FX  ' 'FY  ' 'FAX ' 'FAY ';
  357. 'SINON';
  358.  lmot1 = 'MOTS' 'UX  ' 'UY  ' 'UZ  ' 'AX  ' 'AY  ' 'AZ  ' ;
  359.  lmot2 = 'MOTS' 'FX  ' 'FY  ' 'FZ  ' 'FAX ' 'FAY ' 'FAZ ' ;
  360. 'FINSI';
  361.  
  362.  wsauv = 'ENLEVER' ('REDU' meshint u1) lmot2;
  363.  fsauv = 'ENLEVER' ('REDU' meshint u1) lmot1;    
  364.  
  365. 'FINPROC' wsauv fsauv;
  366.  
  367.  
  368. *----------------------------------------------------------------------
  369. *
  370. * SPLIT : Différencie les champs locaux sur les lèvres de la fissure
  371. *              à partir des inconnus XFEM issus de l'étape globale
  372. *
  373. *----------------------------------------------------------------------
  374.  
  375. 'DEBPROC' SPLIT MeshInt*MAILLAGE fint1*CHPOINT uint1*CHPOINT;
  376.  
  377.  'SI' (('VALEUR' 'DIME') 'EGA' 2);
  378.    lmot1 = 'MOTS' 'UX  ' 'UY  ';
  379.    lmot2 = 'MOTS' 'AX  ' 'AY  ';
  380.    lmot3 = 'MOTS' 'FX  ' 'FY  ';
  381.    lmot4 = 'MOTS' 'FAX ' 'FAY ';
  382.  'SINON';
  383.    lmot1 = 'MOTS' 'UX  ' 'UY  '  'UZ  ';
  384.    lmot2 = 'MOTS' 'AX  ' 'AY  '  'AZ  ';
  385.    lmot3 = 'MOTS' 'FX  ' 'FY  '  'FZ  ';
  386.    lmot4 = 'MOTS' 'FAX ' 'FAY '  'FAZ ';
  387.  'FINSI'; 
  388.  
  389. *-----------------------------< Déplacement
  390.  umoy = 'EXCO' lmot1 uint1;
  391.  uH = 'EXCO' lmot2 uint1;
  392.  uH = 'NOMC' lmot2 lmot1 uH;
  393.  usup = umoy '+' uH;
  394.  uinf = umoy '-' uH;
  395. *-----------------------------< Effort
  396.  tmoy = 'EXCO' lmot3 fint1;
  397.  tH = 'EXCO' lmot4 fint1;
  398.  tH = 'NOMC' lmot4 lmot3 tH;
  399.  tsup = tmoy '+' tH;
  400.  tinf = tmoy '-' tH;
  401.  
  402.  
  403. 'FINPROC' uinf usup tinf tsup;
  404.  
  405.  
  406. *----------------------------------------------------------------------
  407. *
  408. * LOCABA3D : Définit une lèvre plus et moins pour tous les éléments
  409. *            d'interface à partir de la normale orienté dans le sens
  410. *            de parcours des éléments d'interface
  411. *
  412. *----------------------------------------------------------------------
  413.  
  414. 'DEBPROC' LOCABA3D MeshInt*MAILLAGE mod_int*MMODEL;
  415.  
  416.  chn = vsur mod_int 'NORM';
  417.  mo1 = 'EXTRAIRE' chn 'COMP';
  418.  mo2 = 'MOTS' 'NX  ' 'NY  ' 'NZ  ';
  419.  chn = 'CHANGER' 'CHPO' mod_int chn 'MOYE';
  420.  chn = 'NOMC' chn mo1 mo2;
  421.  cht1 = 0. * chn;
  422. * creation des tangentes arbitrairement orienté à partir de la normale
  423.  xt1 = 'NOMC' 'T1X ' ('EXCO' 'NY  ' chn);
  424.  yt1 = 'NOMC' 'T1Y ' ('EXCO' 'NX  ' chn);
  425.  zt1 = 'NOMC' 'T1Z ' ('EXCO' 'NZ  ' cht1);
  426.  'SI' ((('MAXIMUM' xt1) 'EGA' 0) 'ET' (('MAXIMUM' yt1) 'EGA' 0));
  427.    xt1 = 'MANUEL' chpo ('EXTRAIRE' xt1 mail) 'T1X ' 1.;
  428.  'FINSI';
  429.  cht1 = xt1 '+'   yt1 '+'   zt1;
  430. * on norme cht1
  431.  m1 = 'EXTRAIRE' cht1 'COMP';
  432.  cht1 = cht1 '/' (('PSCAL' cht1 cht1 m1 m1)**0.5);
  433.  
  434. * creation de la dernière tangente comme résultat du produit vectoriel
  435.  xn1 = 'NOMC' 'NX  ' ('EXCO' 'NX  ' chn);
  436.  yn1 = 'NOMC' 'NY  ' ('EXCO' 'NY  ' chn);
  437.  zn1 = 'NOMC' 'NZ  ' ('EXCO' 'NZ  ' chn);
  438.  xt1 = 'NOMC' 'T1X ' ('EXCO' 'T1X ' cht1);
  439.  yt1 = 'NOMC' 'T1Y ' ('EXCO' 'T1Y ' cht1);
  440.  zt1 = 'NOMC' 'T1Z ' ('EXCO' 'T1Z ' cht1);
  441.  mxt1 = 'MOTS' 'T1X '; myt1 = 'MOTS' 'T1Y '; mzt1 = 'MOTS' 'T1Z ';
  442.  mxt2 = 'MOTS' 'T2X '; myt2 = 'MOTS' 'T2Y '; mzt2 = 'MOTS' 'T2Z ';
  443.  mxn1 = 'MOTS' 'NX  '; myn1 = 'MOTS' 'NY  '; mzn1 = 'MOTS' 'NZ  ';
  444.  
  445. * pvec ne fonctionne pas avec les CHPO contrairement à ce qui est dit dans l'aide
  446. * on le fait à la main
  447.  xt2 = (yn1 '*' zt1 myn1 mzt1 mxt2) '-' (zn1 '*' yt1 mzn1 myt1 mxt2);
  448.  yt2 = (zn1 * xt1 mzn1 mxt1 myt2) '-' (xn1 * zt1 mxn1 mzt1 myt2);
  449.  zt2 = (xn1 * yt1 mxn1 myt1 mzt2) '-' (yn1 * xt1 myn1 mxt1 mzt2);
  450.  cht2 = xt2 '+' yt2 '+' zt2;
  451.  
  452.  
  453. 'FINPROC' chn cht1 cht2;
  454.  
  455.  
  456. *----------------------------------------------------------------------
  457. *
  458. * LOCABASE : Définit une lèvre plus et moins pour tous les éléments
  459. *            d'interface à partir de la normale orienté dans le sens
  460. *            de parcours des éléments d'interface
  461. * utiliser chn = vsur modele_interface 'NORM'; pour la création des normales
  462. *----------------------------------------------------------------------
  463.  
  464. 'DEBPROC' LOCABASE MeshInt*MAILLAGE;
  465.  
  466.  
  467.  Nbe1em = 'NBEL' MeshInt;
  468.  
  469. * 'SI' (nbe1em '>EG' 6);
  470.    'SI' (nbe1em '>EG' 6E6);
  471.  cht chn chb = @FRENET MeshInt;
  472.  
  473.  nbn1 = 'NBNO' MeshInt;   
  474.  po1 = MeshInt 'POIN' 1;
  475. * Mise à zéro des normales pour le premier et dernier noeud
  476.  chp1 = 'MANUEL' chpo po1 1 'SCAL' 0. natu discret;
  477.  po2 = MeshInt 'POIN' nbn1;
  478.  chp2 = 'MANUEL' chpo po2 1 'SCAL' 0.  natu discret;
  479.  chp3 = chp1  'ET' chp2;
  480.  
  481.  idnode = 2;
  482.  'REPETER' blo1 (nbn1 '-' 2);
  483.  po1 = MeshInt 'POIN' idnode; 
  484.  chp1 = 'MANUEL' chpo po1 1 'SCAL' 1. natu discret;
  485.  chp3 = chp3 'ET' chp1;
  486.  idnode = idnode '+' 1;
  487.  'FIN' blo1;
  488.  
  489.  lmot1 = 'MOTS' 'SCAL';
  490.  lmot2 = 'MOTS' 'UX  '  'UX  ';
  491.  
  492.  cht = chp3 * cht;
  493.  chn = chp3 * chn;
  494.  
  495. * Cas du premier et dernier noeud de l'interface 
  496.  
  497. * Premier noeud
  498.  po1 = MeshInt 'POIN' 1;
  499.  x1 = 'COORDONNEE' 1 po1;
  500.  y1 = 'COORDONNEE' 2 po1;
  501.  
  502.  po2 = MeshInt 'POIN' 2;
  503.  x2 = 'COORDONNEE' 1 po2;
  504.  y2 = 'COORDONNEE' 2 po2;
  505.  
  506.  pmilieu = 'POIN' (x1 '+' ((x2'-' x1) '/' 2))
  507.  (y1 '+' ((y2 '-' y1) '/' 2));
  508.  
  509.  Lseg = (((x1 '-' x2) '**' 2) '+'  ((y1 '-' y2) '**' 2 )) '**' (0.5);
  510.  a = (x2 '-' x1) '/' Lseg;
  511.  b = (y2 '-' y1) '/' Lseg;
  512.  b2 = -1. * b;
  513.  
  514.  
  515.  chp1 = 'MANUEL' chpo po1 2 'TX  ' a 'TY  ' b ;   
  516.  chp2 = 'MANUEL' chpo po1 2 'NX  ' b2 'NY  ' a;
  517.  
  518.  cht = cht '+' chp1;
  519.  chn = chn '+' chp2;
  520.  
  521.  
  522. * Dernier noeud
  523.  po1 = MeshInt 'POIN' (nbn1 '-' 1);
  524.  x1 = 'COORDONNEE' 1 po1;
  525.  y1 = 'COORDONNEE' 2 po1;
  526.  
  527.  po2 = MeshInt 'POIN' nbn1;
  528.  x2 = 'COORDONNEE' 1 po2;
  529.  y2 = 'COORDONNEE' 2 po2;
  530.  
  531.  pmilieu = 'POIN' (x1 '+' ((x2'-' x1) '/' 2))
  532.  (y1 '+' ((y2 '-' y1) '/' 2));
  533.  
  534.  Lseg = (((x1 '-' x2) '**' 2) '+'  ((y1 '-' y2) '**' 2 )) '**' (0.5);
  535.  a = (x2 '-' x1) '/' Lseg;
  536.  b = (y2 '-' y1) '/' Lseg;
  537.  b2 = -1. * b;
  538.  
  539.  chp1 = 'MANUEL' chpo po2 2 'TX  ' a 'TY  ' b ;   
  540.  chp2 = 'MANUEL' chpo po2 2 'NX  ' b2 'NY  ' a;
  541.  
  542.  cht = cht '+' chp1;
  543.  chn = chn '+' chp2;
  544.  
  545.  co1 = 'EXTRAIRE' cht 'COMP';
  546.  co2 = 'MOTS' 'T1X ' 'T1Y ';
  547.  
  548.  cht = 'NOMC' cht co1 co2;
  549.  
  550. *  cas avec moins de 5 éléments (@frenet ne fonctionne avec moins de 5 éléments)
  551.  'SINON';
  552.  
  553.    cht = 'MANUEL' chpo MeshInt 'T1X ' 0. 'T1Y ' 0.;
  554.    chn = 'MANUEL' chpo MeshInt 'NX  ' 0. 'NY  ' 0.;
  555.    Nbe1em = 'NBEL' MeshInt;
  556.    Nbnode =  'NBNO' MeshInt;
  557.    idnode = 1;
  558.  
  559.    'REPETER' blo1 (Nbnode);
  560.  
  561.       'SI' ((idnode 'EGA' 1) 'OU' (idnode 'EGA' Nbnode));
  562.          'SI' (idnode 'EGA' 1);
  563.             ele1 = 'ELEM' MeshInt 1;
  564.          'SINON';
  565.             ele1 = 'ELEM' MeshInt nbe1em;
  566.          'FINSI';
  567.          po1 = (ele1 point 1);
  568.          x1 = 'COORDONNEE' 1 po1;
  569.          y1 = 'COORDONNEE' 2 po1;
  570.  
  571.          po2 = (ele1 point 2);
  572.          x2 = 'COORDONNEE' 1 po2;
  573.          y2 = 'COORDONNEE' 2 po2;
  574.  
  575.    Lseg = (((x1 '-' x2) '**' 2) '+'  ((y1 '-' y2) '**' 2 )) '**' (0.5);
  576.          a = (x2 '-' x1) '/' Lseg;
  577.          b = (y2 '-' y1) '/' Lseg;
  578.          b2 = -1. * b;
  579.  
  580.          'SI' (idnode 'EGA' 1);
  581.             chp1 = 'MANUEL' chpo po1 2 'T1X ' a 'T1Y ' b ;   
  582.             chp2 = 'MANUEL' chpo po1 2 'NX  ' b2 'NY  ' a;
  583.          'SINON';
  584.             chp1 = 'MANUEL' chpo po2 2 'T1X ' a 'T1Y ' b ;   
  585.             chp2 = 'MANUEL' chpo po2 2 'NX  ' b2 'NY  ' a;
  586.          'FINSI';
  587.  
  588.       'SINON';
  589.          po1 = (MeshInt point (idnode '-' 1));
  590.          x1 = 'COORDONNEE' 1 po1;
  591.          y1 = 'COORDONNEE' 2 po1;
  592.          po2 = (MeshInt point idnode);
  593.          x2 = 'COORDONNEE' 1 po2;
  594.          y2 = 'COORDONNEE' 2 po2;
  595.          po3 = (MeshInt point (idnode '+' 1));
  596.          x3 = 'COORDONNEE' 1 po3;
  597.          y3 = 'COORDONNEE' 2 po3;
  598.  
  599.    Lseg1 = (((x1 '-' x2) '**' 2) '+'  ((y1 '-' y2) '**' 2 )) '**' (0.5);
  600.    Lseg2 = (((x3 '-' x2) '**' 2) '+'  ((y3 '-' y2) '**' 2 )) '**' (0.5);
  601.  
  602.          a1 = (x2 '-' x1) '/' Lseg1;
  603.          b1 = (y2 '-' y1) '/' Lseg1;
  604.          a2 = (x3 '-' x2) '/' Lseg2;
  605.          b2 = (y3 '-' y2) '/' Lseg2;
  606.          Lseg3 = (((a1 '+' a2)**2) '+' ((b1 '+' b2)**2)) '**' (0.5); 
  607.          chp1 = 'MANUEL' chpo po2 2 'T1X ' ((a1 + a2) '/' Lseg3) 'T1Y '
  608.                ((b1 '+' b2) '/' Lseg3) ;   
  609.          chp2 = 'MANUEL' chpo po2 2 'NX  '(-1 * (b2 '+' b1) '/' Lseg3)
  610.                               'NY  '  ((a1 + a2) '/' Lseg3);
  611.       'FINSI' ;
  612.  
  613.       cht = cht '+' chp1;
  614.       chn = chn '+' chp2;
  615.       idnode = idnode '+' 1;
  616.  
  617.    'FIN' blo1;
  618.  
  619.  'FINSI';
  620.  
  621. 'FINPROC' chn cht;
  622.  
  623.  
  624. *----------------------------------------------------------------------
  625. *
  626. * GLO2LOC : Projection des champs de l'interface du repère global
  627. *           vers le repère local
  628. *
  629. *----------------------------------------------------------------------
  630.  
  631. 'DEBPROC' GLO2LOC cht1*CHPOINT chn*CHPOINT Tplus*CHPOINT Tmoins*CHPOINT
  632. Wplus*CHPOINT Wmoins*CHPOINT cht2/CHPOINT;
  633.  
  634.  fg1 = 'EXTRAIRE' tplus 'COMP';
  635.  wg1 = 'EXTRAIRE' wplus 'COMP';
  636.  t1 = 'EXTRAIRE' cht1 'COMP';
  637.  n1 = 'EXTRAIRE' chn 'COMP' ;
  638.  'SI' (dim1 'EGA' 3);
  639.    t2 = 'EXTRAIRE' cht2 'COMP';
  640.  'FINSI' ;
  641.  
  642.  a = 'NOMC' 'FT1 ' (psca cht1 Tplus t1 fg1);
  643.  b = 'NOMC' 'FN  ' (psca chn Tplus n1 fg1);
  644.  TlocP = a '+' b;
  645.  
  646.  a = 'NOMC' 'FT1 ' (psca cht1 Tmoins t1 fg1);
  647.  b = 'NOMC' 'FN  ' (psca chn Tmoins n1 fg1);
  648.  TlocM = a '+' b;
  649.  
  650.  a = 'NOMC' 'UT1 ' (psca cht1 Wplus t1 wg1);
  651.  b = 'NOMC' 'UN  ' (psca chn Wplus n1 wg1);
  652.  WlocP = a '+' b;
  653.  
  654.  a = 'NOMC' 'UT1 ' (psca cht1 Wmoins t1 wg1);
  655.  b = 'NOMC' 'UN  ' (psca chn Wmoins  n1 wg1);
  656.  WlocM = a '+' b;
  657.  
  658.  'SI' (('DIME' fg1) 'EGA' 3);
  659.    c = 'NOMC' 'FT2 ' (psca cht2 Tplus t2 fg1);
  660.    TlocP = TlocP '+' c;
  661.    c = 'NOMC' 'FT2 ' (psca cht2 Tmoins t2 fg1);
  662.    TlocM = TlocM '+' c;
  663.    c = 'NOMC' 'UT2 ' (psca cht2 Wplus t2 wg1);
  664.    WlocP = WlocP '+' c;
  665.    c = 'NOMC' 'UT2 ' (psca cht2 Wmoins t2 wg1);
  666.    WlocM = WlocM '+' c;
  667.  'FINSI';
  668.  
  669.  'FINPROC' TlocP TlocM WlocP WlocM;
  670.  
  671.  
  672. *----------------------------------------------------------------------
  673. *
  674. * LOCSTEP1 : Etape locale pour le premier par de temps
  675. *
  676. *----------------------------------------------------------------------
  677.  
  678.  'DEBPROC' LOCSTEP1 Klatin*FLOTTANT f_crack*FLOTTANT MeshInt*MAILLAGE
  679.     TlocP*CHPOINT TlocM*CHPOINT WlocP*CHPOINT WlocM*CHPOINT it1*ENTIER
  680.     WlocPOld*CHPOINT WlocMOld*CHPOINT;
  681.  
  682.  meshloc = 'CHANGER' POI1 meshint;
  683.  wdiffold= wlocmold '-' wlocpold;
  684.  
  685.  wdiff= wlocm '-' wlocp;
  686.  tdiff= tlocm '-' tlocp;
  687.  k_tan = klatin;
  688.  k_nor = klatin;
  689.  
  690.  TlocPnor = 'EXCO' 'FN  ' TlocP;
  691.  TlocPta1 = 'EXCO' 'FT1 ' TlocP;
  692.  TlocMnor = 'EXCO' 'FN  ' TlocM;
  693.  TlocMta1 = 'EXCO' 'FT1 ' TlocM;
  694.  
  695.  WlocPnor = 'EXCO' 'UN  ' WlocP;
  696.  WlocPta1 = 'EXCO' 'UT1 ' WlocP;
  697.  WlocMnor = 'EXCO' 'UN  ' WlocM;
  698.  WlocMta1 = 'EXCO' 'UT1 ' WlocM;
  699.  
  700.  w1old = 'EXCO' 'UT1 ' wdiffold;
  701.  wm1old = 'EXCO' 'UT1 ' wlocmold;
  702.  wp1old = 'EXCO' 'UT1 ' wlocpold;
  703.  
  704. * initialisation booleen pour détecter fissure complètement ouverte
  705.  ouv1 = vrai;
  706.  
  707.  'SI' (('DIME' ('EXTRAIRE' tlocp 'COMP')) 'EGA' 3);
  708.  TlocPta2 = 'EXCO' 'FT2 ' TlocP; TlocMta2 = 'EXCO' 'FT2 ' TlocM;
  709.  WlocPta2 = 'EXCO' 'UT2 ' WlocP; WlocMta2 = 'EXCO' 'UT2 ' WlocM;
  710.  TMta2 = 0. * TlocMta2;TPta2 = 0. * TlocPta2;
  711.  WMta2 = 0. * WlocMta2;WPta2 = 0. * WlocPta2;
  712.  w2old = 'EXCO' 'UT2 ' wdiffold;
  713.  wm2old = 'EXCO' 'UT2 ' wlocmold;
  714.  wp2old = 'EXCO' 'UT2 ' wlocpold;
  715.  B3D = vrai;
  716.  'SINON';
  717.  B3D = faux;
  718.  'FINSI';
  719. *
  720. * Cacul des indicateur de contact et de frottement initial
  721.  ConIni = 0.5 * (('EXCO' 'UN  'wdiff)'-'(('EXCO' 'FN' tdiff)'/'k_nor));
  722.  
  723. * btrolle le 5 novembre pour coller avec ceux qui est écrit dans NewlocalStepOnInterface dans ELFE
  724.  GliIni1 =(0.5*((K_tan*('EXCO' 'UT1 'wdiff))'-'('EXCO' 'FT1' tdiff)))'-'
  725.  (0.5 * K_tan * w1old);
  726.  
  727.  'SI' B3D;
  728.  GliIni2 =(0.5*((K_tan*('EXCO' 'UT2 'wdiff))'-'('EXCO' 'FT2' tdiff)))'-'
  729.  (0.5 * K_tan * w2old);
  730.  GliIni = ((gliIni1**2)'+'(gliIni2**2))'**'(0.5);
  731.  'SINON';
  732.    GliIni = GliIni1;
  733.  'FINSI' ;
  734.  
  735. * seuil pour le glissement
  736.  GSeuil =  f_crack * ('ABS' (K_nor * ConIni)); 
  737.  
  738. * Initialisation des champs correspondant à l'itération i + 1/2 
  739.  TMnor = 0. * TlocMnor; TMta1 = 0. * TlocMta1; 
  740.  TPnor = 0. * TlocMnor; TPta1 = 0. * TlocPta1; 
  741.  WMnor = 0. * WlocMnor; WMta1 = 0. * WlocMta1; 
  742.  WPnor = 0. * WlocMnor; WPta1 = 0. * WlocPta1;
  743.  
  744. * initialisation des compteurs
  745. Ncon = 0;
  746. Nouv = 0;
  747.  
  748. *
  749. * Loi de comportement local i --> i + 1/2 avec Contact with friction (coulomb)
  750. *---
  751.  Nbn1 = 'NBNO' Meshloc;
  752.  idnode = 1;
  753. * boucle sur les noeuds de l'interface
  754.  'REPETER' blo1 (Nbn1);
  755.  
  756.    ptemp = Meshloc 'POIN' idnode;
  757.    ValCI = 'EXTRAIRE' ConIni 'SCAL' ptemp;
  758.  
  759. *  Valeur à l'itération i
  760.    wpn = 'EXTRAIRE' WlocPnor 'SCAL' ptemp;
  761.    wmn = 'EXTRAIRE' WlocMnor 'SCAL' ptemp;
  762.    wpt1 = 'EXTRAIRE' WlocPta1 'SCAL' ptemp;
  763.    wmt1 = 'EXTRAIRE' WlocMta1 'SCAL' ptemp;
  764.  
  765.    tpn = 'EXTRAIRE' TlocPnor 'SCAL' ptemp;
  766.    tmn = 'EXTRAIRE' TlocMnor 'SCAL' ptemp;
  767.    tpt1 = 'EXTRAIRE' TlocPta1 'SCAL' ptemp;
  768.    tmt1 = 'EXTRAIRE' TlocMta1 'SCAL' ptemp;
  769.  
  770.    wpo1 = 'EXTRAIRE' wp1old 'SCAL' ptemp;
  771.    wmo1 = 'EXTRAIRE' wm1old 'SCAL' ptemp;
  772.  
  773.    'SI' B3D;
  774.       tpt2 = 'EXTRAIRE' TlocPta2 'SCAL' ptemp;
  775.       tmt2 = 'EXTRAIRE' TlocMta2 'SCAL' ptemp;
  776.       wpt2 = 'EXTRAIRE' WlocPta2 'SCAL' ptemp;
  777.       wmt2 = 'EXTRAIRE' WlocMta2 'SCAL' ptemp;
  778.       wpo2 = 'EXTRAIRE' wp2old 'SCAL' ptemp;
  779.       wmo2 = 'EXTRAIRE' wm2old 'SCAL' ptemp;
  780.    'FINSI';
  781.  
  782. *  Si ouverture  : >EG 0. pour le cas ou on bloque les sauts passe par là ;-)
  783.    'SI' (ValCI '>EG' 0.);
  784.       Nouv = Nouv '+' 1;
  785. *     'MESSAGE' '---> Ouvert';
  786. *     calcul des nouveaux déplacemtents  (t= 0)
  787.       chp1 = 'MANUEL' chpo ptemp 'SCAL'  (wpn  '-' (tpn '/' K_nor));
  788.       chp2 = 'MANUEL' chpo ptemp 'SCAL'  (wmn  '-' (tmn '/' K_nor));
  789.       chp3 = 'MANUEL' chpo ptemp 'SCAL'  (wpt1  '-' (tpt1 '/' K_nor));
  790.       chp4 = 'MANUEL' chpo ptemp 'SCAL'  (wmt1  '-' (tmt1 '/' K_nor));      
  791.       WPnor = WPnor '+' chp1; WMnor = WMnor '+' chp2;
  792.       WPta1 = WPta1 '+' chp3; WMta1 = WMta1 '+' chp4;
  793.  
  794.    'SI' B3D;    
  795.       chp5 = 'MANUEL' chpo ptemp 'SCAL'  (wpt2  '-' (tpt2 '/' K_nor));
  796.       chp6 = 'MANUEL' chpo ptemp 'SCAL'  (wmt2  '-' (tmt2 '/' K_nor));
  797.       WPta2 = WPta2 '+' chp5; WMta2 = WMta2 '+' chp6;
  798.    'FINSI';
  799.  
  800. *  Si contact   
  801.    'SINON';
  802.    Ncon = Ncon '+' 1;
  803.  
  804. * Pas de test d'ouverture pour le premier et le dernier noeud.. il faudrait le faire juste pour le ou les fronts pour faire les choses correctements...
  805.    'SI' ((idnode > 1) 'ET' (idnode < Nbn1)); ouv1 = faux; 'FINSI';
  806. *  problème normal
  807. *     w
  808.       chp1 = 'MANUEL' chpo ptemp 'SCAL'  (0.5 * ((wpn '+' wmn )  '-'
  809.       ((tpn  '+' tmn) '/'  K_nor)));
  810.       WPnor = WPnor '+' chp1; WMnor = WMnor '+' chp1;
  811.  
  812. *     t
  813.       chp2 = 'MANUEL' chpo ptemp 'SCAL'  (K_nor * ValCI);   
  814.       TPnor = TPnor '+' chp2; TMnor = TMnor '-' chp2;
  815.  
  816. *     extrait les valeurs pour tester adhérence ou glissement   
  817.       ValGI = 'EXTRAIRE' GliIni 'SCAL' ptemp; 
  818.       G = 'EXTRAIRE' GSeuil 'SCAL' ptemp;
  819.       ValGI1 = 'EXTRAIRE' GliIni1 'SCAL' ptemp; 
  820.       si B3D;
  821.          ValGI2 = 'EXTRAIRE' GliIni2 'SCAL' ptemp; 
  822.       'FINSI';
  823.  
  824. *     Adhérence   
  825.       'SI' (('ABS' ValGI) '<'  G);
  826. *         'MESSAGE' '---> ADHERENCE';
  827. *         'MESSAGE' 'ValGI1 = ' ValGI1 'G =' G 'ValCI =' ValCI;
  828. *        t
  829.          chp31 = 'MANUEL' chpo ptemp 'SCAL' valGI1;
  830.          TPta1 = TPta1 '+' chp31; TMta1 = TMta1 '-' chp31;
  831.  
  832. *        w
  833.          chp41='MANU' chpo ptemp'SCAL'((wpt1 '-' wpo1) '+'((valGI1'-'
  834.          tpt1) '/' K_tan)); 
  835.          WPta1 = WPta1 '+' chp41'+' ('MANU' chpo ptemp'SCAL' wpo1);
  836.          WMta1 = WMta1 '+' chp41'+' ('MANU' chpo ptemp'SCAL' wmo1);
  837.  
  838.          'SI' B3D;
  839. *        t
  840.          chp32 = 'MANUEL' chpo ptemp 'SCAL' valGI2;
  841.          TPta2 = TPta2 '+' chp32; TMta2 = TMta2 '-' chp32;
  842.  
  843. *        w         
  844.          chp42='MANU' chpo ptemp 'SCAL' ((wpt2 '-' wpo2)'+'((valGI2'-'
  845.          tpt2) '/' K_tan));
  846.          WPta2 = WPta2 '+' chp42'+' ('MANU' chpo ptemp'SCAL' wpo2);
  847.          WMta2 = WMta2 '+' chp42'+' ('MANU' chpo ptemp'SCAL' wmo2);
  848.          'FINSI';
  849.  
  850. *     Glissement
  851.       'SINON';
  852. *         'LISTE' ptemp;
  853. *         'MESSAGE' '---> Glissement';
  854. *         'MESSAGE' 'ValGI = ' ValGI 'G =' G 'ValCI =' ValCI;
  855.  
  856. *        t
  857.          'SI' (valGI 'EGA' 0.);
  858. *         chp31 = 'MANUEL' chpo ptemp 'SCAL' 0.;
  859.          'SINON' ;
  860.          chp31 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI1'/'('ABS' valGI));
  861. *      chp31 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI1'/'(2*('ABS' valGI)));
  862.          'FINSI';
  863.          TPta1 = TPta1 '+' chp31; TMta1 = TMta1 '-' chp31;
  864.  
  865. *        w
  866. * newlocalstep oninterface
  867.         chp41 = 'MANUEL' chpo ptemp 'SCAL'
  868. (wpt1 '-' wpo1 '+'( ( (G*valGI1'/'('ABS' valGI)) '-' tpt1) '/' K_tan)); 
  869.         chp51 = 'MANUEL' chpo ptemp 'SCAL'
  870. (wmt1 '-' wmo1 '+'(((-1.*G*valGI1'/'('ABS' valGI))'-'tmt1) '/' K_tan)); 
  871.         WPta1 = WPta1 '+' chp41 '+' ('MANU' chpo ptemp'SCAL' wpo1);
  872.         WMta1 = WMta1 '+' chp51 '+' ('MANU' chpo ptemp'SCAL' wmo1);
  873.  
  874.       'SI' B3D;
  875.          'SI' (valGI 'EGA' 0.);
  876.             chp32 = 'MANUEL' chpo ptemp 'SCAL' 0.;
  877.          'SINON' ;
  878.           chp32 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI2'/'('ABS' valGI)); 
  879.          'FINSI';
  880.          chp42 = 'MANUEL' chpo ptemp 'SCAL'
  881. ((wpt2 '- 'wpo2) '+' ( ( (G*valGI2'/'('ABS' valGI))'-'tpt2)'/'K_tan)); 
  882.          chp52 = 'MANUEL' chpo ptemp 'SCAL'
  883. ((wmt2 '-' wmo2)'+'(((-1.*G*valGI2'/'('ABS' valGI))'-'tmt2)'/'K_tan)); 
  884.          TPta2 = TPta2 '+' chp32; TMta2 = TMta2 '-' chp32;
  885.          WPta2 = WPta2 '+' chp42'+' ('MANU' chpo ptemp'SCAL' wpo2);
  886.          WMta2 = WMta2 '+' chp52'+' ('MANU' chpo ptemp'SCAL' wmo2);
  887.       'FINSI';
  888.  
  889.       'FINSI';
  890.  
  891.    'FINSI';
  892.  
  893.    idnode = idnode '+' 1;
  894.  
  895. 'FIN' blo1;
  896.  
  897.  tpnor = 'NOMC' 'FN  ' tpnor natu discret;
  898.  tpta1 = 'NOMC' 'FT1 ' tpta1 natu discret;
  899.  tmnor = 'NOMC' 'FN  ' tmnor natu discret;
  900.  tmta1 = 'NOMC' 'FT1 ' tmta1 natu discret;
  901.  wpnor = 'NOMC' 'UN  ' wpnor natu diffus;
  902.  wpta1 = 'NOMC' 'UT1 ' wpta1 natu diffus;
  903.  wmnor = 'NOMC' 'UN  ' wmnor natu diffus;
  904.  wmta1 = 'NOMC' 'UT1 ' wmta1 natu diffus;
  905.  
  906.  TlocP_a = TPta1 'ET' TPnor;
  907.  TlocM_a = TMta1 'ET' TMnor;
  908.  WlocP_a = WPta1 'ET' WPnor;
  909.  WlocM_a = WMta1 'ET' WMnor;
  910.  
  911.  'SI' B3D;
  912.    tpta2 = 'NOMC' 'FT2 ' tpta2 natu discret;
  913.    tmta2 = 'NOMC' 'FT2 ' tmta2 natu discret;
  914.    wpta2 = 'NOMC' 'UT2 ' wpta2 natu diffus;
  915.    wmta2 = 'NOMC' 'UT2 ' wmta2 natu diffus;
  916.    TlocP_a = TlocP_a 'ET' Tpta2;
  917.    TlocM_a = TlocM_a 'ET' Tmta2;
  918.    WlocP_a = WlocP_a 'ET' WPta2;
  919.    WlocM_a = WlocM_a 'ET' WPta2;
  920.  'FINSI' ;
  921.  
  922.  
  923.  'MESSAGE' ' ';
  924.  'MESSAGE' 'Nombre de noeud en contact = ' Ncon;
  925.  'MESSAGE' 'Nombre de noeud ouvert     = ' Nouv;
  926.  'MESSAGE' ' ';
  927.  
  928.  'FINPROC' TlocP_a TlocM_a WlocP_a WlocM_a ouv1;
  929.  
  930.  
  931. *----------------------------------------------------------------------
  932. *
  933. * LOC2GLO : Projetion des champs à l'interface du repère local vers
  934. *           le repère global
  935. *
  936. *----------------------------------------------------------------------
  937.  
  938. 'DEBPROC' LOC2GLO cht1*CHPOINT chn*CHPOINT TlocP*CHPOINT TlocM*CHPOINT
  939.  WlocP*CHPOINT WlocM*CHPOINT MeshInt*MAILLAGE cht2/CHPOINT;
  940.  
  941.  fl1 = 'EXTRAIRE' TlocM 'COMP' ;
  942.  wl1 = 'EXTRAIRE' WlocM 'COMP' ;
  943.  t1 = 'EXTRAIRE' cht1 'COMP' ;
  944.  n1 = 'EXTRAIRE' chn 'COMP' ;
  945.  vecx = 'MANUEL' chpo MeshInt 'UX  ' 1. 'UY  ' 0.;
  946.  vecy = 'MANUEL' chpo MeshInt 'UX  ' 0. 'UY  ' 1.;
  947.  
  948.  'SI' (('DIME' fl1) 'EGA' 3);
  949.    vecx = 'MANUEL' chpo MeshInt 'UX  ' 1. 'UY  ' 0. 'UZ  ' 0. ;
  950.    vecy = 'MANUEL' chpo MeshInt 'UX  ' 0. 'UY  ' 1. 'UZ  ' 0. ;
  951.    vecz = 'MANUEL' chpo MeshInt 'UX  ' 0. 'UY  ' 0. 'UZ  ' 1. ;
  952.    t2 = 'EXTRAIRE' cht2 'COMP';
  953.    B3D = vrai;
  954.  'SINON';
  955.    B3D = faux;
  956.  'FINSI';
  957.  
  958.  mx1 = 'EXTRAIRE' vecx 'COMP'; 
  959.  
  960.   a = 'NOMC' 'UT1 ' ('PSCA' vecx cht1 mx1 t1);
  961.   b = 'NOMC' 'UN  ' ('PSCA' vecx chn mx1 n1);
  962.   chx = a '+' b;
  963.   'SI' B3D;
  964.       c = 'NOMC' 'UT2 '('PSCA' vecx cht2 mx1 t2);
  965.       chx = chx '+' c;
  966.   'FINSI';
  967.  
  968.   a = 'NOMC' 'UT1 ' ('PSCA' vecy cht1 mx1 t1);
  969.   b = 'NOMC' 'UN  ' ('PSCA' vecy chn mx1 n1);
  970.   chy = a '+' b;
  971.   'SI' B3D;
  972.       c = 'NOMC' 'UT2 '('PSCA' vecy cht2 mx1 t2);
  973.       chy = chy '+' c;
  974.       a = 'NOMC' 'UT1 ' ('PSCA' vecz cht1 mx1 t1);
  975.       b = 'NOMC' 'UN  ' ('PSCA' vecz chn mx1 n1);
  976.       c = 'NOMC' 'UT2 '('PSCA' vecz cht2 mx1 t2);
  977.       chz = a '+' b '+' c;       
  978.   'FINSI';
  979.  
  980.   a = 'NOMC' 'FY  ' (psca chy TlocP wl1 fl1);
  981.   b = 'NOMC' 'FX  ' (psca chx TlocP wl1 fl1);
  982.   TPglo = b '+' a;
  983.  
  984.   a = 'NOMC' 'FY  ' (psca chy TlocM wl1 fl1);
  985.   b = 'NOMC' 'FX  ' (psca chx TlocM wl1 fl1);
  986.   TMglo = b '+' a;
  987.  
  988.   a = 'NOMC' 'UY  ' (psca chy WlocP wl1 wl1);
  989.   b = 'NOMC' 'UX  ' (psca chx WlocP wl1 wl1);
  990.   WPglo = b '+' a;
  991.  
  992.   a = 'NOMC' 'UY  ' (psca chy WlocM wl1 wl1);
  993.   b = 'NOMC' 'UX  ' (psca chx WlocM wl1 wl1);
  994.   WMglo = b '+'a;
  995.  
  996.   'SI' B3D;
  997.    c = 'NOMC' 'FZ  ' (psca chz TlocP wl1 fl1);
  998.    TPglo = TPglo '+' c;
  999.    c = 'NOMC' 'FZ  ' (psca chz TlocM wl1 fl1);
  1000.    TMglo = TMglo '+' c;
  1001.    c = 'NOMC' 'UZ  ' (psca chz WlocP wl1 wl1);
  1002.    WPglo = WPglo '+' c;
  1003.    c = 'NOMC' 'UZ  ' (psca chz WlocM wl1 wl1);
  1004.    WMglo = WMglo '+' c;
  1005.  'FINSI';
  1006.  
  1007. 'FINPROC'  WMglo WPglo TMglo TPglo;
  1008.  
  1009.  
  1010. *----------------------------------------------------------------------
  1011. *
  1012. * LIPSRECO : Recombine les champs locaux des 2 lèvres de la fissure
  1013. *            issus de l'étape locale
  1014. *            sur les maillages de l'étape globale en différenciant
  1015. *            et en calculant les inconnus XFEM correspondantes
  1016. *
  1017. *----------------------------------------------------------------------
  1018.  
  1019. 'DEBPROC' LIPSRECO  Wmoins*CHPOINT Wplus*CHPOINT Tmoins*CHPOINT
  1020.                                                  Tplus*CHPOINT;
  1021.  
  1022.  Tmoy = (Tmoins '+' Tplus) '/' 2; 
  1023.  Wmoy = (Wmoins '+' Wplus) '/' 2;
  1024.  tmo1 = 'EXTRAIRE' tmoy 'COMP';
  1025.  wmo1 = 'EXTRAIRE' wmoy 'COMP';
  1026.  'SI' (('DIME' wmo1) 'EGA' 2);
  1027.    tmo2 = 'MOTS' 'FAX ' 'FAY ';
  1028.    wmo2 = 'MOTS' 'AX  ' 'AY  ';
  1029.  'SINON';
  1030.    tmo2 = 'MOTS' 'FAX ' 'FAY ' 'FAZ ';
  1031.    wmo2 = 'MOTS' 'AX  ' 'AY  ' 'AZ  ';
  1032.  'FINSI';
  1033.  
  1034.  Tsaut = 'NOMC' tmo1 tmo2 ((Tplus '-' Tmoy));
  1035.  Wsaut = 'NOMC' wmo1 wmo2 ((Wplus '-' Wmoy));
  1036.  Tglo = Tmoy '+' Tsaut;
  1037.  Wglo = Wmoy '+' Wsaut;
  1038.  
  1039. 'FINPROC' Wglo Tglo;
  1040.  
  1041.  
  1042. *----------------------------------------------------------------------
  1043. *
  1044. * UPDATRHS : Mise à jour du second membre du sytème linéaire
  1045. *
  1046. *----------------------------------------------------------------------
  1047.  
  1048. 'DEBPROC' UPDATRHS Kuwt*RIGIDITE Tglo*CHPOINT Wglo*CHPOINT
  1049.          meshint*MAILLAGE it1*ENTIER f1*CHPOINT;
  1050.  
  1051.  lmot1 = 'EXTRAIRE' wglo 'COMP';
  1052.  lmot2 = 'EXTRAIRE' Tglo 'COMP';
  1053.  f2 = 'ENLEVER' ('REDU' (kuwt * wglo) MeshInt) lmot1;
  1054.  f3 = 'ENLEVER' ('REDU' (kuwt * tglo) MeshInt) lmot1;
  1055.  f4 = enle ('REDU' (kuwt * tglo) MeshInt) lmot2;
  1056.  f5 = f3 '+' f2  '+' f4;
  1057.  RHS1 = f1 '+' f5;
  1058.  
  1059. 'FINPROC' rhs1;
  1060.  
  1061.  
  1062. *----------------------------------------------------------------------
  1063. *
  1064. * ERRORINI : Calcul de l'indicateur d'erreur INITIALE
  1065. *            basé sur les quantités w et t sur l'interface
  1066. *            dans les bases locales
  1067. *            entre la solition i + 1/2 et i + 1
  1068. *            fixe les dénominateurs pour les appels suivant -> GETERROR
  1069. *----------------------------------------------------------------------
  1070.  
  1071. 'DEBPROC' ERRORINI TlocP_n*CHPOINT TlocM_n*CHPOINT WlocP_n*CHPOINT
  1072. WlocM_n*CHPOINT TlocP*CHPOINT TlocM*CHPOINT WlocP*CHPOINT WlocM*CHPOINT
  1073. Klatin*FLOTTANT;
  1074.  
  1075.  
  1076. *-> i + 1/2 :   TlocP TlocM WlocP WlocM
  1077.  TPnor = 'EXCO' 'FN  ' TlocP;
  1078.  TPta1 = 'EXCO' 'FT1  ' TlocP;
  1079.  TMnor = 'EXCO' 'FN  ' TlocM;
  1080.  TMta1 = 'EXCO' 'FT1  ' TlocM;
  1081.  
  1082.  WPnor = 'EXCO' 'UN  ' WlocP;
  1083.  WPta1 = 'EXCO' 'UT1 ' WlocP;
  1084.  WMnor = 'EXCO' 'UN  ' WlocM;
  1085.  WMta1 = 'EXCO' 'UT1 ' WlocM;
  1086.  
  1087. *-> i + 1   :   TlocP_n TlocM_n WlocP_n WlocM_n ( _n = new )
  1088.  TPnor_n = 'EXCO' 'FN  ' TlocP_n;
  1089.  TPta1_n = 'EXCO' 'FT1 ' TlocP_n;
  1090.  TMnor_n = 'EXCO' 'FN  ' TlocM_n;
  1091.  TMta1_n = 'EXCO' 'FT1 ' TlocM_n;
  1092.  
  1093.  WPnor_n = 'EXCO' 'UN  ' WlocP_n;
  1094.  WPta1_n = 'EXCO' 'UT1 ' WlocP_n;
  1095.  WMnor_n = 'EXCO' 'UN  ' WlocM_n;
  1096.  WMta1_n = 'EXCO' 'UT1 ' WlocM_n; 
  1097.  
  1098.  
  1099. * ---> NUMERATEUR
  1100.  
  1101. * problème normal 
  1102. numN = 'MAXIMUM' (
  1103. ((((TPnor_n '-' TPnor)**2)'+' ((TMnor_n '-' TMnor)**2)) '/' Klatin) '+'
  1104. ((((WPnor_n '-' WPnor)**2)'+' ((WMnor_n '-' WMnor)**2)) '*'  Klatin));
  1105.  
  1106. * problème tangentiel
  1107. numT= 'MAXIMUM' (
  1108. ((((TPta1_n '-' TPta1)**2)'+' ((TMta1_n '-' TMta1)**2)) '/' Klatin) '+'
  1109. ((((WPta1_n '-' WPta1)**2)'+' ((WMta1_n '-' WMta1)**2)) '*'  Klatin));
  1110.  
  1111.  
  1112. * ---> DENOMINATEUR
  1113.  
  1114. * problème normal
  1115. SNor_n = (((TPnor_n * TPnor_n) '+' (TMnor_n * TMnor_n) ) '/' Klatin) '+'
  1116. (Klatin * ((WPnor_n *  WPnor_n) '+' (WMnor_n *  WMnor_n)));
  1117.  
  1118. SNor = (((TPnor * TPnor) '+' (TMnor * TMnor) ) '/' Klatin) '+'
  1119. (Klatin * ((WPnor *  WPnor) '+' (WMnor *  WMnor)));
  1120.  
  1121. denN = ('MAXIMUM' SNor_n ) '+' ('MAXIMUM' SNor );
  1122.  
  1123.  
  1124. * problème tangentiel
  1125. STa1_n = (((TPta1_n * TPta1_n) '+' (TMta1_n * TMta1_n) ) '/' Klatin) '+'
  1126. (Klatin * ((WPta1_n *  WPta1_n) '+' (WMta1_n *  WMta1_n)));
  1127.  
  1128. STa1 = (((TPta1 * TPta1) '+' (TMta1 * TMta1) ) '/' Klatin) '+'
  1129. (Klatin * ((WPta1 *  WPta1) '+' (WMta1 *  WMta1)));
  1130.  
  1131. denT = ('MAXIMUM' STa1_n ) '+' ('MAXIMUM' STa1 );
  1132.  
  1133. * ---> ERREUR
  1134.  
  1135. *-> erreur normale
  1136. etaN = numN '/' denN;
  1137.  
  1138.  
  1139. *-> erreur tangentielle
  1140. etaT = numT '/' denT;
  1141.  
  1142. *-> erreur totale : ComputeLatinErrorNew dans ELFE-3D 
  1143.  totalErr = ((etaN) '+' (etaT )) '**' 0.5;
  1144.  
  1145. *-> erreur totale : publi récente (= manuscrit Epierres)
  1146. *ler1 = 'PROG' (etaT**0.5) (etaN**0.5);
  1147. *totalerr = 'MAXIMUM' ler1;
  1148. *   'SI' (etaT > etaN);
  1149. *   totalErr = etaT;     
  1150. * 'SINON';
  1151. *   totalErr = etaN
  1152. * 'FINSI';
  1153.  
  1154. 'FINPROC' totalErr (denN) (denT) (etaT**0.5) (etaN**0.5);
  1155.  
  1156.  
  1157. *----------------------------------------------------------------------
  1158. *
  1159. * GETERROR : Calcul de l'indicateur d'erreur 
  1160. *            basé sur les quantités w et t sur l'interface
  1161. *            dans les bases locales
  1162. *            entre la solition i + 1/2 et i + 1
  1163. *            Les DENOMINATEURS DOIVENT ETRE FOURNIT EN ENTREE !
  1164. *----------------------------------------------------------------------
  1165.  
  1166. 'DEBPROC' GETERROR TlocP_n*CHPOINT TlocM_n*CHPOINT WlocP_n*CHPOINT
  1167. WlocM_n*CHPOINT TlocP*CHPOINT TlocM*CHPOINT WlocP*CHPOINT WlocM*CHPOINT 
  1168. den_N*FLOTTANT den_T*FLOTTANT Klatin*FLOTTANT ;
  1169.  
  1170. * exclusion des points anciennement et actuellement extrémités
  1171. * de la fissure pour le calcul de l'erreur
  1172.  
  1173.  
  1174. *-> i + 1/2 :   TlocP TlocM WlocP WlocM
  1175.  TPnor = 'EXCO' 'FN  ' TlocP;
  1176.  TPta1 = 'EXCO' 'FT1 ' TlocP;
  1177.  TMnor = 'EXCO' 'FN  ' TlocM;
  1178.  TMta1 = 'EXCO' 'FT1 ' TlocM;
  1179.  
  1180.  WPnor = 'EXCO' 'UN  ' WlocP;
  1181.  WPta1 = 'EXCO' 'UT1 ' WlocP;
  1182.  WMnor = 'EXCO' 'UN  ' WlocM;
  1183.  WMta1 = 'EXCO' 'UT1 ' WlocM;
  1184.  
  1185. *-> i + 1   :   TlocP_n TlocM_n WlocP_n WlocM_n ( _n = new )
  1186.  TPnor_n = 'EXCO' 'FN  ' TlocP_n;
  1187.  TPta1_n = 'EXCO' 'FT1 ' TlocP_n;
  1188.  TMnor_n = 'EXCO' 'FN  ' TlocM_n;
  1189.  TMta1_n = 'EXCO' 'FT1 ' TlocM_n;
  1190.  
  1191.  WPnor_n = 'EXCO' 'UN  ' WlocP_n;
  1192.  WPta1_n = 'EXCO' 'UT1 ' WlocP_n;
  1193.  WMnor_n = 'EXCO' 'UN  ' WlocM_n;
  1194.  WMta1_n = 'EXCO' 'UT1 ' WlocM_n;
  1195.  
  1196.  
  1197. * ---> NUMERATEUR
  1198.  
  1199. * problème normal 
  1200. numN = 'MAXIMUM' (
  1201. ((((TPnor_n '-' TPnor)**2)'+' ((TMnor_n '-' TMnor)**2)) '/' Klatin) '+'
  1202. ((((WPnor_n '-' WPnor)**2)'+' ((WMnor_n '-' WMnor)**2)) '*' Klatin));
  1203.  
  1204.  
  1205. * problème tangentiel
  1206. numT= 'MAXIMUM' (
  1207. ((((TPta1_n '-' TPta1)**2)'+' ((TMta1_n '-' TMta1)**2)) '/' Klatin) '+'
  1208. ((((WPta1_n '-' WPta1)**2)'+' ((WMta1_n '-' WMta1)**2)) '*' Klatin));
  1209.  
  1210.  
  1211. * ---> ERREUR
  1212.  
  1213. *-> erreur normale
  1214.  etaN = numN '/' den_N;
  1215.  
  1216.  
  1217. *-> erreur tangentielle
  1218.  etaT = numT '/' den_T;
  1219.  
  1220. *-> erreur totale : ComputeLatinErrorNew dans ELFE-3D 
  1221.  totalErr = ((etaN) '+' (etaT )) '**' 0.5;
  1222.  
  1223. *-> erreur totale : publi récente (= manuscrit Epierres)
  1224. *ler1 = 'PROG' (etaT**0.5) (etaN**0.5);
  1225. *totalerr = ('MAXIMUM' ler1)'**'(0.5);
  1226. *'SI' (etaT > etaN);
  1227. *   totalErr = etaT;     
  1228. * 'SINON';
  1229. *   totalErr = etaN
  1230. * 'FINSI';
  1231.  
  1232. 'FINPROC' totalerr (etaT**0.5) (etaN**0.5);
  1233.  
  1234.  
  1235. *----------------------------------------------------------------------
  1236. *
  1237. * LOCSTEPI : Etape locale pour le premier par de temps
  1238. *
  1239. *----------------------------------------------------------------------
  1240.  
  1241.  'DEBPROC' LOCSTEPI Klatin*FLOTTANT f_crack*FLOTTANT MeshInt*MAILLAGE
  1242.     TlocP*CHPOINT TlocM*CHPOINT WlocP*CHPOINT WlocM*CHPOINT;
  1243.  
  1244.  meshloc = 'CHANGER' POI1 meshint;
  1245.  wdiff= wlocm '-' wlocp;
  1246.  tdiff= tlocm '-' tlocp;
  1247.  k_tan = klatin;
  1248.  k_nor = klatin;
  1249.  
  1250.  TlocP = 1.*TlocP;
  1251.  TlocM = 1.*TlocM;
  1252.  WlocP = 1.*WlocP;
  1253.  WlocM = 1.*WlocM;
  1254.  
  1255.  TlocPnor = 'EXCO' 'FN  ' TlocP;
  1256.  TlocPta1 = 'EXCO' 'FT1 ' TlocP;
  1257.  TlocMnor = 'EXCO' 'FN  ' TlocM;
  1258.  TlocMta1 = 'EXCO' 'FT1 ' TlocM;
  1259.  
  1260.  WlocPnor = 'EXCO' 'UN  ' WlocP;
  1261.  WlocPta1 = 'EXCO' 'UT1 ' WlocP;
  1262.  WlocMnor = 'EXCO' 'UN  ' WlocM;
  1263.  WlocMta1 = 'EXCO' 'UT1 ' WlocM;
  1264.  
  1265. * initialisation booleen pour détecter fissure complètement ouverte
  1266.  ouv1 = vrai;
  1267.  
  1268.  'SI' (('DIME' ('EXTRAIRE' tlocp 'COMP')) 'EGA' 3);
  1269.  TlocPta2 = 'EXCO' 'FT2 ' TlocP; TlocMta2 = 'EXCO' 'FT2 ' TlocM;
  1270.  WlocPta2 = 'EXCO' 'UT2 ' WlocP; WlocMta2 = 'EXCO' 'UT2 ' WlocM;
  1271.  TMta2 = 0. * TlocMta2;TPta2 = 0. * TlocPta2;
  1272.  WMta2 = 0. * WlocMta2;WPta2 = 0. * WlocPta2;
  1273.  B3D = vrai;
  1274.  'SINON';
  1275.  B3D = faux;
  1276.  'FINSI';
  1277. *
  1278. * Cacul des indicateur de contact et de frottement initial
  1279.  ConIni = 0.5 * (('EXCO' 'UN  'wdiff)'-'(('EXCO' 'FN' tdiff)'/'k_nor));
  1280.  GliIni1 = 0.5*((K_tan*('EXCO' 'UT1 'wdiff))'-'('EXCO' 'FT1' tdiff));
  1281.  
  1282.  'SI' B3D;
  1283.  GliIni2 = 0.5*((K_tan*('EXCO' 'UT2 'wdiff))'-'('EXCO' 'FT2' tdiff));
  1284.  GliIni = ((gliIni1**2)'+'(gliIni2**2))'**'(0.5);
  1285.  'SINON';
  1286.    GliIni = GliIni1;
  1287.  'FINSI' ;
  1288.  
  1289. * seuil pour le glissement
  1290.  GSeuil =  f_crack * ('ABS' (K_nor * ConIni)); 
  1291.  
  1292. * Initialisation des champs correspondant à l'itération i + 1/2 
  1293.  TMnor = 0. * TlocMnor; TMta1 = 0. * TlocMta1; 
  1294.  TPnor = 0. * TlocMnor; TPta1 = 0. * TlocPta1; 
  1295.  WMnor = 0. * WlocMnor; WMta1 = 0. * WlocMta1; 
  1296.  WPnor = 0. * WlocMnor; WPta1 = 0. * WlocPta1;
  1297.  
  1298. * initialisation des compteurs
  1299. Ncon = 0;
  1300. Nouv = 0;
  1301. *
  1302. * Loi de comportement local i --> i + 1/2 avec Contact with friction (coulomb)
  1303. *---
  1304.  Nbn1 = 'NBNO' Meshloc;
  1305.  idnode = 0;
  1306.  
  1307.  'REPETER' blo1 (Nbn1);
  1308.    idnode = idnode '+' 1;
  1309.    ptemp = Meshloc 'POIN' idnode;
  1310.    ValCI = 'EXTRAIRE' ConIni 'SCAL' ptemp;
  1311.  
  1312. *  Valeur à l'itération i
  1313.    wpn = 'EXTRAIRE' WlocPnor 'SCAL' ptemp;
  1314.    wmn = 'EXTRAIRE' WlocMnor 'SCAL' ptemp;
  1315.    wpt1 = 'EXTRAIRE' WlocPta1 'SCAL' ptemp;
  1316.    wmt1 = 'EXTRAIRE' WlocMta1 'SCAL' ptemp;
  1317.  
  1318.    tpn = 'EXTRAIRE' TlocPnor 'SCAL' ptemp;
  1319.    tmn = 'EXTRAIRE' TlocMnor 'SCAL' ptemp;
  1320.    tpt1 = 'EXTRAIRE' TlocPta1 'SCAL' ptemp;
  1321.    tmt1 = 'EXTRAIRE' TlocMta1 'SCAL' ptemp;
  1322.  
  1323.    'SI' B3D;
  1324.       tpt2 = 'EXTRAIRE' TlocPta2 'SCAL' ptemp;
  1325.       tmt2 = 'EXTRAIRE' TlocMta2 'SCAL' ptemp;
  1326.       wpt2 = 'EXTRAIRE' WlocPta2 'SCAL' ptemp;
  1327.       wmt2 = 'EXTRAIRE' WlocMta2 'SCAL' ptemp;
  1328.    'FINSI';
  1329.  
  1330. *  Si ouverture  : >EG 0. pour le cas ou on bloque les sauts passe par là ;-)
  1331.    'SI' (ValCI '>EG' 0.);
  1332. *      'MESSAGE' 'OUVERTURE pour le noeud' idnode;
  1333.       Nouv = Nouv '+' 1;
  1334. *     calcul des noueaux déplacemtents  (t= 0)
  1335.       chp1 = 'MANUEL' chpo ptemp 'SCAL'  (wpn  '-' (tpn '/' K_nor));
  1336.       chp2 = 'MANUEL' chpo ptemp 'SCAL'  (wmn  '-' (tmn '/' K_nor));
  1337.       chp3 = 'MANUEL' chpo ptemp 'SCAL'  (wpt1  '-' (tpt1 '/' K_nor));
  1338.       chp4 = 'MANUEL' chpo ptemp 'SCAL'  (wmt1  '-' (tmt1 '/' K_nor));      
  1339.       WPnor = WPnor '+' chp1; WMnor = WMnor '+' chp2;
  1340.       WPta1 = WPta1 '+' chp3; WMta1 = WMta1 '+' chp4;
  1341.  
  1342.    'SI' B3D;    
  1343.       chp5 = 'MANUEL' chpo ptemp 'SCAL'  (wpt2  '-' (tpt2 '/' K_nor));
  1344.       chp6 = 'MANUEL' chpo ptemp 'SCAL'  (wmt2  '-' (tmt2 '/' K_nor));
  1345.       WPta2 = WPta2 '+' chp5; WMta2 = WMta2 '+' chp6;
  1346.    'FINSI';
  1347.  
  1348. *  Si contact   
  1349.    'SINON';
  1350.          Ncon = Ncon '+' 1;
  1351. *         'MESSAGE' 'CONTACT pour le noeud' idnode;
  1352.    'SI' ((idnode > 1) 'ET' (idnode < Nbn1)); ouv1 = faux; 'FINSI';
  1353.  
  1354. *  problème normal
  1355. *     w
  1356.       chp1 = 'MANUEL' chpo ptemp 'SCAL'  (0.5 * ((wpn '+' wmn )  '-'
  1357.       ((tpn  '+' tmn) '/'  K_nor)));
  1358.       WPnor = WPnor '+' chp1; WMnor = WMnor '+' chp1;
  1359. *     t
  1360.       chp2 = 'MANUEL' chpo ptemp 'SCAL'  (K_nor * ValCI);   
  1361.       TPnor = TPnor '+' chp2; TMnor = TMnor '-' chp2;
  1362.  
  1363. *     extraction seuil adhérence
  1364.       ValGI = 'EXTRAIRE' GliIni 'SCAL' ptemp; 
  1365.       G = 'EXTRAIRE' GSeuil 'SCAL' ptemp;
  1366.       ValGI1 = 'EXTRAIRE' GliIni1 'SCAL' ptemp; 
  1367.       si B3D;
  1368.          ValGI2 = 'EXTRAIRE' GliIni2 'SCAL' ptemp; 
  1369.       'FINSI';
  1370.  
  1371. *     Adhérence   
  1372.       'SI' (('ABS' ValGI) '<'  G);
  1373.  
  1374. *        t
  1375.          chp31 = 'MANUEL' chpo ptemp 'SCAL' valGI1;
  1376.          TPta1 = TPta1 '+' chp31; TMta1 = TMta1 '-' chp31;
  1377.  
  1378. *        w
  1379.         chp41='MANU' chpo ptemp'SCAL'(wpt1'+'((valGI1'-'tpt1)'/'K_tan)); 
  1380.          WPta1 = WPta1 '+' chp41; WMta1 = WMta1 '+' chp41;
  1381.  
  1382.          'SI' B3D;
  1383.         chp32 = 'MANUEL' chpo ptemp 'SCAL' valGI2;
  1384.         chp42='MANU' chpo ptemp'SCAL'(wpt2'+'((valGI2'-'tpt2)'/'K_tan));
  1385.         TPta2 = TPta2 '+' chp32; TMta2 = TMta2 '-' chp32;
  1386.         WPta2 = WPta2 '+' chp42; WMta2 = WMta2 '+' chp42;
  1387.          'FINSI';
  1388.  
  1389. *     Glissement
  1390.       'SINON';
  1391. *         'LISTE' ptemp;
  1392. *         'MESSAGE' '---> Glissement';
  1393.  
  1394. *        t
  1395.          'SI' (valGI 'EGA' 0.);
  1396. *         chp31 = 'MANUEL' chpo ptemp 'SCAL' 0.;
  1397.          'SINON' ;
  1398.          chp31 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI1'/' ('ABS' valGI));
  1399.          'FINSI';
  1400.          TPta1 = TPta1 '+' chp31; TMta1 = TMta1 '-' chp31;
  1401.  
  1402.  
  1403. *        w
  1404.          chp41 = 'MANUEL' chpo ptemp 'SCAL'
  1405.          (wpt1 '+'( ( (G*valGI1'/'('ABS' valGI)) '-' tpt1) '/' K_tan)); 
  1406.          chp51 = 'MANUEL' chpo ptemp 'SCAL'
  1407.          (wmt1 '+'(((-1.*G*valGI1'/'('ABS' valGI))'-'tmt1) '/' K_tan)); 
  1408.          WPta1 = WPta1 '+' chp41; WMta1 = WMta1 '+' chp51;
  1409.  
  1410.       'SI' B3D;
  1411.          'SI' (valGI 'EGA' 0.);
  1412.             chp32 = 'MANUEL' chpo ptemp 'SCAL' 0.;
  1413.          'SINON' ;
  1414.           chp32 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI2'/'('ABS' valGI)); 
  1415.          'FINSI';
  1416.          chp32 = 'MANUEL' chpo ptemp 'SCAL'(G*valGI2'/'('ABS' valGI));         
  1417.          chp42 = 'MANUEL' chpo ptemp 'SCAL'
  1418.          (wpt2 '+'( ( (G*valGI2'/'('ABS' valGI)) '-' tpt2) '/' K_tan)); 
  1419.          chp52 = 'MANUEL' chpo ptemp 'SCAL'
  1420.          (wmt2 '+'(((-1.*G*valGI2'/'('ABS' valGI))'-'tmt2) '/' K_tan)); 
  1421.          TPta2 = TPta2 '+' chp32; TMta2 = TMta2 '-' chp32;
  1422.          WPta2 = WPta2 '+' chp42; WMta2 = WMta2 '+' chp52;
  1423.       'FINSI';
  1424.  
  1425.       'FINSI';
  1426.  
  1427.    'FINSI';
  1428.  
  1429. 'FIN' blo1;
  1430.  
  1431.  tpnor = 'NOMC' 'FN  ' tpnor natu discret;
  1432.  tpta1 = 'NOMC' 'FT1 ' tpta1 natu discret;
  1433.  tmnor = 'NOMC' 'FN  ' tmnor natu discret;
  1434.  tmta1 = 'NOMC' 'FT1 ' tmta1 natu discret;
  1435.  wpnor = 'NOMC' 'UN  ' wpnor natu diffus;
  1436.  wpta1 = 'NOMC' 'UT1 ' wpta1 natu diffus;
  1437.  wmnor = 'NOMC' 'UN  ' wmnor natu diffus;
  1438.  wmta1 = 'NOMC' 'UT1 ' wmta1 natu diffus;
  1439.  
  1440.  TlocP_a = TPta1 'ET' TPnor;
  1441.  TlocM_a = TMta1 'ET' TMnor;
  1442.  WlocP_a = WPta1 'ET' WPnor;
  1443.  WlocM_a = WMta1 'ET' WMnor;
  1444.  
  1445.  'SI' B3D;
  1446.    tpta2 = 'NOMC' 'FT2 ' tpta2 natu discret;
  1447.    tmta2 = 'NOMC' 'FT2 ' tmta2 natu discret;
  1448.    wpta2 = 'NOMC' 'UT2 ' wpta2 natu diffus;
  1449.    wmta2 = 'NOMC' 'UT2 ' wmta2 natu diffus;
  1450.    TlocP_a = TlocP_a 'ET' Tpta2;
  1451.    TlocM_a = TlocM_a 'ET' Tmta2;
  1452.    WlocP_a = WlocP_a 'ET' WPta2;
  1453.    WlocM_a = WlocM_a 'ET' WPta2;
  1454.  'FINSI' ;
  1455.  
  1456.  'MESSAGE' ' ';
  1457.  'MESSAGE' 'Nombre de noeud en contact = ' Ncon;
  1458.  'MESSAGE' 'Nombre de noeud ouvert     = ' Nouv;
  1459.  'MESSAGE' ' ';
  1460.  
  1461.  'FINPROC' TlocP_a TlocM_a WlocP_a WlocM_a;
  1462.  
  1463.  
  1464. ************************************************************************
  1465. *                               AFFICHAGE                              *
  1466. ************************************************************************
  1467.  
  1468. OPTI TRAC PSC EPTR 6 POTR HELVETICA_16;
  1469. GRAPH = vrai;
  1470. COMPLET = faux;
  1471. * opti debug 1;
  1472.  
  1473. * Gestion de l'affichage au cours d'un calcul*
  1474.   Tdisplay = tabl ;
  1475. * Affichage des levels set en début de cycle
  1476.   Tdisplay . levelset = GRAPH;
  1477. * Affichage du repère locale  
  1478.   Tdisplay . rep_local = GRAPH;
  1479. * Affichage des déplacement, déformés, contraintes, au cours de itérations de la LATIN
  1480. *  (pour chaque pas de temps !)
  1481.   Tdisplay . Iteration = faux;
  1482. * Affichage des pressions de contact sur déforméee locale
  1483.   Tdisplay . PasTemps = GRAPH;
  1484. *  Affichage de la convergence de la LATIN
  1485. *  (pour chaque pas de temps !)
  1486.   Tdisplay . convergence = GRAPH;
  1487. * Affichage des FICS
  1488.   Tdisplay . fics = GRAPH;
  1489. * Affichage du champ de déplacement global  
  1490.   Tdisplay . dep_global = vrai;   
  1491.  
  1492.  
  1493.  
  1494. ************************************************************************
  1495. *                        PARAMETRES DU CALCUL                          *
  1496. ************************************************************************
  1497.  
  1498.  dim1 = 2;
  1499.  
  1500. *----------------- < Précision pour solveur non-linéaire
  1501. si COMPLET;
  1502. *  eps1 = 1E-6;
  1503.   eps1 = 1E-8;
  1504. sino;
  1505.   eps1 = 1E-3;
  1506. fins;
  1507.  
  1508. *----------------- < Nombre d'itéraion max
  1509.  Itmax = 1000;
  1510.  
  1511. *----------------- < Materiau
  1512.  Ey1  = 2.06E11;
  1513.  nu1  = 0.3;
  1514.  rho1 = 7800.0;
  1515.  
  1516. *----------------- < Frottement des levres de la Fissure
  1517.  f_crack = 0.1;
  1518.  
  1519. *----------------- < Chargement
  1520.  f_ext = -1.0E9;
  1521.  
  1522. *----------------- < facteur d'amplification pour les tracés
  1523. * effort
  1524. * fac1 = 0;
  1525.  fac1 =  0.05 /(('ABS' f_ext));
  1526.  
  1527. * déformée
  1528.  fac2 = 0.01*('ABS' (Ey1 '/' f_ext));
  1529.  
  1530. *----------------- < Paramètre géométrique Hammouda 
  1531.  L = 0.4;
  1532.  W = 0.2;
  1533.  a = 0.02;
  1534.  beta = 45.;
  1535. *----------------- < Maillage structure
  1536.  N1 = 10;
  1537. * taille élément
  1538.  den4 = a '/' N1;
  1539. *----------------- < Taille élément interface 
  1540.  den3 = den4 '/' 1;
  1541.  
  1542.  
  1543. ************************************************************************
  1544. *                              MAILLAGE                                *
  1545. ************************************************************************
  1546.  
  1547. *----------------- < Maillage structure 1
  1548.  
  1549. *----< surface fine de propagation
  1550.  p2 = (W '+' (2 *  a)) (W '+' (2*a));
  1551.  p3 = (W '+' (2 * a)) (W '-' (2*a));
  1552.  p4 = (W '-' (2 * a)) (W '+' (2*a));
  1553.  p5 = (W '-' (2 * a)) (W '-' (2*a));
  1554.  
  1555.  l5 = 'DROIT' p2 p3 'DINI' den4 'DFIN' den4;
  1556.  l6 = 'DROIT' p3 p5 'DINI' den4 'DFIN' den4;
  1557.  l7 = 'DROIT' p5 p4 'DINI' den4 'DFIN' den4;
  1558.  l8 = 'DROIT' p4 p2 'DINI' den4 'DFIN' den4;
  1559.  lfi1 = l5 'ET' l6 'ET' l7 'ET' l8; 
  1560.  s2 = 'DALLER' l5 l6 l7 l8 'PLAN'; 
  1561.  
  1562. *----< contour grossier
  1563.  p0 = 0. 0.;
  1564.  p1 = (2*W) 0.;
  1565.  p9 = (2*W) L;
  1566.  p8 = 0. L;
  1567.  l1 = 'DROIT' p0 N1 p1 ;
  1568.  l2 = 'DROIT' p1 (N1) p9 ;
  1569.  l3 = 'DROIT' p9 N1 p8 ;
  1570.  l4 = 'DROIT' p8 (N1) p0 ;
  1571.  lco1 = l1 'ET' l2 'ET' l3 'ET' l4;
  1572.  lto1 = lco1 'ET' lfi1;
  1573.  s1 = 'SURFACE' lto1 'PLANE';
  1574.  linbas1 = l1;
  1575.  linhaut = l3;
  1576.  ling = l2;
  1577.  lind = l4;
  1578.  s3 = s1 'ET' s2;
  1579.  
  1580. *----------------- < Maillage interface (fissure) (si 1 seule front, 
  1581. *                    définir tip1 comme le front)
  1582.  tip2 = (W '+' (a * ('COS' beta))) (W '+' (a * ('SIN' beta)));
  1583.  tip1 = (W '-' (a * ('COS' beta))) (W '-' (a * ('SIN' beta)));
  1584.  crack = 'DROIT' tip2 tip1 'DINI' den3 'DFIN' den3 'COULEUR' viol;
  1585.  stot1 = s3 'ET' crack;
  1586.  si (GRAPH); TRAC stot1 'TITR' 'Maillage'; fins;
  1587.  
  1588. * Creation des Levels Sets 
  1589. *---------------
  1590.  psi0 phi0  = PSIPHI s2 crack 'DEUX' tip1 tip2;
  1591.  si( Tdisplay.levelset);
  1592.    TRAC phi0 s2 'TITR' 'Level Set \f';
  1593.    TRAC psi0 s2 'TITR' 'Level Set \y';
  1594.  finsi;
  1595.  
  1596.  
  1597. ************************************************************************
  1598. *        MODELE, ENRICHISSEMENT, MATÉRIAU, MATRICE, RELATIONS          *
  1599. ************************************************************************
  1600.  
  1601. *----------------- < Modele
  1602.  mod1a = MODE s1 MECANIQUE ELASTIQUE ISOTROPE; 
  1603.  mod1b = MODE s2 MECANIQUE ELASTIQUE ISOTROPE xq4r;
  1604. * Il faut d'abord donner le modèle X-FEM pour qu'il soit défini en 
  1605. * premier dans le modèle totale
  1606.  mod1 = mod1b 'ET' mod1a;
  1607.  
  1608. *----------------- < Enrichissement
  1609.  Che1X = tabl;
  1610. * Che1X . 0 = TRIE mod1 psi0 phi0 saut;
  1611.  Che1X . 0 = TRIE mod1 psi0 phi0; 
  1612.  
  1613. * constructionsion des blocages des ddl X-fem non actifs dans 
  1614. * les éléments de transition.
  1615.  Rel1X = tabl;
  1616.  Rel1X . 0 = RELA mod1;
  1617.  
  1618. *----------------- < Matériau
  1619.  mat1 = MATE mod1 'YOUN' Ey1 'NU' nu1 'RHO' rho1;
  1620.  
  1621. *----------------- < Latin
  1622.  Klatin = 1*(Ey1/('MESURE' crack)) ;
  1623.  
  1624. *----------------- < Stabilisation
  1625.  xi = 1*(-1) '/' (Klatin);
  1626. *'OPTION' donn 5;
  1627. *-----------------------------< Creation des differentes matrices blocs 
  1628. *                               du système
  1629.  kuwt = RELA 'ACCRO' 'FAIBL' 'PENA' Klatin 'STAB' xi 'NOMU' mod1 Crack;
  1630. * multiplication par facteur 2 à cause de l'utilisation d'enrichissement 
  1631. * saut sur l'interface 
  1632.  kuwt = 2 * kuwt;
  1633.  
  1634. *-----------------------------< Rigidité formulation XFEM 1 champ
  1635.  k1 = 'RIGIDITE' mod1 mat1;
  1636.  
  1637. *-----------------------------< Condition aux limites 
  1638.  cl1 = (BLOQ  'DEPL' (linbas1 'POIN' 'INITIAL')) 'ET'
  1639.        (BLOQ  'DEPL' (linbas1 'POIN' 'FINAL'));
  1640.  cl1 = cl1 'ET' ('BLOQUE' 'DEPL' 'DIRECTION' (0. 1.) linbas1);
  1641.  
  1642. *-----------------------------< Matrice globale 1 champ
  1643.  K1tot  = k1  'ET' cl1 'ET' (Rel1X.0);
  1644.  
  1645. *-----------------------------< Matrice globale 3 champs
  1646.  Ktot1  = K1tot  'ET' Kuwt;
  1647.  
  1648. *-----------------------------< Effort imposé
  1649.  RHS0 = PRES 'MASS' mod1 (-1.*f_ext) linhaut;
  1650.  si (GRAPH); TRAC (vect RHS0 FORC roug) stot1 'TITR' 'chargement'; fins;
  1651.  
  1652. *-----------------------------< Construction repère local
  1653.  mod_int = 'MODELISER' crack 'MECANIQUE' 'ZCO2';
  1654.  chn cht1  = LOCABASE crack ;
  1655.  chn = -1.*chn;
  1656.  cht1 = -1.*cht1;
  1657.  
  1658. *-----------------------------< Affichage repère local
  1659.  vn = vect chn 'NX  ' 'NY  '  jaune; 
  1660.  vt1 = vect cht1 'T1X ' 'T1Y '  blanc; 
  1661.  vloc = vn 'ET' vt1;
  1662.  'SI' Tdisplay.rep_local;
  1663.    'TRACER' vloc (crack coul bleu);
  1664.  'FINSI';
  1665.  
  1666.  
  1667. ************************************************************************
  1668. *                               CALCUL                                 *
  1669. ************************************************************************
  1670.  
  1671. *-----------------------------< Calcul avec un unique pas de temps, 
  1672. *     seule la solution du problème de rupture nous interresse ici
  1673.  IdPasT = 1;
  1674.  
  1675.  
  1676. *-----------------------------<  Resolution
  1677.  'SI' (IdPasT 'EGA' 1);
  1678.    u0 Tplus Tmoins Wplus Wmoins err = LATIN Tdisplay crack
  1679.    K1tot mod1 RHS0 IdPasT ktot1 kuwt cht1 chn ItMax eps1 klatin
  1680.    f_crack s3;
  1681.  'SINON' ;
  1682.    u0 Tplus Tmoins Wplus Wmoins err = LATIN Tdisplay crack
  1683.    K1tot mod1 RHS0 IdPasT ktot1 kuwt cht1 chn ItMax eps1 klatin
  1684.    f_crack v1 cht2 Tplus Tmoins Wplus Wmoins;
  1685.  'FINSI';
  1686.  
  1687.  
  1688. ************************************************************************
  1689. *                        POST-TRAITEMENTS                              *
  1690. ************************************************************************
  1691.  
  1692. *-----------------------------< Affichage de la pressions de contact 
  1693. *                               sur la déformée
  1694. 'SI'(Tdisplay.PasTemps);
  1695.    b = xfem reco u0 mod1;
  1696.    def0 =  defo s3 b fac2;
  1697.    VTplus = vect Tplus forc fac1 rouge;
  1698.    VTmoins = vect Tmoins forc fac1 vert;
  1699.    def1 = defo wplus crack fac2  VTplus rouge;                          
  1700.    def2 = defo wmoins crack fac2  VTmoins vert;
  1701.    def3 = def1 'ET' def2 'ET' def0;
  1702.    'TRACER' def3 TITR 'deformee';
  1703.    'TRACER' (def1 et def2) TITR 'deformee';
  1704.  'FINSI';
  1705. *-----------------------------< Convergence
  1706.  'SI' (Tdisplay.convergence);
  1707. *  Réglage de la fenêtre de tracé
  1708.    min1 = ('MINIMUM' err);  
  1709.    max1 = ('MAXI' err);
  1710.    abs2 = 'MAXI' ('PROG' ('DIME' err));
  1711.    abs1 = 'PROG' 1. PAS 1. abs2;
  1712.  
  1713. *  Erreur
  1714.    evo1 = 'EVOL' 'MANUEL' abs1 err;
  1715.    DESS  evo1 
  1716.    'GRIL' 'POIN' 'GRIS'
  1717.    'LOGX' 'XBOR' 1 abs2 'TITX' 'Iterations' POSX CENT
  1718.    'LOGY' 'YBOR' min1 max1 'TITY' 'Erreur' POSY CENT
  1719.    TITR 'Convergence de la LATIN';
  1720.  'FINSI';
  1721.  
  1722.  
  1723. *-----------------------------< Ajout des termes de contacts dans le 
  1724. *                               calcul des FICs
  1725.    wsaut = wplus '-' wmoins;
  1726.    lmow = 'EXTRAIRE' wsaut 'COMP';
  1727.    mow1 = 'MOTS' 'AX  ' 'AY  ';
  1728.    wsaut = 'NOMC' lmow mow1 wsaut;
  1729.  
  1730.    tpglo = EXCO (tplus) (mots 'FX' 'FY') 'NOID';
  1731.    lmo1 = 'EXTRAIRE' tpglo  'COMP';
  1732.    lmo2 = 'MOTS' 'SMX ' 'SMY  ';
  1733.    tpglo = 'NOMC' lmo1 lmo2 tpglo;
  1734.  
  1735. *-----------------------------< Calcul des FICs : front 1
  1736.    'MESSAGE' ' ';'MESSAGE' '----------------------------';
  1737.    'MESSAGE' 'Calcul des FICs : front 1';
  1738.    'MESSAGE' '----------------------------';'MESSAGE' ' ';
  1739.    KTAB1 = TABL;
  1740.    KTAB1 . 'OBJECTIF' = MOT 'DECOUPLAGE';
  1741.    KTAB1 . 'PSI' = psi0;
  1742.    KTAB1 . 'PHI' = phi0;
  1743.    KTAB1 . 'FRONT_FISSURE' = tip1;
  1744.    KTAB1 . 'MODELE' = mod1;
  1745.    KTAB1 . 'CARACTERISTIQUES' = mat1;
  1746.    KTAB1 . 'SOLUTION_RESO' = u0;
  1747.    KTAB1 . 'CHARGEMENTS_MECANIQUES' = RHS0;
  1748.    KTAB1 . 'COUCHE' = 2;
  1749.    KTAB1 . 'MODELE_FISSURE' = mod_int;
  1750.    KTAB1 . 'DEPLACEMENT_FISSURE' = wsaut;
  1751.    KTAB1 . 'PRESSION_FISSURE' = TPglo ;   
  1752.  
  1753.    G_THETA KTAB1;
  1754.  
  1755.    K1G_1 = KTAB1 . 'RESULTATS' . 'I';
  1756.    K2G_1 = KTAB1 . 'RESULTATS' . 'II';
  1757.    GK1K2_1 = Ey1 * ((K1G_1**2) + (K2G_1**2)) ;
  1758. * * pas de valeur négative pour KI 
  1759. *    'SI' (K1G_1 < 0.);
  1760. *       K1G_1 = 0.;
  1761. *    'FINSI' ;
  1762.  
  1763. *-----------------------------< Calcul des FICs : front 2
  1764.    'MESSAGE' ' ';'MESSAGE' '----------------------------';
  1765.    'MESSAGE' 'Calcul des FICs : front 2';
  1766.    'MESSAGE' '----------------------------';'MESSAGE' ' ';
  1767.    KTAB2 = TABL;
  1768.    KTAB2 . 'OBJECTIF' = MOT 'DECOUPLAGE';
  1769.    KTAB2 . 'PSI' = psi0;
  1770.    KTAB2 . 'PHI' = phi0;
  1771.    KTAB2 . 'FRONT_FISSURE' = tip2;
  1772.    KTAB2 . 'MODELE' = mod1;
  1773.    KTAB2 . 'CARACTERISTIQUES' = mat1;
  1774.    KTAB2 . 'SOLUTION_RESO' = u0;
  1775.    KTAB2 . 'CHARGEMENTS_MECANIQUES' = RHS0;
  1776.    KTAB2 . 'COUCHE' = 2;
  1777.    KTAB2 . 'MODELE_FISSURE' = mod_int;
  1778.    KTAB2 . 'DEPLACEMENT_FISSURE' = wsaut;
  1779.    KTAB2 . 'PRESSION_FISSURE' = TPglo ;   
  1780.  
  1781.    G_THETA KTAB2;
  1782.    K1G_2 = KTAB2 . 'RESULTATS' . 'I';
  1783.    K2G_2 = KTAB2 . 'RESULTATS' . 'II';
  1784.    GK1K2_2 = Ey1 * ((K1G_2**2) + (K2G_2**2)) ;
  1785. * * pas de valeur négative pour KI 
  1786. *    'SI' (K1G_2 < 0.);
  1787. *       K1G_2 = 0.;
  1788. *    'FINSI' ;
  1789.  
  1790.  'SI' (Tdisplay.fics);
  1791. * Affichage de la nouvelle fissure
  1792.    q7v1 = VECT 2.E-3 (KTAB1 . CHAMP_THET1) 'VERT';
  1793.    q7v2 = VECT 2.E-3 (KTAB2 . CHAMP_THET1) 'BLEU';
  1794.    'TRACER'  (q7v1 et q7v2) (s2 et crack) TITR 'champs \q';
  1795.    'MESSAGE' ' ';'MESSAGE' ' ';'MESSAGE' ' ';
  1796.    'MESSAGE' '----------------------------';
  1797.    'MESSAGE' 'K1 front1 = ' K1G_1;
  1798.    'MESSAGE' 'K2 front1 = ' K2G_1;
  1799.    'MESSAGE' '----------------------------';
  1800.    'MESSAGE' ' ';'MESSAGE' ' ';'MESSAGE' ' ';
  1801.    'MESSAGE' '----------------------------';
  1802.    'MESSAGE' 'K1 front2 = ' K1G_2;
  1803.    'MESSAGE' 'K2 front2 = ' K2G_2;
  1804.    'MESSAGE' '----------------------------';
  1805.  'FINSI';
  1806.  
  1807.  
  1808. 'SI' (Tdisplay . dep_global);
  1809. *-----------------------------<AFfichage champ de déplacement
  1810.  u2 = enle u0 'LX  ';
  1811.  u2 = 'REDU' u2 s3;
  1812.  u2 = xfem reco u2 mod1;
  1813.  ux2 =  'NOMC' 'SCAL' ('EXCO' 'UX  ' u2);
  1814.  uy2 =  'NOMC' 'SCAL'('EXCO' 'UY  ' u2);
  1815.  uxy2 = 1.*(((ux2 **2)'+'(uy2 **2)));
  1816.  'TRACER' uxy2 s3 10;
  1817. 'FINSI';
  1818.  
  1819. mess ' ';mess ' ';mess ' ';
  1820. 'MESSAGE' '----------------------';
  1821. 'MESSAGE' 'FIN DES POST-TRAITEMENTS';
  1822.  
  1823.  
  1824.  
  1825. ************************************************************************
  1826. ***                    TEST DE BON FONCTIONNEMENT                    ***
  1827. ************************************************************************
  1828.  
  1829. * valeur de reference obtenue en 2013 :
  1830. K2ref = -11.5E7;
  1831. xtol  = 0.10;
  1832. err21 = abs ((K2G_1 - K2ref) / K2ref);
  1833. err22 = abs ((K2G_2 - K2ref) / K2ref);
  1834. mess 'ecart relatif sur K2:' err21 err22;
  1835. si ((err21 >eg xtol) ou (err22 >eg xtol));
  1836.   ERRE 5;
  1837. fins;
  1838.  
  1839. * OPTI donn 5 trac X;
  1840. fin ;
  1841.  
  1842.  
  1843.  
  1844.  
  1845.  
  1846.  
  1847.  
  1848.  
  1849.  
  1850.  
  1851.  
  1852.  
  1853.  
  1854.  
  1855.  
  1856.  

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