Télécharger menism.eso

Retour à la liste

Numérotation des lignes :

  1. C MENISM SOURCE CHAT 11/03/16 21:27:52 6902
  2. SUBROUTINE MENISM
  3. *********************************************************************
  4. *
  5. * T_CHPO1 T_CHML1=MESM MOD1 BLOQ1 (('TOUT') or 'ROTA');
  6. *
  7. * CALCUL DES MECHANISMES ELEMENTAIRES D'UNE STRUCTURE
  8. * COMPOSEE DE POUTRE.
  9. *
  10. * P.PEGON (ISPRA) AOUT 1996
  11. *********************************************************************
  12. IMPLICIT INTEGER(I-N)
  13. IMPLICIT REAL*8 (A-H,O-Z)
  14. -INC CCOPTIO
  15. -INC SMRIGID
  16. -INC SMMODEL
  17. -INC SMCOORD
  18. -INC SMELEME
  19. -INC SMCHAML
  20. -INC SMCHPOI
  21. -INC SMTABLE
  22. SEGMENT ICPR(XCOOR(/1)/(IDIM+1))
  23. SEGMENT WORK0
  24. REAL*8 C(MRIGI,2*NCOMP)
  25. REAL*8 Q(2*NCOMP,2*NCOMP)
  26. REAL*8 Q1(2*NCOMP,2*NCOMP+MECAN)
  27. REAL*8 C1(MRIGI,2*NCOMP)
  28. REAL*8 C2(MRIGI,2*NCOMP+MECAN)
  29. REAL*8 XE(3,2)
  30. REAL*8 V1(2*NCOMP)
  31. REAL*8 TH(MRIGI)
  32. ENDSEGMENT
  33. SEGMENT WORK1
  34. REAL*8 A(NEQUA,NINCO)
  35. ENDSEGMENT
  36. SEGMENT WORK2
  37. REAL*8 B(NINCO,NINCO)
  38. ENDSEGMENT
  39. SEGMENT WORK3
  40. REAL*8 B1(NNNCO,NMECA)
  41. ENDSEGMENT
  42. *
  43. PARAMETER(NCOMP=3,NECAN=3,MRIGI=3)
  44. CHARACTER*4 COMP(NCOMP),CCOMP,CANI(NECAN)
  45. PARAMETER (NMOCLE=2)
  46. CHARACTER*4 MOTCLE(NMOCLE)
  47. *
  48. DATA COMP/'UX ','UY ','RZ '/
  49. DATA CANI/'RZP1','RZP2','UP2 '/
  50. DATA MOTCLE/'TOUT','ROTA'/
  51. *
  52. * Lecture des donnees
  53. *
  54. CALL LIROBJ('MMODEL ',MMODEL,1,IRET)
  55. IF(IRET.EQ.0)RETURN
  56. CALL LIROBJ('RIGIDITE',MRIGID,1,IRET)
  57. IF(IRET.EQ.0)RETURN
  58. CALL LIRMOT(MOTCLE,NMOCLE,IVAL,0)
  59. IF(IVAL.EQ.0.OR.IVAL.EQ.1) THEN
  60. MECAN=3
  61. ELSE IF(IVAL.EQ.2)THEN
  62. MECAN=2
  63. ENDIF
  64. *
  65. * On verifie que l'on est en 2D et que les elements sont de
  66. * poutre ou de timo
  67. *
  68. IF (IDIM.NE.2)THEN
  69. WRITE(IOIMP,*)'MENISM: ne fonctionne que en 2D'
  70. CALL ERREUR(19)
  71. RETURN
  72. ENDIF
  73. *
  74. IRET=0
  75. SEGACT,MMODEL
  76. NZONE=KMODEL(/1)
  77. DO IZONE=1,NZONE
  78. IMODEL=KMODEL(IZONE)
  79. SEGACT,IMODEL
  80. IF(NEFMOD.NE.29.AND.NEFMOD.NE.84)THEN
  81. IRET=IRET+1
  82. ENDIF
  83. ENDDO
  84. IF(IRET.GT.0)THEN
  85. WRITE(IOIMP,*)
  86. > 'MENISM: ne fonctionne que pour des POUT ou des TIMO'
  87. CALL ERREUR(16)
  88. GOTO 9999
  89. ENDIF
  90. *
  91. * On numerote les points support actifs
  92. *
  93. SEGINI,ICPR
  94. DO I=1,XCOOR(/1)/(IDIM+1)
  95. ICPR(I)=0
  96. ENDDO
  97. NPOI=0
  98. NBEL=0
  99. DO IZONE=1,NZONE
  100. IMODEL=KMODEL(IZONE)
  101. MELEME=IMAMOD
  102. SEGACT,MELEME
  103. NBELEM=NUM(/2)
  104. NBEL=NBEL+NBELEM
  105. DO I=1,NUM(/1)
  106. DO J=1,NBELEM
  107. IKI=NUM(I,J)
  108. IF (ICPR(IKI).EQ.0)THEN
  109. NPOI=NPOI+1
  110. ICPR(IKI)=NPOI
  111. ENDIF
  112. ENDDO
  113. ENDDO
  114. ENDDO
  115. *
  116. * Verification que les c.l. s'appliquent a des points
  117. * et a des d.o.f.s du probleme
  118. *
  119. NLIM=0
  120. SEGACT,MRIGID
  121. NRI=IRIGEL(/1)
  122. NR=IRIGEL(/2)
  123. DO I=1,NR
  124. MELEME=IRIGEL(1,I)
  125. DESCR =IRIGEL(3,I)
  126. xMATRI=IRIGEL(4,I)
  127. NEGALI=IRIGEL(6,I)
  128. IF (NEGALI.NE.0)THEN
  129. WRITE(IOIMP,*)
  130. > 'MENISM: on ne veut que des rigidites d"egalites'
  131. CALL ERREUR(21)
  132. GOTO 9997
  133. ENDIF
  134. SEGACT,MELEME
  135. IF (ITYPEL.NE.22)THEN
  136. WRITE(IOIMP,*)
  137. > 'MENISM: on ne veut que des rigidites de blocage'
  138. CALL ERREUR(21)
  139. GOTO 9997
  140. ENDIF
  141. NBELEM=NUM(/2)
  142. NBNN=NUM(/1)
  143. NLIM=NLIM+NBELEM
  144. DO J=1,NBELEM
  145. DO K=2,NBNN
  146. IF(ICPR(NUM(K,J)).EQ.0)THEN
  147. WRITE(IOIMP,*)'MENISM: ',
  148. > 'un point soumis a blocage n"est pas dans le modele'
  149. CALL ERREUR(21)
  150. GOTO 9997
  151. ENDIF
  152. ENDDO
  153. ENDDO
  154. SEGACT,xMATRI
  155. NELRIG=re(/3)
  156. IF(NBELEM.NE.NELRIG)THEN
  157. WRITE(IOIMP,*)
  158. > 'MENISM: il y a un probleme sioux dans une des rigidites'
  159. CALL ERREUR(21)
  160. GOTO 9997
  161. ENDIF
  162. SEGACT,DESCR
  163. NLIGRP=LISINC(/2)
  164. IF(NBNN.NE.NLIGRP)THEN
  165. WRITE(IOIMP,*)
  166. > 'MENISM: il y a un probleme sioux dans une des rigidites'
  167. CALL ERREUR(21)
  168. GOTO 9997
  169. ENDIF
  170. DO K=2,NLIGRP
  171. CCOMP=LISINC(K)
  172. DO 1 J=1,NCOMP
  173. IF(COMP(J).EQ.CCOMP)GOTO 2
  174. 1 CONTINUE
  175. WRITE(IOIMP,*)
  176. > 'MENISM: un des d.o.f. de blocage n"est pas dans le modele'
  177. CALL ERREUR(21)
  178. GOTO 9997
  179. 2 CONTINUE
  180. ENDDO
  181. ENDDO
  182. *
  183. * On est en mesure de commencer le travail en se souvenant que
  184. *
  185. * NCOMP=nb de d.o.f par point
  186. * MECAN=nb de mecanisme par element
  187. * MRIGI=nb de contrainte rigide
  188. * NPOI=nb de point support de NCOMP dof
  189. * NBEL=nb d'elements
  190. * NLIM=nb de condition aux limites
  191. *
  192. * En consequence on peut allouer le pb avec
  193. *
  194. * NINCO=nb d'inconnues
  195. * NEQUA=nb d'equations
  196. *
  197. NINCO=NCOMP*NPOI+MECAN*NBEL
  198. NEQUA=MRIGI*NBEL+NLIM
  199. NNNCO=NCOMP*NPOI
  200. NMECA=NINCO-NEQUA
  201. *
  202. * On commence par la matrice A
  203. *
  204. SEGINI,WORK1
  205. CALL ZERO(A,NEQUA,NINCO)
  206. *
  207. * Boucle sur les elements
  208. *
  209. SEGINI,WORK0
  210. IEL=0
  211. DO IZONE=1,NZONE
  212. IMODEL=KMODEL(IZONE)
  213. MELEME=IMAMOD
  214. DO I=1,NUM(/2)
  215. IEL=IEL+1
  216. *
  217. * Calcul de C2 elementaire
  218. *
  219. CALL DOXE(XCOOR,IDIM,2,NUM,I,XE)
  220. CALL TIPOCQ(XE, C,Q,MRIGI,2*NCOMP)
  221. CALL TIPOQ1(Q,2*NCOMP, Q1,2*NCOMP+MECAN)
  222. CALL MATMAT(C,Q1,MRIGI,2*NCOMP,2*NCOMP+MECAN,C2)
  223. *
  224. * Assemblage de C2 dans A
  225. *
  226. * 1) noeud 1 de l'element
  227. *
  228. ILIGN=MRIGI*(IEL-1)
  229. ICOLO=NCOMP*(ICPR(NUM(1,I))-1)
  230. DO J=1,MRIGI
  231. DO K=1,NCOMP
  232. A(ILIGN+J,ICOLO+K)=C2(J,K)
  233. ENDDO
  234. ENDDO
  235. *
  236. * 2) noeud 2 de l'element
  237. *
  238. ICOLO=NCOMP*(ICPR(NUM(2,I))-1)
  239. DO J=1,MRIGI
  240. DO K=1,NCOMP
  241. A(ILIGN+J,ICOLO+K)=C2(J,K+NCOMP)
  242. ENDDO
  243. ENDDO
  244. *
  245. * 3) degres suplementaires de l'element
  246. *
  247. ICOLO=NCOMP*NPOI+MECAN*(IEL-1)
  248. DO J=1,MRIGI
  249. DO K=1,MECAN
  250. A(ILIGN+J,ICOLO+K)=C2(J,K+2*NCOMP)
  251. ENDDO
  252. ENDDO
  253. ENDDO
  254. ENDDO
  255. *
  256. * Assemblage des conditions aux limites
  257. *
  258. ILIGN=MRIGI*NBEL
  259. NRI=IRIGEL(/1)
  260. NR=IRIGEL(/2)
  261. DO I=1,NR
  262. MELEME=IRIGEL(1,I)
  263. DESCR=IRIGEL(3,I)
  264. xMATRI=IRIGEL(4,I)
  265. DO J=1,NUM(/2)
  266. ILIGN=ILIGN+1
  267. * XMATRI=IMATTT(J)
  268. * SEGACT,XMATRI
  269. DO K=2,NUM(/1)
  270. CCOMP=LISINC(K)
  271. DO 3 L=1,NCOMP
  272. IF(COMP(L).EQ.CCOMP)GOTO 4
  273. 3 CONTINUE
  274. 4 CONTINUE
  275. ICOLO=NCOMP*(ICPR(NUM(K,J))-1)
  276. A(ILIGN,ICOLO+L)=RE(1,K,j)
  277. ENDDO
  278. * SEGDES,XMATRI
  279. ENDDO
  280. ENDDO
  281. *
  282. * On alloue B et on appelle le solveur
  283. *
  284. SEGINI,WORK2
  285. CALL GAJOME(A,NEQUA,NINCO,B,IRET,IOIMP)
  286. SEGSUP,WORK1
  287. IF (IRET.NE.0)THEN
  288. SEGSUP,WORK2
  289. WRITE(IOIMP,*)
  290. > 'MENISM: il y a un mecanisme global'
  291. CALL ERREUR(26)
  292. GOTO 9997
  293. ENDIF
  294. *
  295. * On ne garde que la partie active de B dans B1
  296. *
  297. SEGINI,WORK3
  298. DO J=1,NMECA
  299. DO I=1,NNNCO
  300. B1(I,J)=B(I,J+NEQUA)
  301. ENDDO
  302. ENDDO
  303. SEGSUP,WORK2
  304. *
  305. * On prepare les tables de sortie
  306. *
  307. M=NMECA+1
  308. SEGINI,MTABLE
  309. MLOTAB=M
  310. MTABTI(1)='MOT'
  311. CALL POSCHA('SOUSTYPE',IRET)
  312. MTABII(1)=IRET
  313. MTABTV(1)='MOT'
  314. CALL POSCHA('MECANISMES_PAR_NOEUDS',IRET)
  315. MTABIV(1)=IRET
  316. DO I=1,NMECA
  317. MTABTI(I+1)='ENTIER'
  318. MTABII(I+1)=I
  319. MTABTV(I+1)='CHPOINT '
  320. MTABIV(I+1)=0
  321. ENDDO
  322. MTAB1=MTABLE
  323. *
  324. SEGINI,MTABLE=MTAB1
  325. CALL POSCHA('MECANISMES_PAR_ELEMENTS',IRET)
  326. MTABIV(1)=IRET
  327. DO I=1,NMECA
  328. MTABTV(I+1)='MCHAML '
  329. ENDDO
  330. MTAB2=MTABLE
  331. *
  332. * Maillage support des chanpoint solution
  333. *
  334. NBSOUS=0
  335. NBREF=0
  336. NBELEM=NPOI
  337. NBNN=1
  338. SEGINI IPT1
  339. IPT1.ITYPEL=1
  340. DO I=1,XCOOR(/1)/(IDIM+1)
  341. IF (ICPR(I).NE.0)THEN
  342. IPT1.NUM(1,ICPR(I))=I
  343. ENDIF
  344. ENDDO
  345. SEGDES,IPT1
  346. *
  347. * Generation des chanpoint solution
  348. *
  349. DO I=1,NMECA
  350. IF(I.EQ.1)THEN
  351. * prototype MCHPOI
  352. NAT=2
  353. NSOUPO=1
  354. SEGINI,MCHPOI
  355. MOCHDE=' CHPOINT CREE PAR MENISM'
  356. MTYPOI='MECANISM'
  357. JATTRI(1)=1
  358. IFOPOI=IFOUR
  359. MCHPO1=MCHPOI
  360. * prototype MSOUPO
  361. NC=NCOMP
  362. SEGINI,MSOUPO
  363. DO J=1,NCOMP
  364. NOCOMP(J)=COMP(J)
  365. ENDDO
  366. IGEOC=IPT1
  367. MSOUP1=MSOUPO
  368. ELSE
  369. SEGINI,MCHPOI=MCHPO1
  370. SEGINI,MSOUPO=MSOUP1
  371. ENDIF
  372. IPCHP(1)=MSOUPO
  373. SEGDES,MCHPOI
  374. N=NPOI
  375. SEGINI,MPOVAL
  376. IPOVAL=MPOVAL
  377. SEGDES,MSOUPO
  378. DO J=1,NPOI
  379. DO K=1,NCOMP
  380. VPOCHA(J,K)=B1(NCOMP*(J-1)+K,I)
  381. ENDDO
  382. ENDDO
  383. SEGDES,MPOVAL
  384. MTAB1.MTABIV(I+1)=MCHPOI
  385. ENDDO
  386. SEGDES,MTAB1
  387. *
  388. * Generation des mchaml solution
  389. *
  390. DO I=1,NMECA
  391. *
  392. * Chapeau
  393. *
  394. IF(I.EQ.1)THEN
  395. * prototype MCHAML
  396. N2=MECAN
  397. SEGINI,MCHAML
  398. DO J=1,MECAN
  399. NOMCHE(J)=CANI(J)//' '
  400. TYPCHE(J)='REAL*8 '
  401. ENDDO
  402. MCHAM1=MCHAML
  403. * prototype MCHELM
  404. N1=NZONE
  405. L1=9
  406. N3=6
  407. SEGINI,MCHELM
  408. TITCHE='MECANISME'
  409. IFOCHE=IFOUR
  410. DO IZONE=1,NZONE
  411. IMODEL=KMODEL(IZONE)
  412. IMACHE(IZONE)=IMAMOD
  413. INFCHE(IZONE,1)=1
  414. INFCHE(IZONE,2)=0
  415. INFCHE(IZONE,3)=NIFOUR
  416. INFCHE(IZONE,4)=0
  417. INFCHE(IZONE,5)=0
  418. INFCHE(IZONE,6)=0
  419. ENDDO
  420. MCHEL1=MCHELM
  421. *
  422. N1PTEL=1
  423. N2PTEL=0
  424. N2EL=0
  425. ELSE
  426. SEGINI,MCHELM=MCHEL1
  427. ENDIF
  428. *
  429. * Boucle sur les elements
  430. *
  431. DO IZONE=1,NZONE
  432. SEGINI,MCHAML=MCHAM1
  433. ICHAML(IZONE)=MCHAML
  434. IMODEL=KMODEL(IZONE)
  435. MELEME=IMAMOD
  436. N1EL=NUM(/2)
  437. DO J=1,MECAN
  438. SEGINI,MELVAL
  439. IELVAL(J)=MELVAL
  440. ENDDO
  441. DO J=1,N1EL
  442. *
  443. * Calcul de C1 elementaire
  444. *
  445. CALL DOXE(XCOOR,IDIM,2,NUM,J,XE)
  446. CALL TIPOCQ(XE, C,Q,MRIGI,2*NCOMP)
  447. CALL MATMAT(C,Q,MRIGI,2*NCOMP,2*NCOMP,C1)
  448. *
  449. * On remplit V1
  450. *
  451. * 1) noeud 1 de l'element
  452. *
  453. ICOLO=NCOMP*(ICPR(NUM(1,J))-1)
  454. DO K=1,NCOMP
  455. V1(K)=B1(ICOLO+K,I)
  456. ENDDO
  457. *
  458. * 2) noeud 2 de l'element
  459. *
  460. ICOLO=NCOMP*(ICPR(NUM(2,J))-1)
  461. DO K=1,NCOMP
  462. V1(K+NCOMP)=B1(ICOLO+K,I)
  463. ENDDO
  464. *
  465. * Calcul de THeta elementaire
  466. *
  467. CALL MATVE1(C1,V1,MRIGI,2*NCOMP,TH,2)
  468. *
  469. * Remplissage des MELVAL
  470. *
  471. DO K=1,MECAN
  472. MELVAL=IELVAL(K)
  473. VELCHE(1,J)=TH(K)
  474. ENDDO
  475. ENDDO
  476. DO J=1,MECAN
  477. MELVAL=IELVAL(J)
  478. SEGDES,MELVAL
  479. ENDDO
  480. SEGDES,MCHAML
  481. ENDDO
  482. SEGDES,MCHELM
  483. MTAB2.MTABIV(I+1)=MCHELM
  484. ENDDO
  485. SEGDES,MTAB2
  486. SEGSUP,MCHAM1
  487. *
  488. * A faire ...
  489. *
  490. * WARNING!!!!!!!!! SIGNE DES ROTATIONS!!!!!!! fait en modifiant C
  491. SEGSUP,WORK3
  492. SEGSUP,WORK0
  493. *
  494. * On rend la main a GIBIANE
  495. *
  496. CALL ECROBJ('TABLE',MTAB2)
  497. CALL ECROBJ('TABLE',MTAB1)
  498. *
  499. * Desactivation de la rigidite (imatri, descripteur et maillage)
  500. *
  501. 9997 DO I=1,NR
  502. MELEME=IRIGEL(1,I)
  503. SEGDES,MELEME
  504. DESCR=IRIGEL(3,I)
  505. SEGDES,DESCR
  506. xMATRI=IRIGEL(4,I)
  507. SEGDES,xMATRI
  508. ENDDO
  509. SEGDES,MRIGID
  510. *
  511. * Supression de ICPR
  512. *
  513. SEGSUP,ICPR
  514. *
  515. * Desactivation des maillages du modele
  516. *
  517. DO IZONE=1,NZONE
  518. IMODEL=KMODEL(IZONE)
  519. MELEME=IMAMOD
  520. SEGDES,MELEME
  521. ENDDO
  522. *
  523. * Desactivation du modele
  524. *
  525. 9999 DO IZONE=1,NZONE
  526. IMODEL=KMODEL(IZONE)
  527. SEGDES,IMODEL
  528. ENDDO
  529. SEGDES,MMODEL
  530. *
  531. RETURN
  532. END
  533.  
  534.  
  535.  
  536.  
  537.  

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