Télécharger resou.eso

Retour à la liste

Numérotation des lignes :

  1. C RESOU SOURCE PV 16/11/03 21:15:15 9158
  2. SUBROUTINE RESOU
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C
  6. C **** CET OPERATEUR SERT A RESOUDRE UN SYSTEME D EQUATIONS LINEAIRES
  7. C **** CHPOINT = RESOU RIGIDITE CHPOINT
  8. C
  9. C
  10. -INC SMRIGID
  11. -INC SMTEXTE
  12. -INC CCOPTIO
  13. -INC SMTABLE
  14. -INC SMCHPOI
  15. -INC SMELEME
  16. -INC SMLCHPO
  17. PARAMETER(ZERO=0.D0)
  18. SEGMENT IDEMEM(0)
  19. segment ideme0(idemem(/1))
  20. segment ideme1(idemem(/1))
  21. segment idnote(0)
  22. C
  23. CHARACTER*4 LISM(8)
  24. CHARACTER*(8) ITCHPO,CHARRE1
  25. CHARACTER*72 CHARRE
  26. REAL*8 XVA
  27. LOGICAL ILOG,ILUG,casfimp
  28. DATA LISM/'NOID','NOUN','ENSE','GRAD','CHOL','STAB','NOEL','NOST'/
  29. DATA ITCHPO/'CHPOINT '/
  30. DATA ILOG/.FALSE./
  31. C
  32. C-------------------------------------------------------
  33. c LECTURE ET INITIALISATION
  34.  
  35. c LECTURE DES OPTIONS
  36. XVA=REAL(0.D0)
  37. IOB=0
  38.  
  39. iunil=0
  40. noel=0
  41. IMTVID=0
  42. NOUNIL=0
  43. NOID=0
  44. NOEN=1
  45. IGRADJ = 0
  46. ICHSKI = 0
  47. INSYM = 0
  48. KIKI=0
  49. KSYMET = 0
  50. IPSHPO = 0
  51. ISTAB=0
  52. 5 CONTINUE
  53. CALL LIRMOT(LISM,8,KIKI,0)
  54. IF (KIKI.EQ.1) NOID=1
  55. IF (KIKI.EQ.2) NOUNIL=1
  56. IF (KIKI.EQ.3) NOEN=0
  57. * IF (KIKI.EQ.4) IGRADJ = 1
  58. * IF (KIKI.EQ.5) ICHSKI = 1
  59. * IF (KIKI.EQ.4.OR.KIKI.EQ.5) KSYMET = 1
  60. IF (KIKI.eQ.6) ISTAB=1
  61. IF (KIKI.eQ.7) NOEL=1
  62. IF (KIKI.eQ.8) ISTAB=0
  63. IF (KIKI.NE.0) GOTO 5
  64. IF(NUCROU.EQ.0) THEN
  65. ICHSKI=1
  66. ELSEIF(NUCROU.EQ.1) THEN
  67. IGRADJ=1
  68. KSYMET=1
  69. ENDIF
  70. * WRITE(6,*) ' nucrou', nucrou
  71. * IF ( IGRADJ + ICHSKI .EQ. 0 ) ICHSKI = 1
  72.  
  73. c LECTURE DE LA RIGIDITE
  74. CALL LIROBJ('RIGIDITE',IPOIRI,1,IRETOU)
  75. IF(IERR.NE.0) GO TO 5000
  76. IPRIGO=IPOIRI
  77. C
  78. c LECTURE DE LA PRECISION
  79. PREC=REAL(1.D-18)
  80. CALL LIRREE(PREC,0,IRETOU)
  81. IF(IERR.NE.0) GO TO 5000
  82.  
  83. C REMPLISSAGE DU 2ND MEMBRE IDEMEM(**) A PARTIR DE ...
  84. c ... 'CHPOINT'
  85. SEGINI IDEMEM
  86. 1 CONTINUE
  87. CALL LIROBJ(ITCHPO,ISECO,0,IRETOU)
  88. IF(IRETOU.NE.0) THEN
  89. IDEMEM(**)=ISECO
  90. GO TO 1
  91. ENDIF
  92.  
  93. c ... 'TABLE DE SOUS-TYPE LIAISONS_STATIQUES'
  94. CALL LIRTAB('LIAISONS_STATIQUES',ITBAS,0,IRET)
  95.  
  96. c ... 'LISTCHPO'
  97. CALL LIROBJ('LISTCHPO',ISECO,0,IRETOU)
  98. IF(IRETOU.NE.0) THEN
  99. mlchpo=ISECO
  100. segact mlchpo
  101. n1 = ichpoi(/1)
  102. do iu = 1 , n1
  103. idemem(**) = ichpoi(iu)
  104. enddo
  105. segdes mlchpo
  106. segini mlchpo
  107. ipshpo = mlchpo
  108. ENDIF
  109. IF (IERR.NE.0) RETURN
  110.  
  111. IF (ITBAS.NE.0 .AND. IIMPI.EQ.333) THEN
  112. WRITE(IOIMP,*) 'on a lu la table des conditions aux limites'
  113. ENDIF
  114. if (itbas.ne.0) then
  115. mtab1 = itbas
  116. segact mtab1
  117. ima = mtab1.mlotab - 1
  118. segini idnote
  119. im = 0
  120. segdes mtab1
  121. else
  122. goto 90
  123. endif
  124. * boucle en cas de résolutions successives avec table
  125. 80 continue
  126. im = im + 1
  127. itmod = 0
  128. ichp0 = 0
  129. if (im.gt.ima) then
  130. if (idemem(/1).gt.0) goto 90
  131. * pas de champs de force
  132. call erreur(1)
  133. return
  134. endif
  135. CALL ACCTAB(ITBAS,'ENTIER',IM,0.d0,' ',.true.,IP0,
  136. & 'TABLE',I1,X1,CHARRE,.true.,ITMOD)
  137. if (ierr.ne.0) return
  138. c table itmod trouvee --> on recupere la force
  139. if (itmod.gt.0) then
  140. CALL ACCTAB(ITMOD,'MOT',0,0.d0,'FORCE',.true.,IP0,
  141. & 'CHPOINT',I1,X1,CHARRE,.true.,ICHP0)
  142. if (ierr.ne.0) return
  143. if (ichp0.gt.0) then
  144. idemem(**) = ichp0
  145. idnote(**) = im
  146. else
  147. call erreur(1)
  148. return
  149. endif
  150. c on cree le point repere ici
  151. CALL CREPO1 (ZERO, ZERO, ZERO, IPOIN)
  152. CALL ECCTAB(ITMOD,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0,
  153. & 'POINT',0,0.0D0,' ',.TRUE.,IPOIN)
  154. endif
  155. goto 80
  156. IF (IERR.NE.0) RETURN
  157.  
  158. C-------------------------------------------------------
  159. c DEBUT DU TRAVAIL
  160.  
  161. 90 continue
  162. * verification pas de blocage en double
  163. call verlag(ipoiri)
  164. if (ierr.ne.0) return
  165. * y a t il des matrices de relations non unilaterales
  166. ipoir0 = ipoiri
  167. mrigid=ipoiri
  168. C call prrigi(ipoiri,1)
  169. segact mrigid
  170. nrige= irigel(/1)
  171. idepe=0
  172. nbr = irigel(/2)
  173. do 1000 irig = 1,nbr
  174. meleme=irigel(1,irig)
  175. segact meleme
  176. if ((irigel(6,irig).eq.0.or.nounil.eq.1).and.itypel.eq.22)
  177. > idepe=idepe+1
  178. if (irigel(6,irig).ne.0) iunil=1
  179. segdes meleme
  180. 1000 continue
  181. * idepe=0
  182. lagdua=0
  183.  
  184. if(idepe.ne.0) then
  185. C on va separer les raideurs
  186. * write (6,*) ' nounil jrcond vaut ',nounil,jrcond
  187. if (jrcond.eq.0) then
  188. call separm(mrigid,ri1,ri2,nounil,lagdua,noel)
  189. * if (lagdua.ne.0) call ecmail(lagdua)
  190. * write(6,*) ' resou apres separm -- ri2'
  191. * call prrigi(ri2,1)
  192. * write(6,*) ' resou apres separm -- ri1'
  193. * call prrigi ( ri1,1)
  194. * write(6,*) ' fin impression ri1'
  195. segact mrigid*mod
  196. jrelim=ri1
  197. jrgard=ri2
  198. imlag=lagdua
  199. call fusrig(ri1,ri2,ipoir0)
  200. jrtot=ipoir0
  201. else
  202. ri1=jrelim
  203. ri2=jrgard
  204. ipoir0=jrtot
  205. lagdua=imlag
  206. ipt1=lagdua
  207. if (ipt1.ne.0) segact ipt1
  208. endif
  209. iri1s=ri1
  210. iri2s=ri2
  211. C
  212. 1010 continue
  213. C
  214. * mrigid matrice complete
  215. * ri1 dependance
  216. * ri2 les autres matrices
  217. * ri6 matrice de transfert
  218. * ri3 matrice reduite
  219. * ri5 matrice de transfert transposee
  220. C
  221. C on va proceder a la condensation rigidite puis forces
  222. if(jrcond.eq.0) then
  223.  
  224. CALL DEPEN3(RI1,RI6)
  225. * write (6,*) ' resou apres depen3 --- ri6 '
  226. * call prrigi(ri6,1)
  227. call scnd2 (ri2,ri6,ri3)
  228. * write (6,*) ' resou apres scnd2--- ri3 '
  229. * call prrigi(ri3,1)
  230. segact ri3
  231. if(ierr.ne.0) then
  232. segsup ri1,ri2,ri6,idemem
  233. return
  234. endif
  235.  
  236. segact mrigid*mod
  237. jrcond=ri3
  238. JRDEPP=RI6
  239. C dualisation de la (les) matrice(s) de dependance
  240. call dual00(ri6,ri5)
  241. * write(6,*) ' apres dual0 -- ri5'
  242. * call prrigi( ri5,1)
  243. jrdepd=ri5
  244. ipoiri = ri3
  245. else
  246. ipoiri= jrcond
  247. RI6 = JRDEPP
  248. ri5 = jrdepd
  249. endif
  250. * test si ri3 est vide
  251. ri3=jrcond
  252. segact ri3
  253. * write (6,*) ' dans resou ri3.irigel(/2) ',ri3.irigel(/2)
  254. if (ri3.irigel(/2).eq.0) imtvid=1
  255. C
  256. segdes ri1,ri2,mrigid
  257. C maintenant les seconds membres
  258. C write(6,*) ' ipoiri jrdepp jrdepd',ipoiri,ri6,ri5
  259. C call prrigi(ri3)
  260.  
  261. segact idemem*mod
  262. segini ideme0,ideme1
  263.  
  264. C call prrigi(ri5)
  265. * en cas de contacts on ne dualise pas . Ce sera fait dans unilater
  266. ifrot=0
  267. MRIGID=IPOIRI
  268. SEGACT MRIGID*MOD
  269. DO I=1,IRIGEL(/2)
  270. IF(IRIGEL(6,I).ne.0) ifrot=1
  271. enddo
  272. if (nounil.eq.1) ifrot=0
  273. * if (ifrot.eq.0) write (6,*) ' resou non unilateral '
  274. * if (ifrot.eq.1) write (6,*) ' resou unilateral pas de dualisation'
  275. if (lagdua.ne.0) then
  276. ipt8=lagdua
  277. if (ierr.ne.0) return
  278. segact ipt8
  279. * call ecmail(lagdua,0)
  280. endif
  281. do 1050 ig=1,idemem(/1)
  282. ichp2= idemem(ig)
  283. ideme0(ig)=ichp2
  284. * transferer les valeurs imposees des relations sur les inconnues (à éliminer)
  285. * write (6,*) ' valeurs imposees ichp2'
  286. * call ecchpo(ichp2,0)
  287. call transr(ichp2,ri1,ichp3)
  288. * call prrigi(ri1)
  289. * write (6,*) ' apres transfert ichp3'
  290. * call ecchpo(ichp3,0)
  291. ideme1(ig)=ichp3
  292. call mucpri(ichp3,ri2,ichp4)
  293. * write (6,*) ' apres mucpri ichp4 '
  294. * call ecchpo(ichp4,0)
  295. * ri2 est deja dualise. Il faut donc dedualiser ichp4
  296. call dbbcd(ichp4,lagdua)
  297. * write (6,*) ' apres dbbcd ichp4 '
  298. * call ecchpo(ichp4,0)
  299. call adchpo(ichp2,ichp4,ichp1,1.D0,-1.0D0)
  300. call dtchpo(ichp4)
  301. * write (6,*) ' apres adchpo ichp1'
  302. * call ecchpo(ichp1,0)
  303. call mucpri(ichp1,ri5,ichp2)
  304. * write (6,*) ' impression de ri5 '
  305. * call prrigi(ri5,1)
  306. * write (6,*) ' apres mucpri ichp2'
  307. * call ecchpo(ichp2,0)
  308. mchpo1=ichp1
  309. C segact mchpo1
  310. C write(6,*) 'reso mchpo1 jattri(1)',mchpo1.jattri(1)
  311. C segdes mchpo1
  312. C
  313. mchpo2= ichp1
  314. segact mchpo2*mod
  315. mchpo2.jattri(1)=2
  316. mchpo2= ichp2
  317. segact mchpo2*mod
  318. mchpo2.jattri(1)=2
  319. C write(6,*) 'reso mchpo2 jattri(1)',mchpo2.jattri(1)
  320. segdes mchpo2
  321. C
  322. call fuchpo (ichp1,ichp2,ichp3)
  323. * call dtchpo(ichp1)
  324. call dtchpo(ichp2)
  325. idemem(ig)=ichp3
  326. * vecteur a resoudre
  327. * write (6,*) ' le vecteur '
  328. * call ecchpo(ichp3,0)
  329. * dualisation des mult de lagrange
  330. if (lagdua.ne.0.and.ifrot.eq.0) call dbbch(ichp3,lagdua)
  331. * matrice
  332. * write (6,*) ' la matrice '
  333. * call prrigi(ipoiri)
  334. 1050 continue
  335. noid = 1
  336. endif
  337. 1060 continue
  338. C-------------------------------------------------------
  339.  
  340. *
  341. * Si au moins une des matrices n'est pas symétrique, on passera
  342. * par le solveur non-symétrique LDMT.
  343. *
  344. MRIGID=IPOIRI
  345. SEGACT MRIGID*MOD
  346. NRG = IRIGEL(/1)
  347. NBR = IRIGEL(/2)
  348. C ... Ceci peut arriver si par exemple on extrait la partie
  349. C symétrique d'une matrice purement antisymétrique ...
  350. * IF(NBR.EQ.0) THEN
  351. * SEGDES MRIGID
  352. * CALL ERREUR(727)
  353. * RETURN
  354. * ENDIF
  355. C ... Mais avant on va tester si la normalisation des variables
  356. C primales et duales a été demandée - ceci entraîne la perte
  357. C de la symétrie ...
  358. IF(NORINC.GT.0 .AND. NORIND.GT.0) THEN
  359. IF(KSYMET.EQ.1) THEN
  360. CALL ERREUR(19)
  361. SEGDES,MRIGID
  362. RETURN
  363. ENDIF
  364. INSYM = 1
  365. IGRADJ = 0
  366. ICHSKI = 0
  367. GOTO 15
  368. ENDIF
  369.  
  370. IF (NRG.GE.7) THEN
  371. C ... On teste si la matrice contient des matrices non symétriques ...
  372. C ... Si OUI, ce n'est pas la peine de faire les autres tests ...
  373. DO 9 IN = 1,NBR
  374. IANTI=IRIGEL(7,IN)
  375. IF(IANTI.GT.0) THEN
  376. C ... On vérifie si l'utilisateur n'a pas demandé explicitement
  377. C la résolution par Choleski ou gradient conjugué,
  378. C si OUI on râle puis on s'en va !!! ...
  379. IF(KSYMET.EQ.1) THEN
  380. CALL ERREUR(19)
  381. SEGDES,MRIGID
  382. RETURN
  383. ENDIF
  384. IF(NORINC.NE.0.AND.NORIND.EQ.0) THEN
  385. CALL ERREUR(760)
  386. SEGDES,MRIGID
  387. RETURN
  388. ENDIF
  389. INSYM = 1
  390. IGRADJ = 0
  391. ICHSKI = 0
  392. GOTO 15
  393. ENDIF
  394. 9 CONTINUE
  395.  
  396.  
  397. ENDIF
  398.  
  399. 15 CONTINUE
  400. C
  401. C **** ON S'ASSURE QU'IL N'Y A PAS D'APPUIS UNILATERAUX
  402. C
  403. if (iunil.eq.0) goto 30
  404. IF(IRIGEL(/1).LT.6) GO TO 30
  405. IF (NOUNIL.EQ.1) GOTO 30
  406. 21 CONTINUE
  407. C
  408. C **** EXISTENCE DES APPUIS UNILATERAUX
  409. C **** SI ON EST DEJA PASSE DANS LES APPUIS UNILATERAUX
  410. C ISUPEQ POINTE SUR UNE TABLE CONTENANT LES DONNEES A PASSER
  411. C A LA PROCEDURE UNILATER
  412. C
  413. ISUPLO=ISUPEQ
  414. IF (ISUPLO.NE.0) GOTO 27
  415. NNOR=0
  416. DO 22 I=1,IRIGEL(/2)
  417. IF(IRIGEL(6,I).EQ.0) NNOR=NNOR+1
  418. 22 CONTINUE
  419. IF(NNOR.EQ.0) THEN
  420. CALL ERREUR(312)
  421. RETURN
  422. ENDIF
  423. NRIGE=IRIGEL(/1)
  424. NRIGEL=NNOR
  425. SEGINI RI1
  426. NRIGEL=IRIGEL(/2)-NNOR
  427. SEGINI RI2
  428. II1=0
  429. II2=0
  430. DO 23 I=1,IRIGEL(/2)
  431. IF(IRIGEL(6,I).NE.0) THEN
  432. RI3=RI2
  433. II2=II2+1
  434. II=II2
  435. ELSE
  436. RI3=RI1
  437. II1=II1+1
  438. II=II1
  439. ENDIF
  440. DO 24 J=1,NRIGE
  441. RI3.IRIGEL(J,II) = IRIGEL(J,I)
  442. 24 CONTINUE
  443. RI3.COERIG(II)=COERIG(I)
  444. 23 CONTINUE
  445. * RI1 raideur sans condition unilaterale
  446. * RI2 conditions unilaterales
  447. CALL CRTABL(MTABLE)
  448. ISUPEQ=MTABLE
  449. * il faut aussi mettre isupeq dans la raideur initiale
  450. MRIGID=IPRIGO
  451. SEGACT MRIGID*MOD
  452. ISUPEQ=MTABLE
  453. if (idepe.ne.0) then
  454. * il faut extraire de la matrice initiale (ipoir0) les conditions unilaterales
  455. mrigid=iri2s
  456. segact mrigid
  457. nrigel=0
  458. do 40 i=1,irigel(/2)
  459. if (irigel(6,i).eq.0) nrigel=nrigel+1
  460. 40 continue
  461. if (nrigel.eq.0) call erreur(312)
  462. if (ierr.ne.0) return
  463. nrige=irigel(/1)
  464. segini ri4
  465. ii1=0
  466. nrigel=irigel(/2)-nrigel
  467. segini ri5
  468. ii2=0
  469. do 41 i=1,irigel(/2)
  470. if (irigel(6,i).ne.0) goto 42
  471. ii1=ii1+1
  472. do j=1,nrige
  473. ri4.irigel(j,ii1)=irigel(j,i)
  474. enddo
  475. ri4.coerig(ii1)=coerig(i)
  476. goto 41
  477. 42 continue
  478. ii2=ii2+1
  479. do j=1,nrige
  480. ri5.irigel(j,ii2)=irigel(j,i)
  481. enddo
  482. ri5.coerig(ii2)=coerig(i)
  483. 41 continue
  484. segdes mrigid,ri4
  485. endif
  486. CALL ECCTAB(MTABLE,'ENTIER ',1,XVA,' ',ILOG,IOB,
  487. $ 'RIGIDITE',IOB,XVA,' ',ILOG,RI1)
  488. CALL ECCTAB(MTABLE,'ENTIER ',2,XVA,' ',ILOG,IOB,
  489. $ 'RIGIDITE',IOB,XVA,' ',ILOG,RI2)
  490. CALL ECCTAB(MTABLE,'ENTIER ',3,XVA,' ',ILOG,IOB,
  491. $ 'LOGIQUE ',IOB,XVA,' ',ILOG,IOB)
  492. if(idepe.ne.0) then
  493. CALL ECCTAB(MTABLE,'ENTIER ',8,XVA,' ',ILOG,IOB,
  494. $ 'RIGIDITE',IOB,XVA,' ',ILOG,iri1s)
  495. CALL ECCTAB(MTABLE,'ENTIER ',9,XVA,' ',ILOG,IOB,
  496. $ 'RIGIDITE',IOB,XVA,' ',ILOG,ri4 )
  497. CALL ECCTAB(MTABLE,'ENTIER ',12,XVA,' ',ILOG,IOB,
  498. $ 'RIGIDITE',IOB,XVA,' ',ILOG,ri5 )
  499. endif
  500. if (lagdua.ne.0)
  501. > CALL ECCTAB(MTABLE,'ENTIER ',13,XVA,' ',ILOG,IOB,
  502. $ 'MAILLAGE',IOB,XVA,' ',ILOG,lagdua)
  503. ISUPLO=MTABLE
  504. SEGDES RI1,RI2,MTABLE
  505. 27 CONTINUE
  506. MTABLE=ISUPLO
  507. SEGACT MTABLE
  508. IF(INSYM.EQ.1) THEN
  509. ILUG=.TRUE.
  510. ELSE
  511. ILUG=.FALSE.
  512. ENDIF
  513. CALL ECCTAB(MTABLE,'MOT ',4,XVA,'NSYM',ILOG,IOB,
  514. $ 'LOGIQUE ',IOB,XVA,' ',ILUG,IOB)
  515. if(idepe.ne.0) then
  516. ichp3=ideme0(1)
  517. * call copie2(ich,ichp3)
  518. CALL ECCTAB(MTABLE,'ENTIER ',10,XVA,' ',ILOG,IOB,
  519. $ 'CHPOINT ',IOB,XVA,' ',ILOG,ichp3)
  520. ichp3=ideme1(1)
  521. * call copie2(ich,ichp3)
  522. CALL ECCTAB(MTABLE,'ENTIER ',11,XVA,' ',ILOG,IOB,
  523. $ 'CHPOINT ',IOB,XVA,' ',ILOG,ichp3)
  524. endif
  525. SEGDES MRIGID
  526. DO 26 I=IDEMEM(/1),1,-1
  527. ISECO=IDEMEM(I)
  528. CALL ECROBJ (ITCHPO,ISECO )
  529. 26 CONTINUE
  530. SEGSUP IDEMEM
  531. CALL ECROBJ ('TABLE ',ISUPLO)
  532. SEGINI MTEXTE
  533. LTT=8
  534. MTEXT(1:LTT) ='UNILATER'
  535. NCART=8
  536. SEGDES MTEXTE
  537. CALL ECROBJ('TEXTE',MTEXTE)
  538. mrigid=iprigo
  539. segdes mrigid
  540. RETURN
  541.  
  542. C ... On arrive ici dans le cas où il n'y a pas d'appuis unilatéraux ...
  543. 30 CONTINUE
  544. * il se peut que le dernier chp soit du frottement
  545. * on l'enleve car il ne sert a rien si on n'appele pas unilater
  546. if (idemem(/1).gt.1.and.idepe.ne.0) then
  547. mchpoi=ideme0(idemem(/1))
  548. segact MCHPOI
  549. if (mtypoi.eq.'LX ') idemem(/1)=idemem(/1)-1
  550. endif
  551. * frottement
  552. SEGDES IDEMEM
  553. * write(6,*) ' ichski, igradj,insym ',ichski, igradj,insym
  554. * write (6,*) ' imtvid ',imtvid
  555. if (imtvid.eq.1) then
  556. * matrice vide
  557. nsoupo=0
  558. nat=0
  559. segact idemem*mod
  560. do i=1,idemem(/1)
  561. segini mchpoi
  562. idemem(i)=mchpoi
  563. enddo
  564. if (noen.eq.0) then
  565. call ecrent(0)
  566. nbelem=0
  567. nbnn=1
  568. nbsous=0
  569. nbref=0
  570. segini ipt8
  571. ipt8.itypel=1
  572. call ecrobj('MAILLAGE',ipt8)
  573. endif
  574. else
  575. * write(6,*) ' appel resou1 -- idemem(1)'
  576. * segact idemem
  577. * idesec= idemem(1)
  578. * call ecchpo(idesec,0)
  579. * write(6,*) ' appel resou1 -- ipoiri'
  580. * call prrigi ( ipoiri,1)
  581. * write(6,*) ' ichski insym ', ichski, insym
  582. IF(ICHSKI.EQ.1) CALL RESOU1(IPOIRI,IDEMEM,NOID,NOEN,PREC,ISTAB)
  583. IF(IGRADJ.EQ.1) CALL GRACO0(IPOIRI,IDEMEM,NOID,NOEN)
  584. IF(INSYM .EQ.1) CALL LDMT (IPOIRI,IDEMEM,NOID,NOEN,PREC)
  585. IF(IERR.NE.0) GO TO 5001
  586. endif
  587. C
  588. C-------------------------------------------------------
  589. C LA SOLUTION EST CALCULEE --> ON LA MET EN FORME
  590.  
  591. SEGACT IDEMEM
  592. N=IDEMEM(/1)
  593. c-----boucle sur les solutions
  594. DO 3 I=1,N
  595. IRET=IDEMEM(N+1-I)
  596. *
  597. MCHPOI= IRET
  598. C-----------------------
  599. if(idepe.ne.0) then
  600. * reintroduction des inconnues supprimees
  601. call mucpri(mchpoi,ri6,ichp3)
  602. * write (6,*) ' resou - mchpoi '
  603. * call ecchpo(mchpoi,0)
  604. * write (6,*) ' resou - ichp3'
  605. * call ecchpo(ichp3 ,0)
  606. * write (6,*) ' resou - ri6'
  607. * call prrigi(ri6,1)
  608. call adchpo(mchpoi,ichp3,ichp2,1.D0,1.D0)
  609. mchpo1=ideme1(N+1-I)
  610. * write (6,*) ' resou - iret '
  611. call adchpo(ichp2,mchpo1,iret,1.D0,1D0)
  612. * call ecchpo(mchpoi)
  613. * call ecchpo(iret)
  614. call dtchpo(mchpoi)
  615. call dtchpo(ichp3)
  616. call dtchpo(ichp2)
  617. mchpo1=iret
  618. segact mchpo1*mod
  619. mchpo1.jattri(1)=1
  620. C ------------- deplacements complets puis KU
  621. call mucpri(mchpo1,ri2 ,ichp5)
  622. * write (6,*) ' apres mucpri '
  623. * call ecchpo(ichp5,0)
  624. mchpo4=ichp5
  625. segact mchpo4*mod
  626. mchpo4.jattri(1)=1
  627. segdes mchpo4
  628. ichp6= ideme0(N+1-I)
  629. C ------- write(6,*) ' --------- KU - F '
  630. * write (6,*) ' avant adchpo ichp5 '
  631. * call ecchpo(ichp5)
  632. call adchpo(ichp5,ichp6,IRET,1D0,-1D0)
  633. call dtchpo(ichp5)
  634. * write (6,*) ' apres adchpo '
  635. * call ecchpo(iret)
  636. call remplx(ri1,iret,ichp7)
  637. if (noen.eq.0) then
  638. call lirobj('MAILLAGE',ipt8,1,iretou)
  639. if (ierr.ne.0) return
  640. segact ipt8
  641. call lirent(nben,1,iretou)
  642. endif
  643. *
  644. * on n'appelle plus verlx car je ne vois pas pourquoi on voudrait que les multiplicateurs de lagrange non éliminés soient nuls
  645. *
  646. ** call verlx(ri2,iret,mchpo1,noen,ipt8)
  647. if (noen.eq.0) then
  648. call ecrobj('MAILLAGE',ipt8)
  649. nben=ipt8.num(/2)
  650. call ecrent(nben)
  651. endif
  652. call dtchpo(iret)
  653. * write (6,*) ' apres remplx ichp7 '
  654. * call ecchpo(ichp7,0)
  655. call fuchpo(ichp7,mchpo1,iret)
  656. mchpoi=iret
  657. * supression des multiplicateurs dédoublés
  658. if (lagdua.ne.0) then
  659. call dbbcf(mchpoi,lagdua)
  660. ipt1=lagdua
  661. segdes ipt1
  662. endif
  663. endif
  664. * write (6,*) ' mchpoi en fi de resou '
  665. segact mchpoi
  666. * call ecchpo(mchpoi,0)
  667. * les champs de points qui sortent sont de nature diffuse
  668. SEGACT MCHPOI
  669. NAT = MAX(1,JATTRI(/1))
  670. NSOUPO=IPCHP(/1)
  671. SEGADJ MCHPOI
  672. JATTRI(1)=1
  673. SEGDES MCHPOI
  674. IRET = MCHPOI
  675. *
  676. * cas table de liaisons statiques
  677. if (itbas.ne.0) then
  678. il = n + 1 - i
  679. ilo = idnote(il)
  680. CALL ACCTAB(ITBAS,'ENTIER',ILO,0.d0,' ',.true.,IP0,
  681. & 'TABLE',I1,X1,CHARRE,.true.,ITMOD)
  682. if (ierr.ne.0) return
  683. c CALL CREPO1 (ZERO, ZERO, ZERO, IPOIN)
  684. c CALL ECCTAB(ITMOD,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0,
  685. c & 'POINT',0,0.0D0,' ',.TRUE.,IPOIN)
  686.  
  687. CALL ECCTAB(ITMOD,'MOT',0,0.D0,'DEFORMEE',
  688. & .TRUE.,0,'CHPOINT',0,0.D0,' ',.TRUE.,MCHPOI)
  689.  
  690. else if (ipshpo.gt.0) then
  691. mlchpo = ipshpo
  692. ichpoi(N+1-I) = iret
  693. else
  694. CALL ECROBJ (ITCHPO,IRET)
  695. endif
  696.  
  697. 3 CONTINUE
  698. c-----fin de boucle sur les solutions
  699.  
  700.  
  701. C-------------------------------------------------------
  702. c MENAGE AVANT DE QUITTER
  703.  
  704. 5001 CONTINUE
  705. if (itbas.ne.0) then
  706. segdes mtab1
  707. segsup idnote
  708. CALL ECROBJ ('TABLE ',itbas)
  709. endif
  710.  
  711. if (ipshpo.gt.0) then
  712. mlchpo = ipshpo
  713. segdes mlchpo
  714. CALL ECROBJ ('LISTCHPO ',ipshpo)
  715. endif
  716. SEGSUP IDEMEM
  717. C
  718. 5000 CONTINUE
  719. RETURN
  720. END
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  
  727.  
  728.  
  729.  
  730.  
  731.  
  732.  
  733.  
  734.  
  735.  
  736.  
  737.  

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