Télécharger la plaquette Cast3M


fr en

Télécharger qzcons.eso

Retour à la liste

Numérotation des lignes : 

qzcons

  1. C QZCONS    SOURCE    CB215821  23/01/25    21:15:31     11573          
  2.       SUBROUTINE QZCONS(RI1,RI2,RI3,XMATR1,XMATR2,MELEME,MLMOTS,DOUBLE)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5. **********************************************************************************
  6. *                                                                                *
  7. *          CONSTITUTION DES MATRICES A ET B DE TAILLE 2Nx2N POUR LE QZ           *
  8. * ______________________________________________________________________________ *
  9. *                                                                                *
  10. *                                                                                *
  11. *    CREATION : Nicolas BENECH, 24 Fevrier 1995                                  *
  12. *    MODIF : BP, 2016-01-15 on reecrit tout pour permettre que K, M et C n'aient *
  13. *            pas le meme support, mais on maintient la limitation : 1 seule      *
  14. *            inconnue possible (pas forcément ALFA, mais BETA, LX ...) par noeud *
  15. *    MODIF : BP, 2016-06-15 1ddl <=> (1noeud + 1 nom dinconnue)                  *
  16. * ______________________________________________________________________________ *
  17. *                                                                                *
  18. *    MODULE(S) APPELANT(S) : VIBRAC                                              *
  19. *                                                                                *
  20. *    MODULE(S) APPELE(S) : ASSMAT, POSPT                                         *
  21. * ______________________________________________________________________________ *
  22. *                                                                                *
  23. *   EN ENTREE :                                                                  *
  24. *      - RI1 : matrice masse                                                     *
  25. *      - RI2 : matrice raideur                                                   *
  26. *      - RI3 : matrice amortissement                                             *
  27. *      -DOUBLE : indique si on double le maillage                                *
  28. *                                                                                *
  29. *   EN SORTIE :                                                                  *
  30. *      - XMATR1, XMATR2 : les matrices A et B 2nx2n                              *
  31. *               | -K . 0  |            |  C . M  |                               *
  32. *           A = | ....... |        B = | ....... |                               *
  33. *               |  0 . M  |            |  M . 0  |                               *
  34. *                                                                                *
  35. *      - MELEME : le maillage support                                            *
  36. *      - MLMOTS : les composantes (ALFA suelement meme si non coherent !)        *
  37. **********************************************************************************
  38. *
  39.  
  40. -INC PPARAM
  41. -INC CCOPTIO
  42. -INC SMELEME
  43. -INC SMLMOTS
  44. -INC SMRIGID
  45. -INC SMCOORD
  46. *
  47.       SEGMENT ICPR(nbpts)
  48.       SEGMENT KDDL(NK)
  49.       SEGMENT KINCO(JGM,NINCO)
  50.       POINTEUR XMATRM.XMATRI
  51.  
  52.       LOGICAL AFFICH,DOUBLE
  53.       INTEGER MAT(3)
  54.       CHARACTER*8 TYPMAT(3)
  55.       DATA TYPMAT(1),TYPMAT(2),TYPMAT(3)
  56.      .    / 'MASSE   ','RIGIDITE','AMORTISS' /
  57.       CHARACTER*4 MOPRI,MODUA
  58.       PARAMETER(NDDMAX=400)
  59. *
  60. *     Affichage des messages pour verification
  61.       AFFICH = IIMPI.GE.21
  62.  
  63.       MAT(1)=RI1
  64.       MAT(2)=RI2
  65.       MAT(3)=RI3
  66. *
  67. *======================================================================*
  68. *     INITIALISATIONS
  69. *======================================================================*
  70.  
  71. *     du MELEME
  72.       NBREF=0
  73.       NBSOUS=0
  74.       NBELEM=20
  75.       NBNN=1
  76.       SEGINI,MELEME
  77.  
  78. *     du MLMOTS
  79.       JGN=4
  80.       JGM=NBELEM
  81.       SEGINI,MLMOTS
  82.  
  83. *     de XMATR1 XMATR2 ET XMATRM
  84.       NELRIG=1
  85.       NLIGRD=NBELEM
  86.       NLIGRP=NLIGRD
  87.       SEGINI,XMATR1,XMATR2,XMATRM
  88.  
  89. *     de ICPR = tableau de la numerotation locale
  90.       SEGINI,ICPR
  91.       NLOC=0
  92.  
  93. *     de MLMOT1 = tableau des noms d'inconnue locale
  94.       JGN=4
  95.       JGM=10
  96.       SEGINI,MLMOT1
  97.       JGM1=0
  98.  
  99. *     de KINCO = tableau des ddls
  100.       NINCO=10
  101.       SEGINI KINCO
  102.  
  103. *     nombre total de DDL :      
  104.       NDDL=0
  105.  
  106.  
  107. *======================================================================*
  108. *     BOUCLE SUR LES RIGIDITES MASSE, RAIDEUR ET AMORTISSEMENT
  109. *======================================================================*
  110.  
  111.       DO 100 IMAT=1,3
  112.  
  113.         MRIGID=MAT(IMAT)
  114.         IF(MRIGID.EQ.0) GOTO 100
  115.  
  116.         SEGACT,MRIGID
  117.         NRIGEL=IRIGEL(/2)
  118.  
  119. *======================================================================*
  120. *       BOUCLE SUR LES SOUS-RIGIDITES
  121. *======================================================================*
  122.         DO 200 IRI=1,NRIGEL
  123.  
  124.           IF(AFFICH)
  125.      &    WRITE(IOIMP,*)'MATRICE ',TYPMAT(IMAT),IRI,'IEME SOUS-RIGIDITE'
  126.  
  127.           IPT1   = IRIGEL(1,IRI)
  128.           DES1   = IRIGEL(3,IRI)
  129.           XMATRI = IRIGEL(4,IRI)
  130.           SEGACT,DES1,IPT1,XMATRI
  131.  
  132. *         Verification que la matrice est carrée
  133.           NLIG1P=DES1.NOELEP(/1)
  134.           NLIG1D=DES1.NOELED(/1)
  135.           IF(NLIG1P.NE.NLIG1D) THEN
  136.             CALL ERREUR(756)
  137.             SEGDES,DES1,IPT1,XMATRI,MRIGID
  138.             RETURN
  139.           ENDIF
  140. *         rem : on pourrait aussi tester NOELEP(:) = NOELED(:)
  141. *                                        LISDUA(:) = dual de LISINC(:)
  142.           SEGACT,IPT1
  143.           NBNN1   = IPT1.NUM(/1)
  144.           NBELEM1 = IPT1.NUM(/2)
  145.  
  146.           IF(AFFICH) THEN
  147.             WRITE(IOIMP,*)'MATRICE ',IMAT,' : ',IRI,'IEME SOUS-RIGIDITE'
  148.             WRITE(IOIMP,*) '+INCO=',(DES1.LISINC(iou),iou=1,NLIG1P)
  149.             WRITE(IOIMP,*) '     #',(DES1.NOELEP(iou),iou=1,NLIG1P)
  150.           ENDIF
  151.  
  152. *         creation de KDDL = tableau local des ddls
  153.           NK=NLIG1P
  154.           SEGINI,KDDL
  155.  
  156.           XCOEF=COERIG(IRI)
  157.  
  158.  
  159. *----------------------------------------------------------------------*
  160. *         BOUCLE SUR LES ELEMENTS
  161. *----------------------------------------------------------------------*
  162.           DO 300 JEL=1,NBELEM1
  163.  
  164. *----------------------------------------------------------------------*
  165. *           BOUCLE SUR LES DDLS de cet element de cette sous matrice
  166. *----------------------------------------------------------------------*
  167.             DO 400 I=1,NLIG1P
  168.  
  169. c             recup du noeud + nom d'inconnue
  170.               INONO = DES1.NOELEP(I)
  171.               INONO2= DES1.NOELED(I)
  172.               MOPRI = DES1.LISINC(I)
  173.               MODUA = DES1.LISINC(I) 
  174.               IF(INONO.NE.INONO2) THEN
  175. c               La matrice de rigidite n'est pas carree 
  176.                 CALL ERREUR(756)
  177.                 RETURN
  178.               ENDIF
  179. c             rem : il faut aussi tester l'association primal-dual
  180.  
  181. c         --- NUMEROTATION LOCALE DES NOEUD ---
  182.               IP    = IPT1.NUM(INONO,JEL)
  183.               IPLOC = ICPR(IP)
  184. c             NOUVEAU NOEUD : ON L'AJOUTE 
  185.               IF(IPLOC.EQ.0) THEN
  186.                 NLOC =NLOC+1
  187.                 IF(NLOC.GT.NINCO) THEN
  188.                   NINCO=NINCO+10
  189.                   SEGADJ,KINCO
  190.                 ENDIF
  191.  
  192.                 IPLOC=NLOC
  193.                 ICPR(IP)=IPLOC
  194.               ELSE
  195.  
  196. c             NOEUD DEJA VU : IPLOC^ieme NOEUD LOCAL
  197.  
  198.               ENDIF
  199.  
  200.  
  201. c         --- NUMEROTATION LOCALE DES NOMS D'INCONNUES ---
  202.               DO 410 IILOC=1,JGM1
  203.                 IF(MOPRI.EQ.MLMOT1.MOTS(IILOC)) GOTO 411
  204. c               NOM D'INCONNUE DEJA VU : IILOC^ieme INCONNUE
  205.  410          CONTINUE
  206. c             NOUVEAU NOM D'INCONNUE : ON L'AJOUTE
  207.               JGM1 = JGM1 + 1
  208.               IF(JGM1.GT.MLMOT1.MOTS(/2)) THEN
  209.                 JGM = MLMOT1.MOTS(/2) + 10
  210.                 SEGADJ,MLMOT1,KINCO
  211.               ENDIF
  212.               IILOC= JGM1
  213.               MLMOT1.MOTS(IILOC)=MOPRI
  214.  411          CONTINUE
  215.  
  216. c         --- DDL = COUPLE NOEUD + NOM D'INCONNUE ---
  217.               IDDL=KINCO(IILOC,IPLOC)
  218.               IF(IDDL.EQ.0) THEN
  219.                 NDDL=NDDL+1
  220.                 IDDL=NDDL
  221.                 KINCO(IILOC,IPLOC)=IDDL
  222.               ENDIF
  223.  
  224. c             ON REMPLIT LE MELEME DE POI1 + LE MLMOTS
  225.               IF(IDDL.gt.NBELEM) THEN
  226.                 NBELEM=NBELEM+20
  227.                 SEGADJ,MELEME
  228.                 JGM=NBELEM
  229.                 SEGADJ,MLMOTS
  230.                 NLIGRD=NBELEM
  231.                 NLIGRP=NLIGRD
  232.                 SEGADJ,XMATR1,XMATR2,XMATRM                
  233.               ENDIF
  234.               NUM(1,IDDL)=IP
  235.               MOTS(IDDL)=MOPRI
  236.  
  237. c             ON REMPLIT AUSSI LE TABLEAU INVERSE KDDL
  238.               KDDL(I)=IDDL
  239.  
  240.  
  241.  
  242.  400      CONTINUE
  243. *----------------------------------------------------------------------*
  244. *         FIN DE BOUCLE SUR LES DDLS
  245. *----------------------------------------------------------------------*
  246.  
  247.           IF(AFFICH) THEN
  248.             WRITE(IOIMP,*) '+#DDL=',(KDDL(iou),iou=1,NLIG1P)
  249. c             WRITE(IOIMP,*) 'dim de XMATRI=',RE(/1),RE(/2),RE(/3)
  250.           ENDIF
  251.  
  252. *         DEBUT DE REMPLISSAGE DES XMATR*
  253.           GOTO (501,502,503),IMAT
  254.  
  255. *         MASSE -> XMATRM local
  256.  501      CONTINUE
  257.           DO J=1,NLIG1P
  258.             JDDL = KDDL(J)
  259.             DO I=1,NLIG1P
  260.               IDDL = KDDL(I)
  261.               XMATRM.RE(IDDL,JDDL,1)
  262.      &        = XMATRM.RE(IDDL,JDDL,1) + XCOEF*RE(I,J,JEL)
  263.             ENDDO
  264.           ENDDO
  265.           GOTO 509
  266.  
  267. *         RAIDEUR --> 1er quadrant de XMATR1
  268.  502      CONTINUE
  269.           DO J=1,NLIG1P
  270.             JDDL = KDDL(J)
  271.             DO I=1,NLIG1P
  272.               IDDL = KDDL(I)
  273.               XMATR1.RE(IDDL,JDDL,1)
  274.      &        = XMATR1.RE(IDDL,JDDL,1) - XCOEF*RE(I,J,JEL)
  275.             ENDDO
  276.           ENDDO
  277.           GOTO 509
  278.  
  279. *         AMORTISSEMENT --> 1er quadrant de XMATR2
  280.  503      CONTINUE
  281.           DO J=1,NLIG1P
  282.             JDDL = KDDL(J)
  283.             DO I=1,NLIG1P
  284.               IDDL = KDDL(I)
  285.               XMATR2.RE(IDDL,JDDL,1)
  286.      &        = XMATR2.RE(IDDL,JDDL,1) + XCOEF*RE(I,J,JEL)
  287.             ENDDO
  288.           ENDDO
  289.           GOTO 509
  290.  
  291.  509      CONTINUE
  292.  
  293.  
  294.  300      CONTINUE
  295. *----------------------------------------------------------------------*
  296. *         FIN DE BOUCLE SUR LES ELEMENTS
  297. *----------------------------------------------------------------------*
  298.  
  299.           SEGDES,DES1,IPT1,XMATRI
  300.  
  301.  200    CONTINUE
  302. *======================================================================*
  303. *       FIN DE BOUCLE SUR LES SOUS-RIGIDITES
  304. *======================================================================*
  305.  
  306.         SEGDES,MRIGID
  307.  
  308.  100  CONTINUE
  309. *======================================================================*
  310. *     FIN DE BOUCLE SUR LES RIGIDITES MASSE, RAIDEUR ET AMORTISSEMENT
  311. *======================================================================*
  312.  
  313. *======================================================================*
  314. *     MENAGE
  315. *======================================================================*
  316.       SEGSUP,KINCO,MLMOT1
  317.       do i=1,ICPR(/1)
  318.         ICPR(i)=0
  319.       enddo
  320.  
  321. *======================================================================*
  322. *     FINALISATION DES OBJETS RESULTATS
  323. *======================================================================*
  324.  
  325.       if(NDDL.gt.NDDMAX) then
  326.         WRITE(IOIMP,*) 'Probleme de grande taille (',NDDL,'ddls):'
  327.         WRITE(IOIMP,*) 'L execution risque de prendre du temps...'
  328.       endif
  329.  
  330. *     MISE A LA BONNE DIMENSION DES OBJETS MELEME MLMOTS
  331.       NBELEM=NDDL
  332.       SEGADJ,MELEME
  333.       JGM=NDDL
  334.       SEGADJ,MLMOTS
  335.  
  336. *     Duplication du maillage pour la creation des CHPOINTs
  337.       IF (DOUBLE) THEN
  338.         NBELEM=2*NDDL
  339.         SEGADJ,MELEME
  340.         JGM=2*NDDL
  341.         SEGADJ,MLMOTS
  342.         NP2=nbpts
  343.         DO J=1,NDDL
  344.           MOTS(NDDL+J)=MOTS(J)
  345.           IP=NUM(1,J)
  346.           IP2=ICPR(IP)
  347.           IF(IP2.EQ.0)THEN
  348.             NP2=NP2+1
  349.             IP2=NP2
  350.             ICPR(IP)=IP2
  351.           ENDIF
  352.           NUM(1,NDDL+J)=IP2
  353.         ENDDO
  354.         NBPTS=NP2
  355.         SEGADJ,MCOORD
  356. C       CB215821 : Et on ne leur donne pas de coordonnees aux points crees ?!
  357.       ENDIF
  358.  
  359. *     MISE A LA BONNE DIMENSION DES OBJETS XMATR*
  360.       NLIGRD=NDDL
  361.       NLIGRP=NLIGRD
  362.       SEGADJ,XMATRM
  363.       NLIGRD=2*NDDL
  364.       NLIGRP=NLIGRD
  365.       SEGADJ,XMATR1,XMATR2
  366.  
  367. *     ECRITURE DE M DANS XMATR1 (A)  ET XMATR2 (B)
  368. *               | -K . 0  |            |  C . M  |                     *
  369. *           A = | ....... |        B = | ....... |                     *
  370. *               |  0 . M  |            |  M . 0  |                     *
  371.       DO JDDL=1,NDDL
  372.         DO IDDL=1,NDDL
  373.            XMATR1.RE(NDDL+IDDL,NDDL+JDDL,1) = XMATRM.RE(IDDL,JDDL,1)
  374.            XMATR2.RE(NDDL+IDDL,JDDL,1)      = XMATRM.RE(IDDL,JDDL,1)
  375.            XMATR2.RE(IDDL,NDDL+JDDL,1)      = XMATRM.RE(IDDL,JDDL,1)
  376.         ENDDO
  377.       ENDDO
  378.       SEGSUP,XMATRM
  379.  
  380. *     AFFICHAGE DES OBJETS RESULTATS
  381.       IF (AFFICH) THEN
  382.         WRITE(IOIMP,*) 'MATRICE A ='
  383.         WRITE(IOIMP,*) (NUM(1,iou),iou=1,NBELEM)
  384.         WRITE(IOIMP,*) (MOTS(iou),iou=1,JGM)
  385.         DO iou=1,2*NDDL
  386.            WRITE(IOIMP,*) (XMATR1.RE(iou,jou),jou=1,2*NDDL)
  387.         ENDDO
  388.         WRITE(IOIMP,*) 'MATRICE B ='
  389.         WRITE(IOIMP,*) (NUM(1,iou),iou=1,NBELEM)
  390.         WRITE(IOIMP,*) (MOTS(iou),iou=1,JGM)
  391.         DO iou=1,2*NDDL
  392.            WRITE(IOIMP,*) (XMATR2.RE(iou,jou),jou=1,2*NDDL)
  393.         ENDDO
  394.       ENDIF
  395. *
  396. *======================================================================*
  397. *     MENAGE
  398. *======================================================================*
  399.       SEGSUP,ICPR
  400.  
  401.       RETURN
  402.       END
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  

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