Télécharger la plaquette Cast3M


fr en

Télécharger relami.eso

Retour à la liste

Numérotation des lignes : 

relami

  1. C RELAMI    SOURCE    PV090527  26/04/30    21:16:09     12529          
  2.  
  3.       SUBROUTINE RELAMI(IP0,IPT1,IPRIG)
  4.  
  5. C=======================================================================
  6. C CE SOUS-PROGRAMME CONSTRUIT LA RIGIDITE LIANT LINEAIREMENT LES DDL
  7. C DES NOEUDS MILIEUX D'UN MAILLAGE QUADRATIQUE AUX NOEUDS SOMMETS
  8. C
  9. C SYNTHAXE GIBIANE : RIG1 = RELA MILI (LMOTS1) GEO1 ;
  10. C=======================================================================
  11.  
  12. C======================== ZONE DE DECLARATIONS =========================
  13.  
  14.       IMPLICIT INTEGER(I-N)
  15.       IMPLICIT REAL*8 (A-H,O-Z)
  16.  
  17. -INC PPARAM
  18. -INC CCOPTIO
  19. -INC CCHAMP
  20. -INC CCGEOME
  21. -INC SMLMOTS
  22. -INC SMCOORD
  23. -INC SMELEME
  24. -INC SMRIGID
  25. -INC CCREEL 
  26.  
  27.       CHARACTER*4 LESDDL(10),LESDUA(10)
  28.  
  29.       SEGMENT IMILI(NBDDL)
  30.  
  31. C========================= CORPS DU PROGRAMME ==========================
  32.  
  33. C==== NOMS DES INCONNUES DE LA RIGIDITE
  34.  
  35. C     Si pas de LISTMOTS, je prends le inconnues de la meca.
  36.       IF (IP0.EQ.0) THEN
  37.         IF (IFOUR.LT.0) THEN
  38.           NBDDL=2
  39.           LESDDL(1)='UX  '
  40.           LESDDL(2)='UY  '
  41.           LESDUA(1)='FX  '
  42.           LESDUA(2)='FY  '
  43.         ELSEIF (IFOUR.EQ.0) THEN
  44.           NBDDL=2
  45.           LESDDL(1)='UR  '
  46.           LESDDL(2)='UZ  '
  47.           LESDUA(1)='FR  '
  48.           LESDUA(2)='FZ  '
  49.         ELSEIF (IFOUR.EQ.1) THEN
  50.           NBDDL=3
  51.           LESDDL(1)='UR  '
  52.           LESDDL(2)='UT  '
  53.           LESDDL(3)='UZ  '
  54.           LESDUA(1)='FR  '
  55.           LESDUA(2)='FT  '
  56.           LESDUA(3)='FZ  '
  57.         ELSE
  58.           NBDDL=3
  59.           LESDDL(1)='UX  '
  60.           LESDDL(2)='UY  '
  61.           LESDDL(3)='UZ  '
  62.           LESDUA(1)='FX  '
  63.           LESDUA(2)='FY  '
  64.           LESDUA(3)='FZ  '
  65.         ENDIF
  66.       ELSE
  67. C     Sinon, on prends les DDL specifies et on cherche les duals
  68.         MLMOTS=IP0
  69.         SEGACT,MLMOTS
  70.         NBDDL=MOTS(/2)
  71.         DO I=1,NBDDL
  72.            DO J=1,LNOMDD
  73.              IF (MOTS(I).EQ.NOMDD(J)) THEN
  74.                LESDDL(I)=NOMDD(J)
  75.                LESDUA(I)=NOMDU(J)
  76.              ENDIF
  77.            ENDDO
  78.         ENDDO
  79.         SEGDES,MLMOTS
  80.       ENDIF
  81.  
  82. C==== TRANSFORMATION DU MAILLAGE INI. EN SEGMENTS SI BESOIN
  83.  
  84.       IF (IDIM.GE.2) THEN
  85.         CALL ECROBJ('MAILLAGE',IPT1)
  86.         CALL CHANLG
  87.         IF (IPT1.EQ.0) THEN
  88.           CALL ERREUR(16)
  89.         END IF
  90.         CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
  91.         IF (IERR.NE.0) RETURN
  92.       ENDIF
  93.  
  94. C==== CONSTRUCTION DU MAILLAGE SUPPORT DE LA RIGIDITE
  95.  
  96. C     J'initialise un vecteur que je vais remplir de maillages elem.
  97.       SEGINI,IMILI
  98.  
  99.       IDIMP1=IDIM+1
  100. C     Je parcours le maillage ini. et construis les maillages elem.
  101.       SEGACT,IPT1
  102.       NBSOUS1=IPT1.LISOUS(/1)
  103. C     J'ai un maillage simple
  104.       IF (NBSOUS1.EQ.0) THEN
  105.         IF (IPT1.ITYPEL.EQ.3) THEN
  106.           NBEL1=IPT1.NUM(/2)
  107.           segact mcoord*mod
  108.           NBPTI=nbpts
  109.           NBPTS=NBPTI+NBDDL*NBEL1
  110.           SEGADJ,MCOORD
  111.           NBSOUS=0
  112.           NBREF=0
  113.           NBELEM=NBEL1
  114.           NBNN=4
  115.           ICPT1=0
  116.           DO I=1,NBDDL
  117.             SEGINI,IPT2
  118.             IPT2.ITYPEL=22
  119.             DO J=1,NBEL1
  120.               IP=NBPTI+ICPT1+J
  121.               IREF=(IP-1)*IDIMP1
  122.               IREF2=(IPT1.NUM(2,J)-1)*IDIMP1
  123.               DO K=1,IDIMP1
  124.                 XCOOR(IREF+K)=XCOOR(IREF2+K)
  125.               ENDDO
  126.               IPT2.NUM(1,J)=IP
  127.               IPT2.NUM(2,J)=IPT1.NUM(2,J)
  128.               IPT2.NUM(3,J)=IPT1.NUM(1,J)
  129.               IPT2.NUM(4,J)=IPT1.NUM(3,J)
  130.             ENDDO
  131.             SEGDES,IPT2
  132.             IMILI(I)=IPT2
  133.             ICPT1=ICPT1+NBEL1
  134.           ENDDO
  135.         ELSE
  136. C         Si pas de SEG3, ERREUR
  137.           SEGSUP,IMILI
  138.           SEGDES,IPT1
  139.           CALL ERREUR(16)
  140.           RETURN
  141.         ENDIF
  142. C     J'ai un maillage complexe
  143.       ELSE
  144.         NBEL1=0
  145.         DO ISOUS=1,NBSOUS1
  146.           IPT3=IPT1.LISOUS(ISOUS)
  147.           SEGACT,IPT3
  148.           IF (IPT3.ITYPEL.EQ.3) THEN
  149.             NBEL1=NBEL1+IPT3.NUM(/2)
  150.           ENDIF
  151.         ENDDO
  152. C       Si pas de SEG3, ERREUR
  153.         IF (NBEL1.EQ.0) THEN
  154.           SEGSUP,IMILI
  155.           SEGDES,IPT1
  156.           CALL ERREUR(16)
  157.           RETURN
  158.         ENDIF
  159.         segact mcoord*mod
  160.         NBPTI=nbpts
  161.         NBPTS=NBPTI+NBDDL*NBEL1
  162.         SEGADJ,MCOORD
  163.         NBSOUS=0
  164.         NBREF=0
  165.         NBELEM=NBEL1
  166.         NBNN=4
  167.         ICPT1=0
  168.         DO I=1,NBDDL
  169.           ICEL1=0
  170.           SEGINI,IPT2
  171.           IPT2.ITYPEL=22
  172.           DO ISOUS=1,NBSOUS1
  173.             IPT3=IPT1.LISOUS(ISOUS)
  174.             IF (IPT3.ITYPEL.EQ.3) THEN
  175.               NBEL3=IPT3.NUM(/2)
  176.               DO J=1,NBEL3
  177.                 IP=NBPTI+ICPT1+J
  178.                 IREF=(IP-1)*IDIMP1
  179.                 IREF2=(IPT3.NUM(2,J)-1)*IDIMP1
  180.                 DO K=1,IDIMP1
  181.                   XCOOR(IREF+K)=XCOOR(IREF2+K)
  182.                 ENDDO
  183.                 IPT2.NUM(1,ICEL1+J)=IP
  184.                 IPT2.NUM(2,ICEL1+J)=IPT3.NUM(2,J)
  185.                 IPT2.NUM(3,ICEL1+J)=IPT3.NUM(1,J)
  186.                 IPT2.NUM(4,ICEL1+J)=IPT3.NUM(3,J)
  187.               ENDDO
  188.               ICPT1=ICPT1+NBEL3
  189.               ICEL1=ICEL1+NBEL3
  190.             ENDIF
  191.             SEGDES,IPT3
  192.           ENDDO
  193.           SEGDES,IPT2
  194.           IMILI(I)=IPT2
  195.         ENDDO
  196.       ENDIF
  197.       SEGDES,IPT1
  198.  
  199. C==== CONSTRUCTION DE LA RIGIDITE ASSOCIEE AUX RELA.
  200.  
  201.       NRIGEL=NBDDL
  202.       NRIGE=8
  203.       SEGINI,RI1
  204.       RI1.MTYMAT='RIGIDITE'
  205.       RI1.IFORIG=IFOUR
  206.       NLIGRP=4
  207.       NLIGRD=NLIGRP
  208.       DO I=1,NRIGEL
  209. C       On a un segment DESCR par DDL
  210.         SEGINI,DES1
  211.         DES1.LISINC(1)='LX'
  212.         DES1.LISDUA(1)='FLX'
  213.         DES1.NOELEP(1)=1
  214.         DES1.NOELED(1)=1
  215.         DO J=2,4
  216.           DES1.LISINC(J)=LESDDL(I)
  217.           DES1.LISDUA(J)=LESDUA(I)
  218.           DES1.NOELEP(J)=J
  219.           DES1.NOELED(J)=J
  220.         ENDDO
  221.         SEGDES,DES1
  222. C       On a un segment XMATRI par DDL
  223.         nelrig=nbel1
  224.         rigrel=0
  225.         SEGINI,XMATR1
  226.         XMATR1.RE(1,2,1)=1.D0
  227.         XMATR1.RE(1,3,1)=-0.5D0
  228.         XMATR1.RE(1,4,1)=-0.5D0
  229.         XMATR1.RE(2,1,1)=XMATR1.RE(1,2,1)
  230.         XMATR1.RE(3,1,1)=XMATR1.RE(1,3,1)
  231.         XMATR1.RE(4,1,1)=XMATR1.RE(1,4,1)
  232. *        SEGDES,XMATR1
  233. C       On a NBEL1 matrice(s) elementaire(s)
  234.         do ioup=2,nelrig
  235.            do io=1,xmatr1.re(/2)
  236.               do iu=1,xmatr1.re(/1)
  237.                  xmatr1.re(iu,io,ioup)=xmatr1.re(iu,io,1)
  238.                enddo
  239.            enddo
  240.         enddo
  241. *  reprendre xmatr1 pour ponderer en fonction des positions reelles des noeuds.
  242.           meleme=imili(i)
  243.          do ioup=1,nelrig
  244.           segact meleme
  245.           ipm=num(2,ioup)
  246.           ip1=num(3,ioup)
  247.           ip2=num(4,ioup)
  248.           xpm=xcoor((ipm-1)*(idim+1)+1)
  249.           ypm=xcoor((ipm-1)*(idim+1)+2)
  250.           if (idim.eq.3) then              
  251.           zpm=xcoor((ipm-1)*(idim+1)+3)
  252.           else
  253.           zpm=0.d0
  254.           endif
  255.           xp1=xcoor((ip1-1)*(idim+1)+1)
  256.           yp1=xcoor((ip1-1)*(idim+1)+2)
  257.           if (idim.eq.3) then              
  258.           zp1=xcoor((ip1-1)*(idim+1)+3)
  259.           else
  260.           zp1=0.d0
  261.           endif
  262.           xp2=xcoor((ip2-1)*(idim+1)+1)
  263.           yp2=xcoor((ip2-1)*(idim+1)+2)
  264.           if (idim.eq.3) then              
  265.           zp2=xcoor((ip2-1)*(idim+1)+3)
  266.           else
  267.           zp2=0.d0
  268.           endif
  269.           xl12=(xp1-xp2)*(xp1-xp2)+(yp1-yp2)*(yp1-yp2)+
  270.      >    (zp1-zp2)*(zp1-zp2)
  271.           xl2=max(xl2,xpetit)
  272.           xlm1=(xp1-xpm)*(xp1-xp2)+(yp1-ypm)*(yp1-yp2)+
  273.      >    (zp1-zpm)*(zp1-zp2)
  274.           xcoef=xlm1/xl12
  275.           xcoef=min(1.d0,max(0.d0,xcoef))
  276.         XMATR1.RE(1,3,ioup)=-1.d0+xcoef
  277.         XMATR1.RE(1,4,ioup)=-xcoef
  278.         XMATR1.RE(3,1,ioup)=XMATR1.RE(1,3,ioup)
  279.         XMATR1.RE(4,1,ioup)=XMATR1.RE(1,4,ioup)
  280.          enddo
  281.          segdes xmatr1
  282.  
  283. C       On remplit la rigidite
  284.         RI1.COERIG(I)=1.
  285.         RI1.IRIGEL(1,I)=IMILI(I)
  286.         RI1.IRIGEL(2,I)=0
  287.         RI1.IRIGEL(3,I)=DES1
  288.         RI1.IRIGEL(4,I)=xMATR1
  289.         RI1.IRIGEL(5,I)=NIFOUR
  290.         RI1.IRIGEL(6,I)=0
  291.         RI1.IRIGEL(7,I)=0
  292.         RI1.IRIGEL(8,I)=0
  293.       ENDDO
  294.       SEGDES,RI1
  295.       IPRIG=RI1
  296.       SEGSUP IMILI
  297.  
  298.       END
  299.  
  300.  
  301.  
  302.  

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