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

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