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

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