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

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