Télécharger la plaquette Cast3M


fr en

Télécharger massx.eso

Retour à la liste

Numérotation des lignes : 

massx

  1. C MASSX     SOURCE    MB234859  25/09/08    21:15:53     12358          
  2.  
  3.       subroutine MASSX (ivamat,ivacar,NMATT,CMATE,
  4.      &   IMODEL,LOCIRI)
  5. c
  6. C     PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM
  7. c     POUR LE CALCUL DE RIGIDITE ELEMENTAIRES
  8. C
  9. C
  10. C*********************************************************
  11. C       PARTIE DECLARATIVE
  12. C*********************************************************
  13.  
  14.       IMPLICIT INTEGER(I-N)
  15.       IMPLICIT REAL*8(A-H,O-Z)
  16.  
  17. -INC PPARAM
  18. -INC CCOPTIO
  19.  
  20. -INC SMCOORD
  21. -INC SMMODEL
  22. -INC SMCHAML
  23. -INC SMINTE
  24. -INC SMELEME
  25. -INC SMRIGID
  26. -INC SMLREEL
  27.  
  28.       POINTEUR    MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE
  29.  
  30. -INC TMPTVAL
  31.  
  32.       CHARACTER*8 CMATE
  33. CYT   PARAMETER  (NDDLMAX=20,NBNIMAX=10)
  34.       PARAMETER  (NDDLMAX=30,NBNIMAX=10)
  35.       CHARACTER*4 MOTINC(NDDLMAX),MOTDUA(NDDLMAX)
  36. CYT   DATA MOTINC/'UX  ','UY  ','AX  ','AY  ',
  37. CYT  >'B1X ','B1Y ','C1X ','C1Y ','D1X ','D1Y ','E1X ','E1Y ',
  38. CYT  >'B2X ','B2Y ','C2X ','C2Y ','D2X ','D2Y ','E2X ','E2Y '/
  39. CYT   DATA MOTDUA/  'FX  ','FY  ','FAX ','FAY ',
  40. CYT  >'FB1X','FB1Y','FC1X','FC1Y','FD1X','FD1Y','FE1X','FE1Y',
  41. CYT  >'FB2X','FB2Y','FC2X','FC2Y','FD2X','FD2Y','FE2X','FE2Y'/
  42.  
  43.       DATA MOTINC/'UX  ','UY  ','UZ  ','AX  ','AY  ','AZ  ',
  44.      >'B1X ','B1Y ','B1Z ','C1X ','C1Y ','C1Z ','D1X ','D1Y ',
  45.      >'D1Z ','E1X ','E1Y ','E1Z ','B2X ','B2Y ','B2Z ','C2X ',
  46.      >'C2Y ','C2Z ','D2X ','D2Y ','D2Z ','E2X ','E2Y ','E2Z '/
  47.       DATA MOTDUA/'FX  ','FY  ', 'FZ  ','FAX ','FAY ','FAZ ',
  48.      >'FB1X','FB1Y','FB1Z','FC1X','FC1Y','FC1Z','FD1X','FD1Y',
  49.      >'FD1Z','FE1X','FE1Y','FE1Z','FB2X','FB2Y','FB2Z','FC2X',
  50.      >'FC2Y','FC2Z','FD2X','FD2Y','FD2Z','FE2X','FE2Y','FE2Z'/
  51. C
  52.       PARAMETER  (NBENRMAX=5)
  53. C      NBENRMAX deja defini dans rigixr.eso
  54.       INTEGER    LOCIRI(10,(1+NBENRMAX))
  55.       DIMENSION  MLRE(NBENRMAX+1)
  56.  
  57. C
  58.       SEGMENT WRK1
  59.        REAL*8 XE(3,NBBB)
  60.        REAL*8 REL(LRE,LRE),RINT(LRE,LRE)
  61.       ENDSEGMENT
  62. C
  63.       SEGMENT WRK2
  64.        REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
  65. c       REAL*8 LV7WRK(NBENRMA2,2,6,NBNO)
  66.        REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
  67.       ENDSEGMENT
  68. C
  69.       SEGMENT MVELCH
  70.          REAL*8 VALMAT(NV1)
  71.       ENDSEGMENT
  72.  
  73.       SEGMENT MRACC
  74.          INTEGER TLREEL(NBENRMA2,NBI)
  75.       ENDSEGMENT
  76. C
  77. c      write (ioimp,*) '############################'
  78. c      write (ioimp,*) '#      DEBUT DE MASSX     #'
  79. c      write (ioimp,*) '############################'
  80. C
  81. C*********************************************************
  82. c
  83. C       RECUPERATION + ACTIVATIONS + VALEURS UTILES
  84. c
  85. C*********************************************************
  86. C     sont deja actifs dans rigi1.eso :
  87. C           SEGACT MMODEL,IMODEL,MELVAL
  88. c           IVAMAT
  89.       MPTVAL=IVAMAT
  90. c
  91. C++++ Activation au cas ou ++++++++++++++++++++++++++++++
  92.       SEGACT MCOORD
  93.  
  94. C++++ Recup + Activation de la geometrie ++++++++++++++++
  95.       MELEME= IMAMOD
  96. c      SEGACT MELEME    deja actif
  97.  
  98. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  99. c   + OBTENUES PAR ELQUOI DANS MASSE1
  100. C     segment INFO OBTENUES PAR ELQUOI DANS MASSE1 mais supprimé dans MASSE1
  101.       MELE = NEFMOD
  102.       NGAU1 = INFELE(6)
  103.       LRE   = INFELE(9)
  104.       LHOOK = INFELE(10)
  105.       MINT1 = INFMOD(6)
  106.       segact,MINT1
  107.       MINT2 = INFMOD(10)
  108.       if(MINT2.ne.0)   segact,MINT2
  109.       MFR   = INFELE(13)
  110.       IELE  = INFELE(14)
  111.       NDDL  = INFELE(15)
  112.       NSTRS = INFELE(16)
  113. c      write(ioimp,*)'-> MASSX INFELE',(INFELE(iou),iou=1,16)
  114.  
  115. c   + AUTRES INFOS
  116. C     nbre de noeuds par element
  117.       NBNN1 = NUM(/1)
  118. C     nbre d elements
  119.       NBEL1 = NUM(/2)
  120.  
  121. c REM: pour se passer du dimensionnement du nbre d'enrichissement dans
  122. c      elquoi et le realiser localement , on pourrait ecrire:
  123. c      LRE = NDDLMAX*NBNN1
  124. c      NDDL= NDDLMAX
  125.  
  126. C     sous decoupage et points de Gauss de l element geometrique de base
  127. CYT   if(MELE.eq.263)  then
  128.       if(MELE.eq.263.or.MELE.eq.264)  then
  129. c         NGAU2 = 4
  130.          NGAU2 = MINT2.POIGAU(/1)
  131. c         NBSSEF = (NGAU1/NGAU2)**0.5
  132.       endif
  133. c      write(ioimp,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3))
  134. c      write(ioimp,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU)
  135. c      segdes,MINT2
  136.  
  137.  
  138. c     nbre maxi de fonction de forme par noeud (fonction std comprise)
  139. c      NBNI = NDDL/IDIM    inutile!
  140.  
  141.  
  142. C++++ Recup des infos d enrichissement +++++++++++++++++++
  143. c     recup du MCHEX1 d'enrichissement
  144.       NOBMO1 = IVAMOD(/1)
  145.       if(NOBMO1.ne.0) then
  146.        do iobmo1=1,NOBMO1
  147.         if((TYMODE(iobmo1)).eq.'MCHAML') then
  148.          MCHEX1 = IVAMOD(iobmo1)
  149.          segact,MCHEX1
  150.           if((MCHEX1.TITCHE).eq.'ENRICHIS') then
  151.             MCHAM1 = MCHEX1.ICHAML(1)
  152.             segact,MCHAM1
  153.             segdes,MCHEX1
  154.             goto 1000
  155.           endif
  156.          segdes,MCHEX1
  157.         endif
  158.        enddo
  159.        write(ioimp,*) 'Le modele est vide (absence d enrichissement)'
  160. *       return
  161.       else
  162.       write(ioimp,*) 'Aucun MCHAML enrichissement dans le Modele'
  163. *       return
  164.       endif
  165.  
  166.  1000 continue
  167. c     niveau d enrichissement(s) du modele (ddls std u exclus)
  168. c     NBENR1= 0 si std, 1 si H seul,  2 si H+F1,  3 si H+F1+F2, etc...
  169.       if(NOBMO1.ne.0) then
  170.         NBENR1= MCHAM1.IELVAL(/1)
  171.       else
  172.         NBENR1 = 0
  173.       endif
  174. c      write(ioimp,*) 'niveau d enrichissement(s) du modele',NBENR1
  175.  
  176.  
  177. C*********************************************************
  178. c
  179. c     PHASE 2 :    PREPARATION DES OBJETS RESULTATS
  180. c
  181. C*********************************************************
  182. C
  183. C++++ RECHERCHE DES NOMS D'INCONNUES ET DES DUAUX
  184. c     MOTINC et MODUA en dur pour l instant
  185.  
  186. C++++ REMPLISSAGE DE DESCR de taille maxi
  187. c     (on ajustera en enlevant si besoin)
  188.       NLIGRP = LRE
  189.       NLIGRD = LRE
  190.       SEGINI DESCR
  191.       IPDSCR = DESCR
  192. C
  193.       KINC = 0
  194. C----->> BOUCLE SUR LES fonctions de Forme
  195.       DO  IENR=1,NBNIMAX
  196. C------->> BOUCLE SUR LES NOEUDS
  197.         DO  I=1,NBNN1
  198. C--------->> BOUCLE SUR LA DIMENSION
  199.           DO  KDIM=1,IDIM
  200.             KINC = KINC + 1
  201.             LISINC(KINC) = MOTINC((3*(IENR-1))+KDIM)
  202.             LISDUA(KINC) = MOTDUA((3*(IENR-1))+KDIM)
  203. CYT         LISINC(KINC) = MOTINC((IDIM*(IENR-1))+KDIM)
  204. CYT         LISDUA(KINC) = MOTDUA((IDIM*(IENR-1))+KDIM)
  205.             NOELEP(KINC) = I
  206.             NOELED(KINC) = I
  207.           ENDDO
  208. C--------------|
  209.          ENDDO
  210. C-------------|
  211.         ENDDO
  212. C------------|
  213.  
  214.       IF(KINC.NE.LRE)   THEN
  215.         WRITE(*,*) 'IL Y A UNE ERREUR DANS DESCR'
  216.         WRITE(*,*) 'KINC , LRE = ',KINC,LRE
  217.         RETURN
  218.       ENDIF
  219. C
  220.       SEGDES DESCR
  221. C
  222. C*********************************************************
  223. C       INITIALISATIONS...
  224. C*********************************************************
  225. C
  226. c     preparation des tables avec:
  227.  
  228.       if(NBENR1.ne.0) then
  229.       do ienr=1,NBENR1
  230. c        -le nombre d'inconnues de chaque sous-zone
  231. c         determinee depuis le nombre de fonction de forme
  232. c ienr=  1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  233.          nbniJ = 2 + ((ienr-1)*4)
  234.          MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  235. c        -LOCIRI
  236.          LOCIRI(5,1+ienr)= NIFOUR
  237.       enddo
  238.       endif
  239. C     Tables + longues car 1er indice correspond au fontion de forme std
  240.       MLRE(1)    = IDIM*NBNN1*1
  241.       LOCIRI(5,1)= NIFOUR
  242.  
  243.       if(NBENR1.lt.(NBENRMAX+1)) then
  244.         do ienr=(NBENR1+1),(NBENRMAX)
  245.            MLRE(1+ienr) = 0
  246.         enddo
  247.       endif
  248. c
  249. c      ...DU SEGMENT WRK1
  250.       NBENRMA2 = NBENRMAX
  251.       NBBB = NBNN1
  252.       SEGINI,WRK1
  253.  
  254. c      ...DU SEGMENT WRK2
  255. c      NBNO = NBNI
  256.       NBNO = LRE/IDIM
  257. c     ici, pour le calcul de la masse, BGENE ne designe pas B mais N => lhook=2
  258. CTY   LHOOK=2
  259.       LHOOK=IDIM
  260.       SEGINI,WRK2
  261.  
  262. c      ...DU SEGMENT MVELCH
  263.       NV1 = NMATT
  264.       SEGINI,MVELCH
  265. C
  266. c      ...DU SEGMENT MRACC
  267.       NBENRMA2 = NBENRMAX
  268.       NBI = NBNN1
  269.       segini,MRACC
  270. C
  271. C
  272. C*********************************************************
  273. C
  274. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  BOUCLE SUR LES ELEMENTS
  275. C
  276.       DO 2000 J=1,NBEL1
  277. c      write(ioimp,*) '========= EF',J,' /NBEL1 ========='
  278. C
  279. C
  280. C*********************************************************
  281. C     POUR CHAQUE ELEMENT, ...
  282. C*********************************************************
  283. C
  284. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  285.       CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  286. C
  287. c
  288. c      CALL ZEROI(TLREEL,NBENRMA2,NBI)
  289. C++++ NBENRJ = niveau d enrichissement de l element ++++
  290. C     =0 si EF std   =1 si U+H   =2 si U+H+F1   =3 si U+H+F1+F2
  291.       NBENRJ=0
  292.       if(NBENR1.ne.0) then
  293.       do IENR=1,NBENR1
  294.         MELVA1 = MCHAM1.IELVAL(IENR)
  295.         segact,MELVA1
  296.         do I=1,NBNN1
  297.          mlree1 = MELVA1.IELCHE(I,J)
  298. C        on en profite pour remplir MRACC table de raccourcis pour cet element
  299.          TLREEL(IENR,I) = mlree1
  300.          if(mlree1.ne.0)  then
  301.             NBENRJ = max(NBENRJ,IENR)
  302. C           et on active la listreel
  303.             segact,mlree1
  304.          endif
  305.         enddo
  306.         segdes,MELVA1
  307.       enddo
  308.       endif
  309.       if(NBENRMA2.gt.NBENR1) then
  310.         do IENR2=(NBENR1+1),NBENRMA2
  311.         do I=1,NBNN1
  312.            TLREEL(IENR2,I) = 0
  313.         enddo
  314.         enddo
  315.       endif
  316. C
  317.  
  318. c++++ quelques indices pour dimensionner au plus juste
  319. c     nbre total de ddl de l'élément considéré
  320.       NLIGRD = MLRE(1+NBENRJ)
  321.       NLIGRP = MLRE(1+NBENRJ)
  322. c     nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  323.       NBNI = (MLRE(1+NBENRJ)) / IDIM
  324.  
  325. c      write(ioimp,*) 'EF',J,' NBENRJ=',NBENRJ,
  326. c     &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'
  327. c
  328.  
  329. c++++ Selon le niveau d enrichissement,
  330. c
  331. C   + On copie cet element a la geo correspondante
  332.       IPT1 = LOCIRI(1,1+NBENRJ)
  333. c     si elle n existe pas, il faut la creer
  334.       if(IPT1.eq.0) then
  335.          NBNN  =NBNN1
  336.          NBELEM=1
  337.          NBSOUS=0
  338.          NBREF=0
  339.          segini,IPT1
  340.          IPT1.ITYPEL = ITYPEL
  341.          LOCIRI(1,1+NBENRJ)=IPT1
  342. c     si elle existe, il faut l agrandir
  343.       else
  344.          NBNN  =NBNN1
  345.          NBELEM = (IPT1.NUM(/2)) + 1
  346.          NBSOUS=0
  347.          NBREF=0
  348.          segadj,IPT1
  349.       endif
  350. c     copie en cours
  351.       do I=1,NBNN1
  352.          IPT1.NUM(I,NBELEM) = NUM(I,J)
  353.       enddo
  354.  
  355. C   + On crée le descr qui va bien si il n existe pas
  356.       DES1 = LOCIRI(3,1+NBENRJ)
  357.       if(DES1.eq.0) then
  358.          segini,DES1=DESCR
  359. c      write(ioimp,*)'        on ne garde que les',MLRE(1+NBENRJ),
  360. c     &' 1eres inconnues'
  361. c         NLIGRP = MLRE(1+NBENRJ)   fait + haut (ligne 313)
  362. c         NLIGRD = MLRE(1+NBENRJ)
  363.          segadj,DES1
  364.          LOCIRI(3,1+NBENRJ) = DES1
  365.          segdes,DES1
  366.       endif
  367.  
  368. C   + On en profite pour créer IMATRI si il n existe pas
  369.       xMATR1 = LOCIRI(4,1+NBENRJ)
  370.       if(xMATR1.eq.0) then
  371.          NELRIG = 1
  372.          segini,xMATR1
  373.          LOCIRI(4,1+NBENRJ) = xMATR1
  374. c     si il existe, il faut l agrandir
  375.       else
  376.          NELRIG = xmatr1.RE(/3) + 1
  377.          segadj,xMATR1
  378.       endif
  379.  
  380.  
  381. C++++ ON MET A ZERO LA SOMME QUE L ON VA CALCULER
  382. c      CALL ZERO(RINT,LRE,NLIGRP)
  383. *      CALL ZERO(RINT,NLIGRD,NLIGRP)
  384. c |-> impossible car ZERO considere qu il sagit de la dimension de SHPWRK
  385.       CALL ZEROX(RINT,NLIGRD,NLIGRP,LRE)
  386.  
  387. c
  388. c rem: il serait probablement interessant au niveau du tems cpu
  389. c  d'utiliser moins de pts de Gauss lorsque l element n est pas enrichi
  390. c  car cela n'apprte rien de plus
  391. c  On pourrait par ex. utiliser MINT2 qui contiendrait
  392. c  le segment d'integration de l'EF std (QUA4 par ex.)
  393.       if((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  394.          MINTE = MINT2
  395.          NBPGAU= NGAU2
  396.       else
  397.          MINTE = MINT1
  398.          NBPGAU= NGAU1
  399.       endif
  400. C
  401. C*********************************************************
  402. C
  403. C>>>>>>>>>>  BOUCLE SUR LES POINTS DE GAUSS
  404. C
  405.       DO 2500 KGAU=1,NBPGAU
  406. C
  407. C*********************************************************
  408. C     Initialisation à 0
  409. C*********************************************************
  410.  
  411. c ZERO ne serait-il pas facultatif?
  412. *      CALL ZERO(SHPWRK,6,NBNO)
  413.       CALL ZERO(SHPWRK,6,NBNI)
  414. C
  415.       i6zz = 3
  416.       IF (IDIM.EQ.3) i6zz = 4
  417. C
  418. c      do ienr7=1,NBENRMAX
  419.       do ienr7=1,NBENRJ
  420.        do inod7=1,NBNN1
  421. c        do i6=1,6
  422. CYT     do i6=1,3
  423.         do i6=1,i6zz
  424.            LV7WRK(ienr7,1,i6,inod7) = 0
  425.            LV7WRK(ienr7,2,i6,inod7) = 0
  426.         enddo
  427.        enddo
  428.       enddo
  429.  
  430. c*********************************************************
  431. c     Calcul des fonction de forme std dans repere local
  432. c*********************************************************
  433.  
  434. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  435. c     (et donc sur les Ni std)
  436.       DO 2510 I=1,NBNN1
  437.  
  438. C++++ Calcul des Ni std
  439. c     (rappel: 1:Ni  2:Ni,qsi  3:Ni,eta  avec i=1,4)
  440.       SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  441.       SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  442.       SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  443.       IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  444.  
  445.  2510 CONTINUE
  446. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  447.  
  448.  
  449.  
  450. c*********************************************************
  451. c     Passage des fonctions de forme std dans repere global
  452. c*********************************************************
  453. c
  454. C++++ CALCUL DES  Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  455.       CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  456. c
  457. C     PRISE EN COMPTE DU POIDS D'INTEGRATION
  458.       DJAC = ABS(DJAC) * POIGAU(KGAU)
  459.  
  460. c       CAS DES CONTRAINTES PLANES
  461. C       recuperation de l'epaisseur: DIM3 donnée facultative du materiau
  462.         IF (IFOUR.EQ.-2) THEN
  463.           MPTVAL=IVACAR
  464.           IF (IVACAR.NE.0) THEN
  465.             MELVAL=IVAL(1)
  466.             IF (MELVAL.NE.0) THEN
  467.               IGMN=MIN(KGAU,VELCHE(/1))
  468.               IBMN=MIN(J,VELCHE(/2))
  469.               DIM3=VELCHE(IGMN,IBMN)
  470.             ELSE
  471.               DIM3=1.D0
  472.             ENDIF
  473.           ENDIF
  474.           DJAC=DJAC*DIM3
  475.           MPTVAL=IVAMAT
  476.         ENDIF
  477.  
  478.  
  479. c*********************************************************
  480. c     Si on est pas du tout enrichi on peut sauter une grosse partie
  481. c*********************************************************
  482.       if(NBENRJ.eq.0) goto 2999
  483.  
  484. c*********************************************************
  485. c     Calcul des level set + leurs derivees dans repere global
  486. c*********************************************************
  487.  
  488. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  489.       do 2520 ienr=1,NBENRJ
  490.  
  491. c       MELVA1=MCHAM1.IELVAL(IENR)
  492. c       segact,MELVA1
  493.  
  494. ccccc  BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  495.        DO 2521 I=1,NBNN1
  496.  
  497. C++++ Le I eme noeud est-il ienr-enrichi?
  498.         mlree1= TLREEL(IENR,I)
  499.         if(mlree1.eq.0)  goto 2521
  500.  
  501.  
  502. c       Calcul du repere local de fissure(=PSI,PHI)
  503. c       (rappel: 1,1:psi  1,2:phi  2,1 psi,x  ...etc...)
  504.         do 2522 inode=1,NBNN1
  505. c         pour le H-enrichissement, on n a pas gardé PSI (inutile)
  506.           if(ienr.ne.1) then
  507. c          valeur de PSI au inode^ieme noeud
  508.            xpsi1 = mlree1.PROG(inode)
  509. c          qu on multiplie par la valeur de Ni^std au pt de G considéré
  510.            LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  511.      &      + (SHPWRK(1,inode)*xpsi1)
  512.            LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  513.      &      + (SHPWRK(2,inode)*xpsi1)
  514.            LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  515.      &      + (SHPWRK(3,inode)*xpsi1)
  516.       IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  517.      &      + (SHPWRK(4,inode)*xpsi1)
  518. c          valeur de PHI au inode^ieme noeud
  519.            xphi1 = mlree1.PROG(NBNN1+inode)
  520.           else
  521.            xphi1 = mlree1.PROG(inode)
  522.           endif
  523.           LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  524.      &      + (SHPWRK(1,inode)*xphi1)
  525.           LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  526.      &      + (SHPWRK(2,inode)*xphi1)
  527.           LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  528.      &      + (SHPWRK(3,inode)*xphi1)
  529.       IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  530.      &      + (SHPWRK(4,inode)*xphi1)
  531.  2522   continue
  532.  
  533.  2521  continue
  534. ccccc  fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  535.  
  536.  
  537.  2520 CONTINUE
  538. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  539.  
  540. c     on a construit
  541. C     LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  542.  
  543.  
  544. c*********************************************************
  545. c     Ajout des fonctions de forme d enrichissement
  546. c     + leurs derivees dans repere global
  547. c*********************************************************
  548. CYT   CALL SHAPX(XGAU,YGAU,MELE,LV7WRK,NBENRMAX,NBENRJ,
  549.       CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  550.  
  551.  
  552. c     retour a la partie commune aux EF enrichi et non enrichi
  553.  2999 continue
  554.  
  555. C*********************************************************
  556. C       CALCUL DE [N]
  557. C*********************************************************
  558. c ZERO ne serait-il pas facultatif?
  559. c        call ZERO(BGENE,LHOOK,NLIGRP)
  560.         KB=1
  561. C       boucle sur tous les Ni
  562.         DO 3001 II=1,NBNI
  563.  
  564.           BGENE(1,KB)   = SHPWRK(1,II)
  565.           BGENE(2,KB+1) = SHPWRK(1,II)
  566.  
  567.           IF(IDIM.EQ.3) BGENE(3,KB+2) = SHPWRK(1,II)
  568.  
  569.           KB = KB + IDIM
  570. CTY       KB = KB + 2
  571.  
  572.  3001   CONTINUE
  573. C
  574. c      if(J.eq.1.and.KGAU.eq.1) then
  575. c      write(ioimp,*) 'BGENE(1,..)=',(BGENE(1,iou),iou=1,2*NBNI)
  576. c      write(ioimp,*) 'BGENE(2,..)=',(BGENE(2,iou),iou=1,2*NBNI)
  577. c      endif
  578.  
  579. C*********************************************************
  580. C       CALCUL DE rho
  581. C*********************************************************
  582. c
  583. C     on remplit le segment MVELCH =>  valmat(1)=RHO
  584.           IF (IVAL(1).NE.0) THEN
  585.             MELVAL=IVAL(1)
  586.             IBMN=MIN(J   ,VELCHE(/2))
  587.             IGMN=MIN(KGAU,VELCHE(/1))
  588.             VALMAT(1)=VELCHE(IGMN,IBMN)
  589.           ELSE
  590.             VALMAT(1)=0.D0
  591.           ENDIF
  592.  
  593. C       PRISE EN COMPTE DE RHO
  594.         DJAC=DJAC*VALMAT(1)
  595. c        if(J.eq.1)   write(ioimp,*) 'DJAC=',DJAC
  596. c
  597. c
  598. C*********************************************************
  599. C       CALCUL DE  rho.Nt.N
  600. C*********************************************************
  601.  
  602.         CALL ZEROX(REL,NLIGRD,NLIGRP,LRE)
  603.  
  604.         CALL NTNST(BGENE,DJAC,LRE,LHOOK,REL)
  605. c     pour gagner du temps cpu il faudrait qqchose du type:
  606. c     CALL NTNSTX(BGENE,DJAC,DDHOOK,NLIGRP,2,LRE,REL)
  607. c        if(J.eq.1) then
  608. c         write(ioimp,*) 'REL',(REL(1,iou),iou=1,8)
  609. c         write(ioimp,*) 'REL',(REL(2,iou),iou=1,8)
  610. c         write(ioimp,*) 'REL',(REL(3,iou),iou=1,8)
  611. c         write(ioimp,*) 'REL',(REL(4,iou),iou=1,8)
  612. c        endif
  613.  
  614.  
  615. C*********************************************************
  616. C       CUMUL DANS RINT(II,JJ)
  617. C*********************************************************
  618.       DO 4201 II=1,NLIGRD
  619.          DO 4202 JJ=1,NLIGRP
  620.            RINT(II,JJ) =  RINT(II,JJ) + REL(II,JJ)
  621.  4202    CONTINUE
  622.  4201 CONTINUE
  623. c
  624. c
  625.  2500 CONTINUE
  626. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  627. C
  628.  
  629.  
  630. C*********************************************************
  631. c      write(ioimp,*) 'C       STOCKAGE DANS XMATRI de RE(II,JJ)'
  632. C       et de XMATRI dans IMATRI
  633. C*********************************************************
  634. C
  635. *      SEGINI,XMATRI
  636. *      IMATR1.IMATTT(NELRIG) = XMATRI
  637. c
  638.        DO 6001 II=1,NLIGRD
  639.  
  640.          DO 6002 JJ=1,NLIGRP
  641.  
  642.            xmatr1.RE(II,JJ,nelrig) =  RINT(II,JJ)
  643.  
  644. c     si NON enrichi, on met tout a 0 sauf la diago
  645. c          if(mlree1.eq.0)  RE(II,JJ) = 0.      inutile
  646.  
  647.  6002    CONTINUE
  648. c         if(J.eq.1) write(ioimp,*) 'RE_',II,'=',(RE(II,jou),jou=1,NLIGRP)
  649.  
  650.  6001 CONTINUE
  651. c
  652. c     quelques desactivations
  653. CYT   SEGDES,XMATRI
  654.       do IENR=1,NBENR1
  655.         do I=1,NBNN1
  656.          mlree1 = TLREEL(IENR,I)
  657.          if(mlree1.ne.0)  segdes,mlree1
  658.         enddo
  659.       enddo
  660. c
  661. c
  662. c
  663.  2000 CONTINUE
  664. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  665. c
  666. c
  667.  
  668. C*********************************************************
  669. C       SUPPRESSION ET DESACTIVATION DE SEGMENTS
  670. C*********************************************************
  671. C
  672.       SEGSUP,WRK1,WRK2,MVELCH
  673.       SEGSUP,MRACC
  674.  
  675.       segdes,MELEME
  676.       segdes,MINT1
  677.       if(MINT2.ne.0)  segdes,MINT2
  678. c     desactivation des maillages et segment imatri
  679.       do ienr=1,(1+NBENR1)
  680.          IPT1 = LOCIRI(1,ienr)
  681.          if(IPT1.ne.0)  segdes,IPT1
  682.          xMATR1 = LOCIRI(4,ienr)
  683.          if(xMATR1.ne.0)  segdes,xMATR1
  684.       enddo
  685.  
  686. c      SEGDES dans RIGI1 =   IMODEL
  687.  
  688.       RETURN
  689.       END
  690.  
  691.  
  692.  
  693.  

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