Télécharger la plaquette Cast3M


fr en

Télécharger bsigmx.eso

Retour à la liste

Numérotation des lignes : 

bsigmx

  1. C BSIGMX    SOURCE    MB234859  25/09/08    21:15:05     12358          
  2.       subroutine BSIGMX (IMODEL,IVACAR,IVASTR,ncar1,NFOR,
  3.      &         IVAFOR,ADPG,BDPG,CDPG,IIPDPG,IRETER)
  4.  
  5. c
  6. C     PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM
  7. c     POUR LE CALCUL DE la force interne (=B^T sigma)
  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.       PARAMETER  (NDDLMAX=30,NBNIMAX=10)
  18.       PARAMETER  (NBENRMAX=5)
  19.       DIMENSION  MLRE(NBENRMAX+1),NCOMPCUM(NBENRMAX+1)
  20.  
  21. -INC PPARAM
  22. -INC CCOPTIO
  23. -INC CCREEL
  24.  
  25. -INC SMCOORD
  26. -INC SMMODEL
  27. -INC SMCHAML
  28. -INC SMINTE
  29. -INC SMELEME
  30. -INC SMLREEL
  31.  
  32. -INC TMPTVAL
  33.  
  34.       POINTEUR    MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE
  35. C
  36.       SEGMENT WRK1
  37.        REAL*8 XE(3,NBBB)
  38.        REAL*8 XSTRS(LHOOK)
  39.        REAL*8 XFORC(LRE),XFINC(LRE)
  40.       ENDSEGMENT
  41. c
  42.       SEGMENT WRK2
  43.        REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
  44.        REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
  45.       ENDSEGMENT
  46.  
  47.       SEGMENT MRACC
  48.          INTEGER TLREEL(NBENRMA2,NBI)
  49.       ENDSEGMENT
  50. c
  51. *      write (*,*) '############################'
  52. *      write (*,*) '#####  DEBUT DE BSIGMX #####'
  53. *      write (*,*) '############################'
  54.       IVADIM = 0
  55. C
  56. C*********************************************************
  57. c  Introduction du point autour duquel se fait le mouvement
  58. c  de la section en defo plane generalisee
  59. C*********************************************************
  60.       ADPG=XZero
  61.       BDPG=XZero
  62.       CDPG=XZero
  63. C* On est en DEFO PLAN GENE :
  64.       IF (IIPDPG.GT.0) THEN
  65.         IREF=(IIPDPG-1)*(IDIM+1)
  66.         XDPGE=XCOOR(IREF+1)
  67.         YDPGE=XCOOR(IREF+2)
  68.       ELSE
  69.         XDPGE=XZero
  70.         YDPGE=XZero
  71.       ENDIF
  72. C
  73. C*********************************************************
  74. c
  75. C       RECUPERATION + ACTIVATIONS + VALEURS UTILES
  76. c
  77. C*********************************************************
  78. C     SEGACT MMODEL,IMODEL    deja actif
  79. c
  80. C++++ Activation au cas ou ++++++++++++++++++++++++++++++
  81.       SEGACT MCOORD
  82.  
  83. C++++ Recup + Activation de la geometrie ++++++++++++++++
  84.       MELEME= IMAMOD
  85. c      SEGACT MELEME    deja actif
  86. C     nbre de noeuds par element
  87.       NBNN1 = NUM(/1)
  88. C     nbre d elements
  89.       NBEL1 = NUM(/2)
  90.  
  91. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  92. c   + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1
  93. C     segment INFO deja actif dans RIGI1
  94. c   + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait
  95. c     mais du coup, on na pas de segment minte
  96. c     (car depend de si pt de g pour rigi, pour sigma....)
  97. c     c'est + simple de rappeler elquoi ici
  98.       MELE  = NEFMOD
  99.       NGAU1 = INFELE(6)
  100.       LRE   = INFELE(9)
  101.       LHOOK = INFELE(10)
  102.       MINT1 = INFMOD(6)
  103.       segact,MINT1
  104.       MINT2 = INFMOD(10)
  105.       if(MINT2.ne.0)   segact,MINT2
  106.       MFR   = INFELE(13)
  107.       IELE  = INFELE(14)
  108.       NDDL  = INFELE(15)
  109.       NSTRS = INFELE(16)
  110.  
  111. c REM: pour se passer du dimensionnement du nbre d'enrichissement dans
  112. c      elquoi et le realiser localement , on pourrait ecrire:
  113. c      LRE = NDDLMAX*NBNN1
  114. c      NDDL= NDDLMAX
  115.  
  116. C     sous decoupage et points de Gauss de l element geometrique de base
  117.       if(MELE.EQ.263.OR.MELE.EQ.264)  then
  118.          NGAU2 = MINT2.POIGAU(/1)
  119.       endif
  120. c      write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3))
  121. c      write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU)
  122.  
  123. c     nbre maxi de fonction de forme par noeud (fonction std comprise)
  124. c      NBNI = NDDL/IDIM    inutile!
  125.  
  126.  
  127. C++++ Recup des infos d enrichissement +++++++++++++++++++
  128. c     recup du MCHEX1 d'enrichissement
  129.       NOBMO1 = IVAMOD(/1)
  130.       if(NOBMO1.ne.0) then
  131.       do iobmo1=1,NOBMO1
  132.         if((TYMODE(iobmo1)).eq.'MCHAML') then
  133.           MCHEX1 = IVAMOD(iobmo1)
  134.           segact,MCHEX1
  135.           if((MCHEX1.TITCHE).eq.'ENRICHIS') then
  136.             MCHAM1 = MCHEX1.ICHAML(1)
  137.             segact,MCHAM1
  138.             goto 1000
  139.           endif
  140.         endif
  141.       enddo
  142.       write(*,*) 'Le modele est vide (absence d enrichissement)'
  143. *       return
  144.       else
  145.       write(*,*) 'Il n y a pas de MCHEML enrichissement dans le Modele'
  146. *       return
  147.       endif
  148.  
  149.  1000 continue
  150. c     niveau d enrichissement(s) du modele (ddls std u exclus)
  151. c     NBENR1= 0 si std, 1 si H seul,  2 si H+F1,  3 si H+F1+F2, etc...
  152.       if(NOBMO1.ne.0) then
  153.         NBENR1= MCHAM1.IELVAL(/1)
  154.       else
  155.         NBENR1= 0
  156.       endif
  157. c      write(*,*) 'niveau d enrichissement(s) du modele',NBENR1
  158. C
  159. C*********************************************************
  160. C       INITIALISATIONS...
  161. C*********************************************************
  162.       IRETER = 0
  163.       DIM3 = 1.D0
  164.       NCOSOR = 0
  165. c
  166. c     preparation des tables avec:
  167.  
  168.       if(NBENR1.ne.0) then
  169.       do ienr=1,NBENR1
  170. c        -le nombre d'inconnues de chaque sous-zone
  171. c         determinee depuis le nombre de fonction de forme
  172. c ienr=  1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  173.          nbniJ = 2 + ((ienr-1)*4)
  174.          MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  175.          NCOSOR = NCOSOR + (IDIM*nbniJ)
  176.       enddo
  177.       endif
  178. C     Tables + longues car 1er indice correspond au fontion de forme std
  179.       MLRE(1)    = IDIM*NBNN1*1
  180.       NCOSOR = NCOSOR + (IDIM)
  181.  
  182.       if(NBENR1.lt.(NBENRMAX+1)) then
  183.         do ienr=(NBENR1+1),(NBENRMAX)
  184.            MLRE(1+ienr) = 0
  185.         enddo
  186.       endif
  187. c
  188. c      ...DU SEGMENT WRK1
  189.       NBENRMA2 = NBENRMAX
  190.       NBBB = NBNN1
  191.       SEGINI,WRK1
  192.  
  193. c      ...DU SEGMENT WRK2
  194. c      NBNO = NBNI
  195.       NBNO = LRE/IDIM
  196.       SEGINI,WRK2
  197. C
  198. c      ...DU SEGMENT MRACC
  199.       NBENRMA2 = NBENRMAX
  200.       NBI = NBNN1
  201.       segini,MRACC
  202. C
  203. c      ...DU SEGMENT IVAFOR de type MPTVAL
  204.       MPTVAL = IVAFOR
  205. *     on change  NSR = IPOS(/1) et NCOSOR = IVAL(/1) (calculé + haut)
  206.       NSR = 1 + NBENR1
  207. *      write(*,*) 'segadj IVAFOR aux dim',NSR,NCOSOR
  208.       segadj,MPTVAL
  209. c     on remplit NSOF avec le nombre de composante par sous-zone
  210. c     et NCOMPCUM avec le nombre de composante cumulée
  211.       NCUMUL1 = 0
  212.       do ienr1=0,NBENR1
  213.          NCOMPCUM(1+ienr1) = NCUMUL1
  214.          NSOF1 = (MLRE(1+ienr1)) / NBNN1
  215.          NSOF(1+ienr1) = NSOF1
  216.          NCUMUL1 = NCUMUL1 + NSOF1
  217.       enddo
  218.       if(NBENR1.lt.(NBENRMAX+1)) then
  219.       do ienr=(NBENR1+1),(NBENRMAX)
  220.          NCOMPCUM(1+ienr1) = NCUMUL1
  221.       enddo
  222.       endif
  223. c     on recupere les dimensions (MAXI)
  224.       MELVAL = IVAL(1)
  225.       N1PTEL = VELCHE(/1)
  226.       N1EL   = VELCHE(/2)
  227.       N2PTEL = IELCHE(/1)
  228.       N2EL   = IELCHE(/2)
  229. c      write(*,*) 'MELVAL de taille maxi=',N1PTEL,N1EL,N2PTEL,N2EL
  230.  
  231. C*********************************************************
  232. C
  233. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  BOUCLE SUR LES ELEMENTS
  234. C
  235.       DO 2000 J=1,NBEL1
  236. c      write(*,*) '========= EF',J,' /',NBEL1,'========='
  237. C
  238. C*********************************************************
  239. C     POUR CHAQUE ELEMENT, ...
  240. C*********************************************************
  241. C
  242. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  243.       CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  244. C
  245. C++++ NBENRJ = niveau d enrichissement de l element ++++
  246. C     =0 si EF std   =1 si U+H   =2 si U+H+F1   =3 si U+H+F1+F2
  247.       NBENRJ=0
  248.       if(NBENR1.ne.0) then
  249.       do IENR=1,NBENR1
  250.         MELVA1 = MCHAM1.IELVAL(IENR)
  251.         segact,MELVA1
  252.         do I=1,NBNN1
  253.           mlree1 = MELVA1.IELCHE(I,J)
  254. C         on en profite pour remplir MRACC table de raccourcis pour cet element
  255.           TLREEL(IENR,I) = mlree1
  256.           if(mlree1.ne.0)  then
  257.             NBENRJ = max(NBENRJ,IENR)
  258. C           et on active la listreel
  259.             segact,mlree1
  260.           endif
  261.         enddo
  262.       enddo
  263.       endif
  264.       if(NBENRMA2.gt.NBENR1) then
  265.         do IENR2=(NBENR1+1),NBENRMA2
  266.         do I=1,NBNN1
  267.            TLREEL(IENR2,I) = 0
  268.         enddo
  269.         enddo
  270.       endif
  271. c
  272. c++++ quelques indices pour dimensionner au plus juste
  273. c     nbre total de ddl de l'élément considéré
  274.       NLIGRD = MLRE(1+NBENRJ)
  275.       NLIGRP = MLRE(1+NBENRJ)
  276. c     nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  277.       NBNI = (MLRE(1+NBENRJ)) / IDIM
  278. c      write(*,*) 'EF',J,' NBENRJ=',NBENRJ,
  279. c     &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'
  280. c
  281. c++++ Selon le niveau d enrichissement,
  282. c
  283. C   + On copie cet element a la geo correspondante
  284.       IPT1 = IPOS(1+NBENRJ)
  285. c     si elle n existe pas, il faut la creer
  286.       if(IPT1.eq.0) then
  287.          NBNN  =NBNN1
  288.          NBELEM=1
  289.          NBSOUS=0
  290.          NBREF=0
  291.          segini,IPT1
  292.          IPT1.ITYPEL = ITYPEL
  293.          IPOS(1+NBENRJ)=IPT1
  294. c     si elle existe, il faut l agrandir
  295.       else
  296.          NBNN  =NBNN1
  297.          NBELEM = (IPT1.NUM(/2)) + 1
  298.          NBSOUS=0
  299.          NBREF=0
  300.          segadj,IPT1
  301.       endif
  302.       J1 = NBELEM
  303. c     copie en cours
  304.       do I=1,NBNN1
  305.          IPT1.NUM(I,NBELEM) = NUM(I,J)
  306.       enddo
  307.  
  308. C   + On en profite pour créer les MELVA2 si ils n existent pas
  309.       MPTVAL = IVAFOR
  310.       NCUMUL1 = NCOMPCUM(1+NBENRJ)
  311.       NCOMP2  = NSOF(1+NBENRJ)
  312.       do icomp2=1,NCOMP2
  313.          MELVA2 = IVAL(NCUMUL1+icomp2)
  314.          if(MELVA2.eq.0) then
  315. c          N1PTEL N1EL,N2PTEL,N2EL = dim MAX : récupérés + haut
  316.            segini,MELVA2
  317.            IVAL(NCUMUL1+icomp2) = MELVA2
  318.          endif
  319.       enddo
  320.  
  321. C++++ MISE A 0 DES FORCES  ++++
  322. C
  323.       CALL ZERO(XFINC,1,NLIGRP)
  324. c
  325. c rem: il serait probablement interessant au niveau du tems cpu
  326. c  d'utiliser moins de pts de Gauss lorsque l element est élastique
  327. c  On pourrait par ex. utiliser MINT2 = infell(12) qui contient
  328. c  le segment d'integration de l'EF std (QUA4 par ex.)
  329. *      if((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  330. *         MINTE = MINT2
  331. *         NBPGAU= NGAU2
  332. *      else
  333.          MINTE = MINT1
  334.          NBPGAU= NGAU1
  335. *      endif
  336. c
  337.         NBPTEL=NBPGAU
  338. C
  339. C*********************************************************
  340. C
  341.       IDIM1 = IDIM + 1
  342. C>>>>>>>>>>  BOUCLE SUR LES POINTS DE GAUSS
  343. C
  344.       DO 2500 KGAU=1,NBPGAU
  345. C
  346. C*********************************************************
  347. C     Initialisation à 0
  348. C*********************************************************
  349.  
  350. c ZERO ne serait-il pas facultatif?
  351.       CALL ZERO(SHPWRK,6,NBNO)
  352. C
  353.       i6zz = 3
  354.       IF (IDIM.EQ.3) i6zz = 4
  355. C
  356. c      do ienr7=1,NBENRMAX
  357.       do ienr7=1,NBENRJ
  358.        do inod7=1,NBNN1
  359. c        do i6=1,6
  360.         do i6=1,i6zz
  361.            LV7WRK(ienr7,1,i6,inod7) = 0.D0
  362.            LV7WRK(ienr7,2,i6,inod7) = 0.D0
  363.         enddo
  364.        enddo
  365.       enddo
  366.  
  367.  
  368. c*********************************************************
  369. c     Calcul des fonction de forme std dans repere local
  370. c*********************************************************
  371.  
  372. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  373. c     (et donc sur les Ni std)
  374.       DO 2510 I=1,NBNN1
  375.       DO 2511 JJ=1,IDIM1
  376.  
  377. C++++ Calcul des Ni std
  378. c     (rappel: 1:Ni  2:Ni,qsi  3:Ni,eta  avec i=1,4)
  379.       SHPWRK(JJ,I) = SHPTOT(JJ,I,KGAU)
  380.  2511 CONTINUE
  381.  2510 CONTINUE
  382. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  383.  
  384.  
  385.  
  386. c*********************************************************
  387. c     Passage des fonctions de forme std dans repere global
  388. c*********************************************************
  389.  
  390. C++++ CALCUL DES  Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  391.       CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  392. c
  393. C     PRISE EN COMPTE DU POIDS D'INTEGRATION
  394.       DJAC = ABS(DJAC) * POIGAU(KGAU)
  395. c
  396. C     PRISE EN COMPTE DE L EPAISSEUR  (cas Contrainte Planes)
  397.       IF (IFOUR.EQ.-2)THEN
  398.         MPTVAL=IVACAR
  399.         IF (IVACAR.NE.0) THEN
  400.           MELVAL=IVAL(1)
  401.           IF (MELVAL.NE.0) THEN
  402.             IGMN=MIN(KGAU,VELCHE(/1))
  403.             IBMN=MIN(J,VELCHE(/2))
  404.             DIM3=VELCHE(IGMN,IBMN)
  405.             DJAC=DJAC*DIM3
  406. c         ELSE
  407. c           DIM3=1.D0
  408.           ENDIF
  409.         ENDIF
  410.       ENDIF
  411. c
  412. c      if(J.eq.1.and.KGAU.eq.1)
  413. c     &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  414.  
  415. c*********************************************************
  416. c     Si on est pas du tout enrichi on peut sauter une grosse partie
  417. c*********************************************************
  418.       if(NBENRJ.eq.0) goto 2999
  419.  
  420. c*********************************************************
  421. c     Calcul des level set + leurs derivees dans repere global
  422. c*********************************************************
  423.  
  424. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  425.       do 2520 ienr=1,NBENRJ
  426.  
  427. c       MELVA1=MCHAM1.IELVAL(IENR)
  428. c       segact,MELVA1
  429.  
  430. ccccc  BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  431.        DO 2521 I=1,NBNN1
  432.  
  433. C++++ Le I eme noeud est-il ienr-enrichi?
  434.         mlree1= TLREEL(IENR,I)
  435.         if(mlree1.eq.0)  goto 2521
  436.  
  437.  
  438. c       Calcul du repere local de fissure(=PSI,PHI)
  439. c       (rappel: 1,1:psi  1,2:phi  2,1 psi,x  ...etc...)
  440.         do 2522 inode=1,NBNN1
  441. c         pour le H-enrichissement, on n a pas gardé PSI (inutile)
  442.           if(ienr.ne.1) then
  443. c          valeur de PSI au inode^ieme noeud
  444.            xpsi1 = mlree1.PROG(inode)
  445. c          qu on multiplie par la valeur de Ni^std au pt de G considéré
  446.         do 2523 jj=1,IDIM1
  447.            LV7WRK(ienr,1,JJ,I)= LV7WRK(ienr,1,JJ,I)
  448.      &      + (SHPWRK(JJ,inode)*xpsi1)
  449.  2523   continue
  450.  
  451. C     IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  452. C    &      + (SHPWRK(4,inode)*xpsi1)
  453. c          valeur de PHI au inode^ieme noeud
  454.            xphi1 = mlree1.PROG(NBNN1+inode)
  455.           else
  456.            xphi1 = mlree1.PROG(inode)
  457.           endif
  458.         do 2524 jj=1,IDIM1
  459.           LV7WRK(ienr,2,jj,I)= LV7WRK(ienr,2,jj,I)
  460.      &      + (SHPWRK(jj,inode)*xphi1)
  461. C     IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  462. C    &      + (SHPWRK(4,inode)*xphi1)
  463.  2524   continue
  464.  2522   continue
  465.  
  466.  2521  continue
  467. ccccc  fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  468.  
  469.  
  470.  2520 CONTINUE
  471. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  472.  
  473. c     on a construit
  474. C     LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  475.  
  476.  
  477. c*********************************************************
  478. c     Ajout des fonctions de forme d enrichissement
  479. c     + leurs derivees dans repere global
  480. c*********************************************************
  481.       CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  482.  
  483.  
  484. c     retour a la partie commune aux EF enrichi et non enrichi
  485.  2999 continue
  486.  
  487. C*********************************************************
  488. C       CALCUL DE B
  489. C*********************************************************
  490. c ZERO ne serait-il pas facultatif?
  491. c        call ZERO(BGENE,LHOOK,NLIGRP)
  492.         KB=1
  493. C       boucle sur tous les Ni
  494.         DO 3001 II=1,NBNI
  495.  
  496.           BGENE(1,KB)   = SHPWRK(2,II)
  497.           BGENE(2,KB+1) = SHPWRK(3,II)
  498.           BGENE(4,KB)   = SHPWRK(3,II)
  499.           BGENE(4,KB+1) = SHPWRK(2,II)
  500.  
  501.         IF(IDIM.EQ.3) THEN
  502.         BGENE(3,KB+2)=SHPWRK(4,II)
  503.         BGENE(5,KB)=SHPWRK(4,II)
  504.         BGENE(5,KB+2)=SHPWRK(2,II)
  505.         BGENE(6,KB+1)=SHPWRK(4,II)
  506.         BGENE(6,KB+2)=SHPWRK(3,II)
  507.         ENDIF
  508.  
  509.           KB = KB + IDIM
  510.  
  511.  3001   CONTINUE
  512. C
  513. c
  514. C*********************************************************
  515. C       CALCUL DE sigma
  516. C*********************************************************
  517.         MPTVAL = IVASTR
  518.         DO 3200 ICOMP=1,LHOOK
  519.           MELVAL = IVAL(ICOMP)
  520.           IGMN = MIN(KGAU,VELCHE(/1))
  521.           IBMN = MIN(J  ,VELCHE(/2))
  522.           XSTRS(ICOMP) = VELCHE(IGMN,IBMN)
  523. *          XSTRS(ICOMP) = VELCHE(KGAU,J)
  524.  3200   CONTINUE
  525.  
  526.  
  527. C*********************************************************
  528. C       CALCUL DE Fint = B^T * sigma
  529. C*********************************************************
  530. c
  531.         CALL ZERO(XFORC,1,NLIGRP)
  532.         CALL BSIG(BGENE,XSTRS,LHOOK,NLIGRP,DJAC,XFORC)
  533. c
  534. * rem: matrice d'efficacite  presente dans bsigm1.eso  non copiée ici
  535. c
  536. c       cumul
  537.         do ii = 1,NLIGRP
  538.            XFINC(ii) = XFINC(ii) + XFORC(ii)
  539.         enddo
  540. c
  541. c
  542.  2500 CONTINUE
  543. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  544. C
  545. C
  546. C   EXTRACTION DES FORCES AU NOEUD SUPPORT DE LA DEF PLAN GENE
  547. C   ON CALCULE LES RESULTANTES DES FORCES SUR CHAQUE ELEMENT
  548. C
  549.       NFOFO=NFOR
  550.       IF (IIPDPG.GT.0) THEN
  551.       IF (IFOUR.EQ.-3) THEN
  552.         NFOFO=NFOR-3
  553.         ADPG=ADPG+XFINC(NBNN1*NFOFO+1)
  554.         BDPG=BDPG+XFINC(NBNN1*NFOFO+2)
  555.         CDPG=CDPG+XFINC(NBNN1*NFOFO+3)
  556.       ELSE IF (IFOUR.EQ. 7.OR.IFOUR.EQ. 8.OR.IFOUR.EQ.9.OR.
  557.      .         IFOUR.EQ.10.OR.IFOUR.EQ.14) THEN
  558.         NFOFO=NFOR-1
  559.         ADPG=ADPG+XFINC(NBNN1*NFOFO+1)
  560.       ELSE IF (IFOUR.EQ.11) THEN
  561.         NFOFO=NFOR-2
  562.         ADPG=ADPG+XFINC(NBNN1*NFOFO+1)
  563.         BDPG=BDPG+XFINC(NBNN1*NFOFO+2)
  564.       ENDIF
  565.       ENDIF
  566. C
  567. C     ON RANGE XFORC DANS le bon MELVAL
  568. c car plusieurs fois la meme composante dans ivafor pour créer des zones
  569. C
  570.       MPTVAL = IVAFOR
  571.       INITOT = 0
  572. C-->> BOUCLE SUR LES niveaux d'ENRICHISSEMENTS (0:U 1:A 2:BCDE1 3:BCDE2)
  573.       DO IENR=0,NBENRJ
  574. *nbre de fonction(s) de ce niveau d'enrichissement (=1 si std ou H, =4 pour F1,2,...)
  575.         if(IENR .le. 1) then
  576.          NBNIENR = 1
  577.         else
  578.          NBNIENR = 4
  579.         endif
  580. C---->> BOUCLE SUR LES fonctions de forme de ce type d'enrichissement
  581.         DO INI=1,NBNIENR
  582.           INITOT = INITOT + 1
  583. C------>> BOUCLE SUR LA DIMENSION
  584.           DO KDIM=1,IDIM
  585.              ICOMP = (INITOT-1)*IDIM + KDIM
  586. c            MELVAL = IVAL(ICOMP)
  587.             MELVAL = IVAL(NCUMUL1+ICOMP)
  588. C---------->> BOUCLE SUR LES NOEUDS
  589.             DO I=1,NBNN1
  590.               IQ = ((INITOT-1)*NBNN1*IDIM) + ((I-1)*IDIM) + KDIM
  591. *              IBMN=MIN(IB  ,VELCHE(/2))
  592. *              VELCHE(KGAU,IBMN) = XFINC(IQ)
  593.               VELCHE(I,J1) = XFINC(IQ)
  594.             ENDDO
  595.           ENDDO
  596.         ENDDO
  597.       ENDDO
  598. c
  599.  2000 CONTINUE
  600. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  601.  
  602. C*********************************************************
  603. C       SUPPRESSION ET DESACTIVATION DE SEGMENTS
  604. C*********************************************************
  605.       SEGSUP,WRK1,WRK2
  606.       SEGSUP,MRACC
  607.  
  608. c     ajustement avant retour des MELVAL contenus dans IVAFOR
  609. c     en focntion de la géométrie créée
  610.       MPTVAL = IVAFOR
  611.       do 3000 ienr1=0,NBENR1
  612.          IPT1   = IPOS(1+ienr1)
  613.          if(IPT1.eq.0) goto 3000
  614.          N1PTEL = IPT1.NUM(/1)
  615.          N1EL   = IPT1.NUM(/2)
  616.          N2PTEL = 1
  617.          N2EL   = 1
  618.          NCUMUL1= NCOMPCUM(1+ienr1)
  619.          NCOMP2 = NSOF(1+ienr1)
  620. c         write(6,*) '=======ienr1=',ienr1,'/',NBENR1,'======'
  621. c         write(6,*) 'de dim geo=',N1PTEL,N1EL
  622. c         write(6,*) 'avec',NCOMP2,'composantes apres lindice',NCUMUL1
  623.          do icomp2=1,NCOMP2
  624.             MELVA2 = IVAL(NCUMUL1+icomp2)
  625. c         write(6,*) MELVA2,'=',(MELVA2.VELCHE(iou,N1EL),iou=1,4)
  626.             if(MELVA2.ne.0)  segadj,MELVA2
  627.          enddo
  628.  3000 continue
  629.  
  630.       END
  631.  
  632.  
  633.  
  634.  

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