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

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