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    MB234859  25/09/08    21:16:07     12358          
  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.            segini,XMATRI
  360.            LOCIRI(4,IENR)=XMATRI
  361. *  on remet en ro ipt1
  362. C++++ ...ON DESACTIVE IPT1
  363.            IPT1 = LOCIRI(1,IENR)
  364.            segact ipt1
  365.          endif
  366.       enddo
  367.  
  368.  
  369.  
  370. C*********************************************************
  371. C
  372. C>>>>>>>>>>>>>>>>>>>>>>>>>>> 2eme BOUCLE SUR LES ELEMENTS
  373. C
  374.       DO 2001 J=1,NBEL1
  375. c       write(*,*) '===boucle 2=== EF',J,' /NBEL1 ================'
  376. C
  377. C*********************************************************
  378. C     POUR CHAQUE ELEMENT, ...
  379. C*********************************************************
  380. C
  381. C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT
  382.       CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE)
  383.  
  384. C++++ NBENRJ = niveau d enrichissement de l element ++++
  385. C     =0 si EF std   =1 si U+H   =2 si U+H+F1   =3 si U+H+F1+F2
  386.       NBENRJ=0
  387.       if (NBENR1.ne.0) then
  388.       do IENR=1,(NBENR1)
  389.         MELVA1 = MCHAM1.IELVAL(IENR)
  390.         segact,MELVA1
  391.         do I=1,NBNN1
  392.          mlree1 = MELVA1.IELCHE(I,J)
  393. C        on en profite pour remplir MRACC table de raccourcis pour cet element
  394.          TLREEL(IENR,I) = mlree1
  395.          if (mlree1.ne.0)  then
  396.             NBENRJ = max(NBENRJ,IENR)
  397. C           et on active la listreel
  398.             segact,mlree1
  399.          endif
  400.         enddo
  401.       enddo
  402.       endif
  403.       if (NBENRMA2.gt.NBENR1) then
  404.         do IENR2=(NBENR1+1),NBENRMA2
  405.         do I=1,NBNN1
  406.            TLREEL(IENR2,I) = 0
  407.         enddo
  408.         enddo
  409.       endif
  410. C
  411. c++++ on recupere XMATRI
  412.       XMATRI = LOCIRI(4,1+NBENRJ)
  413.  
  414. c++++ quelques indices pour dimensionner au plus juste
  415. c     nbre total de ddl de l'élément considéré
  416.       NLIGRD = MLRE(1+NBENRJ)
  417.       NLIGRP = MLRE(1+NBENRJ)
  418. c     nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim
  419.       NBNI = (MLRE(1+NBENRJ)) / IDIM
  420. c     position de cet element pour cet enrichissement
  421.       IELRIG = MELRIG(J)
  422. c     write(*,*) 'EF',J,' NBENRJ=',NBENRJ,
  423. c    &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme'
  424. c
  425.  
  426. C++++ ON MET A ZERO LA SOMME QUE L ON VA CALCULER
  427. c      CALL ZERO(RINT,LRE,NLIGRP)
  428. *      CALL ZERO(RINT,NLIGRD,NLIGRP)
  429. c |-> impossible car ZERO considere qu il sagit de la dimension de SHPWRK
  430.        CALL ZEROX(RINT,NLIGRD,NLIGRP,LRE)
  431. c
  432. c REM : il est interessant au niveau du tems cpu d'utiliser moins de pts
  433. c       de Gauss lorsque l element n est pas enrichi car cela n'apporte
  434. c       rien de plus.
  435. c       On utilise MINT2 = INFELE(12) qui contient le segment d'integration
  436. C       de l'EF std (QUA4 par ex. pour XQ4R)
  437.       if ((NBENRJ.eq.0).and.(MINT2.ne.0)) then
  438.          MINTE = MINT2
  439.          NBPGAU= NGAU2
  440.       else
  441.          MINTE = MINT1
  442.          NBPGAU= NGAU1
  443.       endif
  444.  
  445. C*********************************************************
  446. C
  447. C>>>>>>>>>>  BOUCLE SUR LES POINTS DE GAUSS
  448. C
  449.       DO 2500 KGAU=1,NBPGAU
  450. c
  451. c          write(ioimp,*) '-pt de G ',KGAU,' / ',NBPGAU,'----'
  452. C
  453. C*********************************************************
  454. C     Initialisation à 0
  455. C*********************************************************
  456.  
  457. c      CALL ZERO(SHPWRK,6,NBNO)
  458.       CALL ZERO(SHPWRK,6,NBNI)
  459. C
  460.       i6zz = 3
  461.       IF (IDIM.EQ.3) i6zz = 4
  462. C
  463. c      do ienr7=1,NBENRMAX
  464.       do ienr7=1,NBENRJ
  465.        do inod7=1,NBNN1
  466. c        do i6=1,6
  467. CYT     do i6=1,3
  468.         do i6=1,i6zz
  469.            LV7WRK(ienr7,1,i6,inod7) = 0.D0
  470.            LV7WRK(ienr7,2,i6,inod7) = 0.D0
  471.         enddo
  472.        enddo
  473.       enddo
  474.  
  475.  
  476. c*********************************************************
  477. c     Calcul des fonction de forme std dans repere local
  478. c*********************************************************
  479.  
  480. ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc
  481. c     (et donc sur les Ni std)
  482.       DO 2510 I=1,NBNN1
  483.  
  484. C++++ Calcul des Ni std
  485. c     (rappel: 1:Ni  2:Ni,qsi  3:Ni,eta  avec i=1,4)
  486.       SHPWRK(1,I) = SHPTOT(1,I,KGAU)
  487.       SHPWRK(2,I) = SHPTOT(2,I,KGAU)
  488.       SHPWRK(3,I) = SHPTOT(3,I,KGAU)
  489.       IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU)
  490.  
  491.  2510 CONTINUE
  492. ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  493.  
  494.  
  495.  
  496. c*********************************************************
  497. c     Passage des fonctions de forme std dans repere global
  498. c*********************************************************
  499.  
  500. C++++ CALCUL DES  Ni,x Ni,y (i=1,4) + CALCUL DE det(J)
  501.       CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC)
  502. c      if(J.eq.1.and.KGAU.eq.1)
  503. c     &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1)
  504.  
  505. c*********************************************************
  506. c     Si on est pas du tout enrichi on peut sauter une grosse partie
  507. c*********************************************************
  508.       if(NBENRJ.eq.0) goto 2999
  509.  
  510. c*********************************************************
  511. c     Calcul des level set + leurs derivees dans repere global
  512. c*********************************************************
  513.  
  514. ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc
  515.       do 2520 ienr=1,NBENRJ
  516.  
  517. c       MELVA1=MCHAM1.IELVAL(IENR)
  518. c       segact,MELVA1
  519.  
  520. ccccc  BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc
  521.        DO 2521 I=1,NBNN1
  522.  
  523. C++++ Le I eme noeud est-il ienr-enrichi?
  524.         mlree1= TLREEL(IENR,I)
  525.         if(mlree1.eq.0)  goto 2521
  526.  
  527.  
  528. c       Calcul du repere local de fissure(=PSI,PHI)
  529. c       (rappel: 1,1:psi  1,2:phi  2,1 psi,x  ...etc...)
  530.         do 2522 inode=1,NBNN1
  531. c         pour le H-enrichissement, on n a pas gardé PSI (inutile)
  532.           if (ienr.ne.1) then
  533. c          valeur de PSI au inode^ieme noeud
  534.            xpsi1 = mlree1.PROG(inode)
  535. c          qu on multiplie par la valeur de Ni^std au pt de G considéré
  536.            LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I)
  537.      &      + (SHPWRK(1,inode)*xpsi1)
  538.            LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I)
  539.      &      + (SHPWRK(2,inode)*xpsi1)
  540.            LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I)
  541.      &      + (SHPWRK(3,inode)*xpsi1)
  542.       IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I)
  543.      &      + (SHPWRK(4,inode)*xpsi1)
  544. c          valeur de PHI au inode^ieme noeud
  545.            xphi1 = mlree1.PROG(NBNN1+inode)
  546.           else
  547.            xphi1 = mlree1.PROG(inode)
  548.           endif
  549.           LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I)
  550.      &      + (SHPWRK(1,inode)*xphi1)
  551.           LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I)
  552.      &      + (SHPWRK(2,inode)*xphi1)
  553.           LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I)
  554.      &      + (SHPWRK(3,inode)*xphi1)
  555.       IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I)
  556.      &      + (SHPWRK(4,inode)*xphi1)
  557.  2522   continue
  558.  
  559.  2521  continue
  560. ccccc  fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc
  561.  
  562.  
  563.  2520 CONTINUE
  564. ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc
  565.  
  566. c     on a construit
  567. C     LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
  568.  
  569.  
  570. c*********************************************************
  571. c     Ajout des fonctions de forme d enrichissement
  572. c     + leurs derivees dans repere global
  573. c*********************************************************
  574.       CALL SHAPX(MELE,LV7WRK,NBENRMAX,NBENRJ,TLREEL,SHPWRK,IRET)
  575.  
  576. c     retour a la partie commune aux EF enrichi et non enrichi
  577.  2999 continue
  578.  
  579. C*********************************************************
  580. c      write(*,*) 'C       CALCUL DE B'
  581. C*********************************************************
  582.         call ZERO(BGENE,LHOOK,LRE)
  583.         KB=1
  584. C       boucle sur tous les Ni
  585.         DO 3001 II=1,NBNI
  586.  
  587.           BGENE(1,KB)   = SHPWRK(2,II)
  588.           BGENE(2,KB+1) = SHPWRK(3,II)
  589.           BGENE(4,KB)   = SHPWRK(3,II)
  590.           BGENE(4,KB+1) = SHPWRK(2,II)
  591.  
  592.           IF (IDIM.EQ.3) THEN
  593.             BGENE(3,KB+2)=SHPWRK(4,II)
  594.             BGENE(5,KB)  =SHPWRK(4,II)
  595.             BGENE(5,KB+2)=SHPWRK(2,II)
  596.             BGENE(6,KB+1)=SHPWRK(4,II)
  597.             BGENE(6,KB+2)=SHPWRK(3,II)
  598.           ENDIF
  599.           KB = KB + IDIM
  600.  3001   CONTINUE
  601.  
  602.  
  603. C*********************************************************
  604. c      write(*,*) 'C       CALCUL DE D (=MATRICE DE HOOKE)'
  605. C*********************************************************
  606. c
  607. c+++++cas DOHOOO
  608.       IF (IMAT.EQ.2) THEN
  609.         MELVAL=IVAL(1)
  610.         IBMN=MIN(J  ,IELCHE(/2))
  611.         IGMN=MIN(KGAU,IELCHE(/1))
  612.         MLREEL=IELCHE(IGMN,IBMN)
  613.         SEGACT MLREEL
  614.         IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1))
  615.      1  CALL DOHOOO(PROG,LHOOK,DDHOOK)
  616.  
  617. c+++++cas DOHMAS
  618.       ELSE IF (IMAT.EQ.1) THEN
  619.         DO 4001 IM=1,NMATT
  620.           IF (IVAL(IM).NE.0) THEN
  621.             MELVAL=IVAL(IM)
  622.             IBMN=MIN(J  ,VELCHE(/2))
  623.             IGMN=MIN(KGAU,VELCHE(/1))
  624.             if (ibmn.gt.0.and.igmn.gt.0) then
  625.               VALMAT(IM)=VELCHE(IGMN,IBMN)
  626.             else
  627.               VALMAT(IM)=0.D0
  628.             endif
  629.           ELSE
  630.             VALMAT(IM)=0.D0
  631.           ENDIF
  632.  4001   CONTINUE
  633. C
  634.          CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRET)
  635.          IF (IRET.EQ.0) THEN
  636.           MOTERR(1:8)=CMATE
  637.           MOTERR(9:16)=NOMFR(MFR/2+1)
  638.           INTERR(1)=IFOUR
  639.           CALL ERREUR(81)
  640.           RETURN
  641.          ENDIF
  642.        ENDIF
  643. C
  644. C*********************************************************
  645. c      write(*,*) 'C       CALCUL DE  Bt.D.B'
  646. C*********************************************************
  647.  
  648. C       PRISE EN COMPTE DU POIDS D'INTEGRATION
  649.         DJAC = ABS(DJAC) * POIGAU(KGAU)
  650.  
  651. c       CAS DES CONTRAINTES PLANES
  652. C       recuperation de l'epaisseur: DIM3 donnée facultative du materiau
  653.         IF (IFOUR.EQ.-2) THEN
  654.           MPTVAL=IVACAR
  655.           IF (IVACAR.NE.0) THEN
  656.             MELVAL=IVAL(1)
  657.             IF (MELVAL.NE.0) THEN
  658.               IGMN=MIN(KGAU,VELCHE(/1))
  659.               IBMN=MIN(J,VELCHE(/2))
  660.               DIM3=VELCHE(IGMN,IBMN)
  661.             ELSE
  662.               DIM3=1.D0
  663.             ENDIF
  664.           ENDIF
  665.           DJAC=DJAC*DIM3
  666.           MPTVAL=IVAMAT
  667.         ENDIF
  668.  
  669.         CALL ZEROX(REL,NLIGRD,NLIGRP,LRE)
  670.  
  671.         CALL BDBST(BGENE,DJAC,DDHOOK,LRE,LHOOK,REL)
  672. *        CALL BDBST(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,REL)
  673. c |-> impossible car BDBST considere qu il sagit de la dimension de bgene
  674. c     pour gagner du temps cpu il faudrait qqchose du type:
  675. c     CALL BDBSTX(BGENE,DJAC,DDHOOK,NLIGRP,LHOOK,LRE,REL)
  676.  
  677. C*********************************************************
  678. C       CUMUL DANS RINT(II,JJ)
  679. C*********************************************************
  680.       DO  II=1,NLIGRD
  681.          DO  JJ=1,NLIGRP
  682.            RINT(II,JJ) =  RINT(II,JJ) + REL(II,JJ)
  683.       enddo
  684.       enddo
  685. c
  686. c
  687.  2500 CONTINUE
  688. C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<<
  689. C
  690. C*********************************************************
  691. C      write(*,*) 'C       STOCKAGE DANS XMATRI de RE(II,JJ)'
  692. c       et de XMATRI dans IMATRI
  693. C*********************************************************
  694. C
  695.        DO 6001 II=1,NLIGRD
  696.  
  697.          DO 6002 JJ=1,NLIGRP
  698.  
  699. c            XMATRI.RE(II,JJ,nelrig) =  RINT(II,JJ)
  700. c            RE(II,JJ,NELRIG) =  RINT(II,JJ)
  701.            RE(II,JJ,IELRIG) =  RINT(II,JJ)
  702.  
  703. c     si NON enrichi, on met tout a 0 sauf la diago
  704. c          if(mlree1.eq.0)  RE(II,JJ) = 0.      inutile
  705.  
  706.  6002    CONTINUE
  707.  
  708.  6001 CONTINUE
  709. c
  710. c
  711.  2001 CONTINUE
  712. C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<<
  713.  
  714.       SEGSUP,WRK1,WRK2,MVELCH
  715.       SEGSUP,MRACC
  716. c
  717. c     desactivation des maillages et segment imatri
  718.       do ienr=1,(1+NBENR1)
  719.          XMATRI = LOCIRI(4,ienr)
  720.          if(XMATRI.ne.0)  segdes,XMATRI
  721.       enddo
  722.  
  723. c      SEGDES dans RIGI1 =   IMODEL
  724.  
  725. C      RETURN
  726.       END
  727.  
  728.  
  729.  
  730.  

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