Télécharger pjmode.eso

Retour à la liste

Numérotation des lignes :

  1. C PJMODE SOURCE PV 17/09/29 21:15:42 9578
  2. SUBROUTINE PJMODE(ipmode)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C=======================================================================
  6. C OPERATEUR PJBA :
  7. C PROJECTION D'UN CHPOINT, D'UN CHARGEMENT OU D'UNE RIGIDITE
  8. C SUR LES ELEMENTS D'UNE BASE MODALE B.
  9. C LE RESULTAT EST DU MEME TYPE.
  10. C
  11. C SYNTAXE :
  12. C * FN = PJBA B OBJET ; SI BASE ELEMENTAIRE
  13. C * FN = PJBA B STR1 (N) OBJET ; SI BASE COMPLEXE
  14. C
  15. C OBJET POUVANT ETRE UNE FORCE OU UN CHARGEMENT,
  16. C OU UNE RIGIDITE DANS LE PREMIER CAS.
  17. C=======================================================================
  18. ***********************************************************
  19. * PROJECTION D'UNE MATRICE SUR UNE BASE DE MODES *
  20. * _______________________________________________________ *
  21. * *
  22. * DATE : le 11 Avril 1995 *
  23. * AUTEUR : Nicolas BENECH *
  24. * _______________________________________________________ *
  25. * *
  26. * MODULE(S) APPELANT(S) : PJBA *
  27. * *
  28. * MODULE(S) APPELE(S) : ACCTAB, YTMX *
  29. * _______________________________________________________ *
  30. * *
  31. * EN ENTREE : *
  32. * MRIGID : Matrice a projeter *
  33. * MTAB1 : Base de modes, reels ou complexes *
  34. * 'REEL' : indique que l'on utilise le produit *
  35. * scalaire reel (pas de conjugaison) *
  36. * *
  37. * EN SORTIE : *
  38. * RI1 : Matrice projetee (partie reelle) *
  39. * RI2 : Matrice projetee (partie imaginaire) *
  40. * _______________________________________________________ *
  41. * *
  42. * REMARQUE : *
  43. * L'operation realisee est : *
  44. * (MTAB1)t . MRIGID. MTAB1 *
  45. * Dans le cas complexe, la transposition est accompagnee *
  46. * de la conjugaison (si REEL n'est pas mentionne). *
  47. *
  48. * voir aussi PROJTA
  49. ***********************************************************
  50. *
  51. -INC SMCHPOI
  52. -INC SMCHARG
  53. -INC SMLCHPO
  54. -INC CCOPTIO
  55. -INC CCGEOME
  56. -INC CCREEL
  57. -INC SMELEME
  58. -INC CCHAMP
  59. -INC SMCHAML
  60. -INC SMMODEL
  61. -INC SMRIGID
  62. -INC SMLMOTS
  63. -INC SMLENTI
  64.  
  65. C
  66. * Declarations
  67. *
  68. PARAMETER(ZERO=0.D0)
  69. REAL*8 XVAL, RMAX
  70. CHARACTER*8 LETYPE
  71. CHARACTER*8 TYPMOD,cmate
  72. LOGICAL MODCOM,dedans,dchpo,l3,lr2,lirl
  73. INTEGER I, J, NBMOD, POS, IREEL, IVALRE, IOBRE
  74. REAL*8 XVALRE
  75. LOGICAL LOGRE
  76. segment plcf
  77. integer lpref(ldepl),ldefo(ldepl),lmade(ldepl)
  78. real*8 prmas(ldepl)
  79. endsegment
  80. segment pmod
  81. integer kdefo(nbmod)
  82. endsegment
  83. segment prigmat
  84. integer lrigmat(nrigmat,2+9)
  85. endsegment
  86. segment pmapmo
  87. integer lmapmo(nmapmo),defpmo(nmapmo),dimpmo(nmapmo)
  88. character*4 compmo(nmapmo)
  89. real*8 coepmo(nmapmo)
  90. endsegment
  91. segment pcompo
  92. character*4 mcol
  93. real*8 valmod(nipmod)
  94. endsegment
  95. LOGICAL L0,L1,lcf
  96. PARAMETER (ncod=8)
  97. CHARACTER*4 IDDL,lcod(ncod),lcof(ncod),motinc
  98. CHARACTER*8 TYPRET,CHARRE
  99. data xlopre/1.d-11/
  100. DATA KZERO/0/
  101. data lcod/'UX','UY','UZ','RX','RY','RZ','UR','UT'/
  102. data lcof/'FX','FY','FZ','MX','MY','MZ','FR','FT'/
  103.  
  104. plcf = 0
  105. jgn = 4
  106. jgm = ncod
  107. segini mlmot5
  108. segini mlmot6
  109. do io = 1,ncod
  110. mlmot5.mots(io) = lcod(io)
  111. mlmot6.mots(io) = lcof(io)
  112. enddo
  113.  
  114. modcom = .false.
  115. dchpo = .false.
  116. iriout = 0
  117. iriout1 = 0
  118. iriout2 = 0
  119. mmodel = ipmode
  120. segact mmodel
  121. n1 = kmodel(/1)
  122. segini mmode1
  123. jn = 0
  124. do im = 1, n1
  125. imodel = kmodel(im)
  126. segact imodel
  127. if (formod(1).eq.'MECANIQUE'.and.MATMOD(1).eq.'ELASTIQUE'
  128. &.and.(MATMOD(2).eq.'MODAL'.or.MATMOD(2).eq.'STATIQUE')) then
  129. jn = jn + 1
  130. mmode1.kmodel(jn) = imodel
  131. endif
  132. segdes imodel
  133. enddo
  134. segdes mmodel
  135. if (jn.ne.0) then
  136. n1 = jn
  137. segadj mmode1
  138. ipmode = mmode1
  139. else
  140. segsup mmode1
  141. * cas de projection non pr�vue
  142. call erreur(5)
  143. return
  144. endif
  145.  
  146. call lirobj('MCHAML',IPCAR1,1,iretou)
  147. if (ierr.ne.0) return
  148.  
  149. ipchpo = 0
  150. iprigi = 0
  151. call lirobj('CHARGEME',IPCHAR,0,iretou)
  152. if (iretou.eq.0) call lirobj('CHPOINT ',IPCHPO,0,iretou)
  153. if (iretou.eq.0) call lirobj('RIGIDITE',IPRIGI,0,iretou)
  154.  
  155. if (iretou.eq.0) then
  156. * manque un op�rande
  157. call erreur(5)
  158. return
  159. endif
  160.  
  161. C CB215821 : Faire appel directement a REDUAF plutot que faire du GIBIANE en ESOPE...
  162. call reduaf ( ipcar1,ipmode,IPCARA,1,iretr,kerre)
  163. if (ierr.ne.0) return
  164. if( iretr.ne.1) then
  165. call erreur (kerre)
  166. return
  167. endif
  168.  
  169. lcf = .false.
  170. mmodel = ipmode
  171. segact mmodel
  172. mchelm = ipcara
  173. segact mchelm
  174. if (ipchar.ne.0) goto 100
  175. if (iprigi.ne.0) goto 200
  176. if (ipchpo.ne.0) then
  177. n = 1
  178. segini mcharg
  179. ipchar = mcharg
  180. segini icharg
  181. kcharg(1) = icharg
  182. ichpo1 = ipchpo
  183. goto 100
  184. endif
  185.  
  186.  
  187. 100 continue
  188. MCHAR1=IPCHAR
  189. SEGINI,MCHARG=MCHAR1
  190. NBCHG=KCHARG(/1)
  191. DO 10 INCHA=1,NBCHG
  192. ICHAR1=KCHARG(INCHA)
  193. SEGINI,ICHARG=ICHAR1
  194. KCHARG(INCHA)=ICHARG
  195. IP1=ICHPO1
  196. c
  197. IRET = 0
  198. c
  199. c deplacement impose => idepi=1
  200. c force imposee => idepi=0
  201. c
  202. IDEPI = 0
  203. c idepi = -1
  204. KDEPI = 0
  205. MCHPOI = IP1
  206. SEGACT MCHPOI
  207. IF (MTYPOI.EQ.'FLX ') IDEPI = 1
  208. SEGDES MCHPOI
  209. c if (idepi.lt.0) then
  210. c moterr(1:8) = 'chpoint'
  211. c call erreur(302)
  212. c return
  213. c endif
  214. c
  215. NBNN = 1
  216. NBREF = 0
  217. NBSOUS = 0
  218. *
  219. LDEPL = kmodel(/1)
  220. if (.not.lcf) segini plcf
  221. c
  222. c
  223. c **** on initialise le chpoint
  224. c
  225. NSOUPO = 1
  226. NAT=1
  227. SEGINI,MCHPOI
  228. IRET = MCHPOI
  229. MTYPOI = ' '
  230. MOCHDE=' J''AI ETE FABRIQUE PAR L''OPERATEUR PJBA'
  231. IFOPOI = IFOUR
  232. * champ de force nodal: nature discrete
  233. JATTRI(1)=2
  234. NC = 1
  235. SEGINI,MSOUPO
  236. IPCHP(1) = MSOUPO
  237. NOHARM(1) = NIFOUR
  238. NOCOMP(1) = 'FALF'
  239.  
  240. do 101 inocomp=1,2
  241.  
  242. N = LDEPL
  243. SEGINI MPOVAL
  244. IPOVAL = MPOVAL
  245. *
  246. NBNN = 1
  247. NBELEM = LDEPL
  248. NBSOUS = 0
  249. NBREF = 0
  250. SEGINI MELEME
  251. IGEOC = MELEME
  252. ITYPEL = 1
  253.  
  254. knum = 0
  255. c
  256. c ****boucle sur les chpoints de depl
  257. c
  258. DO 11 IM = 1,kmodel(/1)
  259.  
  260. imodel = kmodel(im)
  261. segact imodel
  262. nomid = lnomid(2)
  263. segact nomid
  264. if (.not.lcf) then
  265. ipt1 = imamod
  266. segact ipt1
  267. iptr = ipt1.num(1,1)
  268. lpref(im) = iptr
  269.  
  270. indc = 1
  271. 34 if (imache(indc).ne.imamod.or.conche(indc).ne.conmod) then
  272. indc = indc + 1
  273. if (indc.gt.imache(/1)) then
  274. * champ de caracteristiques incomplet
  275. goto 99
  276. endif
  277. goto 34
  278. endif
  279.  
  280. mchaml = ichaml(indc)
  281. segact mchaml
  282. do iij = 1, nomche(/2)
  283. if (nomche(iij).eq.'DEFO') then
  284. melval = ielval(iij)
  285. segact melval
  286. ipp1 = ielche(1,1)
  287. ldefo(im) = ipp1
  288. segdes melval
  289. endif
  290. if (nomche(iij).eq.'MADE') then
  291. melval = ielval(iij)
  292. segact melval
  293. ipp2 = ielche(1,1)
  294. lmade(im) = ipp2
  295. segdes melval
  296. endif
  297. if (nomche(iij).eq.'MASS') then
  298. melval = ielval(iij)
  299. segact melval
  300. ymass = velche(1,1)
  301. prmas(im) = ymass
  302. segdes melval
  303. endif
  304. if(ldefo(im).gt.0.and.lmade(im).gt.0.and.
  305. &prmas(im).gt.0) goto 35
  306. enddo
  307. 35 continue
  308. segdes mchaml
  309. if (ldefo(im).eq.0) goto 99
  310. if (prmas(im).le.0.and.cmatee(1:5).eq.'MODAL') goto 99
  311. if (lmade(im).eq.0.and.cmatee(1:8).eq.'STATIQUE') goto 99
  312.  
  313. endif
  314.  
  315. if (NOCOMP(1).ne.lesobl(1)(1:4)) goto 11
  316. knum = knum + 1
  317.  
  318. iptr = lpref(im)
  319. ipp1 = ldefo(im)
  320. NUM(1,knum) = IPTR
  321. ICOLOR(knum) = IDCOUL
  322. XRET = 0.D0
  323. call xty1(ipp1,ip1,mlmot5,mlmot6,xret)
  324. IF (IDEPI.NE.1) THEN
  325. ELSE
  326. * ??
  327. indn = 1
  328. 45 if (nomche(indn).ne.'FREQ') then
  329. indn = indn + 1
  330. if (indn.gt.nomche(/2)) then
  331. * pas la composante FREQ
  332. goto 99
  333. endif
  334. goto 45
  335. endif
  336.  
  337. melval = ielval(indn)
  338. segact melval
  339. x1 = velche(1,1)
  340. OM = X1
  341. OM = 2.D0 * XPI * OM
  342. OM = OM * OM
  343. XRET = -XRET / OM
  344. ENDIF
  345. IF (IFOUR .EQ. 1) THEN
  346. IF (NIFOUR .NE. 0) THEN
  347. XRET = XRET*XPI
  348. ELSE
  349. XRET = XRET*2.D0*XPI
  350. ENDIF
  351. ENDIF
  352. VPOCHA(knum,1) = XRET
  353. *
  354. if (cmatee(1:5).eq.'MODAL') then
  355. ymass = prmas(im)
  356. elseif (cmatee(1:8).eq.'STATIQUE') then
  357. ipp2 = lmade(im)
  358. call xty1(ipp1,ipp2,mlmot5,mlmot6,ymass)
  359. else
  360. endif
  361. if (lmade(im).gt.0.and.ABS(XRET).gt.(1.d-10*ymass).and.
  362. & ymass.gt.0.and.cmatee(1:5).eq.'MODAL') then
  363. * kich : on enleve la projection sur base modale - a creuser pour statique !
  364. CALL ADCHPO(IP1,IPP2,IP2,1.d0,(XRET/ymass*(-1.d0)))
  365. IP1 = IP2
  366. endif
  367. *
  368. 11 CONTINUE
  369. *
  370. lcf = .true.
  371. *
  372. *
  373. if (knum.eq.kmodel(/1)) goto 102
  374. if (inocomp.eq.1) then
  375. if (knum.eq.0) then
  376. NOCOMP(1) = 'FBET'
  377. else
  378. N = knum
  379. NBELEM = knum
  380. segadj MPOVAL,MELEME
  381. SEGDES MPOVAL,MELEME
  382. NSOUPO = 2
  383. segadj MCHPOI
  384. SEGINI,MSOUPO
  385. IPCHP(2) = MSOUPO
  386. NOCOMP(1) = 'FBET'
  387. endif
  388. endif
  389. 101 continue
  390.  
  391. 102 continue
  392. N = knum
  393. NBELEM = knum
  394. segadj MPOVAL,MELEME
  395. SEGDES MPOVAL,MELEME
  396. SEGDES MSOUPO
  397.  
  398. SEGDES MCHPOI
  399. IF(IERR.NE.0) RETURN
  400. ICHPO1=IRET
  401. SEGDES,ICHARG
  402. 10 CONTINUE
  403. segsup mlmot5,mlmot6,plcf
  404. if (ipchpo.gt.0) then
  405. segsup icharg,mcharg
  406. call ecrobj('CHPOINT ',iret)
  407. goto 999
  408. endif
  409. SEGDES,MCHARG
  410. CALL ECROBJ('CHARGEME',MCHARG)
  411.  
  412. goto 999
  413. 99 segsup mpoval,msoupo,mchpoi
  414. call erreur(26)
  415. return
  416.  
  417.  
  418. 200 continue
  419. ipri1 = iprigi
  420. call SEPA(ipri1,1)
  421. ipri2 = iprigi
  422. call SEPA(ipri2,2)
  423. *
  424. *
  425. *
  426. *
  427. nmapmo = 100
  428. kpmo = 0
  429. segini pmapmo
  430. do isous = 1,kmodel(/1)
  431. imodel = kmodel(isous)
  432. segact imodel
  433. cmate = cmatee
  434. meleme = imamod
  435. segact meleme
  436. if (itypel.ne.1) call erreur(5)
  437. if (num(/1).ne.1) call erreur(5)
  438. if (cmate.eq.'STATIQUE'.or.cmate.EQ.'MODAL') then
  439. do ilp = 1,num(/2)
  440. kpmo = kpmo + 1
  441. if (kpmo.gt.nmapmo) then
  442. nmapmo = nmapmo + 100
  443. segadj pmapmo
  444. endif
  445. lmapmo(kpmo) = num(1,ilp)
  446. if (cmate.eq.'STATIQUE') then
  447. compmo(kpmo) = 'BETA'
  448. elseif (cmate.eq.'MODAL') then
  449. compmo(kpmo) = 'ALFA'
  450. endif
  451. do im = 1 , imache(/1)
  452. if (imache(im).eq.imamod) then
  453. if (conche(im).eq.conmod) then
  454. mchaml = ichaml(im)
  455. segact mchaml
  456. do iv = 1,ielval(/1)
  457. if (nomche(iv).eq.'DEFO') then
  458. melval = ielval(iv)
  459. segact melval
  460. ibmn = min(ilp,ielche(/2))
  461. defpmo(kpmo) = ielche(1,ibmn)
  462. endif
  463. if (nomche(iv).eq.'IDEF') then
  464. melval = ielval(iv)
  465. segact melval
  466. ibmn = min(ilp,ielche(/2))
  467. dimpmo(kpmo) = ielche(1,ibmn)
  468. endif
  469. enddo
  470. endif
  471. endif
  472. enddo
  473.  
  474. enddo
  475. endif
  476. segdes meleme,imodel
  477. enddo
  478.  
  479. nmapmo = kpmo
  480. segadj pmapmo
  481. nbmod = nmapmo
  482. *
  483. N1 = NBMOD
  484. nbcod = 8
  485. segini pmod
  486. SEGINI, MLCHP1
  487. SEGINI, MLCHP2
  488. jgm = 1
  489. jgn = 4
  490. segini mlmot4
  491. *
  492. * Constitution du maillage support et du segment descriptif
  493. *
  494. NBNN = NBMOD
  495. NBELEM = 1
  496. NBSOUS = 0
  497. NBREF = 0
  498. SEGINI, MELEME
  499. ITYPEL = 1
  500. *
  501. NLIGRD=NBMOD
  502. NLIGRP=NBMOD
  503. SEGINI, DESCR
  504. *
  505. mrigid = ipri1
  506. segact mrigid
  507. nrigel = coerig(/1)
  508. if (nrigel.lt.1) goto 250
  509. typmod = ' '
  510. IREEL = -1
  511. C* POS ? IF (POS.EQ.1) IREEL = 1
  512. *
  513. LETYPE = ' '
  514. DO 210 IM=1,NBMOD
  515. IPT1 = 0
  516. *
  517. imodel = kmodel(im)
  518. segact imodel
  519. ipt1 = imamod
  520. segact ipt1
  521. iptr = ipt1.num(1,1)
  522. * Cas reel ou cas complexe ?
  523. *
  524. if (dimpmo(im).gt.0) TYPMOD = 'MODE_COM'
  525.  
  526. IF (TYPMOD .EQ. 'MODE_COM') THEN
  527. MODCOM=.TRUE.
  528. mchpoi = defpmo(im)
  529. MLCHP1.ICHPOI(IM) = MCHPOI
  530. mchpoi = dimpmo(im)
  531. MLCHP2.ICHPOI(IM) = MCHPOI
  532. ELSE
  533. MODCOM = .FALSE.
  534. mchpoi = defpmo(im)
  535. MLCHP1.ICHPOI(IM) = MCHPOI
  536. ENDIF
  537. *
  538. ipt1 = iptr
  539. *
  540. MELEME.NUM(IM,1)=IPT1
  541. *
  542. if (cmatee.eq.'MODAL') then
  543. DESCR.LISINC(IM) = 'ALFA'
  544. DESCR.LISDUA(IM) = 'FALF'
  545. else if (cmatee.eq.'STATIQUE') then
  546. DESCR.LISINC(IM) = 'BETA'
  547. DESCR.LISDUA(IM) = 'FBET'
  548. endif
  549. DESCR.NOELEP(IM) = IM
  550. DESCR.NOELED(IM) = IM
  551. *
  552. 210 CONTINUE
  553. SEGDES, DESCR
  554. SEGDES, MELEME
  555. *
  556. * Constitution des segments XMATRI
  557. *
  558. NLIGRD=NBMOD
  559. NLIGRP=NBMOD
  560. nelrig=1
  561. *
  562. IF (LETYPE .EQ. 'ANNULE') THEN
  563. SEGINI, XMATR1
  564. IF (MODCOM) THEN
  565. SEGINI, XMATR2
  566. SEGDES, XMATR2
  567. ENDIF
  568. SEGDES, XMATR1
  569. GOTO 55
  570. ENDIF
  571. *
  572. * Cas reel
  573. *
  574. SEGINI, XMATR1
  575. DO 20, I=1, NBMOD
  576. MCHPO1 = MLCHP1.ICHPOI(I)
  577. DO 20, J=1, NBMOD
  578. MCHPO2 = MLCHP1.ICHPOI(J)
  579. CALL YTMX (MCHPO2, MCHPO1, MRIGID, XVAL)
  580. XMATR1.RE(I,J,1)=XVAL
  581. 20 CONTINUE
  582. SEGDES, XMATR1
  583. *
  584. * Cas complexe : calcul de termes complementaires
  585. *
  586. IF (MODCOM) THEN
  587. SEGACT, XMATR1*mod
  588. DO 30, I=1, NBMOD
  589. MCHPO1 = MLCHP2.ICHPOI(I)
  590. DO 30, J=1, NBMOD
  591. MCHPO2 = MLCHP2.ICHPOI(J)
  592. CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  593. XMATR1.RE(I,J,1)=XMATR1.RE(I,J,1)-IREEL*XVAL
  594. 30 CONTINUE
  595. SEGDES, XMATR1
  596. *
  597. SEGINI, XMATR2
  598. DO 40, I=1, NBMOD
  599. MCHPO1 = MLCHP1.ICHPOI(I)
  600. DO 40, J=1, NBMOD
  601. MCHPO2 = MLCHP2.ICHPOI(J)
  602. CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  603. XMATR2.RE(I,J,1)=XVAL
  604. 40 CONTINUE
  605. DO 50, I=1, NBMOD
  606. MCHPO1 = MLCHP2.ICHPOI(I)
  607. DO 50, J=1, NBMOD
  608. MCHPO2 = MLCHP1.ICHPOI(J)
  609. CALL YTMX (MCHPO1, MCHPO2, MRIGID, XVAL)
  610. XMATR2.RE(I,J,1)=XMATR2.RE(I,J,1)+IREEL*XVAL
  611. 50 CONTINUE
  612. SEGDES, XMATR2
  613. ENDIF
  614. SEGDES, MLCHP1
  615. SEGDES, MLCHP2
  616. *
  617. SEGACT, MRIGID
  618. LETYPE = MRIGID.MTYMAT
  619. SEGDES, MRIGID
  620. *
  621. * Creation des segments IMATRI
  622. *
  623. 55 NELRIG = 1
  624. * SEGINI, IMATR1
  625. * IMATR1.IMATTT(1) = XMATR1
  626. SEGDES, xMATR1
  627. IF (MODCOM) THEN
  628. * SEGINI, IMATR2
  629. * IMATR2.IMATTT(1) = XMATR2
  630. SEGDES, xMATR2
  631. ENDIF
  632. *
  633. * Creation des rigidites calculees
  634. *
  635. NRIGE=7
  636. NRIGEL=1
  637. SEGINI, RI1
  638. RI1.MTYMAT = LETYPE
  639. RI1.IFORIG = IFOMOD
  640. RI1.IMGEO1 = 0
  641. RI1.IMGEO2 = 0
  642. RI1.COERIG = 1.D0
  643. RI1.IRIGEL(1,1) = MELEME
  644. RI1.IRIGEL(2,1) = 0
  645. RI1.IRIGEL(3,1) = DESCR
  646. RI1.IRIGEL(4,1) = xMATR1
  647. RI1.IRIGEL(5,1) = NIFOUR
  648. RI1.IRIGEL(6,1) = 0
  649. RI1.IRIGEL(7,1) = 2
  650. segact xmatr1*mod
  651. xmatr1.symre=2
  652. segdes xmatr1
  653. SEGDES, RI1
  654. IF (MODCOM) THEN
  655. SEGINI, RI2 = RI1
  656. RI2.IRIGEL(4,1) = xMATR2
  657. SEGDES, RI2
  658. ELSE
  659. RI2 = 0
  660. SEGSUP, MLCHP2
  661. ENDIF
  662. *
  663. iriout1 = ri1
  664. iriout2 = ri2
  665.  
  666. 250 continue
  667. mrigid = ipri2
  668. segact mrigid
  669. nrigel = coerig(/1)
  670. if (nrigel.lt.1) goto 290
  671. typmod = ' '
  672.  
  673. nrigmat =100
  674. kgmat = 0
  675. segini prigmat
  676.  
  677. KRIGEL = 0
  678. nrigel = irigel(/2)
  679. nrige = irigel(/1)
  680. segini ri1
  681. ri1.mtymat = mtymat
  682. nrige0 = nrigel
  683.  
  684. kige = 0
  685. kige1 = 100
  686. nrigel = kige1
  687. segini ri2
  688. ri2.mtymat = mtymat
  689.  
  690. DO ire = 1,nrige0
  691. meleme = irigel (1,ire)
  692. segact meleme
  693. if (itypel.ne.22) then
  694. call erreur(977)
  695. return
  696. endif
  697. nbelem = num(/2)
  698. nbele0 = nbelem
  699. descr = irigel(3,ire)
  700. segact descr
  701. nligrp0 = noelep(/1)
  702. nligrd0 = noeled(/1)
  703. nligrp = nligrp0 + nmapmo
  704. nligrd = nligrd0 + nmapmo
  705.  
  706. nbnn = num(/1)
  707. nbsous = 0
  708. nbref = 0
  709. segini ipt2
  710. ipt2.itypel = itypel
  711. nbelem = 1
  712. nbnn = nligrd
  713. segini ipt1
  714. ipt1.itypel = itypel
  715. ri1.coerig(ire) = coerig(ire)
  716. kele = 0
  717.  
  718. xmatr1 = irigel(4,ire)
  719. segact xmatr1
  720. nelrig0 = xmatr1.re(/3)
  721. nelrig = nelrig0 + nmapmo
  722. segini xmatr2
  723. DO iele = 1,nbele0
  724. ie2 = min(iele,nelrig0)
  725. * xmatr1 = imatr1.imattt(ie2)
  726. * segact xmatr1
  727. nligrp = nligrp0 + nmapmo
  728. nligrd = nligrd0 + nmapmo
  729. nelrig=1
  730. segini des2,xmatri
  731. des2.lisinc(1) = 'LX'
  732. des2.lisdua(1) = 'FLX'
  733. des2.noelep(1) = 1
  734. des2.noeled(1) = 1
  735. * le premier point correspond aux multiplicateurs
  736. CALL CREPO1 (ZERO, ZERO, ZERO, IPTS)
  737. ipt1.num(1,1) = ipts
  738. kgrp = 1
  739. kirp = 1
  740. do ipmo = 1,nmapmo
  741. coepmo(ipmo) = 0.d0
  742. enddo
  743. do igrp = 2,nligrp0
  744. jno = noelep(igrp)
  745. motinc = lisinc(igrp)
  746. IP1 = num(jno,iele)
  747. * recherche association noeud physique - points support déformée
  748. do ilmat = 1,kgmat
  749. if(lrigmat(ilmat,1).eq.ip1) goto 315
  750. enddo
  751.  
  752. kgmat = kgmat+1
  753. ilmat = kgmat
  754. if (kgmat.gt.nrigmat) then
  755. nrigmat = nrigmat + 100
  756. segadj prigmat
  757. endif
  758. kpb = 0
  759. jg = 100
  760. segini mlent3
  761. lrigmat(kgmat,1) = ip1
  762. do ikmo = 1, nmapmo
  763. ichp1 = defpmo(ikmo)
  764. call ecrcha('NOMU')
  765. call ecrcha('MAIL')
  766. call ecrobj('CHPOINT ',ichp1)
  767. call extrai
  768. call lirobj('MAILLAGE',iuu1,1,iretou)
  769. call ecrobj('MAILLAGE',iuu1)
  770. call ecrobj('POINT ',IP1)
  771. call DANS
  772. call lirlog(l3,1,iretou)
  773. if(l3) then
  774. kpb = kpb + 1
  775. if (kpb.gt.jg) then
  776. jg = jg + 100
  777. segadj mlent3
  778. endif
  779. mlent3.lect(kpb) = ikmo
  780. endif
  781. enddo
  782. jg = kpb
  783. segadj mlent3
  784. if (kpb.gt.0) then
  785. lrigmat(ilmat,2) = mlent3
  786. else
  787. lrigmat(ilmat,2) = 0
  788. segsup mlent3
  789. endif
  790.  
  791. 315 continue
  792. ilr3 = lrigmat(ilmat,2)
  793. if (ilr3.eq.0) goto 253
  794. mlent3 = ilr3
  795. segact mlent3
  796. * selection selon nom composante
  797. mlmat = 0
  798. do lmo = 1,9
  799. if (motinc.eq.lcod(lmo)) mlmat = lmo+2
  800. enddo
  801. if (mlmat.eq.0) then
  802. * WRITE(6,*) 'coefs pour cette composante non trouves'
  803. goto 253
  804. endif
  805. if (lrigmat(ilmat,mlmat).ne.0) then
  806. pcompo = lrigmat(ilmat,mlmat)
  807. segact pcompo
  808. nipmod = valmod(/1)
  809. do ilg = 1,nipmod
  810. lkmo = mlent3.lect(ilg)
  811. coepmo(lkmo) = (valmod(ilg)* xmatr1.re(1,igrp,ie2))
  812. & + coepmo(lkmo)
  813. enddo
  814. else
  815. jg = mlent3.lect(/1)
  816. nipmod = jg
  817. segini pcompo
  818. mcol = motinc
  819. do ilg = 1,nipmod
  820. lkmo = mlent3.lect(ilg)
  821. ichp1 = defpmo(lkmo)
  822. CALL EXTRA9(ICHP1,ip1,motinc,KERRE,XFLOT)
  823. coepmo(lkmo) = (xflot * xmatr1.re(1,igrp,ie2))
  824. & + coepmo(lkmo)
  825. valmod(ilg) = xflot
  826. enddo
  827. lrigmat(ilmat,mlmat) = pcompo
  828. endif
  829.  
  830. 253 continue
  831. enddo
  832.  
  833. xmaut1 = 0.d0
  834. do kpmo = 1,nmapmo
  835. xmaut1 = max(xmaut1,ABS(coepmo(kpmo)))
  836. enddo
  837.  
  838. * synthèse
  839. do igrp = 2,nligrp0
  840. jno = noelep(igrp)
  841. motinc = lisinc(igrp)
  842. IP1 = num(jno,iele)
  843. lr2 = .false.
  844. do jgmat = 1,kgmat
  845. if(lrigmat(jgmat,1).eq.ip1) goto 325
  846. enddo
  847. c WRITE(6,*) 'bizarre, point dans l element non repertorie'
  848. call erreur(5)
  849. return
  850. 325 continue
  851. mlmat = 0
  852. do lmo = 1,9
  853. if (motinc.eq.lcod(lmo)) mlmat = lmo+2
  854. enddo
  855. if (mlmat.eq.0) lr2 = .true.
  856. if (lrigmat(jgmat,mlmat).eq.0) lr2 = .true.
  857. if(lr2) then
  858. jirp = 0
  859. do iirp = 1,kgrp
  860. if (ipt1.num(iirp,1).eq.ip1) jirp = iirp
  861. enddo
  862. c recopie
  863. kgrp = kgrp + 1
  864. if (jirp.ne.0) then
  865. des2.noelep(kgrp) = des2.noelep(jirp)
  866. des2.noeled(kgrp) = des2.noeled(jirp)
  867. else
  868. kirp = kirp + 1
  869. ipt1.num(kirp,1) = ip1
  870. des2.noelep(kgrp) = kirp
  871. des2.noeled(kgrp) = kirp
  872. endif
  873. des2.lisinc(kgrp) = lisinc(igrp)
  874. des2.lisdua(kgrp) = lisdua(igrp)
  875. re(1,kgrp,1) = xmatr1.re(1,igrp,ie2)
  876. re(kgrp,1,1) = re(1,kgrp,1)
  877. endif
  878. *
  879. enddo
  880.  
  881. do kpmo = 1,nmapmo
  882. if (ABS(coepmo(kpmo)).gt.xlopre*xmaut1) then
  883. kirp = kirp + 1
  884. kgrp = kgrp + 1
  885. ipt1.num(kirp,1) = lmapmo(kpmo)
  886. des2.noelep(kgrp) = kirp
  887. des2.noeled(kgrp) = kirp
  888. des2.lisinc(kgrp) = compmo(kpmo)
  889. if (compmo(kpmo).eq.'ALFA') des2.lisdua(kgrp) = 'FALF'
  890. if (compmo(kpmo).eq.'BETA') des2.lisdua(kgrp) = 'FBET'
  891. re(1,kgrp,1) = coepmo(kpmo)
  892. re(kgrp,1,1) = re(1,kgrp,1)
  893. endif
  894. enddo
  895. *
  896. lirl = .false.
  897. if (kirp.ne.num(/1)) then
  898. lirl = .true.
  899. else
  900. do io = 1,kirp
  901. if (num(io,iele).ne.ipt1.num(io,1)) lirl=.true.
  902. enddo
  903. endif
  904. c creation d'un irigel
  905. if (lirl) then
  906. kige = kige + 1
  907. if (kige.gt.kige1) then
  908. nrigel = kige1 + 100
  909. segadj ri2
  910. kige1 = nrigel
  911. endif
  912. nbelem = 1
  913. nbnn = kirp
  914. segini ipt3
  915. ipt3.itypel = itypel
  916. do io =1,nbnn
  917. ipt3.num(io,1) = ipt1.num(io,1)
  918. enddo
  919. nligrp = kgrp
  920. nligrd = kgrp
  921. nelrig=1
  922. segadj xmatri,des2
  923. * segini imatr3
  924. * imatr3.imattt(1) = xmatri
  925. segdes ipt3,des2,xmatri
  926. RI2.IRIGEL(1,kige) = IPT3
  927. RI2.IRIGEL(3,kige) = DES2
  928. RI2.IRIGEL(4,kige) = xmatri
  929. RI2.IRIGEL(2,kige) = 0
  930. RI2.IRIGEL(5,kige) = irigel(5,ire)
  931. RI2.IRIGEL(6,kige) = irigel(6,ire)
  932. ri2.coerig(kige) = coerig(ire)
  933. else
  934. * relation non modifiee pour cet element
  935. kele = kele + 1
  936. do ig = 1,nligrp0
  937. ipt2.num(ig,kele) = ipt1.num(ig,1)
  938. enddo
  939. * imatr2.imattt(kele) = xmatr1
  940. * kich : a tester
  941. do ju = 1,kgrp
  942. xmatr2.re(1,ju,kele) = re(1,ju,1)
  943. xmatr2.re(ju,1,kele) = re(ju,1,1)
  944. enddo
  945. segsup xmatri,des2
  946. endif
  947. ENDDO
  948.  
  949. nbelem = kele
  950. nelrig = kele
  951. nligrd=xmatr2.re(/1)
  952. nligrp=xmatr2.re(/2)
  953. if (nbelem.gt.0) then
  954. segadj ipt2
  955. segadj xmatr2
  956. krigel = krigel + 1
  957. RI1.IRIGEL(1,krigel) = IPT2
  958. RI1.IRIGEL(3,krigel) = irigel(3,ire)
  959. RI1.IRIGEL(4,krigel) = xmatr2
  960. RI1.IRIGEL(2,krigel) = 0
  961. RI1.IRIGEL(5,krigel) = irigel(5,ire)
  962. RI1.IRIGEL(6,krigel) = irigel(6,ire)
  963. segdes ipt2,xmatr2
  964. else
  965. segsup ipt2
  966. endif
  967. segsup ipt1
  968. ENDDO
  969.  
  970. iriout = 0
  971. nrigel = krigel
  972. segadj ri1
  973. nrigel = kige
  974. segadj ri2
  975. segdes mrigid,ri1,ri2
  976. if (kige.eq.0) segsup ri2
  977. if (krigel.eq.0) segsup ri1
  978. if (kige.gt.0.and.krigel.gt.0) then
  979. c WRITE(6,*) 'fus', ri1,ri2,kige,krigel
  980. call fusrig(ri1,ri2,iriout)
  981. segsup ri1, ri2
  982. return
  983. endif
  984. if (kige.gt.0) iriout = ri2
  985. if (krigel.gt.0) iriout = ri1
  986. if (iriout.eq.0) call erreur(-5)
  987. c WRITE(6,*) 'iriout', iriout
  988.  
  989. 290 continue
  990. if (iriout.ne.0) iriout3 = iriout
  991. if (iriout1.ne.0) iriout3 = iriout1
  992. if (iriout.ne.0.and.iriout1.ne.0) then
  993. call fusrig(iriout, iriout1,iriout3)
  994. ri1 = iriout
  995. ri2 = iriout1
  996. segsup ri1,ri2
  997. endif
  998.  
  999. call ecrobj('RIGIDITE',iriout3)
  1000. if (modcom) call ecrobj('RIGIDITE',iriout2)
  1001.  
  1002. goto 999
  1003.  
  1004. 199 continue
  1005. segsup descr,meleme,mlchp1,mlchp2
  1006. call erreur(5)
  1007. return
  1008.  
  1009. 999 continue
  1010.  
  1011. do ich = 1, ichaml(/1)
  1012. mchaml = ichaml(ich)
  1013. segact mchaml
  1014. do iel =1,ielval(/1)
  1015. melval = ielval(iel)
  1016. segdes melval
  1017. enddo
  1018. segdes mchaml
  1019. enddo
  1020. do k1 = 1,kmodel(/1)
  1021. imodel = kmodel(k1)
  1022. segdes imodel
  1023. enddo
  1024. C segsup mmode1, mchelm
  1025. if (plcf.ne.0) segsup plcf
  1026. return
  1027. END
  1028.  
  1029.  
  1030.  
  1031.  
  1032.  

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