Download ch_thetx.procedur

Back to the list

   1 : * CH_THETX  PROCEDUR  BP208322  12/10/04    21:15:18     7519           
   2 : *    
   3 : DEBPROC CH_THETX SUPTAB*TABLE ;
   4 : *|=====================================================================|
   5 : *|                                                                     |
   6 : *|  << OBJET >> :                                                      |
   7 : *|                                                                     |
   8 : *|  Procedure analogue a CH_THETA mais utilisee avec des elements XFEM |
   9 : *|  Procedure determinant un champ/point de type THETA, c'est-a-dire   |
  10 : *|  un champ/point dont le norme est constant a l'interieur d'une      |
  11 : *|  courronne entourant le front d'une fissure, zero a l'exterieur de  |
  12 : *|  cette couronne. Le vecteur represente par le champ THETA indique   |
  13 : *|  la direction de propagation eventuelle de la fissure.              |
  14 : *|                                                                     |
  15 : *|  << ENTREE >> :                                                     |
  16 : *|                                                                     |
  17 : *|  SUPTAB = Objet de type TABLE dont les indices sont des             |
  18 : *|           objets de type MOT (a ecrire en toutes lettres) :         |
  19 : *|                                                                     |
  20 : *|  ARGUMENTS OBLIGATOIRES                                             |
  21 : *|  °°°°°°°°°°°°°°°°°°°°°°                                             |
  22 : *|                                                                     |
  23 : *|  SUPTAB.'MAILLAGE' = Objet de type MAILLAGE representant soit       |
  24 : *|                      la structure totale etudiee (maillage          |
  25 : *|                      utilise dans l'analyse par elements finis,     |
  26 : *|                      soit, pour reduire le temps de calcul, le      |
  27 : *|                      maillage entourant le plus grand des contours  |
  28 : *|                      qu'on a defini pour calculer le champ THETA.   |
  29 : *|                                                                     |
  30 : *|  SUPTAB.'PSI'  = Objet de type CHPOINT representant la 1ere level   |
  31 : *|                  decrivant le repere local de la fissure            |
  32 : *|  SUPTAB.'PHI'  = Objet de type CHPOINT representant la 2eme level   |
  33 : *|                  decrivant le repere local de la fissure            |
  34 : *|                                                                     |
  35 : *|  SUPTAB.'FRONT_FISSURE' = Objet representant le front de fissure    |
  36 : *|                     - facultatif et de type POINT en 2D             |
  37 : *|                     - obligatoire et de type MAILLAGE (ligne) en 3D |
  38 : *|                                                                     |
  39 : *|  SUPTAB.'COUCHE'   = Objet de type ENTIER representant le nombre    |
  40 : *|                      de couches d'elements (autour du point de      |
  41 : *|                      fissure) qui se deplacent pour simuler la      |
  42 : *|                      propagtion de la fissure.                      |
  43 : *|                                                                     |
  44 : *|                                                                     |
  45 : *|                                                                     |
  46 : *|  << SORTIE >> :                                                     |
  47 : *|                                                                     |
  48 : *|  TETA = Objet de type :                                             |
  49 : *|                                                                     |
  50 : *|       - TABLE INDICEE PAR DES OBJETS DE TYPE POINT CONTENANT        |
  51 : *|         DES ELEMENTS DE TYPE CHPOINT DANS LE CAS 3 DIMENSIONS.      |
  52 : *|         CHAQUE ELEMENT CONTIENT LE CHAMP THETA AU NOEUD DU          |
  53 : *|         FRONT DE COORDONNEES CELLES DU POINT P : TETA.P. ELLE       |
  54 : *|         EST EGALEMENT INDICEE PAR LE MOT 'GLOBAL' POUR DONNER       |
  55 : *|         LE CHAMP THETA GLOBAL LE LONG DE TOUT FRONT DE LA FISSURE   |
  56 : *|       - ELEMENT DE TYPE CHPOINT CONTENANT LE CHAMP THETA EN 2       |
  57 : *|         DIMENSIONS (OU EN 3 DIMENSIONS AVEC DES ELEMENTS DE         |
  58 : *|         COQUE MINCE) A LA POINTE DE FISSURE                         |
  59 : *|                                                                     |
  60 : *|  TABUTIL = Table avec la direction                                  |
  61 : *|                                                                     |
  62 : *|=====================================================================|
  63 : &DIME = VALE DIME ; 
  64 : &MODE = VALE MODE ;  
  65 : 
  66 : *---------------------------------------------------------*
  67 : *-------- RECUP  + TEST DE COMPABILITE DES DONNEES -------*                     
  68 : *---------------------------------------------------------*  
  69 : 
  70 : SI (NON (EXIS SUPTAB 'MAILLAGE'));                               
  71 :    MESS 'ERREUR : ON N A PAS TROUVE DANS LA' 
  72 :    MESS '         TABLE L OBJET MAILLAGE';         
  73 :    QUIT CH_THETX;                                                  
  74 : SINON;                                                                
  75 :    MAILLAGE = SUPTAB.'MAILLAGE' ;                                       
  76 : *   NB1 = NBNO (CHAN MAILLAGE 'POI1');                               
  77 : FINSI; 
  78 : SI (NON (EXIS SUPTAB 'FRONT_FISSURE'));                               
  79 :    MESS 'ERREUR : LE FRONT DE LA FISSURE N EST PAS DONNE';                   
  80 :    QUIT CH_THETX;                                                  
  81 : SINON;                                                                
  82 :    PFISS = SUPTAB.'FRONT_FISSURE';
  83 :    si(ega (type PFISS) 'POINT');
  84 :      PFISS1 = MANU 'POI1' PFISS;
  85 :    sino; 
  86 :      PFISS1 = CHAN 'POI1' PFISS;
  87 :    fins;
  88 : FINSI;                                                                
  89 : SI ((NON (EXIS SUPTAB 'PSI')) ou (NON (EXIS SUPTAB 'PHI')));                    
  90 :    MESS 'ERREUR : CHPOINT PSI et PHI NON FOURNIS';                  
  91 :    QUIT CH_THETX;                                                  
  92 : SINON;       
  93 :    psy1 = SUPTAB . 'PSI';
  94 :    phy1 = SUPTAB . 'PHI';
  95 : FINSI; 
  96 : SI (NON (EXIS SUPTAB 'COUCHE'));                                 
  97 :    MESS 'ERREUR : ON VEUT LE NOMBRE DE COUCHES D ELEMENTS';           
  98 :    MESS '         AUTOUR DE LA FISSURE QUI SE DEPLACE';               
  99 :    MESS '         POUR SIMULER LA PROPAGATION DE LA FISSURE';         
 100 :    QUIT CH_THETX;                                                  
 101 : SINON;                                                                
 102 :    COUCHE = SUPTAB.'COUCHE' ;                                           
 103 : FINSI;
 104 : 
 105 : *--------------------------------------------------*
 106 : * On veut savoir si une seule ou toutes les 2 levres
 107 : * de la fissure ont été modelisees (bp, 2012-10-04)
 108 : *--------------------------------------------------*
 109 : XMULT = 1.;
 110 : * si une demi-eprouvette est modelisee, 
 111 : * on espere que phi(x)>0 ou <0 qqsoit x
 112 : miny1 = mini phy1;  maxy1 = maxi phy1;
 113 : si ((maxy1*miny1) ' 0);                                                  
 153 :   REPE BCOUCH COUCHE ;                                            
 154 :        MBOUGER = ELEM MAILLAGE 'APPUYE' 'LARG' MBOUGER ;              
 155 :   FIN  BCOUCH ;                                                       
 156 :   FINSI;                                                              
 157 :   MAIL = ELEM MAILLAGE 'APPU' 'LARG' MBOUGER ;
 158 :   mod7 = MODE MAIL 'MECANIQUE' 'ELASTIQUE';
 159 :   
 160 :   
 161 : *--- Vecteur direction unitaire ------*
 162 : 
 163 :   lv7 = (NOMC 'UX' psy1 'NATU' 'DIFFUS') ET
 164 :         (NOMC 'UY' phy1 'NATU' 'DIFFUS') ;
 165 :   glv7 = GRAD mod7 lv7 ;
 166 :   
 167 : * calcul basé sur grad de psi ou sur grad de phi ou les 2 ou autre...?
 168 : * question delicate puisqu'elle met en jeu le pb des fissures courbes 
 169 : * et qui branchent.
 170 : * grad(psi) est + facile a metter en oeuvre, notamment si presence de 2
 171 : * pointes car le repere grad(psi), grad(phi) n'est pas direct des 2 cotés,
 172 : * mais le vecteur direction ne semble pas tourner facilement
 173 :   teta7a= EXCO glv7 (mots 'UX,X' 'UX,Y') (mots 'UX' 'UY');
 174 :   teta7a= CHAN teta7a 'TYPE' 'SCALAIRE';
 175 : *  vteta7a = vect (chan chpo mod7 teta7a) 1.E-2 ROSE;
 176 : * avec grad(phi) il y a 2 produits vectoriel a faire
 177 :   teta7b= EXCO glv7 (mots 'UY,X' 'UY,Y') (mots 'UX' 'UY');
 178 :   teta7b= CHAN teta7b 'TYPE' 'SCALAIRE';
 179 :   taw1 = ((exco teta7a 'UX' 'UZ') * (exco teta7b 'UY' 'UZ'))
 180 :        - ((exco teta7a 'UY' 'UZ') * (exco teta7b 'UX' 'UZ'));
 181 :   taw1 = CHAN  taw1 'TYPE' 'SCALAIRE';
 182 :   teta7c=((exco teta7b 'UY' 'UX') * (exco taw1   'UZ' 'UX'))
 183 :       et (-1.* ((exco teta7b 'UX' 'UY') * (exco taw1   'UZ' 'UY'))); 
 184 : *        - ((exco teta7b 'UX' 'UY') * (exco taw1   'UZ' 'UY'));  bug !!!
 185 : *  vteta7c = vect (chan chpo mod7 teta7c) (mots 'UX' 'UY') 1.E-2 BLEU;
 186 : 
 187 :   teta7 = teta7c;
 188 : * passage au noeud + normalisation du chamelem
 189 :   teta7 = CHAN teta7 mod7 'NOEUD';
 190 :   nteta7 = (PSCA teta7 teta7 (mots 'UX' 'UY') (mots 'UX' 'UY'))**(-0.5);
 191 : *  mess (maxi nteta7) (mini nteta7);
 192 :   teta7 = teta7 * nteta7;
 193 :   TETA = CHAN 'CHPO' mod7 teta7 'MOYE';
 194 :   TETA = REDU TETA MBOUGER;
 195 :   
 196 : * normalisation du chpoint + elargissement à MAIL
 197 :   nTETA = (PSCA TETA TETA (mots 'UX' 'UY') (mots 'UX' 'UY'))**(-0.5);
 198 :   TETA = (TETA * (nTETA * XMULT))
 199 :        + (MANU 'CHPO' MAIL 2 'UX' 0. 'UY' 0. 'NATURE' 'DIFFUS') ; 
 200 : *  vteta = vect TETA 1.E-2 BLEU;  trac vteta MAIL;
 201 :            
 202 :            
 203 : *------ DIRECTIONs a garder dans la table TABUTIL ------*
 204 : 
 205 :   VECTEUR1 = PROI PFISS1 teta7 ;
 206 : *on utilise desormais un CHPOINT et plus un POINT
 207 :     TABUTIL .'DIRECTION1' = VECTEUR1; 
 208 : *on calule systematiquement DIRCISA  = DIRTANG PVEC DIRTETA  
 209 :     VECTEUR2 =  (-1.*(EXCO VECTEUR1 'UY' 'UX')) 
 210 :                   et (EXCO VECTEUR1 'UX' 'UY');
 211 :     TABUTIL .'DIRECTION2' =  VECTEUR2;
 212 :     
 213 : FINSI;
 214 : * FIN DU CAS 2D ---------------------------------*
 215 : *------------------------------------------------*
 216 : 
 217 : 
 218 : 
 219 : *------------------------------------------------*
 220 : * CAS 3D ----------------------------------------*
 221 : SI (&DIME EGA 3);
 222 : *    MESS 'ERREUR : CH_THETX NON PREVU EN 3D POUR L INSTANT';           
 223 : *    MESS 'CH_THETX EN PHASE DE TEST EN 3D POUR L INSTANT';           
 224 : 
 225 : *------ EFISS0 = ELEMENTS CONTENANT LE FRONT DE FISSURE ------*
 226 : *
 227 :   MAILPSIP = EXTR psy1 'MAILLAGE';
 228 :   MAILPSIP = ELEM MAILLAGE 'APPUYE' 'STRICTEMENT' MAILPSIP;
 229 :   psy1e = CHAN 'CHAM' psy1 MAILPSIP;
 230 :   phy1e = CHAN 'CHAM' phy1 MAILPSIP;
 231 :   EFISS0 = ((ELEM psy1e 'EGSUPE' 0.) INTE (ELEM psy1e 'EGINFE' 0.))
 232 :       INTE ((ELEM phy1e 'EGSUPE' 0.) INTE (ELEM phy1e 'EGINFE' 0.));
 233 :   DX0 = ((MESU EFISS0 'VOLUME') / (nbel EFISS0)) ** (1./3.);
 234 :   EPS0P = 1.E-4 * DX0;
 235 :   EPS0N = -1.E-4 * DX0;
 236 :   EFISS0= ((ELEM psy1e 'EGSUP' EPS0N) INTE (ELEM psy1e 'EGINF' EPS0P))
 237 :      INTE ((ELEM phy1e 'EGSUP' EPS0N) INTE (ELEM phy1e 'EGINF' EPS0P));
 238 :           
 239 : *------ Recup de MBOUGER et de MAIL ------*
 240 : *
 241 :   MBOUGER = EFISS0;
 242 : * boucle sur les couches
 243 :   SI (COUCHE '>' 1);                                                  
 244 :   REPE BCOUCH (COUCHE - 1) ;                                            
 245 :        MBOUGER = ELEM MAILLAGE 'APPUYE' 'LARG' MBOUGER ;              
 246 :   FIN  BCOUCH ;
 247 :   FINSI;      
 248 :       
 249 : * MAIL = maillage de definition de TETA 
 250 :   MAIL = ELEM MAILLAGE 'APPU' 'LARG' MBOUGER ;
 251 : * inutile de chercher a calculer + loin que le domaine de def de psi et phi  
 252 :   MAIL = MAIL inte MAILPSIP;
 253 :   mod7 = MODE MAIL 'MECANIQUE' 'ELASTIQUE' ;
 254 : 
 255 : *--- Travail sur les level set 
 256 : *    pour definir le repere local du front de fissure ------*
 257 : *--- Creation du TETA global ------*
 258 : *
 259 : * Vecteur direction unitaire
 260 :   lv7 = (NOMC 'UX' psy1 'NATU' 'DIFFUS') ET
 261 :         (NOMC 'UY' phy1 'NATU' 'DIFFUS') ET
 262 :         (MANU 'CHPO' MAIL 1 'UZ' 0. 'NATU' 'DIFFUS');
 263 :   glv7 = GRAD mod7 lv7 ;
 264 :   gpsy1= EXCO glv7 (mots 'UX,X' 'UX,Y' 'UX,Z') (mots 'UX' 'UY' 'UZ');
 265 :   gpsy1= CHAN gpsy1 'TYPE' 'SCALAIRE';
 266 : *   gpsy1= CHAN gpsy1 'TYPE' 'DEPLACEMENT';
 267 :   gphy1s= EXCO glv7 (mots 'UY,X' 'UY,Y' 'UY,Z') (mots 'UX' 'UY' 'UZ');
 268 :   gphy1s= CHAN gphy1s 'TYPE' 'SCALAIRE';
 269 : *   gphy1= CHAN gpsy1 'TYPE' 'DEPLACEMENT';
 270 : *  mess '* taw1 = (grad PSI) pvec (grad PHI)';
 271 : * ~(grad PSI) mais pas tout a fait...
 272 :   taw1s = ( ((exco gpsy1 'UY' 'UX') * (exco gphy1s 'UZ' 'UX'))
 273 :           - ((exco gpsy1 'UZ' 'UX') * (exco gphy1s 'UY' 'UX')) )
 274 :        et ( ((exco gpsy1 'UZ' 'UY') * (exco gphy1s 'UX' 'UY'))
 275 :           - ((exco gpsy1 'UX' 'UY') * (exco gphy1s 'UZ' 'UY')) )
 276 :        et ( ((exco gpsy1 'UX' 'UZ') * (exco gphy1s 'UY' 'UZ'))
 277 :           - ((exco gpsy1 'UY' 'UZ') * (exco gphy1s 'UX' 'UZ')) );
 278 :   taw1s = CHAN  taw1s 'TYPE' 'SCALAIRE';
 279 :   teta1s = ( ((exco gphy1s 'UY' 'UX') * (exco taw1s 'UZ' 'UX'))
 280 :            - ((exco gphy1s 'UZ' 'UX') * (exco taw1s 'UY' 'UX')) )
 281 :         et ( ((exco gphy1s 'UZ' 'UY') * (exco taw1s 'UX' 'UY'))
 282 :            - ((exco gphy1s 'UX' 'UY') * (exco taw1s 'UZ' 'UY')) )
 283 :         et ( ((exco gphy1s 'UX' 'UZ') * (exco taw1s 'UY' 'UZ'))
 284 :            - ((exco gphy1s 'UY' 'UZ') * (exco taw1s 'UX' 'UZ')) );
 285 :            
 286 : * PASSAGE AU NOEUD DES CHAMELEMS (normalisation inutile?)
 287 :   teta1 = CHAN teta1s mod7 'NOEUD';
 288 : *   nteta1 = (PSCA teta1 teta1 
 289 : *   (mots 'UX' 'UY' 'UZ') (mots 'UX' 'UY' 'UZ'))**(-0.5);
 290 : *   teta1 = teta1 * nteta1;
 291 :   gphy1= CHAN gphy1s mod7 'NOEUD';
 292 :   taw1 = CHAN taw1s mod7 'NOEUD';
 293 : *   ntaw1 = (PSCA taw1  taw1
 294 : *   (mots 'UX' 'UY' 'UZ') (mots 'UX' 'UY' 'UZ'))**(-0.5);
 295 : *   taw1 = taw1 * ntaw1;
 296 : 
 297 : * CREATION + NORMALISATION DES CHPOINTS 
 298 : * + REDU a  MBOUGER + (re)-elargissement à MAIL
 299 : * TETA
 300 :   TETA = CHAN 'CHPO' mod7 teta1 'MOYE';
 301 :   nTETA = (PSCA TETA TETA 
 302 :   (mots 'UX' 'UY' 'UZ') (mots 'UX' 'UY' 'UZ'))**(-0.5);
 303 :   TETA = chan (TETA * nTETA) 'ATTRIBUT' 'NATURE' 'DIFFUS';
 304 :   TABUTIL . 'V1' = TETA;
 305 :   TETA = (REDU TETA MBOUGER);
 306 : *   + (MANU 'CHPO' MAIL 3 'UX' 0. 'UY' 0. 'UZ' 0. 'NATURE' 'DIFFUS') ; 
 307 : * TAW
 308 :   GPHY = CHAN 'CHPO' mod7 gphy1 'MOYE';
 309 :   nGPHY = (PSCA  GPHY GPHY 
 310 :       (mots 'UX' 'UY' 'UZ') (mots 'UX' 'UY' 'UZ'))**(-0.5);
 311 :   GPHY = chan (GPHY * nGPHY) 'ATTRIBUT' 'NATURE' 'DIFFUS';
 312 :   TABUTIL . 'V2' = GPHY;
 313 :   GPHY = (REDU GPHY MBOUGER);
 314 : *   + (MANU 'CHPO' MAIL 3 'UX' 0. 'UY' 0. 'UZ' 0. 'NATURE' 'DIFFUS') ; 
 315 : * TAW
 316 :   TAW = CHAN 'CHPO' mod7 taw1 'MOYE';
 317 :   nTAW = (PSCA  TAW TAW 
 318 :       (mots 'UX' 'UY' 'UZ') (mots 'UX' 'UY' 'UZ'))**(-0.5);
 319 :   TAW = chan (TAW * nTAW) 'ATTRIBUT' 'NATURE' 'DIFFUS';
 320 :   TABUTIL . 'V3' = TAW;
 321 :   TAW = (REDU TAW MBOUGER);
 322 : *   + (MANU 'CHPO' MAIL 3 'UX' 0. 'UY' 0. 'UZ' 0. 'NATURE' 'DIFFUS') ; 
 323 : *
 324 : *  trac ((vect TETA 0.05 'BLEU') et (vect TAW  0.05 'VERT')) 
 325 : *       ((aret MAIL) et PFISS);
 326 :       
 327 :       
 328 : *------ DIRECTIONs a garder dans la table TABUTIL ------*
 329 : 
 330 : * repere local du front
 331 : *   gphy7 = PROI gphy1 mod7 PFISS1;
 332 : *   teta7 = PROI teta1 mod7 PFISS1;
 333 : *   taw7  = PROI taw1  mod7 PFISS1;
 334 : * => pas assez regulier => provoque des oscillations non symetrique 
 335 : * si 1 point du front (symetrique) appartient a la frontiere entre 2 elements
 336 :   gphy7 = INT_COMP MBOUGER GPHY PFISS1;
 337 :   teta7 = INT_COMP MBOUGER TETA PFISS1;
 338 :   taw7  = INT_COMP MBOUGER TAW  PFISS1;
 339 : * normalisation
 340 :   ngphy7 = (PSCA gphy7  gphy7
 341 :   (mots 'UX' 'UY' 'UZ') (mots 'UX' 'UY' 'UZ'))**(-0.5);
 342 :   gphy7 = gphy7 * ngphy7;
 343 :   nteta7 = (PSCA teta7  teta7
 344 :   (mots 'UX' 'UY' 'UZ') (mots 'UX' 'UY' 'UZ'))**(-0.5);
 345 :   teta7 = teta7 * nteta7;
 346 :   ntaw7 = (PSCA taw7  taw7
 347 :   (mots 'UX' 'UY' 'UZ') (mots 'UX' 'UY' 'UZ'))**(-0.5);
 348 :   taw7 = taw7 * ntaw7;  
 349 : *on utilise desormais des CHPOINT et plus des POINTs
 350 :   TABUTIL .'DIRECTION1' = chan teta7 'ATTRIBUT' 'NATURE' 'DIFFUS'; 
 351 :   TABUTIL .'DIRECTION2' = chan gphy7 'ATTRIBUT' 'NATURE' 'DIFFUS';
 352 :   TABUTIL .'DIRECTION3' = chan taw7 'ATTRIBUT' 'NATURE' 'DIFFUS';
 353 : 
 354 : * mess '*------ identification de chaque element du front ------*';
 355 : *------ identification de chaque element du front ------*
 356 : *
 357 : * pour identifier les elem a bouger, on boucle sur les points du front  
 358 :   nefiss1 = 0;
 359 :   EFISSi = tabl;
 360 :   nPFISS1 = NBEL PFISS1;
 361 :   ipfiss1 = 0;
 362 : * boucle sur les points du front
 363 :   repe BPFISS1 nPFISS1;
 364 :     ipfiss1 = ipfiss1 + 1;
 365 :     PFISS1i = POIN PFISS1 ipfiss1;
 366 :     
 367 : * mess 'ch_thetx: ' ipfiss1 ' -> noeud' (noeu PFISS1i) 
 368 : *      'on cherche le(s) element(s)' (nefiss1 + 1);
 369 :     
 370 : *   on recupere le(s) element(s) contenant le ieme point 
 371 : *     Etesti = MBOUGER ELEM 'CONTENANT' PFISS1i;
 372 :     Etesti = MBOUGER ELEM 'CONTENANT' PFISS1i 'TOUS';
 373 : *     list Etesti;
 374 :     
 375 :     si(ega ipfiss1 1); 
 376 :       nefiss1 = nefiss1 + 1;
 377 :       EFISSi . nefiss1 = Etesti;
 378 : * mess '1ere fois ->' nefiss1;
 379 :     sino;
 380 :     
 381 : *       Etesti = diff Etesti (inte Etesti (EFISSi . nefiss1) 'NOVERIF');
 382 : *       si(neg (nbel Etesti) 0);
 383 : *         nefiss1 = nefiss1 + 1;
 384 : *         EFISSi . nefiss1 = Etesti;
 385 : *       fins;
 386 : 
 387 : *     de base on enleve l eventuel element n-2
 388 :       si(nefiss1 >eg 2);
 389 :         Einte2 = inte Etesti (EFISSi . (nefiss1-1)) 'NOVERIF';
 390 :         si(neg (nbel Einte2) 0);
 391 :           Etesti = diff Etesti Einte2;
 392 :           si((nbel Etesti) ega 0); 
 393 : * mess 'tous les elements sont deja dans' (nefiss1-1) 'on itere' ;
 394 :             iter BPFISS1; 
 395 :           fins;
 396 :         fins;
 397 :       fins;
 398 :       
 399 : *     puis on fait le travail      
 400 :       Einte = inte Etesti (EFISSi . nefiss1) 'NOVERIF';
 401 : * mess 'on a' (nbel Einte) 'elements en commun avec' (nefiss1);
 402 :       si(ega (nbel Einte) 0);
 403 :         nefiss1 = nefiss1 + 1;
 404 :         EFISSi . nefiss1 = Etesti;
 405 : * mess '                ->' nefiss1;
 406 :       sino;
 407 : *     il existe des elements en commun
 408 : *        on a peut etre trop pris d element dans l ancien... 
 409 :          Etemp = diff (EFISSi . nefiss1) Einte;
 410 :          si(neg (nbel Etemp) 0);
 411 : * mess ' on parvient a enlever des elements de lancien' nefiss1 
 412 : *      ' qui existe encore' (nbel Etemp);
 413 : *          on enleve ce qu il y a en trop et on met ce qui reste dans le nouveau
 414 :            EFISSi . nefiss1 = Etemp;
 415 : *            Etesti = diff Etesti Etemp ;
 416 :              Einte = inte Etesti Etemp 'NOVERIF';
 417 :              Etesti = diff Etesti Einte;
 418 :            si(neg (nbel Etesti) 0);
 419 :              nefiss1 = nefiss1 + 1;
 420 :              EFISSi . nefiss1 = Etesti;
 421 : * mess ' on stocke ce qui reste (soit ' (nbel Etesti) ' -> ' nefiss1;
 422 :            finsi;
 423 :          sino;
 424 : *          on a peut etre trop pris d element dans le nouveau...
 425 :            Etemp = diff Etesti Einte;         
 426 :            si(neg (nbel Etemp) 0);
 427 :              nefiss1 = nefiss1 + 1;
 428 :              EFISSi . nefiss1 = Etemp;           
 429 : * mess ' on parvient a enlever des elements du nouveau'  
 430 : *      ' qui existe encore' (nbel Etemp) ' -> ' nefiss1;
 431 :            fins;
 432 :          fins;
 433 :          
 434 :       fins;
 435 :     fins;
 436 : *     list EFISSi;
 437 :   fin BPFISS1;
 438 : 
 439 : 
 440 : *------ creation des TETA de chaque element du front ------*
 441 : *
 442 :   TTETA = TABL;
 443 :   coef2 = 1.E-2 * DX0;
 444 : * boucle sur les elements du front 
 445 :   iefiss1 = 0;
 446 :   REPE BEFISS2 nefiss1; 
 447 :     iefiss1 = iefiss1 + 1;
 448 :     
 449 : *   chpoint TETA2 de chaque element contenant le front (=tranche)
 450 :     TETA2 = REDU TETA  EFISSi . iefiss1;
 451 :     
 452 : *   on replace les COUCHES ici en faisant une suite
 453 :     MBOUGERi = EFISSi . iefiss1;
 454 :     SI (COUCHE '>' 1);
 455 :     XCOUCH = 1.;
 456 :     REPE BCOUCH (COUCHE - 1) ;
 457 :          XCOUCH = XCOUCH + 1.;
 458 :          MBOUGERi = ELEM MAILLAGE 'APPUYE' 'LARG' MBOUGERi;
 459 :          TETA2 = TETA2 + (REDU TETA MBOUGERi) ;
 460 : *bp_exp         TETA2 = TETA2 + ((1./XCOUCH) * (REDU TETA MBOUGERi)) ;
 461 :     FIN  BCOUCH ;
 462 : *bp_lin     TETA2 = TETA2 / (COUCHE - 1);
 463 :     FINSI;  
 464 : *bp_+loin    TETA2 = REDU TETA MBOUGERi;     
 465 :     
 466 : *      trac (vect TETA2 DEPL BLEU) EFISS0;
 467 : *     TTETA . (EFISSi . iefiss1) = TETA2;
 468 : *   avancee VECTEUR2 de chaque element du front (=segment)
 469 : *    teta72 = PROI PFISS (redu teta7 EFISSi . iefiss1);
 470 : *    teta72 = INT_COMP EFISSi . iefiss1 TETA2 PFISS;
 471 :     teta72 = INT_COMP EFISS0 TETA2 PFISS;
 472 : *      trac (vect teta72 DEPL BLEU 0.03) ( EFISS0 et PFISS);
 473 : *   calcul de l'aire fracturee par cette avancee virtuelle
 474 : *rem: si discretisation non conforme erreur/XAFISS2 peut etre importante
 475 :     PFISS2 = PFISS PLUS (teta72 * coef2);
 476 :     AFISS2 = PFISS2 REGL PFISS 1;
 477 : *     trac (vect teta72 DEPL BLEU 0.03) (EFISS0 et PFISS et AFISS2);
 478 :     XAFISS2 = (MESU AFISS2 'SURF') / coef2; 
 479 : *     mess iefiss1 'ieme element de surface = ' XAFISS2;
 480 : *     trac (vect teta72 'DEPL') (AFISS2 et PFISS);
 481 :     TTETA . (EFISSi . iefiss1) = TETA2 * (XMULT / XAFISS2);
 482 :     
 483 : *     nom0 = mots 'UX' 'UY' 'UZ';
 484 : *     nom1 = mots (chai 'UX' iefiss1) 
 485 : *                 (chai 'UY' iefiss1) (chai 'UZ' iefiss1);
 486 : *     si(iefiss1 ega 1);
 487 : *       toto = (TTETA . (EFISSi . iefiss1)) nomc nom0 nom1;
 488 : *     sinon;
 489 : *       toto = toto et 
 490 : *       ((TTETA . (EFISSi . iefiss1)) nomc nom0 nom1);
 491 : *     finsi;
 492 : *     
 493 :     
 494 :   FIN  BEFISS2; 
 495 :   
 496 : * opti sort  'THETA.inp';
 497 : * sort 'AVS' MBOUGER toto;
 498 : 
 499 : 
 500 :  
 501 : * champ TETA global
 502 :   teta72 = INT_COMP EFISS0 TETA PFISS;
 503 :   PFISS2 = PFISS PLUS (teta72 * coef2);
 504 :   AFISS2 = PFISS2 REGL PFISS 1;
 505 :   XAFISS2 = (MESU AFISS2 'SURF') / coef2; 
 506 :   TTETA . 'GLOBAL' = TETA * (XMULT / XAFISS2);
 507 : 
 508 : * pour renvoyer la table TTETA
 509 :   TETA = TTETA;
 510 :   
 511 : FINSI;
 512 : * FIN DU CAS 3D ---------------------------------*
 513 : *------------------------------------------------*
 514 : 
 515 : 
 516 : FINPROC TETA TABUTIL;
 517 : *     Fin de la procedure CH_THETX               
 518 : *--------------------------------------------*
 519 : 
 520 :  
 521 :  
 522 :  
 523 :  
 524 :  

© Cast3M 2003 - All rights reserved.
Disclaimer