Télécharger gaunew.eso

Retour à la liste

Numérotation des lignes :

  1. C GAUNEW SOURCE PV 20/03/24 21:17:39 10554
  2. SUBROUTINE GAUNEW
  3. C
  4. C Determination de la direction de descente pour un probleme
  5. C de moindre carre par la methode de Gauss-Newton
  6. C
  7. C Possiblilite d'introduire un parametre de viscosite (facultatif)
  8. C
  9. C TAB2=GANE TAB1 ('AMOR' FLOT1);
  10. C
  11. C CHPO1 RIGI1=GANE TAB1 'MATR' ('AMOR' FLOT1);
  12. C
  13. IMPLICIT INTEGER(I-N)
  14. IMPLICIT REAL*8 (A-H,O-Z)
  15.  
  16. -INC PPARAM
  17. -INC CCOPTIO
  18. -INC SMRIGID
  19. -INC SMCHPOI
  20. -INC SMELEME
  21. -INC CCHAMP
  22. -INC SMLENTI
  23. -INC SMLREEL
  24. -INC SMCOORD
  25. C
  26. CHARACTER*8 TYPOBJ,CHARR0
  27. CHARACTER*30 CHARRE
  28. REAL*8 X0,X1
  29. LOGICAL L0,L1
  30. *
  31. PARAMETER(NOPTIO=2)
  32. CHARACTER*4 OPTION(NOPTIO)
  33. LOGICAL LMATR
  34. DATA OPTION/'AMOR','MATR'/
  35. C
  36. C LECTURE DES OPTIONS
  37. C
  38. LMATR=.FALSE.
  39. 1 CALL LIRMOT(OPTION,NOPTIO,IRETOU,0)
  40. IF(IRETOU.EQ.0)GOTO 9
  41. GOTO(2,3),IRETOU
  42. C
  43. C OPTION AMOR: On introduit un amortissement
  44. C
  45. 2 CALL LIRREE(AMOR,1,IRETOU)
  46. IF(IERR.NE.0)RETURN
  47. GOTO 1
  48. C
  49. C OPTION MATR: Osort avec la matrice et le second membre
  50. C
  51. 3 LMATR=.TRUE.
  52. GOTO 1
  53. C
  54. 9 CONTINUE
  55. C
  56. C LECTURE DE LA TABLE
  57. C
  58. CALL LIROBJ('TABLE',MTABLE,1,IRETOU)
  59. IF(IERR.NE.0)RETURN
  60. C
  61. C INDEX TAB1.SOUSTYPE
  62. C
  63. CALL ACCTAB(MTABLE,'MOT',I0,X0,'SOUSTYPE',L0,IPO0,
  64. $ 'MOT',I1,X1,CHARRE,L1,IPO1)
  65. IF(IERR.NE.0) THEN
  66. IRET=1
  67. RETURN
  68. ENDIF
  69. IF(CHARRE(1:7).NE.'VECTEUR')THEN
  70. CALL ERREUR(648)
  71. IRET=1
  72. RETURN
  73. ENDIF
  74. C
  75. C CHPOINT INDEX 0 (valeur de la fonction)
  76. C
  77. CALL ACCTAB(MTABLE,'ENTIER',0,X0,'SOUSTYPE',L0,IPO0,
  78. $ 'CHPOINT',I1,X1,CHARRE,L1,MCHPOI)
  79. IF(IERR.NE.0) THEN
  80. IRET=1
  81. RETURN
  82. ENDIF
  83. C
  84. C VERIFICATION DU CARACTERE ELEMENTAIRE DU CHPOINT,
  85. C SANS SERIE DE FOURIER
  86. C
  87. SEGACT,MCHPOI
  88. IF(IPCHP(/1).NE.1)THEN
  89. WRITE(IOIMP,*)'GAuss NEwton: work only with elementary CHPO'
  90. MOTERR(1:8)='CHPOINT '
  91. INTERR(1)=MCHPOI
  92. CALL ERREUR(110)
  93. GOTO 9999
  94. ENDIF
  95. IF(IFOPOI.EQ.1)THEN
  96. MOTERR(1:8)='GANE '
  97. CALL ERREUR(333)
  98. GOTO 9999
  99. ENDIF
  100. C
  101. C ACTIVATIONS DIVERSES pour la fonction
  102. C
  103. MSOUPO=IPCHP(1)
  104. SEGACT,MSOUPO
  105. MELEME=IGEOC
  106. MPOVAL=IPOVAL
  107. SEGACT,MELEME,MPOVAL
  108. N1=VPOCHA(/1)
  109. NC1=VPOCHA(/2)
  110. C
  111. C ON PASSE AUX AUTRES INDEX ENTIERS DE TYPE CHPO CONSECUTIFS DE LA TABLE
  112. C
  113. JG=0
  114. SEGINI,MLENTI
  115. 10 CONTINUE
  116. TYPOBJ=' '
  117. JG=JG+1
  118. CALL ACCTAB(MTABLE,'ENTIER',JG,X0,CHARR0,L0,IPO0,
  119. $ TYPOBJ,I1,X1,CHARRE,L1,IPO1)
  120. IF(TYPOBJ.EQ.'CHPOINT')THEN
  121. SEGADJ,MLENTI
  122. LECT(JG)=IPO1
  123. GOTO 10
  124. ENDIF
  125. CONTINUE
  126. C
  127. C ... ET ON VERIFIE QUE TOUT BAIGNE
  128. C
  129. JG=LECT(/1)
  130. IF(JG.EQ.0)THEN
  131. WRITE(IOIMP,*)'GAuss NEwton: The gradient is empty'
  132. CALL ERREUR(314)
  133. GOTO 9998
  134. ENDIF
  135. C
  136. DO IE1=1,JG
  137. MCHPO1=LECT(IE1)
  138. C
  139. C --> VERIFICATION DU CARACTERE ELEMENTAIRE DU CHPOINT,
  140. C SANS SERIE DE FOURIER
  141. C
  142. SEGACT,MCHPO1
  143. IF(MCHPO1.IPCHP(/1).NE.1)THEN
  144. WRITE(IOIMP,*)'GAuss NEwton: work only with elementary CHPO'
  145. MOTERR(1:8)='CHPOINT '
  146. INTERR(1)=MCHPOI
  147. CALL ERREUR(110)
  148. SEGDES,MCHPO1
  149. GOTO 9998
  150. ENDIF
  151. IF(MCHPO1.IFOPOI.EQ.1)THEN
  152. MOTERR(1:8)='GANE '
  153. CALL ERREUR(333)
  154. SEGDES,MCHPO1
  155. GOTO 9998
  156. ENDIF
  157. C
  158. C --> MEME COMPOSANTES, MEME MAILLAGE
  159. C
  160. MSOUP1=MCHPO1.IPCHP(1)
  161. SEGACT,MSOUP1
  162. C
  163. IF (NC1.NE.MSOUP1.NOHARM(/1))THEN
  164. WRITE(IOIMP,*)
  165. > 'GAuss NEwton: all CHPO should have the same nb of components'
  166. CALL ERREUR(665)
  167. SEGDES,MSOUP1,MCHPO1
  168. GOTO 9998
  169. ENDIF
  170. C
  171. DO IE2=1,NC1
  172. IF(NOCOMP(IE2).NE.MSOUP1.NOCOMP(IE2))THEN
  173. WRITE(IOIMP,*)
  174. > 'GAuss NEwton: all CHPO should have the same components',
  175. > ' in the same order'
  176. CALL ERREUR(488)
  177. SEGDES,MSOUP1,MCHPO1
  178. GOTO 9998
  179. ENDIF
  180. ENDDO
  181. C
  182. IPT1=MSOUP1.IGEOC
  183. SEGACT,IPT1
  184. IF(N1.NE.IPT1.ICOLOR(/1))THEN
  185. WRITE(IOIMP,*)
  186. > 'GAuss NEwton: all CHPO should have the same number',
  187. > ' of points'
  188. MOTERR(1:8) ='MAILLAGE'
  189. MOTERR(9:16)='NOEUD '
  190. CALL ERREUR(403)
  191. SEGDES,MSOUP1,MCHPO1,IPT1
  192. GOTO 9998
  193. ENDIF
  194. C
  195. DO IE2=1,N1
  196. IF(NUM(1,IE2).NE.IPT1.NUM(1,IE2))THEN
  197. WRITE(IOIMP,*)
  198. > 'GAuss NEwton: all CHPO should have the same',
  199. > ' points in the same order'
  200. CALL ERREUR(403)
  201. SEGDES,MSOUP1,MCHPO1,IPT1
  202. GOTO 9998
  203. ENDIF
  204. ENDDO
  205. C
  206. LECT(IE1)=MSOUP1.IPOVAL
  207. SEGDES,MSOUP1,MCHPO1,IPT1
  208. C
  209. ENDDO
  210. C
  211. C ON PREPARE LA MATRICE ET LE SECOND MEMBRE
  212. C
  213. C 1) NOEUDS DU SUPER ELEMENT ET DU CHPO
  214. C
  215. segact mcoord*mod
  216. NBPTT=nbpts
  217. NBPTS=NBPTT+JG
  218. SEGADJ,MCOORD
  219. C
  220. C 2) SUPER ELEMENT
  221. C
  222. NBSOUS=0
  223. NBELEM=1
  224. NBNN=JG
  225. NBREF=0
  226. SEGINI,IPT1
  227. IPT1.ITYPEL=28
  228. DO IE1=1,JG
  229. IPT1.NUM(IE1,1)=NBPTT+ IE1
  230. ENDDO
  231. SEGDES,IPT1
  232. C
  233. C 4) DECRIPTEUR POUR LA RIGIDITE
  234. C
  235. NLIGRP=JG
  236. NLIGRD=JG
  237. SEGINI,DESCR
  238. DO IE1=1,JG
  239. LISINC(IE1)=NOMDD(11)
  240. LISDUA(IE1)=NOMDU(11)
  241. NOELEP(IE1)=IE1
  242. NOELED(IE1)=IE1
  243. ENDDO
  244. SEGDES,DESCR
  245. C
  246. C 5) CONTENU DE LA RIGIDITE
  247. C
  248. nelrig=1
  249. SEGINI,XMATRI
  250. DO IE1=1,NLIGRP
  251. MPOVA1=LECT(IE1)
  252. SEGACT,MPOVA1
  253. DO IE2=1,IE1
  254. MPOVA2=LECT(IE2)
  255. SEGACT,MPOVA2
  256. C
  257. XVAL=0.D0
  258. DO IE3=1,NC1
  259. DO IE4=1,N1
  260. XVAL=XVAL+MPOVA1.VPOCHA(IE4,IE3)*MPOVA2.VPOCHA(IE4,IE3)
  261. ENDDO
  262. ENDDO
  263. RE(IE1,IE2,1)=XVAL
  264. RE(IE2,IE1,1)=XVAL
  265. C
  266. SEGDES,MPOVA2
  267. ENDDO
  268. ENDDO
  269. C
  270. C 6) CORRECTION EVENTUELLE DE RIGIDITE SUR LA DIAGONALE
  271. C
  272. IF(AMOR.NE.0.D0)THEN
  273. DO IE1=1,NLIGRP
  274. RE(IE1,IE1,1)=RE(IE1,IE1,1)+AMOR
  275. ENDDO
  276. ENDIF
  277. SEGDES,XMATRI
  278. C
  279. C 7) IMATRI
  280. C
  281. * NELRIG=1
  282. * SEGINI,IMATRI
  283. * IMATTT(1)=XMATRI
  284. * SEGDES,IMATRI
  285. C
  286. C 8) CHAPEAU MRIGID DE LA RIGIDITE
  287. C
  288. NRIGEL=1
  289. NRIGE=8
  290. SEGINI,MRIGID
  291. MTYMAT='RIGIGITE'
  292. COERIG(1)=1.D0
  293. IRIGEL(1,1)=IPT1
  294. IRIGEL(2,1)=0
  295. IRIGEL(3,1)=DESCR
  296. IRIGEL(4,1)=xMATRI
  297. IRIGEL(5,1)=0
  298. IRIGEL(6,1)=0
  299. IRIGEL(7,1)=0
  300. IRIGEL(8,1)=0
  301. ICHOLE=0
  302. IMGEO1=0
  303. IMGEO2=0
  304. IFORIG=IFOPOI
  305. ISUPEQ=0
  306. SEGDES,MRIGID
  307. C
  308. C 9) SUPPORT DU CHPO
  309. C
  310. NBSOUS=0
  311. NBELEM=JG
  312. NBNN=1
  313. NBREF=0
  314. SEGINI,IPT1
  315. IPT1.ITYPEL=1
  316. DO IE1=1,JG
  317. IPT1.NUM(1,IE1)=NBPTT+ IE1
  318. ENDDO
  319. SEGDES,IPT1
  320. C
  321. C 10) VALEUR DU CHPO
  322. C
  323. N=JG
  324. NC=1
  325. SEGINI,MPOVA2
  326. C
  327. DO IE1=1,N
  328. MPOVA1=LECT(IE1)
  329. SEGACT,MPOVA1
  330. C
  331. XVAL=0.D0
  332. DO IE2=1,NC1
  333. DO IE3=1,N1
  334. XVAL=XVAL-VPOCHA(IE3,IE2)*MPOVA1.VPOCHA(IE3,IE2)
  335. ENDDO
  336. ENDDO
  337. MPOVA2.VPOCHA(IE1,1)=XVAL
  338. C
  339. SEGDES,MPOVA1
  340. ENDDO
  341. C
  342. SEGDES,MPOVA2
  343. C
  344. C 11) MSOUPO
  345. C
  346. SEGINI,MSOUP1
  347. MSOUP1.NOCOMP(1)=NOMDU(11)
  348. MSOUP1.IGEOC=IPT1
  349. MSOUP1.IPOVAL=MPOVA2
  350. SEGDES,MSOUP1
  351. C
  352. C 12) MCHPOI
  353. C
  354. NAT=2
  355. NSOUPO=1
  356. SEGINI,MCHPO1
  357. MCHPO1.JATTRI(1)=JATTRI(1)
  358. MCHPO1.IFOPOI=IFOPOI
  359. MCHPO1.IPCHP(1)=MSOUP1
  360. SEGDES,MCHPO1
  361. C
  362. C RETOUR A GIBIANE
  363. C
  364. CALL ECROBJ('RIGIDITE',MRIGID)
  365. CALL ECROBJ('CHPOINT ',MCHPO1)
  366. C
  367. C ON SORT SI LAMATR
  368. C
  369. IF(LMATR)GOTO 9998
  370. C
  371. C DESACTIVATIONS
  372. C
  373. SEGDES,MELEME,MPOVAL,MSOUPO
  374. SEGSUP,MLENTI
  375. SEGDES,MCHPOI
  376. C
  377. C ON APPELLE RESO
  378. C
  379. CALL RESOU
  380. C
  381. C ... ET ON REGARDE LE RESULTAT
  382. C
  383. CALL LIROBJ('CHPOINT',MCHPOI,1,IRETOU)
  384. IF(IERR.NE.0)RETURN
  385. C
  386. SEGACT,MCHPOI
  387. MSOUPO=IPCHP(1)
  388. SEGACT,MSOUPO
  389. MELEME=IGEOC
  390. MPOVAL=IPOVAL
  391. SEGACT,MELEME,MPOVAL
  392. C
  393. SEGINI,MLENTI,MLENT1
  394. DO IE1=1,JG
  395. LECT(IE1)=NUM(1,IE1)
  396. MLENT1.LECT(IE1)=IE1
  397. ENDDO
  398. CALL GENOR2(LECT(1),MLENT1.LECT(1),JG)
  399. C
  400. SEGINI,MLREEL
  401. DO IE1=1,JG
  402. PROG(IE1)=VPOCHA(MLENT1.LECT(IE1),1)
  403. ENDDO
  404. CALL VECTAB(MLREEL,JG,MTABLE)
  405. C
  406. CALL ECROBJ('TABLE ',MTABLE)
  407. C
  408. C DESTRUCTION
  409. C
  410. SEGSUP,MLENTI,MLENT1,MLREEL
  411. C
  412. SEGSUP,MELEME,MPOVAL,MSOUPO,MCHPOI
  413. C
  414. MCHPOI=MCHPO1
  415. SEGACT,MCHPOI
  416. MSOUPO=IPCHP(1)
  417. SEGACT,MSOUPO
  418. MELEME=IGEOC
  419. MPOVAL=IPOVAL
  420. SEGSUP,MELEME,MPOVAL,MSOUPO,MCHPOI
  421. C
  422. SEGACT,MRIGID
  423. MELEME=IRIGEL(1,1)
  424. DESCR=IRIGEL(3,1)
  425. xMATRI=IRIGEL(4,1)
  426. * SEGACT,IMATRI
  427. * XMATRI=IMATTT(1)
  428. SEGSUP,MELEME,DESCR,XMATRI,MRIGID
  429. C
  430. NBPTS=NBPTT
  431. SEGADJ,MCOORD
  432. C
  433. RETURN
  434. C
  435. C ERREURS
  436. C
  437. 9998 SEGDES,MELEME,MPOVAL,MSOUPO
  438. SEGSUP,MLENTI
  439. 9999 SEGDES,MCHPOI
  440. RETURN
  441. END
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  

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