Télécharger qzcon3.eso

Retour à la liste

Numérotation des lignes :

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

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