Télécharger la plaquette Cast3M


fr en

Télécharger relxfe.eso

Retour à la liste

Numérotation des lignes : 

relxfe

  1. C RELXFE    SOURCE    MB234859  25/09/08    21:16:06     12358          
  2. C RELXFE    SOURCE    BP208322  17/04/18    21:15:13     9395           
  3.       SUBROUTINE RELXFE(IMODEL, MRIGID) 
  4. C***********************************************************************
  5. C   cet operateur créé d'une matrice élémentaire de rigidité 
  6. c   pour imposer à zéro les ddl xfem crée mais non actif 
  7. c   (dans les éléments de transition)
  8. C                                                                     
  9. C   ENTREES :                                                         
  10. C   ________                                                          
  11. C                                                                     
  12. C           IMODEL   pointeur sur le modele élémentaire
  13. C                              
  14. C   ENTREES/SORTIE :                                                  
  15. C   ________                                                          
  16. C                                                                     
  17. c           MRIGID  rigidité chapeu dans laquelle on va écrire 
  18. c           (dans la dernière sous matrice) la rigidité voulue                                          
  19. C---------------------------------------------------------------------   
  20. C***********************************************************************
  21.  
  22.       IMPLICIT INTEGER(I-N)
  23.       IMPLICIT REAL*8 (A-H,O-Z)
  24.       CHARACTER*8 CMATE
  25.       PARAMETER  (NDDLMAX=30,NBNIMAX=10)
  26.       INTEGER    LOCIRI(10,(1+NDDLMAX))
  27.       DIMENSION  MLRE(NDDLMAX+1)
  28.  
  29.       CHARACTER*4 MOTINC(NDDLMAX),MOTDUA(NDDLMAX)
  30.       DATA MOTINC/'UX  ','UY  ','UZ  ','AX  ','AY  ','AZ  ',
  31.      >'B1X ','B1Y ','B1Z ','C1X ','C1Y ','C1Z ','D1X ','D1Y ',
  32.      >'D1Z ','E1X ','E1Y ','E1Z ','B2X ','B2Y ','B2Z ','C2X ',
  33.      >'C2Y ','C2Z ','D2X ','D2Y ','D2Z ','E2X ','E2Y ','E2Z '/
  34.       DATA MOTDUA/'FX  ','FY  ', 'FZ  ','FAX ','FAY ','FAZ ',
  35.      >'FB1X','FB1Y','FB1Z','FC1X','FC1Y','FC1Z','FD1X','FD1Y',
  36.      >'FD1Z','FE1X','FE1Y','FE1Z','FB2X','FB2Y','FB2Z','FC2X',
  37.      >'FC2Y','FC2Z','FD2X','FD2Y','FD2Z','FE2X','FE2Y','FE2Z'/
  38.  
  39. -INC PPARAM
  40. -INC CCOPTIO
  41. -INC CCGEOME
  42. -INC SMELEME
  43. -INC SMCOORD
  44. -INC SMRIGID
  45. -INC SMMODEL
  46. -INC CCHAMP
  47. -INC SMCHAML
  48. -INC SMLREEL
  49. -INC SMINTE
  50.  
  51.  
  52.       POINTEUR    MCHEX1.MCHELM
  53.  
  54. c sement raccourcis par éléments
  55.       SEGMENT MRACC
  56.          INTEGER TLREEL(NBENRMA2,NBI)
  57.          INTEGER MELRIG(NBELEM)
  58.       ENDSEGMENT
  59.  
  60. c Segment contenant l'info DDL a mettre à 0
  61.       SEGMENT TBLOQ
  62.          INTEGER   MBLOQ(3,NBPTB)
  63.          INTEGER   NBLOQ(3)
  64.       ENDSEGMENT
  65. C++++ Recup + Activation de la geometrie ++++++++++++++++
  66.       IPT1= IMODEL.IMAMOD
  67.       SEGACT, IPT1
  68. C     nbre de noeuds par element
  69.       NBNN1 = IPT1.NUM(/1)
  70. C     nbre d elements
  71.       NBEL1 = IPT1.NUM(/2)
  72.  
  73. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  74.       MELE  = IMODEL.NEFMOD
  75.       MFR   = INFELE(13)
  76.       IELE  = INFELE(14)
  77.       NDDL  = INFELE(15)
  78.       NSTRS = INFELE(16)
  79.  
  80. C+++++nombre de sous matrice de mrigid (va être ammené a changé)
  81.       NRIGEL = MRIGID.IRIGEL(/2)
  82.  
  83.  
  84. C++++ Recup des infos d enrichissement +++++++++++++++++++
  85. c     recup du MCHEX1 d'enrichissement
  86.       NBENR1=0
  87.       MCHAM1=0
  88.       NOBMO1=IMODEL.IVAMOD(/1)
  89.       if (NOBMO1.ne.0) then
  90.        do iobmo1=1,NOBMO1
  91.         if ((TYMODE(iobmo1)).eq.'MCHAML') then
  92.          MCHEX1 = IVAMOD(iobmo1)
  93.          segact,MCHEX1
  94.           if ((MCHEX1.TITCHE).eq.'ENRICHIS') then
  95.             MCHAM1 = MCHEX1.ICHAML(1)
  96.             segact,MCHAM1
  97.             goto 1000
  98.           endif
  99.         endif
  100.        enddo
  101.        write(ioimp,*) 'Le modele est vide (absence d enrichissement)'
  102. *       return
  103.       else
  104.       write(ioimp,*) 'Aucun MCHAML enrichissement dans le Modele'
  105. *       return
  106.       endif
  107.  
  108.  1000 continue
  109.  
  110. c     niveau d enrichissement(s) du modele (ddls std u exclus)
  111. c     NBENR1= 0 si std, 1 si H seul,  2 si H+F1,  3 si H+F1+F2
  112.       if (MCHAM1.ne.0) NBENR1=MCHAM1.IELVAL(/1)
  113. C++++ INITIALISATIONS...
  114. C
  115. c     ... des tables LOCIRI et MLRE
  116. c     MLRE contient le nombre d'inconnues de chaque sous-zone
  117. c     determinee depuis le nombre de fonctions de forme
  118. c ienr=  1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  119.       if (NBENR1.ne.0) then
  120.       do ienr=1,NBENR1
  121.          nbniJ = 2 + ((ienr-1)*4)
  122.          MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  123. c         write(*,*) 'ienr', ienr, 'mlre', MLRE(1+ienr)
  124. c        -LOCIRI
  125.          LOCIRI(5,1+ienr)= NIFOUR
  126.       enddo
  127.       endif
  128. C     Tables + longues car 1er indice correspond au fontion de forme std
  129.       MLRE(1)    = IDIM*NBNN1*1
  130.       LOCIRI(5,1)= NIFOUR
  131. c     on complete avec des 0
  132.       if (NBENR1.lt.(NDDLMAX+1)) then
  133.         do ienr=(NBENR1+1),(NDDLMAX)
  134.            MLRE(1+ienr) = 0
  135.         enddo
  136.       endif
  137. c
  138. c      ...DU SEGMENT MRACC
  139.  
  140.  
  141.       NBENRMA2 = NDDLMAX
  142.       NBI = NBNN1
  143.       NBELEM = NBEL1
  144.       segini , MRACC
  145.  
  146. C initialisation du tableau des ddl a mettre à zéro
  147.       SEGACT,MCOORD*MOD
  148.       NBPTB= nbpts
  149.       SEGINI,TBLOQ
  150.  
  151. C++++ TBLOQ.MBLOQ(NBENRJ, INUM) = faut il mettre a 0 les ddl ++++
  152. C     =0 si pas encore passé dans les ddl
  153. C     =1 si déja passé dans les ddl pas de mise à 0
  154. c     =2 si déja passé dans les ddl mise à 0 nécéssaire
  155. C ++++ TBLOQ.NBLOQ(NBENRJ) Compteur de ddl de type NBENRJ à mettre à 0
  156.  
  157.  
  158. C*********************************************************
  159. C
  160. C>>>>>>>>>>>>>>>>>>>>>>>>>>> 1ere BOUCLE SUR LES ELEMENTS
  161. C
  162.       NBENR = 0
  163.       DO 2000 J=1,NBEL1
  164.  
  165. C
  166. C++++ NBENRJ = niveau d enrichissement de l element ++++
  167. C     =0 si EF std   =1 si U+H   =2 si U+H+F1   =3 si U+H+F1+F2
  168.           NBENRJ=0
  169.           if (NBENR1.ne.0) then
  170.              do IENR=1,NBENR1
  171.                 MELVA1 = MCHAM1.IELVAL(IENR)
  172.                 segact,MELVA1
  173.                 do I=1,NBNN1
  174.                   mlree1 = MELVA1.IELCHE(I,J)
  175. C        on en profite pour remplir MRACC table de raccourcis pour cet element
  176.                  TLREEL(IENR,I) = mlree1
  177.                   if(mlree1.ne.0)  NBENRJ=max(NBENRJ,IENR)
  178.                 enddo
  179.             enddo
  180.           endif
  181.           NBENR=max(NBENRJ,NBENR)
  182.           NDDLE = MLRE(NBENRJ+1)
  183. c          if (NBENRJ.ne.0) then
  184. c          write(*,*) '***********************************************'
  185. c          write(*,*) 'ELEMENT', J
  186. c          write(*,* ) 'Niveau d enrichssement', NBENRJ            
  187. c          write(*,* ) 'Nb de ddl', NDDLE
  188. c          endif
  189. C*********************************************************
  190. C
  191. C>>>>>>>>>>>>>>>>>>>>>>>>>>>  BOUCLE SUR LES DDL De L' ELEMENT
  192. C
  193.           DO 3000 II=1,NDDLE
  194.  
  195.  
  196. C**********************************************************************
  197. C On cherche les noeuds qui ne sont pas effectivement enrichis 
  198. C pour forcer à 0 les DDL correspondants
  199. C**********************************************************************
  200.  
  201. c     on trouve le type de fonction de ce ddl:  IENR=0 si U, =1 si A,
  202. c     =2 si B1, =3 si C1, = 4 si D1, =5 si E1, =6 si B2, ... =9 si E2
  203.            IENR = ((II-1)/IDIM) / NBNN1
  204.  
  205. c     on trouve le niveau d enrichissement KENR de ce ddl si nonstd
  206.            if (IENR.eq.0) then
  207.              go to 3001
  208.            elseif(IENR.ge.2) then
  209. C           write(*,*) 'DDL' , II, 'Niveu d enr', IENR
  210.              KENR = ((IENR - 2) / 4) + 2
  211. c ci dessus: 4 represente le nombre de fonction de la base d'enrichissement
  212. c            et 2 est le decalage du a U et H
  213.            else
  214.              KENR = IENR
  215.            endif
  216. c     on trouve le noeud correspondant au ddl
  217. CYT        INODE = ((II+1)/IDIM) - ((IENR)*NBNN1)
  218.            INODE = 1 + ((II-1)/IDIM) - ((IENR)*NBNN1)
  219. c     numero global du noeud
  220.            INUM = IPT1.NUM(INODE , J)
  221. c     est ont déja passé dans ce ddl ? 
  222. c           write(*,*) 'INUM', INUM, 'Kenr', KENR 
  223. c           write(*,*) 'mlree1', mlree1,'Mbloq',Tbloq.Mbloq(KENR, INUM)
  224.            if (Tbloq.Mbloq(KENR, INUM).gt.0) GOTO 3001
  225. c     est-ce un noeud vraiment enrichi?
  226.            mlree1 = TLREEL(KENR,INODE)
  227.            Tbloq.Mbloq(KENR, INUM)=Tbloq.Mbloq(KENR, INUM)+1
  228.            if (mlree1.eq.0)  then
  229.               Tbloq.Mbloq(KENR, INUM)=Tbloq.Mbloq(KENR, INUM)+1
  230.               Tbloq.Nbloq(KENR) = Tbloq.Nbloq(KENR)+1
  231.            endif
  232. c           write(*,*) 'Nouveau Mbloq', Tbloq.Mbloq(KENR, INUM)
  233.  
  234.  
  235.  3001 CONTINUE
  236.  3000 CONTINUE
  237. c
  238. c
  239.  2000 CONTINUE
  240.  
  241. C FIN DE BOUCLE SUR LES ELEMENTS 
  242.  
  243. C*********************************************************
  244. C       Creation des matrices de bloquage 
  245. C       pour les DDL non effectivement enrichis
  246. C*********************************************************
  247.  
  248.       NLIGRD = 2
  249.       NLIGRP = 2
  250.  
  251.       NBNN   = 2
  252.       NBPTS1=NBPTB
  253.  
  254. C*********************************************************
  255. C
  256. C>>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les types d'enrichissement  
  257. C
  258. C*********************************************************
  259.  
  260.       IDDL = 3
  261.       DO 4000 IENR =1, NBENR1
  262. C Maillage des noeuds à bloquer
  263. C nombre de noeuds a bloquer -> nombre d'éléments du maillage de blocage 
  264.           if (IENR.eq.1) then
  265.               NBELEM = TBLOQ.NBLOQ(1)
  266.           elseif (IENR.EQ.2) then
  267.               NBELEM = TBLOQ.NBLOQ(2)
  268.           else
  269.               NBELEM = TBLOQ.NBLOQ(3)
  270.           endif
  271.           NFON=1
  272.           IF (IENR.gt.1) NFON=4
  273. C*********************************************************
  274. C
  275. C>>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les fonctions de formes d'enrichissement
  276. C
  277. C*********************************************************
  278.           DO 4001 IFON = 1 , NFON
  279.  
  280. C*********************************************************
  281. C
  282. C>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les composantes                
  283. C
  284. C*********************************************************
  285.  
  286.               DO 4002 ICOMP = 1, IDIM 
  287.                   IDDL=IDDL+1
  288. C si aucun noeud a bloquer on saute ce type d'enrichissement
  289.                   if (NBELEM.EQ.0) goto 4000 
  290.  
  291. C*********************************************************
  292. C Maillage du  blocage
  293.                   NBREF = 0
  294.                   NBSOUS = 0
  295.                   SEGINI IPT2
  296.                   IPT2.ITYPEL = 22
  297.                   Nelem = 0
  298.  
  299. C ajustement de XMCOORD pour ajouter les noeud des multiplicateurs
  300.  
  301.                   NBPTS = NBPTS1 + NBELEM
  302.                   SEGADJ MCOORD
  303.  
  304.  
  305. C>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les noeuds  
  306.                   DO 4010 INUM = 1, NBPTB 
  307.                       IF (TBLOQ.MBLOQ(IENR, INUM).LT.2) goto 4010
  308.                       Nelem = Nelem + 1
  309.                       IPT2.NUM(2,Nelem)=INUM
  310.                       NBPTS1 = NBPTS1 + 1
  311.                       IPT2.NUM(1,Nelem)=NBPTS1
  312.                       IPT2.icolor(Nelem)=1
  313. 4010              CONTINUE
  314. C Fin de boucle sur les noeuds   
  315.  
  316. C coordonées des nouveaux noeud
  317.  
  318. C>>>>>>>>>>>>>>>>>>>>>>>>>> Boucle sur les éléments de ipt2
  319.                   DO 4011 IELE = 1, NBELEM
  320.                       INUM2 = IPT2.NUM(2,IELE)
  321.                       NBPTS2 = NBPTS1 - NBELEM + IELE
  322.                       INEW = (IDIM+1)*(NBPTS2-1)                      
  323.                       IOLD = (IDIM+1)*(INUM2-1)
  324.                       DO ID = 1, IDIM
  325.                         XCOOR(INEW +ID ) =  XCOOR(IOLD +ID )
  326.                       ENDDO
  327. 4011              CONTINUE
  328.  
  329.  
  330. C********************************************************* 
  331. C Matrice de blocage
  332.                   NELRIG = NELEM             
  333.                   SEGINI XMATR1 
  334.                   DO i=1,NELRIG
  335.                       XMATR1.RE(1,1,i)=0.D0
  336.                       XMATR1.RE(2,1,i)=1.D0
  337.                       XMATR1.RE(2,2,i)=0.D0
  338.                       XMATR1.RE(1,2,i)=1.D0
  339.                   ENDDO  
  340.  
  341. C*********************************************************  
  342. C Segment Descripteur
  343.                   SEGINI DES1
  344.                   DES1.LISINC(1)='LX'
  345.                   DES1.LISDUA(1)='FLX'
  346.                   DES1.LISINC(2)=MOTINC(IDDL)
  347.                   DES1.LISDUA(2)=MOTDUA(IDDL)  
  348.  
  349.                   DES1. NOELEP(1)=1
  350.                   DES1. NOELEP(2)=2
  351.                   DES1. NOELED(1)=1
  352.                   DES1. NOELED(2)=2
  353.  
  354. C*********************************************************  
  355. C stockage de la rigidité construite dans MRIGID
  356.  
  357.  
  358.                   NRIGEL = NRIGEL+1
  359.                   SEGADJ MRIGID
  360.                   SEGACT, MRIGID*MOD 
  361.                   MRIGID.IRIGEL(1, NRIGEL)= IPT2   
  362.                   MRIGID.IRIGEL(3, NRIGEL)=DES1 
  363.                   MRIGID.IRIGEL(4, NRIGEL)=XMATR1
  364.                   MRIGID.COERIG(NRIGEL)=1.D0
  365.  
  366. c desactivations 
  367.                  SEGDES XMATR1, DES1
  368. c                 write(*,*) 'SOUS MATRICE', NRIGEL, MOTINC(IDDL)
  369. c                 WRITE(*,*) 'NB de blocages ', NBELEM   
  370. c                 WRITE(*,*) 'maillage ', IPT2 
  371.  4002         CONTINUE
  372.               IF (IDIM.EQ.2) IDDL=IDDL+1
  373.  4001     CONTINUE     
  374.  4000 CONTINUE
  375. C      write (*,*) '***************NRIGEL*******************', NRIGEL
  376. C Fin des boucles sur les niveau d'enrichissement et composantes 
  377.  
  378.  
  379. c
  380. C*********************************************************
  381. C       SUPPRESSION ET DESACTIVATION DE SEGMENTS
  382. C*********************************************************
  383.  
  384.       SEGSUP,MRACC,TBLOQ
  385.  
  386.       segdes, MCOORD
  387.  
  388.       END
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  

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