Télécharger la plaquette Cast3M


fr en

Télécharger menism.eso

Retour à la liste

Numérotation des lignes : 

menism

  1. C MENISM    SOURCE    CB215821  24/04/12    21:16:41     11897          
  2.       SUBROUTINE MENISM
  3. *********************************************************************
  4. *
  5. *     T_CHPO1 T_CHML1=MESM MOD1 BLOQ1 (('TOUT') or 'ROTA');
  6. *
  7. *     CALCUL DES MECHANISMES ELEMENTAIRES D'UNE STRUCTURE
  8. *     COMPOSEE DE POUTRE.
  9. *
  10. *     P.PEGON (ISPRA) AOUT 1996
  11. *********************************************************************
  12.       IMPLICIT INTEGER(I-N)
  13.       IMPLICIT REAL*8 (A-H,O-Z)
  14. -INC PPARAM
  15. -INC CCOPTIO
  16. -INC SMRIGID
  17. -INC SMMODEL
  18. -INC SMCOORD
  19. -INC SMELEME
  20. -INC SMCHAML
  21. -INC SMCHPOI
  22. -INC SMTABLE
  23.       SEGMENT ICPR(nbpts)
  24.       SEGMENT WORK0
  25.         REAL*8 C(MRIGI,2*NCOMP)
  26.         REAL*8 Q(2*NCOMP,2*NCOMP)
  27.         REAL*8 Q1(2*NCOMP,2*NCOMP+MECAN)
  28.         REAL*8 C1(MRIGI,2*NCOMP)
  29.         REAL*8 C2(MRIGI,2*NCOMP+MECAN)
  30.         REAL*8 XE(3,2)
  31.         REAL*8 V1(2*NCOMP)
  32.         REAL*8 TH(MRIGI)
  33.       ENDSEGMENT
  34.       SEGMENT WORK1
  35.         REAL*8 A(NEQUA,NINCO)
  36.       ENDSEGMENT
  37.       SEGMENT WORK2
  38.         REAL*8 B(NINCO,NINCO)
  39.       ENDSEGMENT
  40.       SEGMENT WORK3
  41.         REAL*8 B1(NNNCO,NMECA)
  42.       ENDSEGMENT
  43. *
  44.       PARAMETER(NCOMP=3,NECAN=3,MRIGI=3)
  45.       CHARACTER*(LOCOMP) COMP(NCOMP),CCOMP,CANI(NECAN)
  46.       PARAMETER (NMOCLE=2)
  47.       CHARACTER*4 MOTCLE(NMOCLE)
  48. *
  49.       DATA COMP/'UX  ','UY  ','RZ  '/
  50.       DATA CANI/'RZP1','RZP2','UP2 '/
  51.       DATA MOTCLE/'TOUT','ROTA'/
  52. *
  53. *     Lecture des donnees
  54. *
  55.       CALL LIROBJ('MMODEL  ',MMODEL,1,IRET)
  56.       IF(IRET.EQ.0)RETURN
  57.       CALL LIROBJ('RIGIDITE',MRIGID,1,IRET)
  58.       IF(IRET.EQ.0)RETURN
  59.       CALL LIRMOT(MOTCLE,NMOCLE,IVAL,0)
  60.       IF(IVAL.EQ.0.OR.IVAL.EQ.1) THEN
  61.          MECAN=3
  62.       ELSE IF(IVAL.EQ.2)THEN
  63.          MECAN=2
  64.       ENDIF
  65. *
  66. *     On verifie que l'on est en 2D et que les elements sont de
  67. *     poutre ou de timo
  68. *
  69.       IF (IDIM.NE.2)THEN
  70.         WRITE(IOIMP,*)'MENISM: ne fonctionne que en 2D'
  71.         CALL ERREUR(19)
  72.         RETURN
  73.       ENDIF
  74. *
  75.       IRET=0
  76.       SEGACT,MMODEL
  77.       NZONE=KMODEL(/1)
  78.       DO IZONE=1,NZONE
  79.         IMODEL=KMODEL(IZONE)
  80.         SEGACT,IMODEL
  81.         IF(NEFMOD.NE.29.AND.NEFMOD.NE.84)THEN
  82.           IRET=IRET+1
  83.         ENDIF
  84.       ENDDO
  85.       IF(IRET.GT.0)THEN
  86.         WRITE(IOIMP,*)
  87.      >   'MENISM: ne fonctionne que pour des POUT ou des TIMO'
  88.         CALL ERREUR(16)
  89.         GOTO 9999
  90.       ENDIF
  91. *
  92. *     On numerote les points support actifs
  93. *
  94.       SEGINI,ICPR
  95.       DO  I=1,nbpts
  96.         ICPR(I)=0
  97.       ENDDO
  98.       NPOI=0
  99.       NBEL=0
  100.       DO IZONE=1,NZONE
  101.         IMODEL=KMODEL(IZONE)
  102.         MELEME=IMAMOD
  103.         SEGACT,MELEME
  104.         NBELEM=NUM(/2)
  105.         NBEL=NBEL+NBELEM
  106.         DO I=1,NUM(/1)
  107.           DO J=1,NBELEM
  108.             IKI=NUM(I,J)
  109.             IF (ICPR(IKI).EQ.0)THEN
  110.               NPOI=NPOI+1
  111.               ICPR(IKI)=NPOI
  112.             ENDIF
  113.           ENDDO
  114.         ENDDO
  115.       ENDDO
  116. *
  117. *     Verification que les c.l. s'appliquent a des points
  118. *     et a des d.o.f.s du probleme
  119. *
  120.       NLIM=0
  121.       SEGACT,MRIGID
  122.       NRI=IRIGEL(/1)
  123.       NR=IRIGEL(/2)
  124.       DO I=1,NR
  125.         MELEME=IRIGEL(1,I)
  126.         DESCR =IRIGEL(3,I)
  127.         xMATRI=IRIGEL(4,I)
  128.         NEGALI=IRIGEL(6,I)
  129.         IF (NEGALI.NE.0)THEN
  130.           WRITE(IOIMP,*)
  131.      >   'MENISM: on ne veut que des rigidites d"egalites'
  132.           CALL ERREUR(21)
  133.           GOTO 9997
  134.         ENDIF
  135.         SEGACT,MELEME
  136.         IF (ITYPEL.NE.22)THEN
  137.           WRITE(IOIMP,*)
  138.      >   'MENISM: on ne veut que des rigidites de blocage'
  139.           CALL ERREUR(21)
  140.           GOTO 9997
  141.         ENDIF
  142.         NBELEM=NUM(/2)
  143.         NBNN=NUM(/1)
  144.         NLIM=NLIM+NBELEM
  145.         DO J=1,NBELEM
  146.           DO K=2,NBNN
  147.             IF(ICPR(NUM(K,J)).EQ.0)THEN
  148.               WRITE(IOIMP,*)'MENISM: ',
  149.      >          'un point soumis a blocage n"est pas dans le modele'
  150.               CALL ERREUR(21)
  151.               GOTO 9997
  152.             ENDIF
  153.           ENDDO
  154.         ENDDO
  155.         SEGACT,xMATRI
  156.         NELRIG=re(/3)
  157.         IF(NBELEM.NE.NELRIG)THEN
  158.           WRITE(IOIMP,*)
  159.      >      'MENISM: il y a un probleme sioux dans une des rigidites'
  160.           CALL ERREUR(21)
  161.           GOTO 9997
  162.         ENDIF
  163.         SEGACT,DESCR
  164.         NLIGRP=LISINC(/2)
  165.         IF(NBNN.NE.NLIGRP)THEN
  166.           WRITE(IOIMP,*)
  167.      >      'MENISM: il y a un probleme sioux dans une des rigidites'
  168.           CALL ERREUR(21)
  169.           GOTO 9997
  170.         ENDIF
  171.         DO K=2,NLIGRP
  172.           CCOMP=LISINC(K)
  173.           DO 1 J=1,NCOMP
  174.             IF(COMP(J).EQ.CCOMP)GOTO 2
  175.  1        CONTINUE
  176.           WRITE(IOIMP,*)
  177.      >      'MENISM: un des d.o.f. de blocage n"est pas dans le modele'
  178.           CALL ERREUR(21)
  179.           GOTO 9997
  180.  2        CONTINUE
  181.         ENDDO
  182.       ENDDO
  183. *
  184. *     On est en mesure de commencer le travail en se souvenant que
  185. *
  186. *       NCOMP=nb de d.o.f par point
  187. *       MECAN=nb de mecanisme par element
  188. *       MRIGI=nb de contrainte rigide
  189. *       NPOI=nb de point support de NCOMP dof
  190. *       NBEL=nb d'elements
  191. *       NLIM=nb de condition aux limites
  192. *
  193. *     En consequence on peut allouer le pb avec
  194. *
  195. *       NINCO=nb d'inconnues
  196. *       NEQUA=nb d'equations
  197. *
  198.       NINCO=NCOMP*NPOI+MECAN*NBEL
  199.       NEQUA=MRIGI*NBEL+NLIM
  200.       NNNCO=NCOMP*NPOI
  201.       NMECA=NINCO-NEQUA
  202. *
  203. *     On commence par la matrice A
  204. *
  205.       SEGINI,WORK1
  206.       CALL ZERO(A,NEQUA,NINCO)
  207. *
  208. *     Boucle sur les elements
  209. *
  210.       SEGINI,WORK0
  211.       IEL=0
  212.       DO IZONE=1,NZONE
  213.         IMODEL=KMODEL(IZONE)
  214.         MELEME=IMAMOD
  215.         DO I=1,NUM(/2)
  216.           IEL=IEL+1
  217. *
  218. *     Calcul de C2 elementaire
  219. *
  220.           CALL DOXE(XCOOR,IDIM,2,NUM,I,XE)
  221.           CALL TIPOCQ(XE, C,Q,MRIGI,2*NCOMP)
  222.           CALL TIPOQ1(Q,2*NCOMP, Q1,2*NCOMP+MECAN)
  223.           CALL MATMAT(C,Q1,MRIGI,2*NCOMP,2*NCOMP+MECAN,C2)
  224. *
  225. *     Assemblage de C2 dans A
  226. *
  227. *         1) noeud 1 de l'element
  228. *
  229.           ILIGN=MRIGI*(IEL-1)
  230.           ICOLO=NCOMP*(ICPR(NUM(1,I))-1)
  231.           DO J=1,MRIGI
  232.             DO K=1,NCOMP
  233.               A(ILIGN+J,ICOLO+K)=C2(J,K)
  234.             ENDDO
  235.           ENDDO
  236. *
  237. *         2) noeud 2 de l'element
  238. *
  239.           ICOLO=NCOMP*(ICPR(NUM(2,I))-1)
  240.           DO J=1,MRIGI
  241.             DO K=1,NCOMP
  242.               A(ILIGN+J,ICOLO+K)=C2(J,K+NCOMP)
  243.             ENDDO
  244.           ENDDO
  245. *
  246. *         3) degres suplementaires de l'element
  247. *
  248.           ICOLO=NCOMP*NPOI+MECAN*(IEL-1)
  249.           DO J=1,MRIGI
  250.             DO K=1,MECAN
  251.               A(ILIGN+J,ICOLO+K)=C2(J,K+2*NCOMP)
  252.             ENDDO
  253.           ENDDO
  254.         ENDDO
  255.       ENDDO
  256. *
  257. *     Assemblage des conditions aux limites
  258. *
  259.       ILIGN=MRIGI*NBEL
  260.       NRI=IRIGEL(/1)
  261.       NR=IRIGEL(/2)
  262.       DO I=1,NR
  263.         MELEME=IRIGEL(1,I)
  264.         DESCR=IRIGEL(3,I)
  265.         xMATRI=IRIGEL(4,I)
  266.         DO J=1,NUM(/2)
  267.           ILIGN=ILIGN+1
  268. *          XMATRI=IMATTT(J)
  269. *          SEGACT,XMATRI
  270.           DO K=2,NUM(/1)
  271.             CCOMP=LISINC(K)
  272.             DO 3 L=1,NCOMP
  273.               IF(COMP(L).EQ.CCOMP)GOTO 4
  274.  3          CONTINUE
  275.  4          CONTINUE
  276.             ICOLO=NCOMP*(ICPR(NUM(K,J))-1)
  277.             A(ILIGN,ICOLO+L)=RE(1,K,j)
  278.           ENDDO
  279. *          SEGDES,XMATRI
  280.         ENDDO
  281.       ENDDO
  282. *
  283. *     On alloue B et on appelle le solveur
  284. *
  285.       SEGINI,WORK2
  286.       CALL GAJOME(A,NEQUA,NINCO,B,IRET,IOIMP)
  287.       SEGSUP,WORK1
  288.       IF (IRET.NE.0)THEN
  289.         SEGSUP,WORK2
  290.         WRITE(IOIMP,*)
  291.      >  'MENISM: il y a un mecanisme global'
  292.         CALL ERREUR(26)
  293.         GOTO 9997
  294.       ENDIF
  295. *
  296. *     On ne garde que la partie active de B dans B1
  297. *
  298.       SEGINI,WORK3
  299.       DO J=1,NMECA
  300.         DO I=1,NNNCO
  301.           B1(I,J)=B(I,J+NEQUA)
  302.         ENDDO
  303.       ENDDO
  304.       SEGSUP,WORK2
  305. *
  306. *     On prepare les tables de sortie
  307. *
  308.       M=NMECA+1
  309.       SEGINI,MTABLE
  310.       MLOTAB=M
  311.       MTABTI(1)='MOT'
  312.       CALL POSCHA('SOUSTYPE',IRET)
  313.       MTABII(1)=IRET
  314.       MTABTV(1)='MOT'
  315.       CALL POSCHA('MECANISMES_PAR_NOEUDS',IRET)
  316.       MTABIV(1)=IRET
  317.       DO I=1,NMECA
  318.         MTABTI(I+1)='ENTIER'
  319.         MTABII(I+1)=I
  320.         MTABTV(I+1)='CHPOINT '
  321.         MTABIV(I+1)=0
  322.       ENDDO
  323.       MTAB1=MTABLE
  324. *
  325.       SEGINI,MTABLE=MTAB1
  326.       CALL POSCHA('MECANISMES_PAR_ELEMENTS',IRET)
  327.       MTABIV(1)=IRET
  328.       DO I=1,NMECA
  329.         MTABTV(I+1)='MCHAML  '
  330.       ENDDO
  331.       MTAB2=MTABLE
  332. *
  333. *     Maillage support des chanpoint solution
  334. *
  335.       NBSOUS=0
  336.       NBREF=0
  337.       NBELEM=NPOI
  338.       NBNN=1
  339.       SEGINI IPT1
  340.       IPT1.ITYPEL=1
  341.       DO I=1,nbpts
  342.         IF (ICPR(I).NE.0)THEN
  343.           IPT1.NUM(1,ICPR(I))=I
  344.         ENDIF
  345.       ENDDO
  346.       SEGDES,IPT1
  347. *
  348. *     Generation des chanpoint solution
  349. *
  350.       DO I=1,NMECA
  351.         IF(I.EQ.1)THEN
  352. *     prototype MCHPOI
  353.           NAT=2
  354.           NSOUPO=1
  355.           SEGINI,MCHPOI
  356.           MOCHDE='             CHPOINT CREE PAR MENISM'
  357.           MTYPOI='MECANISM'
  358.           JATTRI(1)=1
  359.           IFOPOI=IFOUR
  360.           MCHPO1=MCHPOI
  361. *     prototype MSOUPO
  362.           NC=NCOMP
  363.           SEGINI,MSOUPO
  364.           DO J=1,NCOMP
  365.             NOCOMP(J)=COMP(J)
  366.           ENDDO
  367.           IGEOC=IPT1
  368.           MSOUP1=MSOUPO
  369.         ELSE
  370.           SEGINI,MCHPOI=MCHPO1
  371.           SEGINI,MSOUPO=MSOUP1
  372.         ENDIF
  373.         IPCHP(1)=MSOUPO
  374.         SEGDES,MCHPOI
  375.         N=NPOI
  376.         SEGINI,MPOVAL
  377.         IPOVAL=MPOVAL
  378.         SEGDES,MSOUPO
  379.         DO J=1,NPOI
  380.           DO K=1,NCOMP
  381.             VPOCHA(J,K)=B1(NCOMP*(J-1)+K,I)
  382.           ENDDO
  383.         ENDDO
  384.         SEGDES,MPOVAL
  385.         MTAB1.MTABIV(I+1)=MCHPOI
  386.       ENDDO
  387.       SEGDES,MTAB1
  388. *
  389. *     Generation des mchaml solution
  390. *
  391.       DO I=1,NMECA
  392. *
  393. *     Chapeau
  394. *
  395.         IF(I.EQ.1)THEN
  396. *     prototype MCHAML
  397.           N2=MECAN
  398.           SEGINI,MCHAML
  399.           DO J=1,MECAN
  400.             NOMCHE(J)=CANI(J)
  401.             TYPCHE(J)='REAL*8          '
  402.           ENDDO
  403.           MCHAM1=MCHAML
  404. *     prototype MCHELM
  405.           N1=NZONE
  406.           L1=9
  407.           N3=6
  408.           SEGINI,MCHELM
  409.           TITCHE='MECANISME'
  410.           IFOCHE=IFOUR
  411.           DO IZONE=1,NZONE
  412.             IMODEL=KMODEL(IZONE)
  413.             IMACHE(IZONE)=IMAMOD
  414.             INFCHE(IZONE,1)=1
  415.             INFCHE(IZONE,2)=0
  416.             INFCHE(IZONE,3)=NIFOUR
  417.             INFCHE(IZONE,4)=0
  418.             INFCHE(IZONE,5)=0
  419.             INFCHE(IZONE,6)=0
  420.           ENDDO
  421.           MCHEL1=MCHELM
  422. *
  423.           N1PTEL=1
  424.           N2PTEL=0
  425.           N2EL=0
  426.         ELSE
  427.           SEGINI,MCHELM=MCHEL1
  428.         ENDIF
  429. *
  430. *     Boucle sur les elements
  431. *
  432.         DO IZONE=1,NZONE
  433.           SEGINI,MCHAML=MCHAM1
  434.           ICHAML(IZONE)=MCHAML
  435.           IMODEL=KMODEL(IZONE)
  436.           MELEME=IMAMOD
  437.           N1EL=NUM(/2)
  438.           DO J=1,MECAN
  439.             SEGINI,MELVAL
  440.             IELVAL(J)=MELVAL
  441.           ENDDO
  442.           DO J=1,N1EL
  443. *
  444. *     Calcul de C1 elementaire
  445. *
  446.             CALL DOXE(XCOOR,IDIM,2,NUM,J,XE)
  447.             CALL TIPOCQ(XE, C,Q,MRIGI,2*NCOMP)
  448.             CALL MATMAT(C,Q,MRIGI,2*NCOMP,2*NCOMP,C1)
  449. *
  450. *     On remplit V1
  451. *
  452. *           1) noeud 1 de l'element
  453. *
  454.             ICOLO=NCOMP*(ICPR(NUM(1,J))-1)
  455.             DO K=1,NCOMP
  456.               V1(K)=B1(ICOLO+K,I)
  457.             ENDDO
  458. *
  459. *           2) noeud 2 de l'element
  460. *
  461.             ICOLO=NCOMP*(ICPR(NUM(2,J))-1)
  462.             DO K=1,NCOMP
  463.               V1(K+NCOMP)=B1(ICOLO+K,I)
  464.             ENDDO
  465. *
  466. *     Calcul de THeta elementaire
  467. *
  468.             CALL MATVE1(C1,V1,MRIGI,2*NCOMP,TH,2)
  469. *
  470. *     Remplissage des MELVAL
  471. *
  472.             DO K=1,MECAN
  473.               MELVAL=IELVAL(K)
  474.               VELCHE(1,J)=TH(K)
  475.             ENDDO
  476.           ENDDO
  477.           DO J=1,MECAN
  478.             MELVAL=IELVAL(J)
  479.             SEGDES,MELVAL
  480.           ENDDO
  481.           SEGDES,MCHAML
  482.         ENDDO
  483.         SEGDES,MCHELM
  484.         MTAB2.MTABIV(I+1)=MCHELM
  485.       ENDDO
  486.       SEGDES,MTAB2
  487.       SEGSUP,MCHAM1
  488. *
  489. *     A faire ...
  490. *
  491. *     WARNING!!!!!!!!! SIGNE DES ROTATIONS!!!!!!! fait en modifiant C
  492.       SEGSUP,WORK3
  493.       SEGSUP,WORK0
  494. *
  495. *     On rend la main a GIBIANE
  496. *
  497.       CALL ECROBJ('TABLE',MTAB2)
  498.       CALL ECROBJ('TABLE',MTAB1)
  499. *
  500. *     Desactivation de la rigidite (imatri, descripteur et maillage)
  501. *
  502.  9997 DO I=1,NR
  503.         MELEME=IRIGEL(1,I)
  504.         SEGDES,MELEME
  505.         DESCR=IRIGEL(3,I)
  506.         SEGDES,DESCR
  507.         xMATRI=IRIGEL(4,I)
  508.         SEGDES,xMATRI
  509.       ENDDO
  510.       SEGDES,MRIGID
  511. *
  512. *     Supression de ICPR
  513. *
  514.       SEGSUP,ICPR
  515. *
  516. *     Desactivation des maillages du modele
  517. *
  518.       DO IZONE=1,NZONE
  519.         IMODEL=KMODEL(IZONE)
  520.         MELEME=IMAMOD
  521.         SEGDES,MELEME
  522.       ENDDO
  523. *
  524. *     Desactivation du modele
  525. *
  526.  9999 DO IZONE=1,NZONE
  527.         IMODEL=KMODEL(IZONE)
  528.         SEGDES,IMODEL
  529.       ENDDO
  530.       SEGDES,MMODEL
  531. *
  532.       RETURN
  533.       END
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  
  540.  
  541.  
  542.  

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