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

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