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

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