Télécharger qzcons.eso

Retour à la liste

Numérotation des lignes :

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

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