Télécharger la plaquette Cast3M


fr en

Télécharger rigix.eso

Retour à la liste

Numérotation des lignes : 

rigix

  1. C RIGIX     SOURCE    PV090527  26/04/28    21:16:21     12529          
  2.       subroutine RIGIX (ivamat,ivacar,NMATT,CMATE,
  3.      &   IMODEL,IPINF,LOCIRI,NBGMAT,NELMAT,IMAT)
  4. C*********************************************************
  5. c
  6. C     PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM
  7. c     POUR LE CALCUL DE RIGIDITE ELEMENTAIRES
  8. C
  9. C*********************************************************
  10. C Liste des modifications :
  11. C -as 2009/09/03 : ajout de IMAT en entrée de RIGIX pour le traitement
  12. C                  d'une matrice de Hooke en entrée (cas IMAT=2)
  13. C -YT 2010/02/11 : Developpement XFEM 3D
  14. C -BP 2011/01    : modification des boucles pour eviter trop de segadj
  15. C -BP 2013/10    : ajout de DIM3 pour les contraintes planes
  16. C -GG 2017/09    : Modification de mise à 0 les DDL non enrichis
  17. C
  18. C*********************************************************
  19. C       PARTIE DECLARATIVE
  20. C*********************************************************
  21.  
  22.       IMPLICIT INTEGER(I-N)
  23.       IMPLICIT REAL*8(A-H,O-Z)
  24.  
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. -INC SMCOORD
  28. -INC SMMODEL
  29. -INC SMCHAML
  30. -INC SMINTE
  31. -INC SMRIGID
  32. -INC SMELEME
  33. -INC SMLREEL
  34. -INC CCHAMP
  35. C
  36.       POINTEUR    MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE
  37.  
  38. -INC TMPTVAL
  39.  
  40.       CHARACTER*8 CMATE
  41. CYT   PARAMETER  (NDDLMAX=20,NBNIMAX=10)
  42.       PARAMETER  (NDDLMAX=30,NBNIMAX=10)
  43.       CHARACTER*4 MOTINC(NDDLMAX),MOTDUA(NDDLMAX)
  44.       DATA MOTINC/'UX  ','UY  ','UZ  ','AX  ','AY  ','AZ  ',
  45.      >'B1X ','B1Y ','B1Z ','C1X ','C1Y ','C1Z ','D1X ','D1Y ',
  46.      >'D1Z ','E1X ','E1Y ','E1Z ','B2X ','B2Y ','B2Z ','C2X ',
  47.      >'C2Y ','C2Z ','D2X ','D2Y ','D2Z ','E2X ','E2Y ','E2Z '/
  48.       DATA MOTDUA/'FX  ','FY  ', 'FZ  ','FAX ','FAY ','FAZ ',
  49.      >'FB1X','FB1Y','FB1Z','FC1X','FC1Y','FC1Z','FD1X','FD1Y',
  50.      >'FD1Z','FE1X','FE1Y','FE1Z','FB2X','FB2Y','FB2Z','FC2X',
  51.      >'FC2Y','FC2Z','FD2X','FD2Y','FD2Z','FE2X','FE2Y','FE2Z'/
  52. C
  53. CYT    DATA MOTINC/'UX  ','UY  ','AX  ','AY  ',
  54. CYT   >'B1X ','B1Y ','C1X ','C1Y ','D1X ','D1Y ','E1X ','E1Y ',
  55. CYT   >'B2X ','B2Y ','C2X ','C2Y ','D2X ','D2Y ','E2X ','E2Y '/
  56. CYT    DATA MOTDUA/  'FX  ','FY  ','FAX ','FAY ',
  57. CYT   >'FB1X','FB1Y','FC1X','FC1Y','FD1X','FD1Y','FE1X','FE1Y',
  58. CYT   >'FB2X','FB2Y','FC2X','FC2Y','FD2X','FD2Y','FE2X','FE2Y'/
  59.  
  60.       PARAMETER  (NBENRMAX=30)
  61. C      NBENRMAX deja defini dans rigixr.eso
  62.       INTEGER    LOCIRI(10,(1+NBENRMAX))
  63.       DIMENSION  MLRE(NBENRMAX+1)
  64.  
  65. C  Segment (type LISTENTI) contenant les informations sur un element
  66.  
  67.       SEGMENT WRK1
  68.        REAL*8 DDHOOK(LHOOK,LHOOK),XE(3,NBBB)
  69.        REAL*8 REL(LRE,LRE),RINT(LRE,LRE)
  70.       ENDSEGMENT
  71. C
  72.       SEGMENT WRK2
  73.        REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE)
  74. c       REAL*8 LV7WRK(NBENRMA2,2,6,NBNO)
  75.        REAL*8 LV7WRK(NBENRMA2,2,6,NBBB)
  76.       ENDSEGMENT
  77. C
  78.       SEGMENT MVELCH
  79.          REAL*8 VALMAT(NV1)
  80.       ENDSEGMENT
  81. C
  82.       SEGMENT MRACC
  83.          INTEGER TLREEL(NBENRMA2,NBI)
  84.          INTEGER MELRIG(NBELEM)
  85.       ENDSEGMENT
  86.  
  87. C*********************************************************
  88. c
  89. C       RECUPERATION + ACTIVATIONS + VALEURS UTILES
  90. c
  91. C*********************************************************
  92. C     sont deja actifs dans rigi1.eso :
  93. C           SEGACT MMODEL,IMODEL,MELVAL
  94. c           IVAMAT
  95.       MPTVAL=IVAMAT
  96.  
  97. C++++ Recup + Activation de la geometrie ++++++++++++++++
  98.       MELEME= IMAMOD
  99.       SEGACT MELEME
  100.  
  101. C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++
  102. c   + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1
  103. C     segment INFO deja actif dans RIGI1
  104. c   + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait
  105. c     mais du coup, on na pas de segment minte
  106. c     (car depend de si pt de g pour rigi, pour sigma....)
  107. c     c'est + simple de rappeler elquoi ici
  108.       MELE = NEFMOD
  109.       NGAU1 = INFELE(6)
  110.       LRE   = INFELE(9)
  111.       LHOOK = INFELE(10)
  112.       MINT1 = INFMOD(6)
  113.       segact,MINT1
  114.       MINT2 = INFMOD(10)
  115.       MFR   = INFELE(13)
  116.       IELE  = INFELE(14)
  117.       NDDL  = INFELE(15)
  118.       NSTRS = INFELE(16)
  119. c      write(ioimp,*)'-> RIGIX INFELE',(INFELE(iou),iou=1,16)
  120.  
  121. c   + AUTRES INFOS
  122. C     nbre de noeuds par element
  123.       NBNN1 = NUM(/1)
  124. C     nbre d elements
  125.       NBEL1 = NUM(/2)
  126.  
  127. c REM: pour se passer du dimensionnement du nbre d'enrichissement dans
  128. c      elquoi et le realiser localement , on pourrait ecrire:
  129. c      LRE = NDDLMAX*NBNN1
  130. c      NDDL= NDDLMAX
  131.  
  132. C     sous decoupage et points de Gauss de l element geometrique de base
  133. CYT     if(MELE.eq.263)  then
  134. cbp      if(MELE.eq.263.or.MELE.eq.264)  then
  135.       if (MINT2.ne.0) then
  136.           segact,MINT2
  137.           NGAU2 = MINT2.POIGAU(/1)
  138.       else
  139.           NGAU2=0
  140.       endif
  141. c      write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3))
  142. c      write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU2)
  143.  
  144.  
  145. c     nbre maxi de fonction de forme par noeud (fonction std comprise)
  146. c      NBNI = NDDL/IDIM    inutile!
  147.  
  148.  
  149. C++++ Recup des infos d enrichissement +++++++++++++++++++
  150. c     recup du MCHEX1 d'enrichissement
  151.       NBENR1=0
  152.       MCHAM1=0
  153.       NOBMO1=IVAMOD(/1)
  154.       if (NOBMO1.ne.0) then
  155.        do iobmo1=1,NOBMO1
  156.         if ((TYMODE(iobmo1)).eq.'MCHAML') then
  157.          MCHEX1 = IVAMOD(iobmo1)
  158.          segact,MCHEX1
  159.           if ((MCHEX1.TITCHE).eq.'ENRICHIS') then
  160.             MCHAM1 = MCHEX1.ICHAML(1)
  161.             segact,MCHAM1
  162.             goto 1000
  163.           endif
  164.         endif
  165.        enddo
  166.        write(ioimp,*) 'Le modele est vide (absence d enrichissement)'
  167. *       return
  168.       else
  169.       write(ioimp,*) 'Aucun MCHAML enrichissement dans le Modele'
  170. *       return
  171.       endif
  172.  
  173.  1000 continue
  174. c     niveau d enrichissement(s) du modele (ddls std u exclus)
  175. c     NBENR1= 0 si std, 1 si H seul,  2 si H+F1,  3 si H+F1+F2, etc...
  176.       if (MCHAM1.ne.0) NBENR1=MCHAM1.IELVAL(/1)
  177. c*       write(ioimp,*) 'niveau d enrichissement(s) du modele',NBENR1
  178.  
  179. C*********************************************************
  180. c
  181. c     PREPARATION DES OBJETS RESULTATS
  182. c
  183. C*********************************************************
  184. C
  185. C++++ RECHERCHE DES NOMS D'INCONNUES ET DES DUAUX
  186. c     MOTINC et MODUA en dur pour l instant
  187.  
  188. C++++ REMPLISSAGE DE DESCR de taille maxi
  189. c     (on ajustera en enlevant si besoin)
  190.       NLIGRP = LRE
  191.       NLIGRD = LRE
  192.       SEGINI DESCR
  193.       IPDSCR = DESCR
  194. C
  195.       KINC = 0
  196. C----->> BOUCLE SUR LES fonctions de Forme
  197.       DO  IENR=1,NBNIMAX
  198. C------->> BOUCLE SUR LES NOEUDS
  199.         DO  I=1,NBNN1
  200. C--------->> BOUCLE SUR LA DIMENSION
  201.           DO  KDIM=1,IDIM
  202.             KINC = KINC + 1
  203. CYT         LISINC(KINC) = MOTINC((IDIM*(IENR-1))+KDIM)
  204. CYT         LISDUA(KINC) = MOTDUA((IDIM*(IENR-1))+KDIM)
  205.             LISINC(KINC) = MOTINC((3*(IENR-1))+KDIM)
  206.             LISDUA(KINC) = MOTDUA((3*(IENR-1))+KDIM)
  207.             NOELEP(KINC) = I
  208.             NOELED(KINC) = I
  209.           ENDDO
  210.         ENDDO
  211.       ENDDO
  212.  
  213.       IF (KINC.NE.LRE)   THEN
  214.         WRITE(ioimp,*) 'IL Y A UNE ERREUR DANS DESCR'
  215.         WRITE(ioimp,*) 'KINC , LRE = ',KINC,LRE
  216.         RETURN
  217.       ENDIF
  218. C
  219.       SEGDES DESCR
  220. C
  221. C
  222. C++++ INITIALISATIONS...
  223. C
  224. c     ... des tables LOCIRI et MLRE
  225. c     MLRE contient le nombre d'inconnues de chaque sous-zone
  226. c     determinee depuis le nombre de fonctions de forme
  227. c ienr=  1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10)
  228.       if (NBENR1.ne.0) then
  229.       do ienr=1,NBENR1
  230.          nbniJ = 2 + ((ienr-1)*4)
  231.          MLRE(1+ienr) = IDIM*NBNN1*nbniJ
  232. c        -LOCIRI
  233.          LOCIRI(5,1+ienr)= NIFOUR
  234.       enddo
  235.       endif
  236. C     Tables + longues car 1er indice correspond au fontion de forme std
  237.       MLRE(1)    = IDIM*NBNN1*1
  238.       LOCIRI(5,1)= NIFOUR
  239. c     on complete avec des 0
  240.       if (NBENR1.lt.(NBENRMAX+1)) then
  241.         do ienr=(NBENR1+1),(NBENRMAX)
  242.            MLRE(1+ienr) = 0
  243.         enddo
  244.       endif
  245. c
  246. c      ...DU SEGMENT WRK1
  247.       NBENRMA2 = NBENRMAX
  248.       NBBB = NBNN1
  249.       SEGINI,WRK1
  250.  
  251. c      ...DU SEGMENT WRK2
  252. c      NBNO = NBNI
  253.       NBNO = LRE/IDIM
  254.       SEGINI,WRK2
  255.  
  256. c      ...DU SEGMENT MVELCH
  257.       NV1 = NMATT
  258.       SEGINI,MVELCH
  259. C
  260. c      ...DU SEGMENT MRACC
  261.       NBENRMA2 = NBENRMAX
  262.       NBI = NBNN1
  263.       NBELEM = NUM(/2)
  264.       segini,MRACC
  265. C
  266. C
  267.  
  268. C*********************************************************
  269. C
  270. C>>>>>>>>>>>>>>>>>>>>>>>>>>> 1ere BOUCLE SUR LES ELEMENTS
  271. C
  272.       NBENR = 0
  273.       DO 2000 J=1,NBEL1
  274. c       write(ioimp,*) '===boucle1=== EF',J,' /NBEL1 ================'
  275. C
  276. C
  277. C*********************************************************
  278. C     POUR CHAQUE ELEMENT, ...
  279. C*********************************************************
  280. C
  281. C++++ NBENRJ = niveau d enrichissement de l element ++++
  282. C     =0 si EF std   =1 si U+H   =2 si U+H+F1   =3 si U+H+F1+F2
  283.       NBENRJ=0
  284.       if (NBENR1.ne.0) then
  285.       do IENR=1,NBENR1
  286.         MELVA1 = MCHAM1.IELVAL(IENR)
  287.         segact,MELVA1
  288.         do I=1,NBNN1
  289.          mlree1 = MELVA1.IELCHE(I,J)
  290.          if(mlree1.ne.0)  NBENRJ=max(NBENRJ,IENR)
  291.         enddo
  292.       enddo
  293.       endif
  294.       NBENR=max(NBENRJ,NBENR)
  295. c
  296. c++++ Selon le niveau d enrichissement NBENRJ...
  297. C
  298. c++++ ...on ajuste la 3eme dimension de XMATRI.RE
  299. c     (=nbre d element de ce niveau d enrichissement)
  300. C     on la stocke temporairement dans LOCIRI(4, ) meme si c est pas propre
  301.       LOCIRI(4,1+NBENRJ) = LOCIRI(4,1+NBENRJ) + 1
  302. c
  303. C++++ ...on copie cet element a la geo correspondante
  304.       IPT1 = LOCIRI(1,1+NBENRJ)
  305. c     si elle n existe pas, il faut la creer
  306.       if (IPT1 .eq. 0) then
  307.          NBNN   = NBNN1
  308.          NBELEM = 1
  309.          NBSOUS = 0
  310.          NBREF  = 0
  311.          segini,IPT1
  312.          IPT1.ITYPEL = ITYPEL
  313.          LOCIRI(1,1+NBENRJ)=IPT1
  314. c     si elle existe, il faut l agrandir
  315.       else
  316.          NBNN   = NBNN1
  317.          NBELEM = (IPT1.NUM(/2)) + 1
  318.          NBSOUS = 0
  319.          NBREF  = 0
  320.          segadj,IPT1
  321.       endif
  322.  
  323. c     Ajout de la connectivite de l'element courant
  324.       do I=1,NBNN1
  325.          IPT1.NUM(I,NBELEM) = NUM(I,J)
  326.       enddo
  327.  
  328. c     on retient la position NBELEM de l element J
  329.       MELRIG(J)= NBELEM
  330.  
  331. C++++ ...on crée le descr qui va bien si il n existe pas
  332.       DES1 = LOCIRI(3,1+NBENRJ)
  333.       if (DES1.eq.0) then
  334.          segini,DES1=DESCR
  335.          NLIGRP = MLRE(1+NBENRJ)
  336.          NLIGRD = MLRE(1+NBENRJ)
  337.          segadj,DES1
  338.          LOCIRI(3,1+NBENRJ) = DES1
  339. *  on remet en ro des1
  340.          segact,DES1
  341.       endif
  342.  
  343.  2000 CONTINUE
  344. C FIN DE 1ere BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  345. C
  346. C*********************************************************
  347.  
  348.  
  349. C*********************************************************
  350. C     POUR CHAQUE NIVEAU D ENRICHISSEMENT, ...
  351. C*********************************************************
  352.       do IENR=1,(NBENR1+1)
  353. c        attention : IENR = IENRvrai + 1
  354.          NELRIG = LOCIRI(4,IENR)
  355. C++++ ...ON CRÉE XMATRI DE LA BONNE TAILLE TOUT DE SUITE
  356.          if (NELRIG.ne.0) then
  357.            NLIGRP = MLRE(IENR)
  358.            NLIGRD = MLRE(IENR)
  359.            rigrel=0
  360.            segini,XMATRI
  361.            LOCIRI(4,IENR)=XMATRI
  362. *  on remet en ro ipt1
  363. C++++ ...ON DESACTIVE IPT1
  364.            IPT1 = LOCIRI(1,IENR)
  365.            segact ipt1
  366.          endif
  367.       enddo
  368.  
  369.  
  370.  
  371. C*********************************************************
  372. C
  373. C>>>>>>>>>>>>>>>>>>>>>>>>>>> 2eme BOUCLE SUR LES ELEMENTS
  374. C
  375.       DO 2001 J=1,NBEL1
  376. c       write(*,*) '===boucle 2=== EF',J,' /NBEL1 ================'
  377. C
  378. C*********************************************************
  379. C     POUR CHAQUE ELEMENT, ...
  380. C*********************************************************
  381. C
  382. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  383.       CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  384.  
  385. C++++ NBENRJ = niveau d enrichissement de l element ++++
  386. C     =0 si EF std   =1 si U+H   =2 si U+H+F1   =3 si U+H+F1+F2
  387.       NBENRJ=0
  388.       if (NBENR1.ne.0) then
  389.       do IENR=1,(NBENR1)
  390.         MELVA1 = MCHAM1.IELVAL(IENR)
  391.         segact,MELVA1
  392.         do I=1,NBNN1
  393.          mlree1 = MELVA1.IELCHE(I,J)
  394. C        on en profite pour remplir MRACC table de raccourcis pour cet element
  395.          TLREEL(IENR,I) = mlree1
  396.          if (mlree1.ne.0)  then
  397.             NBENRJ = max(NBENRJ,IENR)
  398. C           et on active la listreel
  399.             segact,mlree1
  400.          endif
  401.         enddo
  402.       enddo
  403.       endif
  404.       if (NBENRMA2.gt.NBENR1) then
  405.         do IENR2=(NBENR1+1),NBENRMA2
  406.         do I=1,NBNN1
  407.            TLREEL(IENR2,I) = 0
  408.         enddo
  409.         enddo
  410.       endif
  411. C
  412. c++++ on recupere XMATRI
  413.       XMATRI = LOCIRI(4,1+NBENRJ)
  414.  
  415. c++++ quelques indices pour dimensionner au plus juste
  416. c     nbre total de ddl de l'élément considéré
  417.       NLIGRD = MLRE(1+NBENRJ)
  418.       NLIGRP = MLRE(1+NBENRJ)
  419. c     nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  420.       NBNI = (MLRE(1+NBENRJ)) / IDIM
  421. c     position de cet element pour cet enrichissement
  422.       IELRIG = MELRIG(J)
  423. c     write(*,*) 'EF',J,' NBENRJ=',NBENRJ,
  424. c    &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'
  425. c
  426.  
  427. C++++ ON MET A ZERO LA SOMME QUE L ON VA CALCULER
  428. c      CALL ZERO(RINT,LRE,NLIGRP)
  429. *      CALL ZERO(RINT,NLIGRD,NLIGRP)
  430. c |-> impossible car ZERO considere qu il sagit de la dimension de SHPWRK
  431.        CALL ZEROX(RINT,NLIGRD,NLIGRP,LRE)
  432. c
  433. c REM : il est interessant au niveau du tems cpu d'utiliser moins de pts
  434. c       de Gauss lorsque l element n est pas enrichi car cela n'apporte
  435. c       rien de plus.
  436. c       On utilise MINT2 = INFELE(12) qui contient le segment d'integration
  437. C       de l'EF std (QUA4 par ex. pour XQ4R)
  438.       if ((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  439.          MINTE = MINT2
  440.          NBPGAU= NGAU2
  441.       else
  442.          MINTE = MINT1
  443.          NBPGAU= NGAU1
  444.       endif
  445.  
  446. C*********************************************************
  447. C
  448. C>>>>>>>>>>  BOUCLE SUR LES POINTS DE GAUSS
  449. C
  450.       DO 2500 KGAU=1,NBPGAU
  451. c
  452. c          write(ioimp,*) '-pt de G ',KGAU,' / ',NBPGAU,'----'
  453. C
  454. C*********************************************************
  455. C     Initialisation à 0
  456. C*********************************************************
  457.  
  458. c      CALL ZERO(SHPWRK,6,NBNO)
  459.       CALL ZERO(SHPWRK,6,NBNI)
  460. C
  461.       i6zz = 3
  462.       IF (IDIM.EQ.3) i6zz = 4
  463. C
  464. c      do ienr7=1,NBENRMAX
  465.       do ienr7=1,NBENRJ
  466.        do inod7=1,NBNN1
  467. c        do i6=1,6
  468. CYT     do i6=1,3
  469.         do i6=1,i6zz
  470.            LV7WRK(ienr7,1,i6,inod7) = 0.D0
  471.            LV7WRK(ienr7,2,i6,inod7) = 0.D0
  472.         enddo
  473.        enddo
  474.       enddo
  475.  
  476.  
  477. c*********************************************************
  478. c     Calcul des fonction de forme std dans repere local
  479. c*********************************************************
  480.  
  481. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  482. c     (et donc sur les Ni std)
  483.       DO 2510 I=1,NBNN1
  484.  
  485. C++++ Calcul des Ni std
  486. c     (rappel: 1:Ni  2:Ni,qsi  3:Ni,eta  avec i=1,4)
  487.       SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  488.       SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  489.       SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  490.       IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  491.  
  492.  2510 CONTINUE
  493. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  494.  
  495.  
  496.  
  497. c*********************************************************
  498. c     Passage des fonctions de forme std dans repere global
  499. c*********************************************************
  500.  
  501. C++++ CALCUL DES  Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  502.       CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  503. c      if(J.eq.1.and.KGAU.eq.1)
  504. c     &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  505.  
  506. c*********************************************************
  507. c     Si on est pas du tout enrichi on peut sauter une grosse partie
  508. c*********************************************************
  509.       if(NBENRJ.eq.0) goto 2999
  510.  
  511. c*********************************************************
  512. c     Calcul des level set + leurs derivees dans repere global
  513. c*********************************************************
  514.  
  515. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  516.       do 2520 ienr=1,NBENRJ
  517.  
  518. c       MELVA1=MCHAM1.IELVAL(IENR)
  519. c       segact,MELVA1
  520.  
  521. ccccc  BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  522.        DO 2521 I=1,NBNN1
  523.  
  524. C++++ Le I eme noeud est-il ienr-enrichi?
  525.         mlree1= TLREEL(IENR,I)
  526.         if(mlree1.eq.0)  goto 2521
  527.  
  528.  
  529. c       Calcul du repere local de fissure(=PSI,PHI)
  530. c       (rappel: 1,1:psi  1,2:phi  2,1 psi,x  ...etc...)
  531.         do 2522 inode=1,NBNN1
  532. c         pour le H-enrichissement, on n a pas gardé PSI (inutile)
  533.           if (ienr.ne.1) then
  534. c          valeur de PSI au inode^ieme noeud
  535.            xpsi1 = mlree1.PROG(inode)
  536. c          qu on multiplie par la valeur de Ni^std au pt de G considéré
  537.            LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  538.      &      + (SHPWRK(1,inode)*xpsi1)
  539.            LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  540.      &      + (SHPWRK(2,inode)*xpsi1)
  541.            LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  542.      &      + (SHPWRK(3,inode)*xpsi1)
  543.       IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  544.      &      + (SHPWRK(4,inode)*xpsi1)
  545. c          valeur de PHI au inode^ieme noeud
  546.            xphi1 = mlree1.PROG(NBNN1+inode)
  547.           else
  548.            xphi1 = mlree1.PROG(inode)
  549.           endif
  550.           LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  551.      &      + (SHPWRK(1,inode)*xphi1)
  552.           LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  553.      &      + (SHPWRK(2,inode)*xphi1)
  554.           LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  555.      &      + (SHPWRK(3,inode)*xphi1)
  556.       IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  557.      &      + (SHPWRK(4,inode)*xphi1)
  558.  2522   continue
  559.  
  560.  2521  continue
  561. ccccc  fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  562.  
  563.  
  564.  2520 CONTINUE
  565. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  566.  
  567. c     on a construit
  568. C     LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  569.  
  570.  
  571. c*********************************************************
  572. c     Ajout des fonctions de forme d enrichissement
  573. c     + leurs derivees dans repere global
  574. c*********************************************************
  575.       CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  576.  
  577. c     retour a la partie commune aux EF enrichi et non enrichi
  578.  2999 continue
  579.  
  580. C*********************************************************
  581. c      write(*,*) 'C       CALCUL DE B'
  582. C*********************************************************
  583.         call ZERO(BGENE,LHOOK,LRE)
  584.         KB=1
  585. C       boucle sur tous les Ni
  586.         DO 3001 II=1,NBNI
  587.  
  588.           BGENE(1,KB)   = SHPWRK(2,II)
  589.           BGENE(2,KB+1) = SHPWRK(3,II)
  590.           BGENE(4,KB)   = SHPWRK(3,II)
  591.           BGENE(4,KB+1) = SHPWRK(2,II)
  592.  
  593.           IF (IDIM.EQ.3) THEN
  594.             BGENE(3,KB+2)=SHPWRK(4,II)
  595.             BGENE(5,KB)  =SHPWRK(4,II)
  596.             BGENE(5,KB+2)=SHPWRK(2,II)
  597.             BGENE(6,KB+1)=SHPWRK(4,II)
  598.             BGENE(6,KB+2)=SHPWRK(3,II)
  599.           ENDIF
  600.           KB = KB + IDIM
  601.  3001   CONTINUE
  602.  
  603.  
  604. C*********************************************************
  605. c      write(*,*) 'C       CALCUL DE D (=MATRICE DE HOOKE)'
  606. C*********************************************************
  607. c
  608. c+++++cas DOHOOO
  609.       IF (IMAT.EQ.2) THEN
  610.         MELVAL=IVAL(1)
  611.         IBMN=MIN(J  ,IELCHE(/2))
  612.         IGMN=MIN(KGAU,IELCHE(/1))
  613.         MLREEL=IELCHE(IGMN,IBMN)
  614.         SEGACT MLREEL
  615.         IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  616.      1  CALL DOHOOO(PROG,LHOOK,DDHOOK)
  617.  
  618. c+++++cas DOHMAS
  619.       ELSE IF (IMAT.EQ.1) THEN
  620.         DO 4001 IM=1,NMATT
  621.           IF (IVAL(IM).NE.0) THEN
  622.             MELVAL=IVAL(IM)
  623.             IBMN=MIN(J  ,VELCHE(/2))
  624.             IGMN=MIN(KGAU,VELCHE(/1))
  625.             if (ibmn.gt.0.and.igmn.gt.0) then
  626.               VALMAT(IM)=VELCHE(IGMN,IBMN)
  627.             else
  628.               VALMAT(IM)=0.D0
  629.             endif
  630.           ELSE
  631.             VALMAT(IM)=0.D0
  632.           ENDIF
  633.  4001   CONTINUE
  634. C
  635.          CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRET)
  636.          IF (IRET.EQ.0) THEN
  637.           MOTERR(1:8)=CMATE
  638.           MOTERR(9:16)=NOMFR(MFR/2+1)
  639.           INTERR(1)=IFOUR
  640.           CALL ERREUR(81)
  641.           RETURN
  642.          ENDIF
  643.        ENDIF
  644. C
  645. C*********************************************************
  646. c      write(*,*) 'C       CALCUL DE  Bt.D.B'
  647. C*********************************************************
  648.  
  649. C       PRISE EN COMPTE DU POIDS D'INTEGRATION
  650.         DJAC = ABS(DJAC) * POIGAU(KGAU)
  651.  
  652. c       CAS DES CONTRAINTES PLANES
  653. C       recuperation de l'epaisseur: DIM3 donnée facultative du materiau
  654.         IF (IFOUR.EQ.-2) THEN
  655.           MPTVAL=IVACAR
  656.           IF (IVACAR.NE.0) THEN
  657.             MELVAL=IVAL(1)
  658.             IF (MELVAL.NE.0) THEN
  659.               IGMN=MIN(KGAU,VELCHE(/1))
  660.               IBMN=MIN(J,VELCHE(/2))
  661.               DIM3=VELCHE(IGMN,IBMN)
  662.             ELSE
  663.               DIM3=1.D0
  664.             ENDIF
  665.           ENDIF
  666.           DJAC=DJAC*DIM3
  667.           MPTVAL=IVAMAT
  668.         ENDIF
  669.  
  670.         CALL ZEROX(REL,NLIGRD,NLIGRP,LRE)
  671.  
  672.         CALL BDBST(BGENE,DJAC,DDHOOK,LRE,LHOOK,REL)
  673. *        CALL BDBST(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,REL)
  674. c |-> impossible car BDBST considere qu il sagit de la dimension de bgene
  675. c     pour gagner du temps cpu il faudrait qqchose du type:
  676. c     CALL BDBSTX(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,LRE,REL)
  677.  
  678. C*********************************************************
  679. C       CUMUL DANS RINT(II,JJ)
  680. C*********************************************************
  681.       DO  II=1,NLIGRD
  682.          DO  JJ=1,NLIGRP
  683.            RINT(II,JJ) =  RINT(II,JJ) + REL(II,JJ)
  684.       enddo
  685.       enddo
  686. c
  687. c
  688.  2500 CONTINUE
  689. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  690. C
  691. C*********************************************************
  692. C      write(*,*) 'C       STOCKAGE DANS XMATRI de RE(II,JJ)'
  693. c       et de XMATRI dans IMATRI
  694. C*********************************************************
  695. C
  696.        DO 6001 II=1,NLIGRD
  697.  
  698.          DO 6002 JJ=1,NLIGRP
  699.  
  700. c            XMATRI.RE(II,JJ,nelrig) =  RINT(II,JJ)
  701. c            RE(II,JJ,NELRIG) =  RINT(II,JJ)
  702.            RE(II,JJ,IELRIG) =  RINT(II,JJ)
  703.  
  704. c     si NON enrichi, on met tout a 0 sauf la diago
  705. c          if(mlree1.eq.0)  RE(II,JJ) = 0.      inutile
  706.  
  707.  6002    CONTINUE
  708.  
  709.  6001 CONTINUE
  710. c
  711. c
  712.  2001 CONTINUE
  713. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  714.  
  715.       SEGSUP,WRK1,WRK2,MVELCH
  716.       SEGSUP,MRACC
  717. c
  718. c     desactivation des maillages et segment imatri
  719.       do ienr=1,(1+NBENR1)
  720.          XMATRI = LOCIRI(4,ienr)
  721.          if(XMATRI.ne.0)  segdes,XMATRI
  722.       enddo
  723.  
  724. c      SEGDES dans RIGI1 =   IMODEL
  725.  
  726. C      RETURN
  727.       END
  728.  
  729.  
  730.  
  731.  
  732.  

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