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

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