Télécharger qzcon3.eso

Retour à la liste

Numérotation des lignes :

  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