Télécharger la plaquette Cast3M


fr en

Télécharger accro2.eso

Retour à la liste

Numérotation des lignes : 

accro2

  1. C ACCRO2    SOURCE    CB215821  24/04/12    21:15:01     11897          
  2.       SUBROUTINE ACCRO2
  3. C========================================================================
  4. C   Cree la matrice de liaison entre le champ u et w
  5. C   avec une formulation faible a 3 champs : u w et t
  6. C   FAIBLE => On utilise les fonctions de forme
  7. c
  8. C   Option PENA et STAB pour créer les matrices de penalite Kww
  9. C   et de stabilisation Ktt (utile si LATIN utilisée)
  10. c
  11. C   En sortie, on a Kfaibl = [Ktt Ktu Ktw ; Kut 0 0 ; Kwt 0 Kww]
  12. C   [0  Ktu Ktw ; sym. ] = \int t * (u-w) ds
  13. c
  14. c   ajout des enrichissements XFEM
  15. c
  16. C   Creation : Benoit Trolle (SNCF-Lamcos), avril 2012
  17. C   Integration dans Cast3m : BP, decembre 2012
  18. C========================================================================
  19.  
  20.       IMPLICIT INTEGER(I-N)
  21.       IMPLICIT REAL*8 (A-H,O-Z)
  22.  
  23. -INC CCREEL
  24.  
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. -INC SMCOORD
  28. -INC SMELEME
  29. -INC SMRIGID
  30. -INC SMLMOTS
  31. -INC CCHAMP
  32. -INC CCGEOME
  33. -INC SMMODEL
  34. -INC SMINTE
  35. -INC SMCHAML
  36. -INC SMLREEL
  37.  
  38.  
  39. C  Segment contenant les informations sur un element
  40.       SEGMENT INFO
  41.         INTEGER INFELL(JG)
  42.       ENDSEGMENT
  43.  
  44.       POINTEUR INFO1.INFO,INFO2.INFO
  45.       POINTEUR    MCHEX1.MCHELM
  46.       external shape
  47.       DATA EPSI/1.D-9/
  48.       DIMENSION ICOR(6),IMEL(6),imtt(10)
  49.       DIMENSION QSI(3),XPO(3)
  50.  
  51. C  INITIALISATION DES INCONNUES obligatoires et facultatives
  52.       PARAMETER (NOBL=3,NFAC=9,NMOCLE=3)
  53.       CHARACTER*4 DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2
  54.       CHARACTER*4 DUAOBL(NOBL),DUAFAC(NFAC),MOCLE(NMOCLE)
  55.       CHARACTER*4 DDLOB2(NOBL),DDLFA2(NOBL)
  56.       CHARACTER*4 DUAOB2(NOBL),DUAFA2(NOBL)
  57.  
  58.       DATA DDLOBL/'UX  ','UY  ','UZ  '/
  59.       DATA DDLFAC/'AX  ','AY  ','AZ  ',
  60.      >'B1X ','B1Y ','B1Z ',
  61.      >'B2X ','B2Y ','B2Z '/
  62.       DATA DUAOBL/'FX  ','FY  ','FZ  '/
  63.       DATA DUAFAC/'FAX ','FAY ','FAZ ',
  64.      >'FB1X','FB1Y','FB1Z',
  65.      >'FB2X','FB2Y','FB2Z'/
  66.  
  67. c     nom de composantes speciales
  68.       DATA DDLOB2/'UX  ','UY  ','UZ  '/
  69.       DATA DDLFA2/'AX  ','AY  ','AZ  '/
  70.       DATA DUAOB2/'FX  ','FY  ','FZ  '/
  71.       DATA DUAFA2/'FAX ','FAY ','FAZ '/
  72.  
  73.       DATA MOCLE /'PENA','STAB','NOMU'/
  74.  
  75. c-----------------------------------------------------
  76. c     Segment de travail
  77. c-----------------------------------------------------
  78.       SEGMENT MTRAV
  79.         REAL*8 SHPP2(6,NBNN2)
  80.         REAL*8 XE2(3,NBNN2),XEL2(3,NBNN2),NGENE2(1,NBNN2)
  81.         REAL*8 XE3(3,NBNN2)
  82.         REAL*8 SHPP1(6,NBNN1)
  83.         REAL*8 XE1(3,NBNN1),XEL1(3,NBNN1)
  84.         REAL*8 NGENE1(1,NBNN1)
  85.         REAL*8 REL(nbnn2,nbnn2),REL1(nbnn2,nbnn1)
  86.         REAL*8 RINT(nbnn2,nbnn2)
  87.         REAL*8 RINT0(nbnn2,ngau2*nbnn1),RINT1(nbnn2,ngau2*nbnn1)
  88.         REAL*8 RINT2(nbnn2,ngau2*nbnn1),RINT3(nbnn2,ngau2*nbnn1)
  89.         INTEGER ITYMAT(2,4*nbnn1*ngau2)
  90.         INTEGER ITAMYT(ngau2*NBNN1,4)
  91.         INTEGER PRELEM(3*nbnn1*ngau2)
  92.         INTEGER ElemPG(ngau2)
  93.       ENDSEGMENT
  94.  
  95. c      tableaux comptant le nbre d EF de chaque ddl
  96.       PARAMETER(NDDLMAX=6)
  97.       INTEGER NELDDL(NDDLMAX)
  98.  
  99.       if(iimpi.ge.2) then
  100.       write(ioimp,*) '-----------------------------------------------'
  101.       write(ioimp,*) '               ENTREE dans ACCRO2'
  102.       write(ioimp,*) '-----------------------------------------------'
  103.       endif
  104.  
  105.       SEGACT MCOORD
  106.  
  107. c-----------------------------------------------------
  108. c     LECTURE DES MOT CLES PENAlite et STABilistation
  109. c-----------------------------------------------------
  110.       RLATIN=0.d0
  111.       XISTAB=0.d0
  112.       INOMU=0
  113.  700  ICLE=0
  114.       CALL LIRMOT(MOCLE,NMOCLE,ICLE,0)
  115.       IF (ICLE.EQ.1) GOTO 701
  116.       IF (ICLE.EQ.2) GOTO 702
  117.       IF (ICLE.EQ.3) THEN
  118.         INOMU=1
  119.         GOTO 700
  120.       ENDIF
  121.       GOTO 799
  122.  
  123. c     LECTURE DE LA DIRECTION DE RECHERCHE (FLOTTANT)
  124.  701  CALL LIRREE(RLATIN,1,IRETOU)
  125.       IF (IERR.NE.0) RETURN
  126.       if(iimpi.ge.2) write(ioimp,*) 'On a lu la raideur :',RLATIN
  127.       GOTO 700
  128.  
  129. c     LECTURE DE LA STABILISATION (FLOTTANT)
  130.  702  CALL LIRREE(XISTAB,1,IRETOU)
  131.       IF (IERR.NE.0) RETURN
  132.       if(iimpi.ge.2) write(ioimp,*) 'On a lu la STABILISATION :',XISTAB
  133.       GOTO 700
  134.  
  135.  799  CONTINUE
  136.  
  137.  
  138. c-----------------------------------------------------
  139. c     RECUPERATION DU MAILLAGE MASSIF
  140. c-----------------------------------------------------
  141.       CALL LIROBJ('MMODEL  ',IPMODL,1,IRETOU)
  142.       CALL ACTOBJ('MMODEL  ',IPMODL,1)
  143.       IF(IERR.NE.0) RETURN
  144.       MMODE1=IPMODL
  145.       segact,MMODE1
  146. c  Récupération du nombre de zones du modèle
  147.       N1 = MMODE1.KMODEL(/1)
  148.       if(N1.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
  149.       IMODE1 = MMODE1.KMODEL(1)
  150. C Récupération du maillage et du numéro d'élément du modèle
  151.       segact,IMODE1
  152.       nele1 = IMODE1.NEFMOD
  153.       IPT1 = IMODE1.IMAMOD
  154.       SEGACT IPT1
  155. C Récupération du numéro d'élément du maillage, du nombre de noeuds et d'éléments
  156.       iele1 = IPT1.itypel
  157.       nbnn1 = IPT1.num(/1)
  158.       nbel1 = IPT1.num(/2)
  159. c   récupération des caractéristique EF IPT1
  160. cbp : on essaie d'abord avec INFELE
  161.       mele1 = IMODE1.IMAMOD
  162.       MINTE1 = IMODE1.INFELE(11)
  163. cbp : si echec, on retente avec elquoi
  164.       if(mele1.eq.0.or.MINTE1.eq.0) then
  165.         call elquoi(nele1,0,3,IPTR1,IMODE1)
  166.         INFO = IPTR1
  167.         segact,INFO
  168.         mele1 =  INFELL(1)
  169.         MINTE1 = INFELL(11)
  170.         segdes,INFO
  171.       endif
  172. cbp : si toujours pb, erreur...
  173.       if(mele1.eq.0.or.MINTE1.eq.0) then
  174.         write(ioimp,*) 'donnees relatives a l element fini inconnues'
  175.         call erreur(591)
  176.       endif
  177.  
  178.  
  179. c-----------------------------------------------------
  180. c     RECUPERATION DU MAILLAGE INTERFACE
  181. c-----------------------------------------------------
  182.       IPMAI2 = 0
  183.       CALL LIROBJ('MMODEL  ',IPMODL,0,IRETOU)
  184.       IF(IERR.NE.0) RETURN
  185.  
  186. C    -DANS LE CAS OU ON A UN MODELE EN ENTREE
  187.       IF (IRETOU.EQ.1)THEN
  188.          CALL ACTOBJ('MMODEL  ',IPMODL,1)
  189.          MMODE2=IPMODL
  190.          segact,MMODE2
  191.          N2 = MMODE2.KMODEL(/1)
  192.          if(N2.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!'
  193.          IMODE2 = MMODE2.KMODEL(1)
  194.          segact,IMODE2
  195.          nele2 = IMODE2.NEFMOD
  196.          IPT2 = IMODE2.IMAMOD
  197.          SEGACT IPT2
  198. c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
  199.          iele2 = IPT2.itypel
  200.          nbnn2 = IPT2.num(/1)
  201.          nbel2 = IPT2.num(/2)
  202.  
  203. c   recuperation des caracteritiques de l'element
  204. cbp : on essaie d'abord avec INFELE
  205.          MINTE2 = IMODE2.INFELE(11)
  206.          ngau2  = IMODE2.INFELE(6)
  207. cbp : si echec, on retente avec elquoi
  208.          if(MINTE2.eq.0.or.ngau2.eq.0) then
  209. c          write(ioimp,*) 'appel elquoi',iele2,nele2
  210.            call elquoi(nele2,0,3,IPTR2,IMODE2)
  211.            INFO = IPTR2
  212.            segact,INFO
  213.            MINTE2 = INFELL(11)
  214.            ngau2 = INFELL(6)
  215.            segdes,INFO
  216.          endif
  217. cbp : si toujours pb, erreur...
  218.         if(mele1.eq.0.or.MINTE1.eq.0) then
  219.            write(ioimp,*) 'donnees relatives a l element fini INTERFACE'
  220.      &     ,'inconnues'
  221.            call erreur(591)
  222.          endif
  223.  
  224. C    -DANS LE CAS OU ON A UN MAILLAGE EN ENTREE
  225.       ELSE
  226.          CALL LIROBJ('MAILLAGE',IPMAI2,1,IRETOU)
  227.          CALL ACTOBJ('MAILLAGE',IPMAI2,1)
  228.          IF(IERR.NE.0) RETURN
  229.          IPT2 = IPMAI2
  230.          SEGACT IPT2
  231. c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso)
  232.          iele2 = IPT2.itypel
  233.          nele2 = iele2
  234.          if (nele2.lt.2.or.nele2.gt.26) then
  235.          write(ioimp,*)'element geometrique different de l element fini'
  236.             call erreur(16)
  237.          endif
  238.          nbnn2 = IPT2.num(/1)
  239.          nbel2 = IPT2.num(/2)
  240.          ngau2 = -3
  241. c     SEG2
  242.          if (nele2.EQ.2) ngau2 = 2
  243. c     SEG3
  244.          if (nele2.EQ.3) ngau2 = 3
  245. c     TRI3
  246.          if (nele2.EQ.4) ngau2 = 1
  247. c     TRI6
  248.          if (nele2.EQ.6) ngau2 = 4
  249. c     QUA4
  250.          if (nele2.EQ.8) ngau2 = 4
  251. c     QUA8
  252.          if (nele2.EQ.10) ngau2 = 9
  253.          if(ngau2.eq.-3) then
  254.          write(ioimp,*)'nombre de points d integration inconnu pour',
  255.      &         'ce type d element geometrique'
  256.          call erreur(16)
  257.          endif
  258.       ENDIF
  259.  
  260.       call RESHPT(ngau2,nbnn2,iele2,nele2,0,IPTR2,IRET)
  261.       MINTE2 = IPTR2
  262.       segact,MINTE1,MINTE2
  263.  
  264.  
  265. c-----------------------------------------------------
  266. c     RECHERCHE DU MCHAML ISSU MCHEX1 D ENRICHISSEMENT
  267. c-----------------------------------------------------
  268.       MCHAM1=0
  269.       NBENR2=0
  270.       segact,IMODE1
  271.       NOBMOD = IMODE1.IVAMOD(/1)
  272.       IF (NOBMOD.NE.0) THEN
  273.          DO 1002 iobmo1=1,NOBMOD
  274.             if((IMODE1.TYMODE(iobmo1)).ne.'MCHAML')   goto 1002
  275.             MCHEX1 = IMODE1.IVAMOD(iobmo1)
  276.             segact,MCHEX1
  277.             if((MCHEX1.TITCHE).ne.'ENRICHIS')  goto 1003
  278.             MCHAM1 = MCHEX1.ICHAML(1)
  279.             segact,MCHAM1
  280.             NBENR2 = MCHAM1.IELVAL(/1)
  281.             do ienr2=1,NBENR2
  282.                MELVA1=MCHAM1.IELVAL(IENR2)
  283.                if(MELVA1.ne.0) segact,MELVA1
  284.             enddo
  285.  1003       continue
  286.  1002    CONTINUE
  287.       ENDIF
  288.  
  289. c-------------------------------------
  290. c     INITIALISATION DES OBJETS DE TRAVAIL : MTRAV et ITYMAT
  291. c-------------------------------------
  292.  
  293.       segini,MTRAV
  294.  
  295. *bp : NTYMAT = (U ou H ou HB1 ou HB1B2) * dim * nbre de configs possibles
  296. *     NTYMAT = nbre de types de matrices
  297. *      nbre de configs varie entre :
  298. *     -on peut avoir tous les pt de Gauss dans 1 seul elements de structure
  299. *     -jusqu'à chaque pt de Gauss dans un element distinct
  300.       ity=0
  301. c     boucle sur les enrichissements possibles U, H , HB1, HB1B2
  302.       if(iimpi.ge.3) write(ioimp,*) 'ity , ino1 , ienr1'
  303.       do ienr1=1,4
  304. c       boucle sur les nombre de noeuds "structure"
  305.         do ino1=(nbnn1),(ngau2*nbnn1)
  306.           ity=ity+1
  307.           ITYMAT(1,ity)=ino1
  308.           ITYMAT(2,ity)=ienr1
  309.           ITAMYT(ino1,ienr1)=ity
  310.           if(iimpi.ge.3) write(ioimp,*) ity,ino1,ienr1
  311.         enddo
  312.       enddo
  313.       NTYMAT = 4*(nbnn1*(ngau2-1)+1)
  314.       if(iimpi.ge.3) write(ioimp,*) 'ity*idim=',ity,NTYMAT
  315.  
  316.  
  317. c--------------------------------------------------------------------
  318. c     INITIALISATION DU SEGMENT MRIGID
  319. c--------------------------------------------------------------------
  320.       NRIGEL = NTYMAT*idim
  321.       segini,MRIGID
  322.       IFORIG = IFOUR
  323.       MTYMAT ='RIGIDITE'
  324. c    -on prepare le meleme
  325.       NBSOUS = 0
  326.       NBREF = 0
  327.  
  328.       ityty=0
  329. c on initialise la taille matrice en fonction du type de matrice
  330.       do ity=1,NTYMAT
  331.       do iidim=1,idim
  332.  
  333.          ityty=ityty+1
  334.          COERIG(ityty) = 1.D0
  335.  
  336. *     dim de la matrice RE elementaire = nbnoeud*TX + nbnoeud*AX + nbnoeud*AX
  337.          nbno1 = ITYMAT(1,ity)
  338.          nenr1 = ITYMAT(2,ity)
  339.          if(nenr1.le.2) then
  340.            NLIGRP = nbnn2 + nbnn2 + nbno1
  341.          else
  342.            NLIGRP = nbnn2 + nbnn2 + ((nenr1-1)*nbno1)
  343.          endif
  344.          NLIGRD = NLIGRP
  345.  
  346. c       -creation du MELEME
  347.          NBNN = nbnn2 + nbno1
  348.          NBELEM=0
  349.          SEGINI,MELEME
  350.          ITYPEL=28
  351.          IRIGEL(1,ityty) = MELEME
  352.  
  353. c       -remplissage du DESCR
  354.          SEGINI,DESCR
  355.          IRIGEL(3,ityty) = DESCR
  356.          iddl=0
  357. c        remplissage des ddl LX de la fissure (ici appelé FX)
  358. c        on iverse donc inconnue et duale []*FX=UX
  359.          do ino2=1,nbnn2
  360.            iddl=iddl+1
  361.            if (nenr1.eq.1) then
  362.                LISINC(iddl)=DUAOB2(iidim)
  363.                LISDUA(iddl)=DDLOB2(iidim)
  364.            else
  365.              LISINC(iddl)=DUAFA2(iidim)
  366.              LISDUA(iddl)=DDLFA2(iidim)
  367.            endif
  368.            NOELEP(iddl)=ino2
  369.            NOELED(iddl)=ino2
  370.          enddo
  371. c        remplissage des ddl UX de la fissure (ici appelé WX) []*WX=TX
  372.          do ino2=1,nbnn2
  373.            iddl=iddl+1
  374.            if (nenr1.eq.1) then
  375.                LISINC(iddl)=DDLOB2(iidim)
  376.                LISDUA(iddl)=DUAOB2(iidim)
  377.            else
  378.              LISINC(iddl)=DDLFA2(iidim)
  379.              LISDUA(iddl)=DUAFA2(iidim)
  380.            endif
  381.            NOELEP(iddl)=ino2
  382.            NOELED(iddl)=ino2
  383.          enddo
  384. c        remplissage des ddl de la structure
  385.          if (nenr1.eq.1) then
  386.            do ino1=1,nbno1
  387.                iddl=iddl+1
  388.                LISINC(iddl)=DDLOBL(iidim)
  389.                LISDUA(iddl)=DUAOBL(iidim)
  390.                NOELEP(iddl)=nbnn2+ino1
  391.                NOELED(iddl)=nbnn2+ino1
  392.            enddo
  393.          else
  394.            do ini1=1,(nenr1-1)
  395.            do ino1=1,nbno1
  396.                iddl=iddl+1
  397.                LISINC(iddl)=DDLFAC(iidim+(3*(ini1-1)))
  398.                LISDUA(iddl)=DUAFAC(iidim+(3*(ini1-1)))
  399.                NOELEP(iddl)=nbnn2+ino1
  400.                NOELED(iddl)=nbnn2+ino1
  401.            enddo
  402.            enddo
  403.          endif
  404.          if(iimpi.ge.3) write(ioimp,*) ityty,(LISINC(iou),iou=1,NLIGRP)
  405.          SEGDES,DESCR
  406.  
  407. c       -initialisation du XMATRI
  408.          NELRIG=0
  409.          SEGINI,XMATRI
  410.          IRIGEL(4,ityty) = XMATRI
  411.          IRIGEL(5,ityty) = NIFOUR
  412.          IRIGEL(6,ityty) = 0
  413.          IRIGEL(7,ityty) = 0
  414.          IRIGEL(8,ityty) = 0
  415.  
  416.       enddo
  417.       enddo
  418.  
  419.  
  420. c----------------------------------------------------------------------
  421. c     1. RECHERCHE DES ELEMENTS DE STRUCTURE CONTENANT DES POINTS DE GAUSS
  422. c        DES ELEMENTS DE LA FISSURE
  423. c     2. REMPLISSAGE DU MELEME + MRIGID en FONCTION
  424. c----------------------------------------------------------------------
  425.  
  426. C
  427. c==== Boucle sur les elements de fissure ==============================
  428.       DO 1100 iem2=1,nbel2
  429.  
  430.          call doxe(xcoor,idim,nbnn2,ipt2.num,iem2,xe2)
  431.          nbenrj = 0
  432.          nbnno2 = 0
  433.          call ZERO(RINT,nbnn2,nbnn2)
  434.          call ZERO(RINT0,nbnn2,ngau2*nbnn1)
  435.          call ZERO(RINT1,nbnn2,ngau2*nbnn1)
  436.          call ZERO(RINT2,nbnn2,ngau2*nbnn1)
  437.          call ZERO(RINT3,nbnn2,ngau2*nbnn1)
  438.  
  439. c======= Boucle sur les pt de G de fissure ============================
  440.          DO 1132 igau2=1,ngau2
  441.             izok = 0
  442. c   récupération des coordonnees du point de gauss dans le repère global
  443.             XPO(1) = 0.D0
  444.             XPO(2) = 0.D0
  445.             XPO(3) = 0.D0
  446.             DO  ino=1,nbnn2
  447.                XPO(1) = XPO(1) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(1,ino))
  448.                XPO(2) = XPO(2) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(2,ino))
  449.                IF(IDIM.EQ.3) then
  450.                XPO(3) = XPO(3) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(3,ino))
  451.                ENDIF
  452.             ENDDO
  453.  
  454. c            if(iimpi.ge.5) write(ioimp,*) 'Fissure: element',iem2,
  455. c      &     ' point de Gauss',igau2,' x=',XPO(1),XPO(2)
  456.  
  457. c---------- Boucle sur les elements de structure ----------------------
  458.             DO 1131 iem1=1,nbel1
  459.  
  460. c   si pas d'enrichissement, on travaille sur tous les elements
  461.                if(MCHAM1.eq.0) goto 1133
  462. c   on saute les elements non enrichi car a priori ne contiennent pas la fissure
  463.                do ienr2=1,NBENR2
  464.                   MELVA1=MCHAM1.IELVAL(IENR2)
  465.                   if(MELVA1.ne.0) then
  466.                     do inode1=1,nbnn1
  467.                   if(MELVA1.IELCHE(inode1,iem1).ne.0) goto 1133
  468.                     enddo
  469.                   endif
  470.                enddo
  471.                goto 1131
  472.  1133          continue
  473.  
  474. c   recuperation des coordonnées des noeuds de IPT1 : xel1 (dans le repère x,y,z)
  475.                call doxe(xcoor,idim,nbnn1,ipt1.num,iem1,xel1)
  476.  
  477. c   calcul des fonctions de formes de IPT1 au pt de Gauss de IPT2
  478.                call QSIJS(xel1,iele1,nbnn1,idim,XPO,SHPP1,qsi,IRET)
  479. c            if(iimpi.ge.5) write(ioimp,*) 'Structure: element',iem1,
  480. c      &     ' noeuds:',(ipt1.num(iou,iem1),iou=1,nbnn1)
  481. c            if(iimpi.ge.5) write(ioimp,*)
  482. c      &     '  N_i=',(SHPP1(1,iou),iou=1,nbnn1)
  483.  
  484.  
  485. c   test pour savoir si PG est dans EF de IPT1
  486.                DO 1130 ino1=1,NBNN1
  487.                   if (SHPP1(1,ino1).LT.-1.01D-7) then
  488.                      go to 1131
  489.                   endif
  490.  1130          continue
  491.                izok = 1
  492. c              ON a trouvé l'élément de structure !
  493. c           write(ioimp,*) 'ON a trouve l element de structure !'
  494. c     DETECTION DU TYPE D'ENRICHISSEMENT MAX DE CET ELEMENT = nbenrj
  495.                DO 3001 IENR2=1,NBENR2
  496.                  MELVA1=MCHAM1.IELVAL(IENR2)
  497.                  IF(MELVA1.eq.0) GOTO 3001
  498.                  DO 3002 ino1=1,nbnn1
  499.                    MLREEL = MELVA1.IELCHE(ino1,iem1)
  500. c                  Test pour savoir si le noeud est enrichi
  501.                    IF(MLREEL.eq.0) GOTO 3002
  502.                    nbenrj=max(nbenrj,IENR2)
  503.  3002            continue
  504.  3001          continue
  505.  
  506.                if(iimpi.ge.3) write(ioimp,*) 'EF fissure ',iem2,
  507.      &         ' ptdeG ',igau2,' -> EF MASSIF ',iem1,' nbenrj=',nbenrj
  508.  
  509. c        Préparation du Remplissage du MELEME
  510.                do 3003 ino1=1,nbnn1
  511.                   iino1=IPT1.NUM(ino1,iem1)
  512.                   if(nbnno2.gt.0) then
  513.                   do iino2=1,nbnno2
  514.                     if(PRELEM(iino2).eq.iino1) goto 3003
  515.                   enddo
  516.                   endif
  517. c                 nouveau noeud : on l'ajoute
  518.                   nbnno2=nbnno2+1
  519.                   PRELEM(nbnno2)=iino1
  520.  3003          continue
  521.                if(iimpi.ge.3) then
  522.                write(ioimp,*) ' PRELEM=',( PRELEM(iou),iou=1,nbnno2)
  523.                write(ioimp,*) 'EF fissure ',iem2,' ptdeG ',igau2
  524.                endif
  525.  
  526. C        CALCUL DE REL ET REL1
  527.  
  528. c matrices des fonctions de forme de l'element associé à IPT2
  529.  
  530.                do ino2=1,NBNN2
  531.                   xe3(1,ino2) = sqrt( xe2(1,ino2)*xe2(1,ino2)
  532.      1                              + xe2(2,ino2)*xe2(2,ino2) )
  533.                   NGENE2(1,ino2)=MINTE2.SHPTOT(1,ino2,IGAU2)
  534.                   do ipp2=1,6
  535.                      SHPP2(ipp2,ino2)=MINTE2.SHPTOT(ipp2,ino2,IGAU2)
  536.                   enddo
  537.                   if(iimpi.ge.4) then
  538.                   write(ioimp,*) 'SHPP2',(SHPP2(iou,ino2),iou=1,6)
  539.                   write(ioimp,*) 'xe2',(xe2(iou,ino2),iou=1,2)
  540.                   endif
  541.                enddo
  542. c              calcul du jacobien
  543.                if(idim.eq.2) then
  544. c              calcul du jacobien en 2D pour un seg2
  545.                   dXdQsi=0.D0
  546.                   dYdQsi=0.D0
  547.                   DO i=1,NBNN2
  548.                      dXdQsi=dXdQsi+SHPP2(2,i)*XE2(1,i)
  549.                      dYdQsi=dYdQsi+SHPP2(2,i)*XE2(2,i)
  550.                   ENDDO
  551. c           Dans le cas des seg2, poids des pts de gauss = 1
  552.                   DJAC2=SQRT(dXdQsi*dXdQsi+dYdQsi*dYdQsi)
  553.                 else
  554.                       dXdQsi=0.D0
  555.                       dYdQsi=0.D0
  556.                       dZdQsi=0.D0
  557.                       dXdEta=0.D0
  558.                       dYdEta=0.D0
  559.                       dZdEta=0.D0
  560.                       DO i=1,NBNN2
  561.                          dXdQsi=dXdQsi+SHPP2(2,i)*XE2(1,i)
  562.                          dXdEta=dXdEta+SHPP2(3,i)*XE2(1,i)
  563.                          dYdQsi=dYdQsi+SHPP2(2,i)*XE2(2,i)
  564.                          dYdEta=dYdEta+SHPP2(3,i)*XE2(2,i)
  565.                          dZdQsi=dZdQsi+SHPP2(2,i)*XE2(3,i)
  566.                          dZdEta=dZdEta+SHPP2(3,i)*XE2(3,i)
  567.                       ENDDO
  568.                       z = (dXdQsi*dYdEta-dXdEta*dYdQsi)
  569.                       x = (dYdQsi*dZdEta-dYdEta*dZdQsi)
  570.                       y = (dZdQsi*dXdEta-dZdEta*dXdQsi)
  571. c           Dans le cas des tri3, poids du pts de gauss = 0.5
  572.                       DJAC2 = sqrt(x*x+y*y+z*z)*0.5
  573.                endif
  574.  
  575. c matrices des fonctions de forme de l'element associé à IPT1
  576.                do ino1=1,NBNN1
  577.                   NGENE1(1,ino1)=SHPP1(1,ino1)
  578.                enddo
  579.                if(iimpi.ge.4) then
  580.                  write(ioimp,*) 'DJAC2=',DJAC2
  581.                  write(ioimp,*) 'NGENE1=',(NGENE1(1,iou),iou=1,NBNN1)
  582.                  write(ioimp,*) 'NGENE2=',(NGENE2(1,iou),iou=1,NBNN2)
  583.                endif
  584.  
  585. c  calcul des matrices elementaires = produits des fonctions de forme
  586.                DO II=1,nbnn2
  587.                   XGENE2 = NGENE2(1,II)*DJAC2
  588. c calcul de K_ww, Kwt et Ktt
  589.                    DO JJ=1,nbnn2
  590.                       jjj = JJ
  591. C           CC = Jacobien * transposee(N)(II,i) * N(i,JJ) (somme sur i)
  592.                       CC = XGENE2*NGENE2(1,jjj)
  593. C           REL est une matrice symetrique
  594.                       REL(II,JJ) = CC
  595.                       REL(JJ,II) = REL(II,JJ)
  596.                    ENDDO
  597. c calcul de K_tu
  598.                    DO JJ=1,nbnn1
  599. C           CC = Jacobien * transposee(N)(II,i) * N(i,JJ) (somme sur i)
  600.                       CC=XGENE2*NGENE1(1,JJ)
  601. c                     write(ioimp,*) ' CC = ', CC
  602. c                      REL1(II,JJ)=-1.D0*CC
  603.                       REL1(II,JJ)=CC
  604. c                      REL(JJdec,II)=REL(II,JJdec)
  605.                    ENDDO
  606.                   if(iimpi.ge.4) write(ioimp,*)
  607.      1            'REL1(',II,',:)=',(REL1(II,iou),iou=1,NBNN1)
  608.                ENDDO
  609.  
  610. C         CUMUL DANS RINT ET RINT1
  611. C     CUMUL DANS RINT(II,JJ) pour Kww, Kwt, Ktw et Ktt
  612. c     (4 fois la meme matrices = N^{w T} * N^{w} = N^{w T} * N^{t} = ...)
  613.                DO  II=1,nbnn2
  614.                    DO  JJ=1,nbnn2
  615.                          RINT(II,JJ) = RINT(II,JJ) + REL(II,JJ)
  616.                    ENDDO
  617.                ENDDO
  618.  
  619. C     CUMUL DANS RINT1(II,JJ) pour Kut ( = N^{t T} * N^{u})
  620. C    INTEGRATION POSSIBLEMENT INCOMPATIBLE POUR KUT !
  621. c on peut prendre le 1er enrichissement (ou le standard) car on a
  622. c toujours les memes noeuds au debut (sans tenir compte des B1X)
  623.  
  624. c              Stockage de l'element iem1 pour les differents PG
  625.                elemPG(igau2) = iem1
  626.  
  627.                DO 4205 JJ=1,nbnn1
  628.                  iino1=IPT1.NUM(JJ,iem1)
  629.  
  630.                  do JJ1 = 1,nbnno2
  631.                    if(PRELEM(JJ1).eq.iino1) goto 4207
  632.                  enddo
  633.                  write(IOIMP,*) 'On n a pas trouvé le noeud ',iino1
  634.                  call ERREUR(5)
  635.  4207            continue
  636.                  itoto = 0
  637.  
  638. c                UX
  639.                  DO II=1,nbnn2
  640.                      RINT0(II,JJ1) = RINT0(II,JJ1) + REL1(II,JJ)
  641.                  ENDDO
  642. c                AX
  643.                  DO II=1,nbnn2
  644.                    MELVA1 = MCHAM1.IELVAL(1)
  645.                    MLREEL = MELVA1.IELCHE(JJ,iem1)
  646.                    if(MLREEL.ne.0) then
  647.                      RINT1(II,JJ1) = RINT1(II,JJ1) + REL1(II,JJ)
  648. c                     write(ioimp,*) ' RINT1(',II,JJ1,') = ', RINT1(II,JJ1)
  649.                    endif
  650.                  ENDDO
  651.  
  652.  
  653. c                B1X
  654.                  IF(nbenrj.ge.2) THEN
  655.                    MELVA1 = MCHAM1.IELVAL(2)
  656.                    MLREEL = MELVA1.IELCHE(JJ,iem1)
  657.                    if(MLREEL.ne.0) then
  658.                      PSIX = 0.d0
  659.                      SEGACT,MLREEL
  660.                      do iii0 = 1,nbnn1
  661.                        PSIX  = PSIX + (SHPP1(1,iii0) * PROG(iii0))
  662.                      enddo
  663.                      RX05= SQRT(ABS(PSIX))
  664. c                      write(ioimp,*) 'RX05 = ', RX05
  665.                      do II=1,nbnn2
  666.                        RINT2(II,JJ1) = RINT2(II,JJ1) + RX05*REL1(II,JJ)
  667.                      enddo
  668.                   endif
  669.                  ENDIF
  670.  
  671.  
  672. c                B2X
  673.                  IF(nbenrj.ge.3) THEN
  674.                    MELVA1 = MCHAM1.IELVAL(3)
  675.                    MLREEL = MELVA1.IELCHE(JJ,iem1)
  676.                    if(MLREEL.ne.0) then
  677.                      PSIX = 0.d0
  678.                      SEGACT,MLREEL
  679.                      do iii0 = 1,nbnn1
  680.                        PSIX  = PSIX + (SHPP1(1,iii0) * PROG(iii0))
  681.                      enddo
  682.                      RX05= SQRT(ABS(PSIX))
  683.                      do II=1,nbnn2
  684.                         RINT3(II,JJ1) = RINT3(II,JJ1) + RX05*REL1(II,JJ)
  685.                      enddo
  686.                   endif
  687.                  ENDIF
  688.  
  689.  4205          CONTINUE
  690.  
  691.                if(iimpi.ge.4) then
  692.                  DO II=1,nbnn2
  693.           write(ioimp,*) 'RINT1(',II,',:)=',(RINT1(II,iou),iou=1,nbnno2)
  694.                  ENDDO
  695.                endif
  696.  
  697.  
  698. c              on a trouvé le iem1 élément et on a fait le travail :
  699. c              on passe au point de gauss suivant
  700.                goto 1132
  701.  
  702.  1131       CONTINUE
  703. c---------- fin de la Boucle sur les elements de structure ----------------------
  704.             if(iimpi.ge.2) write(ioimp,*)
  705.      1      'Fin de la boucle sur les elements de structure'
  706.             if(izok.eq.0) write(ioimp,*) 'Attention le pt de Gauss '
  707.      &      ,igau2,'de l element ',iem2,' est hors support'
  708.  
  709.  
  710.  1132    CONTINUE
  711. c======= fin de la Boucle sur les pt de G de fissure ============================
  712.  
  713. c    --- Remplissage du MELEME et DU XMATRI ---
  714.          DO ityp1=1,2
  715.  
  716. c         on commence par le MELEME des ddls Obligatoire (U)
  717.           izone = 1
  718. c         on poursuit avec le MELEME des ddls Sauts (H ou HB1 ou HB1B2)
  719.           if(ityp1.eq.2) izone=nbenrj+1
  720.           ity = ITAMYT(nbnno2,izone)
  721.  
  722.           do iidim=1,idim
  723.           ityty = ((ity-1) * idim) + iidim
  724.           if(iimpi.ge.4) write(ioimp,*) ity,nbnno2,izone,' -> ',ityty
  725.  
  726. c     --- Remplissage du MELEME ---
  727.           MELEME = IRIGEL(1,ityty)
  728.           NBNN   = NUM(/1)
  729.           NBELEM = NUM(/2) + 1
  730.           if(iimpi.ge.4) then
  731.           write(ioimp,*) 'on remplit le type de matrice ',ityty,
  732.      1                   ' NBNN,NBELEM=',NBNN,NBELEM
  733.           endif
  734.           segadj,MELEME
  735. c         element de fissure
  736.           do ino2=1,nbnn2
  737.             NUM(ino2,NBELEM) = IPT2.NUM(ino2,iem2)
  738. c             NUM(ino2+nbnn2,NBELEM) = IPT2.NUM(ino2,iem2)
  739.           enddo
  740. c         element de structure
  741.           do iino2=1,nbnno2
  742. c             NUM((2*nbnn2)+iino2,NBELEM) = PRELEM(iino2)
  743.             NUM(nbnn2+iino2,NBELEM) = PRELEM(iino2)
  744.           enddo
  745.  
  746. c     --- Remplissage du XMATRI.RE ---
  747.           XMATRI = IRIGEL(4,ityty)
  748.           NLIGRD = RE(/1)
  749.           NLIGRP = RE(/2)
  750.           NELRIG = RE(/3) + 1
  751.           segadj,XMATRI
  752.  
  753.           DO II=1,nbnn2
  754.  
  755. c            -Ktt
  756. c             on choisit de ne garder que la partie diagonale
  757.              JJ=II
  758.               XMATRI.RE(II,JJ,NBELEM) = RINT(II,JJ)*XISTAB
  759.  
  760. c            -Kwt
  761.               DO JJ=1,nbnn2
  762.                  XMATRI.RE(II,JJ+nbnn2,NBELEM) = RINT(II,JJ)
  763.                  XMATRI.RE(JJ+nbnn2,II,NBELEM) = RINT(II,JJ)
  764. c          write(ioimp,*) 'XMATRI.RE(',II,JJ+nbnn2,NBELEM,') =', RINT(II,JJ)
  765.               ENDDO
  766.  
  767. C            -Kww
  768. c             on choisit de ne garder que la partie diagonale
  769.              JJ=II
  770.                XMATRI.RE(nbnn2+II,nbnn2+JJ,NBELEM) = RINT(II,JJ)*RLATIN
  771.  
  772. C            -Kut (on a deja fait le travail de positionnement avant)
  773.               if(ityp1.eq.1) then
  774. c               UX
  775.                 JJdec = 2*nbnn2
  776.                 DO JJ=1,nbnno2
  777.                   XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT0(II,JJ)
  778.                   XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT0(II,JJ)
  779. c                 write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA
  780. c     1TRI.RE(II,JJdec+JJ,NBELEM)
  781.                 ENDDO
  782.               else
  783. c               AX
  784.                 JJdec = 2*nbnn2
  785.                 DO JJ=1,nbnno2
  786.                   XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT1(II,JJ)
  787.                   XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT1(II,JJ)
  788. c                 write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA
  789. c     1TRI.RE(II,JJdec+JJ,NBELEM)
  790.                 ENDDO
  791.  
  792. c               B1X
  793.                 IF(nbenrj.ge.2) THEN
  794.                    JJdec=JJdec+nbnno2
  795.                    DO JJ=1,nbnno2
  796.                       XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT2(II,JJ)
  797.                       XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT2(II,JJ)
  798.  
  799. c             Cas ou les inconnues aux noeuds ne sont pas présentes
  800. c             mis à 0 de ces inconnues en mettant un 1 sur la diagonale
  801.  
  802. c              Début de la boucle pour tester les noeuds de structure d'ancien point de gauss
  803. cbp                DO 5002 iigau2 = 1,igau2
  804.                 DO 5002 iigau2 = 1,ngau2
  805.                    iiem1 = elemPG(iigau2)
  806.                    DO JJ1=1,nbnn1
  807.                       iino1=IPT1.NUM(JJ1,iiem1)
  808.                       if(PRELEM(JJ).eq.iino1) then
  809.                          MELVA1 = MCHAM1.IELVAL(2)
  810.                          MLREEL = MELVA1.IELCHE(JJ1,iiem1)
  811.                          if(MLREEL.ne.0) goto 5002
  812.                          XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM) = 1
  813. c                         write(ioimp,*) 'XMATRI.RE(',JJd ec+JJ,JJdec+JJ,NBEL
  814. c     1EM,') =', XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM)
  815.                       endif
  816.                    ENDDO
  817. 5002             CONTINUE
  818.              ENDDO
  819.  
  820.  
  821.  
  822.                 ENDIF
  823. c               B2X
  824.                 IF(nbenrj.ge.3) THEN
  825.                 JJdec=JJdec+nbnno2
  826.                 DO JJ=1,nbnno2
  827.                   XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT3(II,JJ)
  828.                   XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT3(II,JJ)
  829.                 ENDDO
  830.                 ENDIF
  831.               endif
  832.  
  833.           ENDDO
  834.  
  835.        enddo
  836.  
  837.  
  838.        ENDDO
  839.  
  840.  
  841.  
  842.  1100 CONTINUE
  843. c==== fin de la Boucle sur les elements de fissure ==============================
  844. c----------------------------------------------------------------------
  845.  
  846.  
  847. c------------------------------------
  848. c     AJUSTEMENT AVANT DE QUITTER
  849. c------------------------------------
  850.  
  851.       if(iimpi.ge.2) write(ioimp,*) 'AJUSTEMENT AVANT DE QUITTER'
  852.  
  853.  
  854. C     BOUCLE SUR LES SOUS RIGIDITÉS
  855.       ityok = 0
  856.       DO 2000 ityty=1,(idim*NTYMAT)
  857.  
  858.          MELEME = IRIGEL(1,ityty)
  859.          DESCR  = IRIGEL(3,ityty)
  860.          XMATRI = IRIGEL(4,ityty)
  861.          NBELEM = NUM(/2)
  862.          if(NBELEM.ne.0) then
  863.            ityok = ityok + 1
  864.            IRIGEL(1,ityok) = MELEME
  865.            IRIGEL(3,ityok) = DESCR
  866.            IRIGEL(4,ityok) = XMATRI
  867.            segdes,XMATRI
  868.          else
  869.            segsup,MELEME,DESCR,XMATRI
  870.          endif
  871.  2000 continue
  872.       NRIGEL = ityok
  873.       segadj,MRIGID
  874.  
  875. c-------------------------------
  876. c     MENAGE AVANT DE QUITTER
  877. c-------------------------------
  878.  
  879.       segsup,MTRAV
  880.  
  881. c-------------------------------
  882. c     ON REMPLACE LES FX PAR DES LX
  883. c-------------------------------
  884.       IF(INOMU.eq.0) THEN
  885.         JGN=4
  886.         JGM=2*idim
  887.         segini,MLMOT1
  888.         do iidim=1,idim
  889.            MLMOT1.MOTS(iidim) = DUAOBL(iidim)
  890.            MLMOT1.MOTS(iidim+idim) = DUAFAC(iidim)
  891.         enddo
  892.         ILMOT1 = MLMOT1
  893.         IPRIG1 = MRIGID
  894.         if(iimpi.ge.2) write(ioimp,*)'       APPEL A FX2LX       '
  895.         CALL FX2LX(IPRIG1,ILMOT1,IPRIG2)
  896.         MRIGID = IPRIG2
  897.         SEGSUP,MLMOT1
  898.       ENDIF
  899.  
  900. c       write(ioimp,*) 'desactivation du MRIGID',MRIGID
  901.       SEGDES,MRIGID
  902.  
  903.       if(iimpi.ge.3) write(ioimp,*) 'ecriture du MRIGID',MRIGID
  904.       CALL ECROBJ('RIGIDITE',MRIGID)
  905.  
  906.       END
  907.  
  908.  
  909.  
  910.  

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