Télécharger gaunew.eso

Retour à la liste

Numérotation des lignes :

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

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