Télécharger la plaquette Cast3M


fr en

Télécharger relsur.eso

Retour à la liste

Numérotation des lignes : 

relsur

  1. C RELSUR    SOURCE    CB215821  24/04/12    21:17:07     11897          
  2. C RELSUR   SOURCE    BP208322  17/04/18    21:15:13     9395           
  3.       SUBROUTINE RELSUR(IMODEL, MRIGID) 
  4. C***********************************************************************
  5. C   cet operateur crée une matrice élémentaire de rigidité 
  6. c   pour imposer les relation portées par un modele de sure
  7. C                                                                     
  8. C   ENTREES :                                                         
  9. C   ________                                                          
  10. C                                                                     
  11. C           IMODEL   pointeur sur le modele élémentaire
  12. C                              
  13. C   ENTREES/SORTIE :                                                  
  14. C   ________                                                          
  15. C                                                                     
  16. c           MRIGID  rigidité chapeau dans laquelle on va écrire 
  17. c           (à la suite des sous matrices déja présentes) les rigidités
  18. c           voulues.                                          
  19. C***********************************************************************
  20.  
  21.       IMPLICIT REAL*8 (A-H,O-Z)
  22.       IMPLICIT INTEGER (I-N)
  23. C
  24.  
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. -INC SMRIGID
  28. C-INC SMINTE
  29. -INC SMMODEL
  30. -INC SMELEME
  31. -INC SMCHAML
  32. -INC CCHAMP
  33. -INC CCGEOME
  34. -INC SMCOORD
  35. C
  36.  
  37. C
  38. C Petit tableau des "couleurs" des relations de conformite
  39.       DIMENSION LCOLOR(6)
  40.       DATA LCOLOR / 1, 3, 6, 10, 16, 24 /
  41. C
  42. C nombre de sous matrice de mrigid (va être ammené a changé)
  43.       NRIGEL = MRIGID.IRIGEL(/2)
  44.  
  45.       IPT1=IMODEL.imamod
  46.       SEGACT, IPT1
  47.       IDEBUT = LCOLOR(IPT1.ICOLOR(1)) - 3
  48.  
  49. c récupérations du nom des composantes obligatoire du modele de sure
  50.       nomid=IMODEL.lnomid(1)
  51.       SEGACT, nomid
  52.       SEGACT,MCOORD*MOD
  53. C**********************************************************************
  54. C Boucle sur les composantes primales obligatoires du sure
  55. C**********************************************************************
  56.       ICOMP=0
  57.  
  58.       DO 2 ICOMP=1,nomid.lesobl(/2)
  59.  
  60. C====================
  61. c creation d'un maillage avec un premier noeud par lélément 
  62. c correspondant à un multiplicateur de lagrange
  63. C====================
  64.           NBNN=IPT1.NUM(/1)+1
  65.           NBELEM=IPT1.NUM(/2)
  66.           NBSOUS=0
  67.           NBREF=0
  68.           SEGINI, IPT4
  69.           IPT4.ITYPEL=22
  70.           DO 1 J=1, NBELEM
  71.               ipt4.icolor(j)=IPT1.icolor(j)
  72.               DO 11 I=1,IPT1.NUM(/1)
  73.                   IPT4.NUM(I+1,J)=IPT1.NUM(I,J)
  74.   11          CONTINUE
  75.   1       CONTINUE
  76. C creation n'un nouveau noeud pour supporter chaque multiplicateur de lagrange
  77.           NBPT1= nbpts
  78.           NBPTS=NBPT1+(IPT4.NUM(/2))
  79.           SEGADJ MCOORD
  80.           DO 12 J=1, NBELEM
  81.              NGLOB=(NBPT1+J-1)*(IDIM+1)
  82. C remplissage des coordonees des nouveux points
  83.              DO 13 ID= 1,IDIM
  84.                  XCOOR(NGLOB+ID)=XCOOR((IPT4.NUM(2,J)-1)*(IDIM+1)+ID)
  85.   13         CONTINUE
  86.              IPT4.NUM(1,J) = NBPT1 + J
  87.   12      CONTINUE
  88.  
  89. C====================
  90. C *** SEGMENT XMATRI
  91. C====================
  92.          NLIGRD=IPT4.NUM(/1)
  93.          NLIGRP=NLIGRD
  94.          NELRIG=IPT4.NUM(/2)
  95.          SEGINI, XMATRI
  96.          DO 3 I=1,NELRIG
  97.             RE(1,2,i)=-1.
  98.             RE(2,1,i)=-1.
  99.             DO 4 ICAZ=3,NLIGRD
  100.                RE(1,ICAZ,i)=XCOEFF(IDEBUT+ICAZ)
  101.                RE(ICAZ,1,i)=RE(1,ICAZ,i)
  102.   4        CONTINUE
  103.   3      CONTINUE
  104. C         write(*,*) 'COMPOSANTE OBLIGATOIRE DU SURE :'
  105. C         write(*,*) (nomid.lesobl(ICOMP))
  106. C         write(*,*) 'MATRICE elementaire :' 
  107. C         DO 5 I=1,NLIGRD
  108. C            write(*,*) (RE(i,iou,1), iou=1,NLIGRD)
  109. C  5      CONTINUE
  110.  
  111. C====================
  112. C *** SEGMENT DESCR
  113. C====================
  114.  
  115.          NEXIST=0
  116.          DO 6 I=1,LNOMDD
  117.            IF (NOMDD(I).EQ.nomid.lesobl(ICOMP)) NEXIST = I
  118.   6      CONTINUE
  119.  
  120.          IF (NEXIST.EQ.0) THEN
  121.            CALL ERREUR(837)
  122.          ENDIF
  123.  
  124.          SEGINI, DESCR
  125.          LISINC(1)='LX  '
  126.          LISDUA(1)='FLX '
  127.          NOELEP(1)=1
  128.          NOELED(1)=1
  129.          DO 7 I=2,NLIGRD
  130.            LISINC(I)=NOMDD(NEXIST)
  131.            LISDUA(I)=NOMDU(NEXIST)
  132.            NOELEP(i)=i
  133.            NOELED(i)=i
  134.   7      CONTINUE
  135.  
  136. C====================
  137. C *** stockage de la rigidité construite dans MRIGID
  138. C====================
  139.  
  140.  
  141.  
  142.          NRIGEL = NRIGEL+1
  143.          SEGADJ, MRIGID
  144.          SEGACT, MRIGID*mod
  145.  
  146.          COERIG(NRIGEL)=1.
  147.          IRIGEL(1,NRIGEL)=IPT4
  148.          IRIGEL(3,NRIGEL)=DESCR
  149.          IRIGEL(4,NRIGEL)=XMATRI
  150.  
  151.  
  152. c         SEGDES, MRIGID
  153.          SEGDES, DESCR
  154.          SEGDES, XMATRI
  155.          SEGDES, IPT4
  156. C**********************************************************************
  157. C FIN Boucle sur les composantes primales obligatoires du sure
  158. C**********************************************************************
  159.   2   CONTINUE
  160. c      write (*,*) 'nomid.lesfac(/2)', nomid.lesfac(/2)
  161.       IF (nomid.lesfac(/2).EQ.0) goto 100      
  162. c récupération du champ d'enrichissement
  163. c      MCHEX1= IMODEL.IVAMOD(1) 
  164. c      SEGACT, MCHEX1    
  165. c      MCHAM1= MCHEX1.ICHAML(1)
  166.       MCHAM1= IMODEL.IVAMOD(1)
  167.       SEGACT, MCHAM1
  168.  
  169.  
  170. C**********************************************************************
  171. C Boucle sur les composantes primales facultatives du sure
  172. C**********************************************************************
  173.       ICOMP=0
  174.       DO 101 ICOMP=1,nomid.lesfac(/2)
  175. C++++    choix du type d'enrichisement de la composante ICOMP
  176.          IF (nomid.lesfac(ICOMP).EQ.'AX') MELVA1=MCHAM1.IELVAL(1)
  177.          IF (nomid.lesfac(ICOMP).EQ.'AY') MELVA1=MCHAM1.IELVAL(1)
  178.          IF (nomid.lesfac(ICOMP).EQ.'AZ') MELVA1=MCHAM1.IELVAL(1)
  179.  
  180.          IF (nomid.lesfac(ICOMP).EQ.'B1X') MELVA1=MCHAM1.IELVAL(2)
  181.          IF (nomid.lesfac(ICOMP).EQ.'B1Y') MELVA1=MCHAM1.IELVAL(2)
  182.          IF (nomid.lesfac(ICOMP).EQ.'B1Z') MELVA1=MCHAM1.IELVAL(2)
  183.          IF (nomid.lesfac(ICOMP).EQ.'C1X') MELVA1=MCHAM1.IELVAL(2)
  184.          IF (nomid.lesfac(ICOMP).EQ.'C1Y') MELVA1=MCHAM1.IELVAL(2)
  185.          IF (nomid.lesfac(ICOMP).EQ.'C1Z') MELVA1=MCHAM1.IELVAL(2)
  186.          IF (nomid.lesfac(ICOMP).EQ.'D1X') MELVA1=MCHAM1.IELVAL(2)
  187.          IF (nomid.lesfac(ICOMP).EQ.'D1Y') MELVA1=MCHAM1.IELVAL(2)
  188.          IF (nomid.lesfac(ICOMP).EQ.'D1Z') MELVA1=MCHAM1.IELVAL(2)
  189.          IF (nomid.lesfac(ICOMP).EQ.'E1X') MELVA1=MCHAM1.IELVAL(2)
  190.          IF (nomid.lesfac(ICOMP).EQ.'E1Y') MELVA1=MCHAM1.IELVAL(2)
  191.          IF (nomid.lesfac(ICOMP).EQ.'E1Z') MELVA1=MCHAM1.IELVAL(2)
  192.  
  193.          IF (nomid.lesfac(ICOMP).EQ.'B2X') MELVA1=MCHAM1.IELVAL(3)
  194.          IF (nomid.lesfac(ICOMP).EQ.'B2Y') MELVA1=MCHAM1.IELVAL(3)
  195.          IF (nomid.lesfac(ICOMP).EQ.'B2Z') MELVA1=MCHAM1.IELVAL(3)
  196.          IF (nomid.lesfac(ICOMP).EQ.'C2X') MELVA1=MCHAM1.IELVAL(3)
  197.          IF (nomid.lesfac(ICOMP).EQ.'C2Y') MELVA1=MCHAM1.IELVAL(3)
  198.          IF (nomid.lesfac(ICOMP).EQ.'C2Z') MELVA1=MCHAM1.IELVAL(3)
  199.          IF (nomid.lesfac(ICOMP).EQ.'D2X') MELVA1=MCHAM1.IELVAL(3)
  200.          IF (nomid.lesfac(ICOMP).EQ.'D2Y') MELVA1=MCHAM1.IELVAL(3)
  201.          IF (nomid.lesfac(ICOMP).EQ.'D2Z') MELVA1=MCHAM1.IELVAL(3)
  202.          IF (nomid.lesfac(ICOMP).EQ.'E2X') MELVA1=MCHAM1.IELVAL(3)
  203.          IF (nomid.lesfac(ICOMP).EQ.'E2Y') MELVA1=MCHAM1.IELVAL(3)
  204.          IF (nomid.lesfac(ICOMP).EQ.'E2Z') MELVA1=MCHAM1.IELVAL(3)
  205.  
  206.          SEGACT MELVA1
  207.  
  208. C====================
  209. c creation d'un maillage de multiplicateurs de lagranges enrichis
  210. C====================
  211.           NBNN=IPT1.NUM(/1)+1
  212.           NBELEM=IPT1.NUM(/2)
  213.           NBSOUS=0
  214.           NBREF=0
  215.           SEGINI, IPT4
  216.           IPT4.ITYPEL=22
  217.           IELENR=0
  218. c++++ BOUCLE sur les éléments de ipt1
  219.           DO 102 JELEM=1,IPT1.NUM(/2)
  220.              NEXIST=0
  221.              ipt4.icolor(jelem)=IPT1.icolor(jelem)
  222.              JNUM = IPT1.NUM(1,JELEM)
  223. c+++ Recherche d'une valeur non nulle du champ d'enrichissement 
  224.              VENR1 = MELVA1.VELCHE(1,JELEM)
  225.  
  226. C  On prend les elements dont le hanging node est enrichi
  227.              IF (VENR1.GT.0) THEN
  228.                  NEXIST=NEXIST+1
  229. C  On prend les element dont tout les autres noeuds sont enrichis
  230.              ELSE 
  231.                  DO 121 JNOEUD= 2 , IPT1.NUM(/1)
  232.                      VENR1 = MELVA1.VELCHE(JNOEUD,JELEM)
  233.                      IF (VENR1.eq.0) GOTO 121
  234.                          NEXIST=NEXIST+1
  235.   121            CONTINUE 
  236.              ENDIF
  237.              IF (nexist.eq.0) GOTO 102
  238.              IELENR= IELENR+1
  239. C  On recopie dans IPT4 les elements de ipt1 sur lequel on veux
  240. c  imposer une relation de compatibilité 
  241.              DO 122 I=1,IPT1.NUM(/1)
  242.                   IPT4.NUM(I+1,IELENR)=IPT1.NUM(I,JELEM)
  243.   122        CONTINUE
  244.  
  245.   102     CONTINUE
  246.  
  247.           NBELEM=IELENR
  248.           SEGADJ IPT4 
  249.           IF (ielenr.eq.0) then 
  250.             segsup ipt4          
  251.             goto 101
  252.           endif
  253.  
  254.  
  255.  
  256. C=======================================================================
  257. C creation n'un nouveau noeud pour supporter chaque multiplicateur de lagrange
  258. C=======================================================================
  259.           NBPT1=nbpts
  260.           NBPTS=NBPT1+IELENR
  261.           SEGADJ,MCOORD
  262.           DO 103 J=1,IPT4.NUM(/2)
  263.              NGLOB=(NBPT1+J-1)*(IDIM+1)
  264. C remplissage des coordonees des nouveux points
  265.              DO 131 ID= 1,IDIM
  266.                  XCOOR(NGLOB+ID)=XCOOR((IPT4.NUM(2,J)-1)*(IDIM+1)+ID)
  267.   131        CONTINUE
  268.              IPT4.NUM(1,J) = NBPT1 + J
  269.   103     CONTINUE
  270.  
  271. C====================
  272. C *** SEGMENT XMATRI
  273. C====================
  274.           NLIGRD=IPT4.NUM(/1)
  275.           NLIGRP=NLIGRD
  276.           NELRIG=IPT4.NUM(/2)
  277.           SEGINI, XMATRI
  278. c++++ BOUCLE sur les éléments de ipt4
  279.           DO 104 IELEM=1,NELRIG              
  280.              RE(1,2,IELEM)=-1.
  281.              RE(2,1,IELEM)=-1.
  282.              DO 141 ICAZ=3,NLIGRD
  283.                 RE(1,ICAZ,IELEM)=XCOEFF(IDEBUT+ICAZ)
  284.                 RE(ICAZ,1,IELEM)=RE(1,ICAZ,IELEM)
  285.   141        CONTINUE
  286.   104     CONTINUE
  287. C         write(*,*) 'COMPOSANTE FACULTATIVE DU SURE :'
  288. C         write(*,*) (nomid.lesfac(ICOMP))
  289. C         write(*,*) 'MATRICE elementaire :' 
  290. C         DO 6666 I=1,NLIGRD
  291. C            write(*,*) (RE(i,iou,1), iou=1,NLIGRD)
  292. C  6666   CONTINUE
  293.  
  294. C====================
  295. C *** SEGMENT DESCR
  296. C====================
  297.  
  298.          NEXIST=0
  299.          DO 105 ICO1=1,LNOMDD
  300.            IF (NOMDD(ICO1).EQ.nomid.lesfac(ICOMP)) NEXIST = ICO1
  301.   105    CONTINUE
  302.  
  303.          IF (NEXIST.EQ.0) THEN
  304.            CALL ERREUR(837)
  305.          ENDIF
  306.  
  307.          SEGINI, DESCR
  308.          LISINC(1)='LX  '
  309.          LISDUA(1)='FLX '
  310.          NOELEP(1)=1
  311.          NOELED(1)=1
  312.          DO 106 ICO2=2,NLIGRD
  313.            LISINC(ICO2)=NOMDD(NEXIST)
  314.            LISDUA(ICO2)=NOMDU(NEXIST)
  315.            NOELEP(ico2)=ico2
  316.            NOELED(ico2)=ico2
  317.   106    CONTINUE
  318.  
  319. C====================
  320. C *** stockage de la rigidité construite dans MRIGID
  321. C====================
  322. C       Ajustement du segment rigidite 
  323.          NRIGEL=NRIGEL +1
  324.          SEGADJ, MRIGID
  325.  
  326.  
  327. C*    ICHOLE=0
  328. C*    IMGEO1=0
  329. C*    IMGEO2=0
  330. C*    IFORIG=0
  331. C*    ISUPEQ=0
  332.          COERIG(NRIGEL)=1.
  333.          IRIGEL(1,NRIGEL)=IPT4
  334.          IRIGEL(3,NRIGEL)=DESCR
  335.          IRIGEL(4,NRIGEL)=XMATRI
  336.  
  337.          SEGDES, DESCR
  338.          SEGDES, XMATRI
  339.          SEGDES, IPT4
  340. C**********************************************************************
  341. C FIN Boucle sur les composantes primales  facultatives  du sure
  342. C**********************************************************************
  343.   101   CONTINUE
  344.   100   CONTINUE
  345.       SEGDES, nomid
  346.       segdes,IPT1, MCOORD
  347.  
  348. C      RETURN
  349.       END
  350.  
  351.  
  352.  
  353.  

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