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    PV090527  26/04/30    21:15:53     12529          
  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.          rigrel=0
  373.          segini,xMATR1
  374.          LOCIRI(4,1+NBENRJ) = xMATR1
  375. c     si il existe, il faut l agrandir
  376.       else
  377.          NELRIG = xmatr1.RE(/3) + 1
  378.          rigrel=0
  379.          segadj,xMATR1
  380.       endif
  381.  
  382.  
  383. C++++ ON MET A ZERO LA SOMME QUE L ON VA CALCULER
  384. c      CALL ZERO(RINT,LRE,NLIGRP)
  385. *      CALL ZERO(RINT,NLIGRD,NLIGRP)
  386. c |-> impossible car ZERO considere qu il sagit de la dimension de SHPWRK
  387.       CALL ZEROX(RINT,NLIGRD,NLIGRP,LRE)
  388.  
  389. c
  390. c rem: il serait probablement interessant au niveau du tems cpu
  391. c  d'utiliser moins de pts de Gauss lorsque l element n est pas enrichi
  392. c  car cela n'apprte rien de plus
  393. c  On pourrait par ex. utiliser MINT2 qui contiendrait
  394. c  le segment d'integration de l'EF std (QUA4 par ex.)
  395.       if((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  396.          MINTE = MINT2
  397.          NBPGAU= NGAU2
  398.       else
  399.          MINTE = MINT1
  400.          NBPGAU= NGAU1
  401.       endif
  402. C
  403. C*********************************************************
  404. C
  405. C>>>>>>>>>>  BOUCLE SUR LES POINTS DE GAUSS
  406. C
  407.       DO 2500 KGAU=1,NBPGAU
  408. C
  409. C*********************************************************
  410. C     Initialisation à 0
  411. C*********************************************************
  412.  
  413. c ZERO ne serait-il pas facultatif?
  414. *      CALL ZERO(SHPWRK,6,NBNO)
  415.       CALL ZERO(SHPWRK,6,NBNI)
  416. C
  417.       i6zz = 3
  418.       IF (IDIM.EQ.3) i6zz = 4
  419. C
  420. c      do ienr7=1,NBENRMAX
  421.       do ienr7=1,NBENRJ
  422.        do inod7=1,NBNN1
  423. c        do i6=1,6
  424. CYT     do i6=1,3
  425.         do i6=1,i6zz
  426.            LV7WRK(ienr7,1,i6,inod7) = 0
  427.            LV7WRK(ienr7,2,i6,inod7) = 0
  428.         enddo
  429.        enddo
  430.       enddo
  431.  
  432. c*********************************************************
  433. c     Calcul des fonction de forme std dans repere local
  434. c*********************************************************
  435.  
  436. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  437. c     (et donc sur les Ni std)
  438.       DO 2510 I=1,NBNN1
  439.  
  440. C++++ Calcul des Ni std
  441. c     (rappel: 1:Ni  2:Ni,qsi  3:Ni,eta  avec i=1,4)
  442.       SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  443.       SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  444.       SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  445.       IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  446.  
  447.  2510 CONTINUE
  448. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  449.  
  450.  
  451.  
  452. c*********************************************************
  453. c     Passage des fonctions de forme std dans repere global
  454. c*********************************************************
  455. c
  456. C++++ CALCUL DES  Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  457.       CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  458. c
  459. C     PRISE EN COMPTE DU POIDS D'INTEGRATION
  460.       DJAC = ABS(DJAC) * POIGAU(KGAU)
  461.  
  462. c       CAS DES CONTRAINTES PLANES
  463. C       recuperation de l'epaisseur: DIM3 donnée facultative du materiau
  464.         IF (IFOUR.EQ.-2) THEN
  465.           MPTVAL=IVACAR
  466.           IF (IVACAR.NE.0) THEN
  467.             MELVAL=IVAL(1)
  468.             IF (MELVAL.NE.0) THEN
  469.               IGMN=MIN(KGAU,VELCHE(/1))
  470.               IBMN=MIN(J,VELCHE(/2))
  471.               DIM3=VELCHE(IGMN,IBMN)
  472.             ELSE
  473.               DIM3=1.D0
  474.             ENDIF
  475.           ENDIF
  476.           DJAC=DJAC*DIM3
  477.           MPTVAL=IVAMAT
  478.         ENDIF
  479.  
  480.  
  481. c*********************************************************
  482. c     Si on est pas du tout enrichi on peut sauter une grosse partie
  483. c*********************************************************
  484.       if(NBENRJ.eq.0) goto 2999
  485.  
  486. c*********************************************************
  487. c     Calcul des level set + leurs derivees dans repere global
  488. c*********************************************************
  489.  
  490. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  491.       do 2520 ienr=1,NBENRJ
  492.  
  493. c       MELVA1=MCHAM1.IELVAL(IENR)
  494. c       segact,MELVA1
  495.  
  496. ccccc  BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  497.        DO 2521 I=1,NBNN1
  498.  
  499. C++++ Le I eme noeud est-il ienr-enrichi?
  500.         mlree1= TLREEL(IENR,I)
  501.         if(mlree1.eq.0)  goto 2521
  502.  
  503.  
  504. c       Calcul du repere local de fissure(=PSI,PHI)
  505. c       (rappel: 1,1:psi  1,2:phi  2,1 psi,x  ...etc...)
  506.         do 2522 inode=1,NBNN1
  507. c         pour le H-enrichissement, on n a pas gardé PSI (inutile)
  508.           if(ienr.ne.1) then
  509. c          valeur de PSI au inode^ieme noeud
  510.            xpsi1 = mlree1.PROG(inode)
  511. c          qu on multiplie par la valeur de Ni^std au pt de G considéré
  512.            LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  513.      &      + (SHPWRK(1,inode)*xpsi1)
  514.            LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  515.      &      + (SHPWRK(2,inode)*xpsi1)
  516.            LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  517.      &      + (SHPWRK(3,inode)*xpsi1)
  518.       IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  519.      &      + (SHPWRK(4,inode)*xpsi1)
  520. c          valeur de PHI au inode^ieme noeud
  521.            xphi1 = mlree1.PROG(NBNN1+inode)
  522.           else
  523.            xphi1 = mlree1.PROG(inode)
  524.           endif
  525.           LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  526.      &      + (SHPWRK(1,inode)*xphi1)
  527.           LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  528.      &      + (SHPWRK(2,inode)*xphi1)
  529.           LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  530.      &      + (SHPWRK(3,inode)*xphi1)
  531.       IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  532.      &      + (SHPWRK(4,inode)*xphi1)
  533.  2522   continue
  534.  
  535.  2521  continue
  536. ccccc  fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  537.  
  538.  
  539.  2520 CONTINUE
  540. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  541.  
  542. c     on a construit
  543. C     LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  544.  
  545.  
  546. c*********************************************************
  547. c     Ajout des fonctions de forme d enrichissement
  548. c     + leurs derivees dans repere global
  549. c*********************************************************
  550. CYT   CALL SHAPX(XGAU,YGAU,MELE,LV7WRK,NBENRMAX,NBENRJ,
  551.       CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  552.  
  553.  
  554. c     retour a la partie commune aux EF enrichi et non enrichi
  555.  2999 continue
  556.  
  557. C*********************************************************
  558. C       CALCUL DE [N]
  559. C*********************************************************
  560. c ZERO ne serait-il pas facultatif?
  561. c        call ZERO(BGENE,LHOOK,NLIGRP)
  562.         KB=1
  563. C       boucle sur tous les Ni
  564.         DO 3001 II=1,NBNI
  565.  
  566.           BGENE(1,KB)   = SHPWRK(1,II)
  567.           BGENE(2,KB+1) = SHPWRK(1,II)
  568.  
  569.           IF(IDIM.EQ.3) BGENE(3,KB+2) = SHPWRK(1,II)
  570.  
  571.           KB = KB + IDIM
  572. CTY       KB = KB + 2
  573.  
  574.  3001   CONTINUE
  575. C
  576. c      if(J.eq.1.and.KGAU.eq.1) then
  577. c      write(ioimp,*) 'BGENE(1,..)=',(BGENE(1,iou),iou=1,2*NBNI)
  578. c      write(ioimp,*) 'BGENE(2,..)=',(BGENE(2,iou),iou=1,2*NBNI)
  579. c      endif
  580.  
  581. C*********************************************************
  582. C       CALCUL DE rho
  583. C*********************************************************
  584. c
  585. C     on remplit le segment MVELCH =>  valmat(1)=RHO
  586.           IF (IVAL(1).NE.0) THEN
  587.             MELVAL=IVAL(1)
  588.             IBMN=MIN(J   ,VELCHE(/2))
  589.             IGMN=MIN(KGAU,VELCHE(/1))
  590.             VALMAT(1)=VELCHE(IGMN,IBMN)
  591.           ELSE
  592.             VALMAT(1)=0.D0
  593.           ENDIF
  594.  
  595. C       PRISE EN COMPTE DE RHO
  596.         DJAC=DJAC*VALMAT(1)
  597. c        if(J.eq.1)   write(ioimp,*) 'DJAC=',DJAC
  598. c
  599. c
  600. C*********************************************************
  601. C       CALCUL DE  rho.Nt.N
  602. C*********************************************************
  603.  
  604.         CALL ZEROX(REL,NLIGRD,NLIGRP,LRE)
  605.  
  606.         CALL NTNST(BGENE,DJAC,LRE,LHOOK,REL)
  607. c     pour gagner du temps cpu il faudrait qqchose du type:
  608. c     CALL NTNSTX(BGENE,DJAC,DDHOOK,NLIGRP,2,LRE,REL)
  609. c        if(J.eq.1) then
  610. c         write(ioimp,*) 'REL',(REL(1,iou),iou=1,8)
  611. c         write(ioimp,*) 'REL',(REL(2,iou),iou=1,8)
  612. c         write(ioimp,*) 'REL',(REL(3,iou),iou=1,8)
  613. c         write(ioimp,*) 'REL',(REL(4,iou),iou=1,8)
  614. c        endif
  615.  
  616.  
  617. C*********************************************************
  618. C       CUMUL DANS RINT(II,JJ)
  619. C*********************************************************
  620.       DO 4201 II=1,NLIGRD
  621.          DO 4202 JJ=1,NLIGRP
  622.            RINT(II,JJ) =  RINT(II,JJ) + REL(II,JJ)
  623.  4202    CONTINUE
  624.  4201 CONTINUE
  625. c
  626. c
  627.  2500 CONTINUE
  628. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  629. C
  630.  
  631.  
  632. C*********************************************************
  633. c      write(ioimp,*) 'C       STOCKAGE DANS XMATRI de RE(II,JJ)'
  634. C       et de XMATRI dans IMATRI
  635. C*********************************************************
  636. C
  637. *      SEGINI,XMATRI
  638. *      IMATR1.IMATTT(NELRIG) = XMATRI
  639. c
  640.        DO 6001 II=1,NLIGRD
  641.  
  642.          DO 6002 JJ=1,NLIGRP
  643.  
  644.            xmatr1.RE(II,JJ,nelrig) =  RINT(II,JJ)
  645.  
  646. c     si NON enrichi, on met tout a 0 sauf la diago
  647. c          if(mlree1.eq.0)  RE(II,JJ) = 0.      inutile
  648.  
  649.  6002    CONTINUE
  650. c         if(J.eq.1) write(ioimp,*) 'RE_',II,'=',(RE(II,jou),jou=1,NLIGRP)
  651.  
  652.  6001 CONTINUE
  653. c
  654. c     quelques desactivations
  655. CYT   SEGDES,XMATRI
  656.       do IENR=1,NBENR1
  657.         do I=1,NBNN1
  658.          mlree1 = TLREEL(IENR,I)
  659.          if(mlree1.ne.0)  segdes,mlree1
  660.         enddo
  661.       enddo
  662. c
  663. c
  664. c
  665.  2000 CONTINUE
  666. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  667. c
  668. c
  669.  
  670. C*********************************************************
  671. C       SUPPRESSION ET DESACTIVATION DE SEGMENTS
  672. C*********************************************************
  673. C
  674.       SEGSUP,WRK1,WRK2,MVELCH
  675.       SEGSUP,MRACC
  676.  
  677.       segdes,MELEME
  678.       segdes,MINT1
  679.       if(MINT2.ne.0)  segdes,MINT2
  680. c     desactivation des maillages et segment imatri
  681.       do ienr=1,(1+NBENR1)
  682.          IPT1 = LOCIRI(1,ienr)
  683.          if(IPT1.ne.0)  segdes,IPT1
  684.          xMATR1 = LOCIRI(4,ienr)
  685.          if(xMATR1.ne.0)  segdes,xMATR1
  686.       enddo
  687.  
  688. c      SEGDES dans RIGI1 =   IMODEL
  689.  
  690.       RETURN
  691.       END
  692.  
  693.  
  694.  
  695.  
  696.  
  697.  

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