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    PV090527  26/04/30    21:16:03     12529          
  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.       rigrel=0
  88.       SEGINI,XMATR1,XMATR2,XMATRM
  89.  
  90. *     de ICPR = tableau de la numerotation locale
  91.       SEGINI,ICPR
  92.       NLOC=0
  93.  
  94. *     de MLMOT1 = tableau des noms d'inconnue locale
  95.       JGN=4
  96.       JGM=10
  97.       SEGINI,MLMOT1
  98.       JGM1=0
  99.  
  100. *     de KINCO = tableau des ddls
  101.       NINCO=10
  102.       SEGINI KINCO
  103.  
  104. *     nombre total de DDL :      
  105.       NDDL=0
  106.  
  107.  
  108. *======================================================================*
  109. *     BOUCLE SUR LES RIGIDITES MASSE, RAIDEUR ET AMORTISSEMENT
  110. *======================================================================*
  111.  
  112.       DO 100 IMAT=1,3
  113.  
  114.         MRIGID=MAT(IMAT)
  115.         IF(MRIGID.EQ.0) GOTO 100
  116.  
  117.         SEGACT,MRIGID
  118.         NRIGEL=IRIGEL(/2)
  119.  
  120. *======================================================================*
  121. *       BOUCLE SUR LES SOUS-RIGIDITES
  122. *======================================================================*
  123.         DO 200 IRI=1,NRIGEL
  124.  
  125.           IF(AFFICH)
  126.      &    WRITE(IOIMP,*)'MATRICE ',TYPMAT(IMAT),IRI,'IEME SOUS-RIGIDITE'
  127.  
  128.           IPT1   = IRIGEL(1,IRI)
  129.           DES1   = IRIGEL(3,IRI)
  130.           XMATRI = IRIGEL(4,IRI)
  131.           SEGACT,DES1,IPT1,XMATRI
  132.  
  133. *         Verification que la matrice est carrée
  134.           NLIG1P=DES1.NOELEP(/1)
  135.           NLIG1D=DES1.NOELED(/1)
  136.           IF(NLIG1P.NE.NLIG1D) THEN
  137.             CALL ERREUR(756)
  138.             SEGDES,DES1,IPT1,XMATRI,MRIGID
  139.             RETURN
  140.           ENDIF
  141. *         rem : on pourrait aussi tester NOELEP(:) = NOELED(:)
  142. *                                        LISDUA(:) = dual de LISINC(:)
  143.           SEGACT,IPT1
  144.           NBNN1   = IPT1.NUM(/1)
  145.           NBELEM1 = IPT1.NUM(/2)
  146.  
  147.           IF(AFFICH) THEN
  148.             WRITE(IOIMP,*)'MATRICE ',IMAT,' : ',IRI,'IEME SOUS-RIGIDITE'
  149.             WRITE(IOIMP,*) '+INCO=',(DES1.LISINC(iou),iou=1,NLIG1P)
  150.             WRITE(IOIMP,*) '     #',(DES1.NOELEP(iou),iou=1,NLIG1P)
  151.           ENDIF
  152.  
  153. *         creation de KDDL = tableau local des ddls
  154.           NK=NLIG1P
  155.           SEGINI,KDDL
  156.  
  157.           XCOEF=COERIG(IRI)
  158.  
  159.  
  160. *----------------------------------------------------------------------*
  161. *         BOUCLE SUR LES ELEMENTS
  162. *----------------------------------------------------------------------*
  163.           DO 300 JEL=1,NBELEM1
  164.  
  165. *----------------------------------------------------------------------*
  166. *           BOUCLE SUR LES DDLS de cet element de cette sous matrice
  167. *----------------------------------------------------------------------*
  168.             DO 400 I=1,NLIG1P
  169.  
  170. c             recup du noeud + nom d'inconnue
  171.               INONO = DES1.NOELEP(I)
  172.               INONO2= DES1.NOELED(I)
  173.               MOPRI = DES1.LISINC(I)
  174.               MODUA = DES1.LISINC(I) 
  175.               IF(INONO.NE.INONO2) THEN
  176. c               La matrice de rigidite n'est pas carree 
  177.                 CALL ERREUR(756)
  178.                 RETURN
  179.               ENDIF
  180. c             rem : il faut aussi tester l'association primal-dual
  181.  
  182. c         --- NUMEROTATION LOCALE DES NOEUD ---
  183.               IP    = IPT1.NUM(INONO,JEL)
  184.               IPLOC = ICPR(IP)
  185. c             NOUVEAU NOEUD : ON L'AJOUTE 
  186.               IF(IPLOC.EQ.0) THEN
  187.                 NLOC =NLOC+1
  188.                 IF(NLOC.GT.NINCO) THEN
  189.                   NINCO=NINCO+10
  190.                   SEGADJ,KINCO
  191.                 ENDIF
  192.  
  193.                 IPLOC=NLOC
  194.                 ICPR(IP)=IPLOC
  195.               ELSE
  196.  
  197. c             NOEUD DEJA VU : IPLOC^ieme NOEUD LOCAL
  198.  
  199.               ENDIF
  200.  
  201.  
  202. c         --- NUMEROTATION LOCALE DES NOMS D'INCONNUES ---
  203.               DO 410 IILOC=1,JGM1
  204.                 IF(MOPRI.EQ.MLMOT1.MOTS(IILOC)) GOTO 411
  205. c               NOM D'INCONNUE DEJA VU : IILOC^ieme INCONNUE
  206.  410          CONTINUE
  207. c             NOUVEAU NOM D'INCONNUE : ON L'AJOUTE
  208.               JGM1 = JGM1 + 1
  209.               IF(JGM1.GT.MLMOT1.MOTS(/2)) THEN
  210.                 JGM = MLMOT1.MOTS(/2) + 10
  211.                 SEGADJ,MLMOT1,KINCO
  212.               ENDIF
  213.               IILOC= JGM1
  214.               MLMOT1.MOTS(IILOC)=MOPRI
  215.  411          CONTINUE
  216.  
  217. c         --- DDL = COUPLE NOEUD + NOM D'INCONNUE ---
  218.               IDDL=KINCO(IILOC,IPLOC)
  219.               IF(IDDL.EQ.0) THEN
  220.                 NDDL=NDDL+1
  221.                 IDDL=NDDL
  222.                 KINCO(IILOC,IPLOC)=IDDL
  223.               ENDIF
  224.  
  225. c             ON REMPLIT LE MELEME DE POI1 + LE MLMOTS
  226.               IF(IDDL.gt.NBELEM) THEN
  227.                 NBELEM=NBELEM+20
  228.                 SEGADJ,MELEME
  229.                 JGM=NBELEM
  230.                 SEGADJ,MLMOTS
  231.                 NLIGRD=NBELEM
  232.                 NLIGRP=NLIGRD
  233.                 rigrel=0
  234.                 SEGADJ,XMATR1,XMATR2,XMATRM                
  235.               ENDIF
  236.               NUM(1,IDDL)=IP
  237.               MOTS(IDDL)=MOPRI
  238.  
  239. c             ON REMPLIT AUSSI LE TABLEAU INVERSE KDDL
  240.               KDDL(I)=IDDL
  241.  
  242.  
  243.  
  244.  400      CONTINUE
  245. *----------------------------------------------------------------------*
  246. *         FIN DE BOUCLE SUR LES DDLS
  247. *----------------------------------------------------------------------*
  248.  
  249.           IF(AFFICH) THEN
  250.             WRITE(IOIMP,*) '+#DDL=',(KDDL(iou),iou=1,NLIG1P)
  251. c             WRITE(IOIMP,*) 'dim de XMATRI=',RE(/1),RE(/2),RE(/3)
  252.           ENDIF
  253.  
  254. *         DEBUT DE REMPLISSAGE DES XMATR*
  255.           GOTO (501,502,503),IMAT
  256.  
  257. *         MASSE -> XMATRM local
  258.  501      CONTINUE
  259.           DO J=1,NLIG1P
  260.             JDDL = KDDL(J)
  261.             DO I=1,NLIG1P
  262.               IDDL = KDDL(I)
  263.               XMATRM.RE(IDDL,JDDL,1)
  264.      &        = XMATRM.RE(IDDL,JDDL,1) + XCOEF*RE(I,J,JEL)
  265.             ENDDO
  266.           ENDDO
  267.           GOTO 509
  268.  
  269. *         RAIDEUR --> 1er quadrant de XMATR1
  270.  502      CONTINUE
  271.           DO J=1,NLIG1P
  272.             JDDL = KDDL(J)
  273.             DO I=1,NLIG1P
  274.               IDDL = KDDL(I)
  275.               XMATR1.RE(IDDL,JDDL,1)
  276.      &        = XMATR1.RE(IDDL,JDDL,1) - XCOEF*RE(I,J,JEL)
  277.             ENDDO
  278.           ENDDO
  279.           GOTO 509
  280.  
  281. *         AMORTISSEMENT --> 1er quadrant de XMATR2
  282.  503      CONTINUE
  283.           DO J=1,NLIG1P
  284.             JDDL = KDDL(J)
  285.             DO I=1,NLIG1P
  286.               IDDL = KDDL(I)
  287.               XMATR2.RE(IDDL,JDDL,1)
  288.      &        = XMATR2.RE(IDDL,JDDL,1) + XCOEF*RE(I,J,JEL)
  289.             ENDDO
  290.           ENDDO
  291.           GOTO 509
  292.  
  293.  509      CONTINUE
  294.  
  295.  
  296.  300      CONTINUE
  297. *----------------------------------------------------------------------*
  298. *         FIN DE BOUCLE SUR LES ELEMENTS
  299. *----------------------------------------------------------------------*
  300.  
  301.           SEGDES,DES1,IPT1,XMATRI
  302.  
  303.  200    CONTINUE
  304. *======================================================================*
  305. *       FIN DE BOUCLE SUR LES SOUS-RIGIDITES
  306. *======================================================================*
  307.  
  308.         SEGDES,MRIGID
  309.  
  310.  100  CONTINUE
  311. *======================================================================*
  312. *     FIN DE BOUCLE SUR LES RIGIDITES MASSE, RAIDEUR ET AMORTISSEMENT
  313. *======================================================================*
  314.  
  315. *======================================================================*
  316. *     MENAGE
  317. *======================================================================*
  318.       SEGSUP,KINCO,MLMOT1
  319.       do i=1,ICPR(/1)
  320.         ICPR(i)=0
  321.       enddo
  322.  
  323. *======================================================================*
  324. *     FINALISATION DES OBJETS RESULTATS
  325. *======================================================================*
  326.  
  327.       if(NDDL.gt.NDDMAX) then
  328.         WRITE(IOIMP,*) 'Probleme de grande taille (',NDDL,'ddls):'
  329.         WRITE(IOIMP,*) 'L execution risque de prendre du temps...'
  330.       endif
  331.  
  332. *     MISE A LA BONNE DIMENSION DES OBJETS MELEME MLMOTS
  333.       NBELEM=NDDL
  334.       SEGADJ,MELEME
  335.       JGM=NDDL
  336.       SEGADJ,MLMOTS
  337.  
  338. *     Duplication du maillage pour la creation des CHPOINTs
  339.       IF (DOUBLE) THEN
  340.         NBELEM=2*NDDL
  341.         SEGADJ,MELEME
  342.         JGM=2*NDDL
  343.         SEGADJ,MLMOTS
  344.         NP2=nbpts
  345.         DO J=1,NDDL
  346.           MOTS(NDDL+J)=MOTS(J)
  347.           IP=NUM(1,J)
  348.           IP2=ICPR(IP)
  349.           IF(IP2.EQ.0)THEN
  350.             NP2=NP2+1
  351.             IP2=NP2
  352.             ICPR(IP)=IP2
  353.           ENDIF
  354.           NUM(1,NDDL+J)=IP2
  355.         ENDDO
  356.         NBPTS=NP2
  357.         SEGADJ,MCOORD
  358. C       CB215821 : Et on ne leur donne pas de coordonnees aux points crees ?!
  359.       ENDIF
  360.  
  361. *     MISE A LA BONNE DIMENSION DES OBJETS XMATR*
  362.       NLIGRD=NDDL
  363.       NLIGRP=NLIGRD
  364.       rigrel=0
  365.       SEGADJ,XMATRM
  366.       NLIGRD=2*NDDL
  367.       NLIGRP=NLIGRD
  368.       SEGADJ,XMATR1,XMATR2
  369.  
  370. *     ECRITURE DE M DANS XMATR1 (A)  ET XMATR2 (B)
  371. *               | -K . 0  |            |  C . M  |                     *
  372. *           A = | ....... |        B = | ....... |                     *
  373. *               |  0 . M  |            |  M . 0  |                     *
  374.       DO JDDL=1,NDDL
  375.         DO IDDL=1,NDDL
  376.            XMATR1.RE(NDDL+IDDL,NDDL+JDDL,1) = XMATRM.RE(IDDL,JDDL,1)
  377.            XMATR2.RE(NDDL+IDDL,JDDL,1)      = XMATRM.RE(IDDL,JDDL,1)
  378.            XMATR2.RE(IDDL,NDDL+JDDL,1)      = XMATRM.RE(IDDL,JDDL,1)
  379.         ENDDO
  380.       ENDDO
  381.       SEGSUP,XMATRM
  382.  
  383. *     AFFICHAGE DES OBJETS RESULTATS
  384.       IF (AFFICH) THEN
  385.         WRITE(IOIMP,*) 'MATRICE A ='
  386.         WRITE(IOIMP,*) (NUM(1,iou),iou=1,NBELEM)
  387.         WRITE(IOIMP,*) (MOTS(iou),iou=1,JGM)
  388.         DO iou=1,2*NDDL
  389.            WRITE(IOIMP,*) (XMATR1.RE(iou,jou),jou=1,2*NDDL)
  390.         ENDDO
  391.         WRITE(IOIMP,*) 'MATRICE B ='
  392.         WRITE(IOIMP,*) (NUM(1,iou),iou=1,NBELEM)
  393.         WRITE(IOIMP,*) (MOTS(iou),iou=1,JGM)
  394.         DO iou=1,2*NDDL
  395.            WRITE(IOIMP,*) (XMATR2.RE(iou,jou),jou=1,2*NDDL)
  396.         ENDDO
  397.       ENDIF
  398. *
  399. *======================================================================*
  400. *     MENAGE
  401. *======================================================================*
  402.       SEGSUP,ICPR
  403.  
  404.       RETURN
  405.       END
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  

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