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

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