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

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