Télécharger la plaquette Cast3M


fr en

Télécharger sigmax.eso

Retour à la liste

Numérotation des lignes : 

sigmax

  1. C SIGMAX    SOURCE    MB234859  25/09/08    21:16:11     12358          
  2.  
  3.       subroutine SIGMAX (MATE,IMAT,NBGMAT,NELMAT,NMATT,CMATE,
  4.      &                  IVAMAT,IMODEL,IREPS2,IVADEP,
  5.      &                  IVASTR,UZDPG,RYDPG,RXDPG,IIPDPG,IRETER)
  6.  
  7. c
  8. C     PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM
  9. c     POUR LE CALCUL DE la contrainte (élastique)
  10. C
  11. C
  12. C*********************************************************
  13. C       PARTIE DECLARATIVE
  14. C*********************************************************
  15.  
  16.       IMPLICIT INTEGER(I-N)
  17.       IMPLICIT REAL*8(A-H,O-Z)
  18.  
  19. -INC PPARAM
  20. -INC CCOPTIO
  21.  
  22. -INC SMCOORD
  23. -INC SMMODEL
  24. -INC SMCHAML
  25. -INC SMINTE
  26. -INC SMELEME
  27. -INC SMLREEL
  28. C
  29.       POINTEUR    MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE
  30.  
  31. -INC TMPTVAL
  32.  
  33. c
  34.       SEGMENT WRK1
  35.        REAL*8 XE(3,NBBB)
  36.        REAL*8 DDHOOK(LHOOK,LHOOK)
  37.        REAL*8 XDDL(LRE),XSTRS(LHOOK)
  38.       ENDSEGMENT
  39. c
  40.       SEGMENT WRK2
  41.        REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
  42. c       REAL*8 LV7WRK(NBENRMA2,2,6,NBNO)
  43.        REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
  44.       ENDSEGMENT
  45. C
  46.       SEGMENT,MVELCH
  47.        REAL*8 VALMAT(NV1)
  48.       ENDSEGMENT
  49. C
  50.       SEGMENT MRACC
  51.          INTEGER TLREEL(NBENRMA2,NBI)
  52.       ENDSEGMENT
  53. C
  54.       SEGMENT MTRACE
  55.         REAL*8 TRACE(3,NBPTEL)
  56.       ENDSEGMENT
  57. C
  58.       PARAMETER  (NDDLMAX=30,NBNIMAX=10)
  59. CTY   PARAMETER  (NDDLMAX=20,NBNIMAX=10)
  60.  
  61.       PARAMETER  (NBENRMAX=5)
  62.       DIMENSION  MLRE(NBENRMAX+1)
  63. C
  64.       CHARACTER*8 CMATE
  65.       DIMENSION UDPGE(3)
  66.       LOGICAL BDPGE
  67. C
  68. c      write (*,*) '#############################'
  69. c      write (*,*) '#####  DEBUT DE SIGMAX  #####'
  70. c      write (*,*) '#############################'
  71. C
  72. C*********************************************************
  73. c  Introduction du point autour duquel se fait le mouvement
  74. c  de la section en defo plane generalisee
  75. C*********************************************************
  76. C  En 1D : pas de rotation
  77.       BDPGE=.FALSE.
  78.       NDPGE=0
  79.       XDPGE=0.D0
  80.       YDPGE=0.D0
  81.       IF (IFOUR.EQ.-3) THEN
  82.         BDPGE=.TRUE.
  83.         NDPGE=3
  84.         UDPGE(1)=UZDPG
  85.         UDPGE(2)=RYDPG
  86.         UDPGE(3)=RXDPG
  87.         SEGACT,MCOORD
  88.         IREF=(IIPDPG-1)*(IDIM+1)
  89.         XDPGE=XCOOR(IREF+1)
  90.         YDPGE=XCOOR(IREF+2)
  91.       ELSE IF (IDIM.EQ.1) THEN
  92.         IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR.
  93.      .      IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN
  94.           BDPGE=.TRUE.
  95.           NDPGE=1
  96.           UDPGE(1)=UZDPG
  97.         ELSE IF (IFOUR.EQ.11) THEN
  98.           BDPGE=.TRUE.
  99.           NDPGE=2
  100.           UDPGE(1)=UZDPG
  101.           UDPGE(2)=RXDPG
  102.         ENDIF
  103.       ENDIF
  104. C
  105.       NV1=NMATT
  106.       SEGINI,MVELCH
  107. C
  108. C
  109. C*********************************************************
  110. c
  111. C       RECUPERATION + ACTIVATIONS + VALEURS UTILES
  112. c
  113. C*********************************************************
  114. C     SEGACT MMODEL,IMODEL    deja activé dans epsi1
  115. c
  116. C++++ Activation au cas ou ++++++++++++++++++++++++++++++
  117.       SEGACT MCOORD
  118.  
  119. C++++ Recup + Activation de la geometrie ++++++++++++++++
  120.       MELEME= IMAMOD
  121. c      SEGACT MELEME    deja activé dans epsi1
  122.  
  123.  
  124. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  125. c   + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1
  126. C     segment INFO deja actif dans RIGI1
  127. c   + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait
  128. c     mais du coup, on na pas de segment minte
  129. c     (car depend de si pt de g pour rigi, pour sigma....)
  130. c     c'est + simple de rappeler elquoi ici
  131.       MELE = NEFMOD
  132.       NGAU1 = INFELE(6)
  133.       LRE   = INFELE(9)
  134.       LHOOK = INFELE(10)
  135.       MINT1 = INFMOD(6)
  136.       segact,MINT1
  137.       MINT2 = INFMOD(10)
  138.       if(MINT2.ne.0)   segact,MINT2
  139.       MFR   = INFELE(13)
  140.       IELE  = INFELE(14)
  141.       NDDL  = INFELE(15)
  142.       NSTRS = INFELE(16)
  143. c      write(6,*)'-> EPSIX INFELE',(INFELE(iou),iou=1,16)
  144.  
  145. c   + AUTRES INFOS
  146. C     nbre de noeuds par element
  147.       NBNN1 = NUM(/1)
  148. C     nbre d elements
  149.       NBEL1 = NUM(/2)
  150.  
  151. c REM: pour se passer du dimensionnement du nbre d'enrichissement dans
  152. c      elquoi et le realiser localement , on pourrait ecrire:
  153. c      LRE = NDDLMAX*NBNN1
  154. c      NDDL= NDDLMAX
  155.  
  156. C     sous decoupage et points de Gauss de l element geometrique de base
  157. CTY   if(MELE.eq.263)  then
  158.       if(MELE.EQ.263.OR.MELE.EQ.264)  then
  159.          NGAU2 = MINT2.POIGAU(/1)
  160.       endif
  161. c      write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3))
  162. c      write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU)
  163.  
  164.  
  165. c     nbre maxi de fonction de forme par noeud (fonction std comprise)
  166. c      NBNI = NDDL/IDIM    inutile!
  167.  
  168.  
  169. C++++ Recup des infos d enrichissement +++++++++++++++++++
  170. c     recup du MCHEX1 d'enrichissement
  171.       NOBMO1 = IVAMOD(/1)
  172.       if(NOBMO1.ne.0) then
  173.        do iobmo1=1,NOBMO1
  174.         if((TYMODE(iobmo1)).eq.'MCHAML') then
  175.          MCHEX1 = IVAMOD(iobmo1)
  176.          segact,MCHEX1
  177.           if((MCHEX1.TITCHE).eq.'ENRICHIS') then
  178.             MCHAM1 = MCHEX1.ICHAML(1)
  179.             segact,MCHAM1
  180.             goto 1000
  181.           endif
  182.         endif
  183.        enddo
  184.        write(*,*) 'Le modele est vide (absence d enrichissement)'
  185. *       return
  186.       else
  187.       write(*,*) 'Il n y a pas de MCHEML enrichissement dans le Modele'
  188. *       return
  189.       endif
  190.  
  191.  1000 continue
  192. c     niveau d enrichissement(s) du modele (ddls std u exclus)
  193. c     NBENR1= 0 si std, 1 si H seul,  2 si H+F1,  3 si H+F1+F2, etc...
  194.       if(NOBMO1.ne.0) then
  195.         NBENR1= MCHAM1.IELVAL(/1)
  196.       else
  197.         NBENR1 = 0
  198.       endif
  199. c      write(*,*) 'niveau d enrichissement(s) du modele',NBENR1
  200. c
  201. C*********************************************************
  202. C       INITIALISATIONS...
  203. C*********************************************************
  204.       IRETER = 0
  205. c
  206. c     preparation des tables avec:
  207.  
  208.       if(NBENR1.ne.0) then
  209.       do ienr=1,NBENR1
  210. c        -le nombre d'inconnues de chaque sous-zone
  211. c         determinee depuis le nombre de fonction de forme
  212. c ienr=  1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  213.          nbniJ = 2 + ((ienr-1)*4)
  214.          MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  215.       enddo
  216.       endif
  217. C     Tables + longues car 1er indice correspond au fontion de forme std
  218.       MLRE(1)    = IDIM*NBNN1*1
  219.  
  220.       if(NBENR1.lt.(NBENRMAX+1)) then
  221.         do ienr=(NBENR1+1),(NBENRMAX)
  222.            MLRE(1+ienr) = 0
  223.         enddo
  224.       endif
  225. c
  226. c      ...DU SEGMENT WRK1
  227.       NBENRMA2 = NBENRMAX
  228.       NBBB = NBNN1
  229.       SEGINI,WRK1
  230.  
  231. c      ...DU SEGMENT WRK2
  232. c      NBNO = NBNI
  233.       NBNO = LRE/IDIM
  234.       SEGINI,WRK2
  235. C
  236. c      ...DU SEGMENT MRACC
  237.       NBENRMA2 = NBENRMAX
  238.       NBI = NBNN1
  239.       segini,MRACC
  240. C
  241. C     du nombre d erreur sur le noms de composantes
  242.       NBERR1=0
  243.  
  244. C*********************************************************
  245. C
  246. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  BOUCLE SUR LES ELEMENTS
  247. C
  248.       DO 2000 J=1,NBEL1
  249. c      write(*,*) '==================================='
  250. c      write(*,*) '========= EF',J,' /NBEL1 ========='
  251. C
  252. C
  253. C*********************************************************
  254. C     POUR CHAQUE ELEMENT, ...
  255. C*********************************************************
  256. C
  257. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  258.       CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  259. C
  260. c
  261. C++++ NBENRJ = niveau d enrichissement de l element ++++
  262. C     =0 si EF std   =1 si U+H   =2 si U+H+F1   =3 si U+H+F1+F2
  263.       NBENRJ=0
  264.       if(NBENR1.ne.0) then
  265.       do IENR=1,NBENR1
  266.         MELVA1 = MCHAM1.IELVAL(IENR)
  267.         segact,MELVA1
  268.         do I=1,NBNN1
  269.          mlree1 = MELVA1.IELCHE(I,J)
  270. C        on en profite pour remplir MRACC table de raccourcis pour cet element
  271.          TLREEL(IENR,I) = mlree1
  272.          if(mlree1.ne.0)  then
  273.             NBENRJ = max(NBENRJ,IENR)
  274. C           et on active la listreel
  275.             segact,mlree1
  276.          endif
  277.         enddo
  278.       enddo
  279.       endif
  280.       if(NBENRMA2.gt.NBENR1) then
  281.         do IENR2=(NBENR1+1),NBENRMA2
  282.         do I=1,NBNN1
  283.            TLREEL(IENR2,I) = 0
  284.         enddo
  285.         enddo
  286.       endif
  287. C
  288. c
  289. c++++ quelques indices pour dimensionner au plus juste
  290. c     nbre total de ddl de l'élément considéré
  291.       NLIGRD = MLRE(1+NBENRJ)
  292.       NLIGRP = MLRE(1+NBENRJ)
  293. c     nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  294.       NBNI = (MLRE(1+NBENRJ)) / IDIM
  295.  
  296. c      write(*,*) 'EF',J,' NBENRJ=',NBENRJ,
  297. c     &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'
  298. c
  299. C++++ ON CALCULE XDDL ++++
  300.       MPTVAL = IVADEP
  301.       NCOSOU = IVAL(/1)
  302.  
  303.       INITOT = 0
  304. C-->> BOUCLE SUR LES niveaux d'ENRICHISSEMENTS (0:U 1:A 2:BCDE1 3:BCDE2)
  305.       DO IENR=0,NBENRJ
  306. *nbre de fonction(s) de ce niveau d'enrichissement (=1 si std ou H, =4 pour F1,2,...)
  307.         if(IENR .le. 1) then
  308.          NBNIENR = 1
  309.         else
  310.          NBNIENR = 4
  311.         endif
  312. C---->> BOUCLE SUR LES fonctions de forme de ce type d'enrichissement
  313.         DO INI=1,NBNIENR
  314.           INITOT = INITOT + 1
  315. C------>> BOUCLE SUR LA DIMENSION
  316.           DO 2220 KDIM=1,IDIM
  317.             ICOMP = (INITOT-1)*IDIM + KDIM
  318.  
  319. c         --cas ou on n'a pas trouvé assez de composantes--
  320.             if(ICOMP.GT.NCOSOU) GOTO 2221
  321.  
  322.             MELVAL = IVAL(ICOMP)
  323. c         --cas ou on a pas trouvé cette composante dans cette zone du
  324. c           chpoint solution devenu mchaml  --
  325.             if(MELVAL.eq.0)then
  326. c             Avait on besoin de cette composante?
  327. c             oui, si c'est une composante obligatoire
  328.               if(IENR.eq.0) goto 2991
  329. c             oui, si l'un des noeuds est enrichi
  330.               do I=1,NBNN1
  331.                 if(TLREEL(IENR,I).ne.0) goto 2991
  332.               enddo
  333. c             non, si c'est facultatif et qu'on n'est pas enrichi -> on saute
  334.               goto 2220
  335. c           ->AVERTISSEMENT puis on saute
  336.  2991         NBERR1=NBERR1+1
  337.               if(IIMPI.lt.1) goto 2220
  338. c               write(IOIMP,*) 'PB OPERATEUR SIGM :'
  339.               write(IOIMP,991) ICOMP,IENR,INI,KDIM
  340.   991         format(2X,'ABSENCE DANS LE CHPOINT DEPLACEMENT DE LA ',I3,
  341.      $                  'ieme COMPOSANTE (enrichissement',I3,
  342.      $                   ', fonction',I3,', direction ',I3,
  343.      $                   ') NECESSAIRE POUR L UN DES NOEUDS SUIVANTS :')
  344.               write(IOIMP,*)'  noeuds :',(NUM(iou,J),iou=1,NBNN1)
  345.               goto 2220
  346.             endif
  347.  
  348. C---------->> BOUCLE SUR LES NOEUDS
  349.             DO I=1,NBNN1
  350.               IQ = ((INITOT-1)*NBNN1*IDIM) + ((I-1)*IDIM) + KDIM
  351.               XDDL(IQ) = VELCHE(I,J)
  352.             ENDDO
  353.  
  354.  2220     CONTINUE
  355.         ENDDO
  356.       ENDDO
  357.  
  358. c   --cas normal (toutes les composantes souhaitees etaient presentes)--
  359.       GOTO 2223
  360.  
  361. c   --cas ou on n'a pas trouvé assez de composantes--
  362.  2221 CONTINUE
  363.       if (IIMPI.ge.1) then
  364.        WRITE(IOIMP,2222) J,NCOSOU,ICOMP
  365.  2222  FORMAT(2X,'ATTENTION : ELEMENT ',I6,
  366.      & ' LE CHAMP DE DEPLACEMENT CONTIENT ',I3,' COMPOSANTES',
  367.      & ' PAR NOEUD AU LIEU DE ',I3,' ATTENDUES')
  368.       endif
  369.       NDDL=NCOSOU*NBNN1
  370.       NBENRJ=IENR
  371.  
  372.  2223 CONTINUE
  373. c   --cas ou on a une ou des erreurs--
  374.       IF (NBERR1.gt.0.and.J.eq.NBEL1) THEN
  375.         write(IOIMP,*) 'OPERATEUR GRAD : ABSENCE DANS LE CHPOINT ',
  376.      &   'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE'
  377.       ENDIF
  378.  
  379. c
  380. c rem: il serait probablement interessant au niveau du tems cpu
  381. c  d'utiliser moins de pts de Gauss lorsque l element est élastique
  382. c  On pourrait par ex. utiliser MINT2 = infele(12) qui contient
  383. c  le segment d'integration de l'EF std (QUA4 par ex.)
  384. *      if((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  385. *         MINTE = MINT2
  386. *         NBPGAU= NGAU2
  387. *      else
  388.          MINTE = MINT1
  389.          NBPGAU= NGAU1
  390. *      endif
  391. c
  392. c       pour les def quadratiques
  393.         ISDJC=0
  394.         NBPTEL=NBPGAU
  395.         IF (IREPS2.EQ.1) SEGINI MTRACE
  396. c
  397. C
  398. C*********************************************************
  399. C
  400. C>>>>>>>>>>  BOUCLE SUR LES POINTS DE GAUSS
  401. C
  402.       DO 2500 KGAU=1,NBPGAU
  403. C
  404. C*********************************************************
  405. C     Initialisation à 0
  406. C*********************************************************
  407.  
  408. c ZERO ne serait-il pas facultatif?
  409.       CALL ZERO(SHPWRK,6,NBNO)
  410. C
  411.       i6zz = 3
  412.       IF (IDIM.EQ.3) i6zz = 4
  413. c      do ienr7=1,NBENRMAX
  414.       do ienr7=1,NBENRJ
  415.        do inod7=1,NBNN1
  416. c        do i6=1,6
  417.         do i6=1,i6zz
  418.            LV7WRK(ienr7,1,i6,inod7) = 0
  419.            LV7WRK(ienr7,2,i6,inod7) = 0
  420.         enddo
  421.        enddo
  422.       enddo
  423.  
  424.  
  425. c*********************************************************
  426. c     Calcul des fonction de forme std dans repere local
  427. c*********************************************************
  428.  
  429. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  430. c     (et donc sur les Ni std)
  431.       DO 2510 I=1,NBNN1
  432.  
  433. C++++ Calcul des Ni std
  434. c     (rappel: 1:Ni  2:Ni,qsi  3:Ni,eta  avec i=1,4)
  435.       SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  436.       SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  437.       SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  438.       IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  439.  
  440.  2510 CONTINUE
  441. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  442.  
  443.  
  444.  
  445. c*********************************************************
  446. c     Passage des fonctions de forme std dans repere global
  447. c*********************************************************
  448.  
  449. C++++ CALCUL DES  Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  450.       CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  451. c      if(J.eq.1.and.KGAU.eq.1)
  452. c     &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  453.  
  454. c*********************************************************
  455. c     Si on est pas du tout enrichi on peut sauter une grosse partie
  456. c*********************************************************
  457.       if(NBENRJ.eq.0) goto 2999
  458.  
  459. c*********************************************************
  460. c     Calcul des level set + leurs derivees dans repere global
  461. c*********************************************************
  462.  
  463. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  464.       do 2520 ienr=1,NBENRJ
  465.  
  466. c       MELVA1=MCHAM1.IELVAL(IENR)
  467. c       segact,MELVA1
  468.  
  469. ccccc  BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  470.        DO 2521 I=1,NBNN1
  471.  
  472. C++++ Le I eme noeud est-il ienr-enrichi?
  473.         mlree1= TLREEL(IENR,I)
  474.         if(mlree1.eq.0)  goto 2521
  475.  
  476.  
  477. c       Calcul du repere local de fissure(=PSI,PHI)
  478. c       (rappel: 1,1:psi  1,2:phi  2,1 psi,x  ...etc...)
  479.         do 2522 inode=1,NBNN1
  480. c         pour le H-enrichissement, on n a pas gardé PSI (inutile)
  481.           if(ienr.ne.1) then
  482. c          valeur de PSI au inode^ieme noeud
  483.            xpsi1 = mlree1.PROG(inode)
  484. c          qu on multiplie par la valeur de Ni^std au pt de G considéré
  485.            LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  486.      &      + (SHPWRK(1,inode)*xpsi1)
  487.            LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  488.      &      + (SHPWRK(2,inode)*xpsi1)
  489.            LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  490.      &      + (SHPWRK(3,inode)*xpsi1)
  491.       IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  492.      &      + (SHPWRK(4,inode)*xpsi1)
  493. c          valeur de PHI au inode^ieme noeud
  494.            xphi1 = mlree1.PROG(NBNN1+inode)
  495.           else
  496.            xphi1 = mlree1.PROG(inode)
  497.           endif
  498.           LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  499.      &      + (SHPWRK(1,inode)*xphi1)
  500.           LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  501.      &      + (SHPWRK(2,inode)*xphi1)
  502.           LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  503.      &      + (SHPWRK(3,inode)*xphi1)
  504.       IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  505.      &      + (SHPWRK(4,inode)*xphi1)
  506.  2522   continue
  507.  
  508.  2521  continue
  509. ccccc  fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  510.  
  511.  
  512.  2520 CONTINUE
  513. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  514.  
  515. c     on a construit
  516. C     LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  517.  
  518. c*********************************************************
  519. c     Ajout des fonctions de forme d enrichissement
  520. c     + leurs derivees dans repere global
  521. c*********************************************************
  522.       CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  523.  
  524. c     retour a la partie commune aux EF enrichi et non enrichi
  525.  2999 continue
  526.  
  527. C*********************************************************
  528. C       CALCUL DE B'
  529. C*********************************************************
  530. c ZERO ne serait-il pas facultatif?
  531. c        call ZERO(BGENE,LHOOK,NLIGRP)
  532.         KB=1
  533. C       boucle sur tous les Ni
  534.         DO 3001 II=1,NBNI
  535.  
  536.           BGENE(1,KB)   = SHPWRK(2,II)
  537.           BGENE(2,KB+1) = SHPWRK(3,II)
  538.           BGENE(4,KB)   = SHPWRK(3,II)
  539.           BGENE(4,KB+1) = SHPWRK(2,II)
  540.  
  541.         IF(IDIM.EQ.3) THEN
  542.         BGENE(3,KB+2)=SHPWRK(4,II)
  543.         BGENE(5,KB)=SHPWRK(4,II)
  544.         BGENE(5,KB+2)=SHPWRK(2,II)
  545.         BGENE(6,KB+1)=SHPWRK(4,II)
  546.         BGENE(6,KB+2)=SHPWRK(3,II)
  547.         ENDIF
  548.  
  549.           KB = KB + IDIM
  550.  
  551.  3001   CONTINUE
  552. C
  553. c      if(J.eq.5.and.KGAU.eq.1) then
  554. c      if(KGAU.eq.1) then
  555. c      write(*,*) 'BGENE(1,..)=',(BGENE(1,iou),iou=1,2*NBNI)
  556. c      write(*,*) 'BGENE(2,..)=',(BGENE(2,iou),iou=1,2*NBNI)
  557. c      write(*,*) 'BGENE(4,..)=',(BGENE(3,iou),iou=1,2*NBNI)
  558. c      endif
  559.  
  560. C*********************************************************
  561. C       CALCUL DE D
  562. C*********************************************************
  563.  
  564. c         on cherche les matrices de Hooke
  565.           MPTVAL=IVAMAT
  566.           IF (IMAT.EQ.2) THEN
  567.             MELVAL=IVAL(1)
  568.             IBMN=MIN(J   ,IELCHE(/2))
  569.             IGMN=MIN(KGAU,IELCHE(/1))
  570.             MLREEL=IELCHE(IGMN,IBMN)
  571.             SEGACT MLREEL
  572.             IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  573.      1        CALL DOHOOO(PROG,LHOOK,DDHOOK)
  574.  
  575.           ELSE IF (IMAT.EQ.1) THEN
  576.             DO 9004 IM=1,NMATT
  577.               IF (IVAL(IM).NE.0) THEN
  578.                 MELVAL=IVAL(IM)
  579.                 IBMN=MIN(J   ,VELCHE(/2))
  580.                 IGMN=MIN(KGAU,VELCHE(/1))
  581.                 VALMAT(IM)=VELCHE(IGMN,IBMN)
  582.               ELSE
  583.                 VALMAT(IM)=0.D0
  584.               ENDIF
  585.  9004       CONTINUE
  586. c
  587.             IF (MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN
  588. *              IF (IGAU.LE.NBGMAT)
  589. *     1          CALL DOHMAO(VALMAT,CMATE,IFOUR,IDIM,TXR,XLOC,XGLOB,D1HO,
  590. *     2                      ROTH,DDHOOK,LHOOK,1,IRTD)
  591.                write(*,*) 'cas non prevu a ce jour pour les EF xfem'
  592.                return
  593.             ELSE
  594.               IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  595.      1          CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRTD)
  596.             ENDIF
  597.           ENDIF
  598. c
  599. c
  600. C*********************************************************
  601. C       CALCUL DE sigma = D*B*u
  602. C*********************************************************
  603. c
  604.           CALL DBST(BGENE,DDHOOK,XDDL,(NBNI*IDIM),LHOOK,XSTRS)
  605. c
  606. c       calcul des eps 2 (termes quadratiques de la deformation)
  607. c
  608.           IF (IREPS2.EQ.1)
  609.      1      CALL DBST2(SHPWRK,DDHOOK,XDDL,XE,NBNI,IFOUR,LHOOK,XSTRS,
  610.      2                 TRACE,KGAU,XDPGE,YDPGE,UDPGE,NIFOUR)
  611. c
  612. c       remplissage du segment contenant les contraintes
  613. c
  614.           MPTVAL=IVASTR
  615.           DO 7004 ICOMP=1,LHOOK
  616.             MELVAL=IVAL(ICOMP)
  617.             IBMN=MIN(J,VELCHE(/2))
  618.             VELCHE(KGAU,IBMN)=XSTRS(ICOMP)
  619.  7004     CONTINUE
  620. c
  621. *      if(J.eq.5 .and. KGAU.eq.1) then
  622. c      if(KGAU.eq.1) then
  623. c      write(*,*) J,KGAU,'sigma(..)=',(XSTRS(iou),iou=1,LHOOK)
  624. c      endif
  625. c
  626.  2500 CONTINUE
  627. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  628. C
  629. c     quelques suppressions
  630. c     (ici element non-incompressible=> pas besoin de MTRACE (cf epsi2)
  631.       IF (IREPS2.EQ.1) THEN
  632.         SEGSUP MTRACE
  633.       ENDIF
  634. c
  635.  2000 CONTINUE
  636. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  637. c
  638. C*********************************************************
  639. C       SUPPRESSION ET DESACTIVATION DE SEGMENTS
  640. C*********************************************************
  641. C
  642.       SEGSUP,WRK1,WRK2,MVELCH
  643.  
  644.       SEGSUP,MRACC
  645.  
  646.       RETURN
  647.       END
  648.  
  649.  
  650.  
  651.  

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