Télécharger qzcons.eso

Retour à la liste

Numérotation des lignes :

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

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