Télécharger la plaquette Cast3M


fr en

Télécharger qzcon3.eso

Retour à la liste

Numérotation des lignes : 

qzcon3

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

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