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