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

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