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

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