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    MB234859  25/09/08    21:15:01     12358          
  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.          SEGINI,XMATRI
  364.  
  365.          IRIGEL(1,ityty) = MELEME
  366.          IRIGEL(3,ityty) = DESCR
  367.          IRIGEL(4,ityty) = XMATRI
  368.          IRIGEL(5,ityty) = NIFOUR
  369.          IRIGEL(6,ityty) = 0
  370.          IRIGEL(7,ityty) = 0
  371.          IRIGEL(8,ityty) = 0
  372.  
  373.       enddo
  374.       enddo
  375.  
  376.  
  377. c----------------------------------------------------------------------
  378. c     1. RECHERCHE DES ELEMENTS DE STRUCTURE CONTENANT DES POINTS DE GAUSS
  379. c        DES ELEMENTS DE LA FISSURE
  380. c     2. REMPLISSAGE DU MELEME + MRIGID en FONCTION
  381. c----------------------------------------------------------------------
  382.  
  383. C
  384. c==== Boucle sur les elements de fissure ==============================
  385.       DO 1100 iem2=1,nbel2
  386.  
  387.          call doxe(xcoor,idim,nbnn2,ipt2.num,iem2,xe2)
  388.          nbenrj = 0
  389.          nbnno2 = 0
  390.          call ZERO(RINT,nbnn2,nbnn2)
  391.          call ZERO(RINT0,nbnn2,ngau2*nbnn1)
  392.          call ZERO(RINT1,nbnn2,ngau2*nbnn1)
  393.          call ZERO(RINT2,nbnn2,ngau2*nbnn1)
  394.          call ZERO(RINT3,nbnn2,ngau2*nbnn1)
  395.  
  396. c======= Boucle sur les pt de G de fissure ============================
  397.          DO 1132 igau2=1,ngau2
  398.             izok = 0
  399. c   récupération des coordonnees du point de gauss dans le repère global
  400.             XPO(1) = 0.D0
  401.             XPO(2) = 0.D0
  402.             XPO(3) = 0.D0
  403.             DO  ino=1,nbnn2
  404.                XPO(1) = XPO(1) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(1,ino))
  405.                XPO(2) = XPO(2) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(2,ino))
  406.                IF(IDIM.EQ.3) then
  407.                XPO(3) = XPO(3) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(3,ino))
  408.                ENDIF
  409.             ENDDO
  410.  
  411. c            if(iimpi.ge.5) write(ioimp,*) 'Fissure: element',iem2,
  412. c      &     ' point de Gauss',igau2,' x=',XPO(1),XPO(2)
  413.  
  414. c---------- Boucle sur les elements de structure ----------------------
  415.             DO 1131 iem1=1,nbel1
  416.  
  417. c   si pas d'enrichissement, on travaille sur tous les elements
  418.                if(MCHAM1.eq.0) goto 1133
  419. c   on saute les elements non enrichi car a priori ne contiennent pas la fissure
  420.                do ienr2=1,NBENR2
  421.                   MELVA1=MCHAM1.IELVAL(IENR2)
  422.                   if(MELVA1.ne.0) then
  423.                     do inode1=1,nbnn1
  424.                   if(MELVA1.IELCHE(inode1,iem1).ne.0) goto 1133
  425.                     enddo
  426.                   endif
  427.                enddo
  428.                goto 1131
  429.  1133          continue
  430.  
  431. c   recuperation des coordonnées des noeuds de IPT1 : xel1 (dans le repère x,y,z)
  432.                call doxe(xcoor,idim,nbnn1,ipt1.num,iem1,xel1)
  433.  
  434. c   calcul des fonctions de formes de IPT1 au pt de Gauss de IPT2
  435.                qsi(1) = 0.D0
  436.                qsi(2) = 0.D0
  437.                qsi(3) = 0.D0
  438.                call QSIJS(xel1,iele1,nbnn1,idim,XPO,SHPP1,qsi,IRET)
  439. c            if(iimpi.ge.5) write(ioimp,*) 'Structure: element',iem1,
  440. c      &     ' noeuds:',(ipt1.num(iou,iem1),iou=1,nbnn1)
  441. c            if(iimpi.ge.5) write(ioimp,*)
  442. c      &     '  N_i=',(SHPP1(1,iou),iou=1,nbnn1)
  443.  
  444.  
  445. c   test pour savoir si PG est dans EF de IPT1
  446.                DO 1130 ino1=1,NBNN1
  447.                   if (SHPP1(1,ino1).LT.-1.01D-7) then
  448.                      go to 1131
  449.                   endif
  450.  1130          continue
  451.                izok = 1
  452. c              ON a trouvé l'élément de structure !
  453. c           write(ioimp,*) 'ON a trouve l element de structure !'
  454. c     DETECTION DU TYPE D'ENRICHISSEMENT MAX DE CET ELEMENT = nbenrj
  455.                DO 3001 IENR2=1,NBENR2
  456.                  MELVA1=MCHAM1.IELVAL(IENR2)
  457.                  IF(MELVA1.eq.0) GOTO 3001
  458.                  DO 3002 ino1=1,nbnn1
  459.                    MLREEL = MELVA1.IELCHE(ino1,iem1)
  460. c                  Test pour savoir si le noeud est enrichi
  461.                    IF(MLREEL.eq.0) GOTO 3002
  462.                    nbenrj=max(nbenrj,IENR2)
  463.  3002            continue
  464.  3001          continue
  465.  
  466.                if(iimpi.ge.3) write(ioimp,*) 'EF fissure ',iem2,
  467.      &         ' ptdeG ',igau2,' -> EF MASSIF ',iem1,' nbenrj=',nbenrj
  468.  
  469. c        Préparation du Remplissage du MELEME
  470.                do 3003 ino1=1,nbnn1
  471.                   iino1=IPT1.NUM(ino1,iem1)
  472.                   if(nbnno2.gt.0) then
  473.                   do iino2=1,nbnno2
  474.                     if(PRELEM(iino2).eq.iino1) goto 3003
  475.                   enddo
  476.                   endif
  477. c                 nouveau noeud : on l'ajoute
  478.                   nbnno2=nbnno2+1
  479.                   PRELEM(nbnno2)=iino1
  480.  3003          continue
  481.                if(iimpi.ge.3) then
  482.                write(ioimp,*) ' PRELEM=',( PRELEM(iou),iou=1,nbnno2)
  483.                write(ioimp,*) 'EF fissure ',iem2,' ptdeG ',igau2
  484.                endif
  485.  
  486. C        CALCUL DE REL ET REL1
  487.  
  488. c matrices des fonctions de forme de l'element associé à IPT2
  489.  
  490.                do ino2=1,NBNN2
  491.                   xe3(1,ino2) = sqrt( xe2(1,ino2)*xe2(1,ino2)
  492.      1                              + xe2(2,ino2)*xe2(2,ino2) )
  493.                   NGENE2(1,ino2)=MINTE2.SHPTOT(1,ino2,IGAU2)
  494.                   do ipp2=1,6
  495.                      SHPP2(ipp2,ino2)=MINTE2.SHPTOT(ipp2,ino2,IGAU2)
  496.                   enddo
  497.                   if(iimpi.ge.4) then
  498.                   write(ioimp,*) 'SHPP2',(SHPP2(iou,ino2),iou=1,6)
  499.                   write(ioimp,*) 'xe2',(xe2(iou,ino2),iou=1,2)
  500.                   endif
  501.                enddo
  502. c              calcul du jacobien
  503.                if(idim.eq.2) then
  504. c              calcul du jacobien en 2D pour un seg2
  505.                   dXdQsi=0.D0
  506.                   dYdQsi=0.D0
  507.                   DO i=1,NBNN2
  508.                      dXdQsi=dXdQsi+SHPP2(2,i)*XE2(1,i)
  509.                      dYdQsi=dYdQsi+SHPP2(2,i)*XE2(2,i)
  510.                   ENDDO
  511. c           Dans le cas des seg2, poids des pts de gauss = 1
  512.                   DJAC2=SQRT(dXdQsi*dXdQsi+dYdQsi*dYdQsi)
  513.                 else
  514.                       dXdQsi=0.D0
  515.                       dYdQsi=0.D0
  516.                       dZdQsi=0.D0
  517.                       dXdEta=0.D0
  518.                       dYdEta=0.D0
  519.                       dZdEta=0.D0
  520.                       DO i=1,NBNN2
  521.                          dXdQsi=dXdQsi+SHPP2(2,i)*XE2(1,i)
  522.                          dXdEta=dXdEta+SHPP2(3,i)*XE2(1,i)
  523.                          dYdQsi=dYdQsi+SHPP2(2,i)*XE2(2,i)
  524.                          dYdEta=dYdEta+SHPP2(3,i)*XE2(2,i)
  525.                          dZdQsi=dZdQsi+SHPP2(2,i)*XE2(3,i)
  526.                          dZdEta=dZdEta+SHPP2(3,i)*XE2(3,i)
  527.                       ENDDO
  528.                       z = (dXdQsi*dYdEta-dXdEta*dYdQsi)
  529.                       x = (dYdQsi*dZdEta-dYdEta*dZdQsi)
  530.                       y = (dZdQsi*dXdEta-dZdEta*dXdQsi)
  531. c           Dans le cas des tri3, poids du pts de gauss = 0.5
  532.                       DJAC2 = sqrt(x*x+y*y+z*z)*0.5
  533.                endif
  534.  
  535. c matrices des fonctions de forme de l'element associé à IPT1
  536.                do ino1=1,NBNN1
  537.                   NGENE1(1,ino1)=SHPP1(1,ino1)
  538.                enddo
  539.                if(iimpi.ge.4) then
  540.                  write(ioimp,*) 'DJAC2=',DJAC2
  541.                  write(ioimp,*) 'NGENE1=',(NGENE1(1,iou),iou=1,NBNN1)
  542.                  write(ioimp,*) 'NGENE2=',(NGENE2(1,iou),iou=1,NBNN2)
  543.                endif
  544.  
  545. c  calcul des matrices elementaires = produits des fonctions de forme
  546.                DO II=1,nbnn2
  547.                   XGENE2 = NGENE2(1,II)*DJAC2
  548. c calcul de K_ww, Kwt et Ktt
  549.                    DO JJ=1,nbnn2
  550. C           CC = Jacobien * transposee(N)(II,i) * N(i,JJ) (somme sur i)
  551.                       CC = XGENE2*NGENE2(1,jj)
  552. C           REL est une matrice symetrique
  553.                       REL(II,JJ) = CC
  554.                       REL(JJ,II) = CC
  555.                    ENDDO
  556. c calcul de K_tu
  557.                    DO JJ=1,nbnn1
  558. C           CC = Jacobien * transposee(N)(II,i) * N(i,JJ) (somme sur i)
  559.                       CC=XGENE2*NGENE1(1,JJ)
  560. c                     write(ioimp,*) ' CC = ', CC
  561.                       REL1(II,JJ)=CC
  562.                    ENDDO
  563.                   if(iimpi.ge.4) write(ioimp,*)
  564.      1            'REL1(',II,',:)=',(REL1(II,iou),iou=1,NBNN1)
  565.                ENDDO
  566.  
  567. C         CUMUL DANS RINT ET RINT1
  568. C     CUMUL DANS RINT(II,JJ) pour Kww, Kwt, Ktw et Ktt
  569. c     (4 fois la meme matrices = N^{w T} * N^{w} = N^{w T} * N^{t} = ...)
  570.                DO  JJ=1,nbnn2
  571.                    DO  II=1,nbnn2
  572.                          RINT(II,JJ) = RINT(II,JJ) + REL(II,JJ)
  573.                    ENDDO
  574.                ENDDO
  575.  
  576. C     CUMUL DANS RINT1(II,JJ) pour Kut ( = N^{t T} * N^{u})
  577. C    INTEGRATION POSSIBLEMENT INCOMPATIBLE POUR KUT !
  578. c on peut prendre le 1er enrichissement (ou le standard) car on a
  579. c toujours les memes noeuds au debut (sans tenir compte des B1X)
  580.  
  581. c              Stockage de l'element iem1 pour les differents PG
  582.                elemPG(igau2) = iem1
  583.  
  584.                DO 4205 JJ=1,nbnn1
  585.                  iino1=IPT1.NUM(JJ,iem1)
  586.  
  587.                  do JJ1 = 1,nbnno2
  588.                    if(PRELEM(JJ1).eq.iino1) goto 4207
  589.                  enddo
  590.                  write(IOIMP,*) 'On n a pas trouvé le noeud ',iino1
  591.                  call ERREUR(5)
  592.  4207            continue
  593.                  itoto = 0
  594.  
  595. c                UX
  596.                  DO II=1,nbnn2
  597.                      RINT0(II,JJ1) = RINT0(II,JJ1) + REL1(II,JJ)
  598.                  ENDDO
  599. c                AX
  600.                  MELVA1 = MCHAM1.IELVAL(1)
  601.                  MLREEL = MELVA1.IELCHE(JJ,iem1)
  602.                  if(MLREEL.ne.0) then
  603.                    DO II=1,nbnn2
  604.                      RINT1(II,JJ1) = RINT1(II,JJ1) + REL1(II,JJ)
  605. c                     write(ioimp,*) ' RINT1(',II,JJ1,') = ', RINT1(II,JJ1)
  606.                    ENDDO
  607.                  endif
  608.  
  609.  
  610. c                B1X
  611.                  IF(nbenrj.ge.2) THEN
  612.                    MELVA1 = MCHAM1.IELVAL(2)
  613.                    MLREEL = MELVA1.IELCHE(JJ,iem1)
  614.                    if(MLREEL.ne.0) then
  615.                      SEGACT,MLREEL
  616.                      PSIX = 0.d0
  617.                      do iii0 = 1,nbnn1
  618.                        PSIX  = PSIX + (SHPP1(1,iii0) * PROG(iii0))
  619.                      enddo
  620.                      RX05= SQRT(ABS(PSIX))
  621. c                      write(ioimp,*) 'RX05 = ', RX05
  622.                      do II=1,nbnn2
  623.                        RINT2(II,JJ1) = RINT2(II,JJ1) + RX05*REL1(II,JJ)
  624.                      enddo
  625.                   endif
  626.                  ENDIF
  627.  
  628.  
  629. c                B2X
  630.                  IF(nbenrj.ge.3) THEN
  631.                    MELVA1 = MCHAM1.IELVAL(3)
  632.                    MLREEL = MELVA1.IELCHE(JJ,iem1)
  633.                    if(MLREEL.ne.0) then
  634.                      SEGACT,MLREEL
  635.                      PSIX = 0.d0
  636.                      do iii0 = 1,nbnn1
  637.                        PSIX  = PSIX + (SHPP1(1,iii0) * PROG(iii0))
  638.                      enddo
  639.                      RX05= SQRT(ABS(PSIX))
  640.                      do II=1,nbnn2
  641.                         RINT3(II,JJ1) = RINT3(II,JJ1) + RX05*REL1(II,JJ)
  642.                      enddo
  643.                   endif
  644.                  ENDIF
  645.  
  646.  4205          CONTINUE
  647.  
  648.                if(iimpi.ge.4) then
  649.                  DO II=1,nbnn2
  650.           write(ioimp,*) 'RINT1(',II,',:)=',(RINT1(II,iou),iou=1,nbnno2)
  651.                  ENDDO
  652.                endif
  653.  
  654.  
  655. c              on a trouvé le iem1 élément et on a fait le travail :
  656. c              on passe au point de gauss suivant
  657.                goto 1132
  658.  
  659.  1131       CONTINUE
  660. c---------- fin de la Boucle sur les elements de structure ----------------------
  661.             if(iimpi.ge.2) write(ioimp,*)
  662.      1      'Fin de la boucle sur les elements de structure'
  663.             if(izok.eq.0) write(ioimp,*) 'Attention le pt de Gauss '
  664.      &      ,igau2,'de l element ',iem2,' est hors support'
  665.  
  666.  
  667.  1132    CONTINUE
  668. c======= fin de la Boucle sur les pt de G de fissure ============================
  669.  
  670. c    --- Remplissage du MELEME et DU XMATRI ---
  671.          DO ityp1=1,2
  672.  
  673. c         on commence par le MELEME des ddls Obligatoire (U)
  674.           izone = 1
  675. c         on poursuit avec le MELEME des ddls Sauts (H ou HB1 ou HB1B2)
  676.           if(ityp1.eq.2) izone=nbenrj+1
  677.           ity = ITAMYT(nbnno2,izone)
  678.  
  679.           do iidim=1,idim
  680.           ityty = ((ity-1) * idim) + iidim
  681.           if(iimpi.ge.4) write(ioimp,*) ity,nbnno2,izone,' -> ',ityty
  682.  
  683. c     --- Remplissage du MELEME ---
  684.           MELEME = IRIGEL(1,ityty)
  685.           NBNN   = NUM(/1)
  686.           NBELEM = NUM(/2) + 1
  687.           if(iimpi.ge.4) then
  688.           write(ioimp,*) 'on remplit le type de matrice ',ityty,
  689.      1                   ' NBNN,NBELEM=',NBNN,NBELEM
  690.           endif
  691.           segadj,MELEME
  692. c         element de fissure
  693.           do ino2=1,nbnn2
  694.             NUM(ino2,NBELEM) = IPT2.NUM(ino2,iem2)
  695. c             NUM(ino2+nbnn2,NBELEM) = IPT2.NUM(ino2,iem2)
  696.           enddo
  697. c         element de structure
  698.           do iino2=1,nbnno2
  699. c             NUM((2*nbnn2)+iino2,NBELEM) = PRELEM(iino2)
  700.             NUM(nbnn2+iino2,NBELEM) = PRELEM(iino2)
  701.           enddo
  702.  
  703. c     --- Remplissage du XMATRI.RE ---
  704.           XMATRI = IRIGEL(4,ityty)
  705.           NLIGRD = RE(/1)
  706.           NLIGRP = RE(/2)
  707.           NELRIG = RE(/3) + 1
  708.           segadj,XMATRI
  709.  
  710.           DO II=1,nbnn2
  711.  
  712. c            -Ktt
  713. c             on choisit de ne garder que la partie diagonale
  714.              JJ=II
  715.               XMATRI.RE(II,JJ,NBELEM) = RINT(II,JJ)*XISTAB
  716.  
  717. c            -Kwt
  718.               DO JJ=1,nbnn2
  719.                  XMATRI.RE(II,JJ+nbnn2,NBELEM) = RINT(II,JJ)
  720.                  XMATRI.RE(JJ+nbnn2,II,NBELEM) = RINT(II,JJ)
  721. c          write(ioimp,*) 'XMATRI.RE(',II,JJ+nbnn2,NBELEM,') =', RINT(II,JJ)
  722.               ENDDO
  723.  
  724. C            -Kww
  725. c             on choisit de ne garder que la partie diagonale
  726.              JJ=II
  727.                XMATRI.RE(nbnn2+II,nbnn2+JJ,NBELEM) = RINT(II,JJ)*RLATIN
  728.  
  729. C            -Kut (on a deja fait le travail de positionnement avant)
  730.               if(ityp1.eq.1) then
  731. c               UX
  732.                 JJdec = 2*nbnn2
  733.                 DO JJ=1,nbnno2
  734.                   XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT0(II,JJ)
  735.                   XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT0(II,JJ)
  736. c                 write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA
  737. c     1TRI.RE(II,JJdec+JJ,NBELEM)
  738.                 ENDDO
  739.               else
  740. c               AX
  741.                 JJdec = 2*nbnn2
  742.                 DO JJ=1,nbnno2
  743.                   XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT1(II,JJ)
  744.                   XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT1(II,JJ)
  745. c                 write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA
  746. c     1TRI.RE(II,JJdec+JJ,NBELEM)
  747.                 ENDDO
  748.  
  749. c               B1X
  750.                 IF(nbenrj.ge.2) THEN
  751.                    JJdec=JJdec+nbnno2
  752.                    DO JJ=1,nbnno2
  753.                       XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT2(II,JJ)
  754.                       XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT2(II,JJ)
  755.  
  756. c             Cas ou les inconnues aux noeuds ne sont pas présentes
  757. c             mis à 0 de ces inconnues en mettant un 1 sur la diagonale
  758.  
  759. c              Début de la boucle pour tester les noeuds de structure d'ancien point de gauss
  760. cbp                DO 5002 iigau2 = 1,igau2
  761.                 DO 5002 iigau2 = 1,ngau2
  762.                    iiem1 = elemPG(iigau2)
  763.                    DO JJ1=1,nbnn1
  764.                       iino1=IPT1.NUM(JJ1,iiem1)
  765.                       if(PRELEM(JJ).eq.iino1) then
  766.                          MELVA1 = MCHAM1.IELVAL(2)
  767.                          MLREEL = MELVA1.IELCHE(JJ1,iiem1)
  768.                          if(MLREEL.ne.0) goto 5002
  769.                          XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM) = 1
  770. c                         write(ioimp,*) 'XMATRI.RE(',JJd ec+JJ,JJdec+JJ,NBEL
  771. c     1EM,') =', XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM)
  772.                       endif
  773.                    ENDDO
  774. 5002             CONTINUE
  775.              ENDDO
  776.  
  777.                 ENDIF
  778. c               B2X
  779.                 IF(nbenrj.ge.3) THEN
  780.                 JJdec=JJdec+nbnno2
  781.                 DO JJ=1,nbnno2
  782.                   XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT3(II,JJ)
  783.                   XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT3(II,JJ)
  784.                 ENDDO
  785.                 ENDIF
  786.               endif
  787.  
  788.           ENDDO
  789.  
  790.        enddo
  791.  
  792.        ENDDO
  793.  
  794.  1100 CONTINUE
  795. c==== fin de la Boucle sur les elements de fissure ==============================
  796. c----------------------------------------------------------------------
  797.  
  798.  
  799. c------------------------------------
  800. c     AJUSTEMENT AVANT DE QUITTER
  801. c------------------------------------
  802.  
  803.       if(iimpi.ge.2) write(ioimp,*) 'AJUSTEMENT AVANT DE QUITTER'
  804.  
  805.  
  806. C     BOUCLE SUR LES SOUS RIGIDITÉS
  807.       ityok = 0
  808.       DO 2000 ityty=1,(idim*NTYMAT)
  809.  
  810.          MELEME = IRIGEL(1,ityty)
  811.          DESCR  = IRIGEL(3,ityty)
  812.          XMATRI = IRIGEL(4,ityty)
  813.          NBELEM = NUM(/2)
  814.          if(NBELEM.ne.0) then
  815.            ityok = ityok + 1
  816.            IRIGEL(1,ityok) = MELEME
  817.            IRIGEL(3,ityok) = DESCR
  818.            IRIGEL(4,ityok) = XMATRI
  819.            segdes,XMATRI
  820.          else
  821.            segsup,MELEME,DESCR,XMATRI
  822.          endif
  823.  2000 continue
  824.       NRIGEL = ityok
  825.       segadj,MRIGID
  826.  
  827. c-------------------------------
  828. c     MENAGE AVANT DE QUITTER
  829. c-------------------------------
  830.  
  831.       segsup,MTRAV
  832.  
  833. c-------------------------------
  834. c     ON REMPLACE LES FX PAR DES LX
  835. c-------------------------------
  836.       IF(INOMU.eq.0) THEN
  837.         JGN=4
  838.         JGM=2*idim
  839.         segini,MLMOT1
  840.         do iidim=1,idim
  841.            MLMOT1.MOTS(iidim) = DUAOBL(iidim)
  842.            MLMOT1.MOTS(iidim+idim) = DUAFAC(iidim)
  843.         enddo
  844.         ILMOT1 = MLMOT1
  845.         IPRIG1 = MRIGID
  846.         if(iimpi.ge.2) write(ioimp,*)'       APPEL A FX2LX       '
  847.         CALL FX2LX(IPRIG1,ILMOT1,IPRIG2)
  848.         MRIGID = IPRIG2
  849.         SEGSUP,MLMOT1
  850.       ENDIF
  851.  
  852. c       write(ioimp,*) 'desactivation du MRIGID',MRIGID
  853.       SEGDES,MRIGID
  854.  
  855.       if(iimpi.ge.3) write(ioimp,*) 'ecriture du MRIGID',MRIGID
  856.       CALL ECROBJ('RIGIDITE',MRIGID)
  857.  
  858.       END
  859.  
  860.  
  861.  
  862.  

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