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

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