1 : C ACT3 SOURCE GOUNAND 10/02/11 21:15:00 6625 2 : SUBROUTINE ACT3 3 : C-------------------------------------------------------------------- 4 : C ACCELERATION POUR DES CHAMPOINTS 5 : C 6 : C APPELE DANS INCREME 7 : C 8 : C LIT DANS L'ORDRE 3 CHAMPS PRIMAUX 9 : C ET 4 CHAMPS DE DUAUX 10 : C 11 : C REND UN CHAMP PRIMAL EXTRAPOLE QUI MINIMISE LE CHAMP DUAL. 12 : C METHODE : ON ESTIME L'APPLICATION TANGENTE ET ON L'UTILISE 13 : C POUR RESOUDRE LE PROBLEME 14 : C 15 : C-------------------------------------------------------------------- 16 : C 17 : C AOUT 97 : RAJOUTE L'ACCELERATION AVEC 2 OU 3 CHAMPS DE DEPLACEMENTS 18 : C 19 : C CHAT 05/10/2007 : suppression de la limitation à certaines composantes 20 : C 21 : C PM 19/07/2007 : cas à 2 points mettait IDCH0 à 0 au lieu de IDCH1 22 : C 23 : C-------------------------------------------------------------------- 24 : IMPLICIT INTEGER(I-N) 25 : IMPLICIT REAL*8(A-H,O-Z) 26 : -INC CCREEL 27 : -INC SMCHPOI 28 : -INC SMELEME 29 : -INC SMLMOTS 30 : -INC CCOPTIO 31 : C 32 : C FABRICATION DES LISTMOTS 33 : C 34 : JGN=4 35 : JGM=100 36 : SEGINI MLMOTS 37 : JGM=0 38 : 39 : MLMOTF=MLMOTS 40 : * 40 : CALL LIROBJ('CHPOINT',ICHD2,1,IRETOU) 40 : CALL LIROBJ('CHPOINT',ICHD1,1,IRETOU) 40 : CALL LIROBJ('CHPOINT',ICHD0,1,IRETOU) 40 : CALL LIROBJ('CHPOINT',IFO3 ,1,IRETOU) 40 : CALL LIROBJ('CHPOINT',IFO2 ,1,IRETOU) 40 : CALL LIROBJ('CHPOINT',IFO1 ,1,IRETOU) 40 : CALL LIROBJ('CHPOINT',IFO0 ,1,IRETOU) 41 : IF (IERR.NE.0) RETURN 42 : 43 : * Acceleration avec 3 points 44 : IF (ICHD2.EQ.ICHD1) THEN 45 : write(6,*) ' acceleration sur 3 points ' 46 : ICHD2 = 0 47 : IFO3 = 0 48 : ENDIF 49 : 50 : * Acceleration avec 2 points 51 : IF (ICHD1.EQ.ICHD0) THEN 52 : write(6,*) ' acceleration sur 2 points ' 53 : ICHD2 = 0 54 : *PM ICHD0 = 0 55 : ICHD1 = 0 56 : IFO3 = 0 57 : IFO2 = 0 58 : ENDIF 59 : 60 : * 61 : * Recherche des noms de composante des forces 62 : * 63 : MCHPOI=IFO3 64 : IPA=1 65 : 1 CONTINUE 66 : IF (MCHPOI.NE.0) THEN 67 : SEGACT MCHPOI 68 : DO IA=1,IPCHP(/1) 69 : MSOUPO=IPCHP(IA) 70 : SEGACT MSOUPO 71 : DO IB=1,NOCOMP(/2) 72 : DO IC=1,JGM 73 : IF (MOTS(IC).EQ.NOCOMP(IB)) GO TO 2 74 : ENDDO 75 : * on ajoute la nouvelle composante à la liste 76 : JGM=JGM+1 77 : IF (JGM.GT.MOTS(/2)) THEN 78 : SEGADJ MLMOTS 79 : ENDIF 80 : MOTS(JGM)=NOCOMP(IB) 81 : 2 CONTINUE 82 : ENDDO 83 : ENDDO 84 : ENDIF 85 : * On fait de même pour les trois autres champs en entrée 86 : IF (IPA.EQ.1) THEN 87 : IPA=2 88 : MCHPOI=IFO2 89 : GO TO 1 90 : ELSEIF(IPA.EQ.2) THEN 91 : IPA=3 92 : MCHPOI=IFO1 93 : GO TO 1 94 : ELSEIF(IPA.EQ.3) THEN 95 : IPA=4 96 : MCHPOI=IFO0 97 : GO TO 1 98 : ENDIF 99 : SEGADJ MLMOTS 100 : * 101 : * On utilise un procédé d'orthogonalisation de Schmidt 102 : * pour obtenir une base de l'espace engendrée par les quatre forces 103 : * puis on effectue la projection (de coordonnées XL, XM et XN) du 104 : * champ de force nulle sur cet espace 105 : * 106 : IF (IFO3.NE.0) THEN 107 : 108 : * CAS A 4 POINTS : 109 : * ----------------- 110 : * famille des 3 vecteurs variations du champ dual 110 : CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0) 110 : CALL ADCHPO(IFO2,IFO0,IFO20,1D0,-1D0) 110 : CALL ADCHPO(IFO3,IFO0,IFO30,1D0,-1D0) 111 : * 111 : CALL XTX1(IFO10,F12) 111 : CALL XTX1(IFO20,F22) 111 : CALL XTX1(IFO30,F32) 112 : if(f12.LT.XPETIT*1d2) f12=XPETIT*1d2 113 : if(f22.LT.XPETIT*1d2) f22=XPETIT*1d2 114 : if(f32.LT.XPETIT*1d2) f32=XPETIT*1d2 115 : * write (6,*) ' dans act3 f12/f22 ',f12/f22,f12/f32 116 : * 117 : * Détection de flip-flop ? 118 : * (1 et 3 sont trop proches, malgré 2 différent de 1) 119 : if (f12/f22.gt.1.2d0 .and. f12/f32.gt.0.707d0 .and. 120 : > f12/f32.lt.1.414d0) then 121 : * act3 : plan B 121 : CALL ERREUR(858) 122 : XL = (sqrt(f12)+sqrt(f22))/(2*sqrt(f12)) 123 : XM = 0. 124 : XN = 0. 125 : GOTO 1010 126 : endif 127 : * 128 : * Construction d'une famille libre orthogonale 128 : CALL XTY1(IFO10,IFO20,MLMOTF,MLMOTF,F1F2) 128 : CALL XTY1(IFO10,IFO30,MLMOTF,MLMOTF,F1F3) 128 : CALL XTY1(IFO20,IFO30,MLMOTF,MLMOTF,F2F3) 128 : CALL XTY1(IFO0 ,IFO10,MLMOTF,MLMOTF,F0F1) 128 : CALL XTY1(IFO0 ,IFO20,MLMOTF,MLMOTF,F0F2) 128 : CALL XTY1(IFO0 ,IFO30,MLMOTF,MLMOTF,F0F3) 129 : * 130 : A1 = F12 131 : B1 = F1F2 132 : C1 = F1F3 133 : XJ1 = -F0F1 134 : D1 = F1F2 135 : E1 = F22 136 : F1 = F2F3 137 : XK1 = -F0F2 138 : G1 = F1F3 139 : H1 = F2F3 140 : XI1 = F32 141 : XL1 = -F0F3 142 : * write (6,*) ' f12 f22 f32 ',f12,f22,f32 143 : 144 : * Calcul des coordonnées de IFO0 dans la famille libre 145 : * et élimination des vecteurs liés 146 : IF (ABS(XI1) .LT. xpetit) GOTO 2000 147 : * 148 : A2 = A1 - ( G1*C1/XI1) 149 : B2 = B1 - ( H1*C1/XI1) 150 : XJ2 = XJ1 - (XL1*C1/XI1) 151 : D2 = D1 - ( G1*F1/XI1) 152 : E2 = E1 - ( H1*F1/XI1) 153 : XK2 = XK1 - (XL1*F1/XI1) 154 : * 155 : IF (ABS(E2) .LE. ABS(E1)*1D-6) THEN 156 : * write (6,*) ' act3 e2 e1 ' 157 : GOTO 2000 158 : ENDIF 159 : * 160 : A3 = A2 -( D2*B2/E2) 161 : XJ3 = XJ2-(XK2*B2/E2) 162 : * 163 : IF (ABS(A3) .LE. ABS(A1)*1D-6) THEN 164 : * write (6,*) ' act3 a3 a2 ' 165 : GOTO 2000 166 : ENDIF 167 : IF (ABS(A2) .LE. ABS(A1)*1D-6) THEN 168 : * write (6,*) ' act3 a2 a1 ' 169 : GOTO 2000 170 : ENDIF 171 : 172 : XL=XJ3/A3 173 : * 174 : B4 = B2 175 : XJ4 = XJ2 - (A2*XL) 176 : * 177 : IF (ABS(B4) .LE. ABS(B1)*1D-6) THEN 178 : * write (6,*) ' act3 b4 b1 ' 179 : GOTO 2000 180 : ENDIF 181 : * 182 : XM = XJ4/B4 183 : * 184 : C5 = C1 185 : XJ5 = XJ1 - (A1*XL) - (B1*XM) 186 : * 187 : IF (ABS(C5) .LT. xpetit) GOTO 2000 188 : * 189 : XN=XJ5/C5 190 : * 191 : ref=sqrt(a1**2+b1**2+c1**2) 192 : * write (6,*) ' act3 ',ref,XI1,A3,b4,c5 193 : 194 : GOTO 1010 195 : 196 : 2000 CONTINUE 197 : * le systeme est singulier on essaye avec un vecteur de moins 198 : * 199 : * write (6,*) ' act3 ',xi1,e2,e1,a3,a1,b4,b1,c5 200 : IF (ABS(E1) .LT. xpetit) GOTO 2100 201 : * 202 : A2 = A1 - (D1 *B1/E1) 203 : XJ2 = XJ1 - (XK1*B1/E1) 204 : * 205 : IF (ABS(A2) .LE. ABS(A1)*1D-6) THEN 206 : * write (6,*) ' act3 a2 a1 ' 207 : GOTO 2100 208 : ENDIF 209 : * 210 : XL = XJ2/A2 211 : * 212 : B3 = B1 213 : XJ3 = XJ1 - (A1*XL) 214 : * 215 : IF (ABS(B3) .LT. xpetit) GOTO 2100 216 : * 217 : XM = XJ3/B3 218 : * 219 : XN = 0 220 : 221 : * act3 : reduction a %i1 dimensions 222 : INTERR(1)=2 222 : * CALL ERREUR(859) 223 : GOTO 1010 224 : 225 : 2100 CONTINUE 226 : * le systeme est singulier on essaye avec un vecteur de moins 227 : * 228 : IF (ABS(A1) .LT. xpetit ) GOTO 1000 229 : * 230 : XL = XJ1/A1 231 : * 232 : XM = 0 233 : XN = 0 234 : 235 : * act3 : reduction a %i1 dimensions 236 : INTERR(1)=1 237 : 237 : CALL ERREUR(859) 238 : GOTO 1010 239 : 240 : 1000 CONTINUE 241 : * act3 : pas d accélération 241 : call erreur(860) 242 : XL = 0 243 : XM = 0 244 : XN = 0 245 : 246 : 1010 CONTINUE 247 : * 248 : * WRITE (6,*) ' f12*f22/(f1f2*2) ',f12*f22/(f1f2**2) 249 : * WRITE (6,*) ' f12*f32/(f1f3*2) ',f12*f32/(f1f3**2) 250 : * WRITE (6,*) ' f22*f32/(f2f3*2) ',f22*f32/(f2f3**2) 251 : * write (6,*) ' xi1,e2,a3,b4,c5 ',xi1,e2,a3,b4,c5 252 : * write (6,*) ' coefficients acc ',-xl,-xm,-xn 253 : * xl=sign(min(abs(xl),1d1),xl) 254 : * xm=sign(min(abs(xm),1d1),xm) 255 : * xn=sign(min(abs(xn),1d1),xn) 256 : * write (6,*) ' coefficients cor ',-xl,-xm,-xn 257 : 258 : * Vérification qu'on ne retombe pas sur l'un des duals entrés 259 : * Sinon, on réessaye en tenant compte de moins de vecteurs 260 : if (abs(xl-1).lt.1e-1 .and. abs(xm).lt.1e-1 261 : > .and. abs(xn).lt.1e-1) goto 1000 262 : if (abs(xm-1).lt.1e-1 .and. abs(xn).lt.1e-1 263 : > .and. abs(xl).lt.1e-1) goto 2100 264 : if (abs(xn-1).lt.1e-1 .and. abs(xl).lt.1e-1 265 : > .and. abs(xm).lt.1e-1) goto 2000 266 : * 266 : CALL ADCHPO(ICHD0,ICHD1,ITEMP4,XL+XM+XN,XM+XN) 266 : CALL ADCHPO(ICHD2,ITEMP4,ITEMP5,XN,1D0) 267 : * 268 : * ajustement de la nature du champ 269 : * on prend la meme nature que ICHDO 270 : MCHPOI = ICHD0 271 : MCHPO1 = ITEMP5 272 : SEGACT MCHPOI,MCHPO1 273 : NAT = JATTRI(/1) 274 : NSOUPO = MCHPO1.IPCHP(/1) 275 : SEGADJ MCHPO1 276 : DO 1020 I=1,NAT 277 : MCHPO1.JATTRI(I)=JATTRI(I) 278 : 1020 CONTINUE 279 : SEGDES MCHPOI,MCHPO1 280 : * 281 : * Affectation 282 : * 282 : CALL ECROBJ('CHPOINT',ITEMP5) 283 : * 284 : * Quelques destructions 285 : * 286 : MLMOTS=MLMOTF 287 : SEGSUP MLMOTS 287 : CALL DTCHPO(IFO10) 287 : CALL DTCHPO(IFO20) 287 : CALL DTCHPO(IFO30) 287 : CALL DTCHPO(ITEMP4) 288 : RETURN 289 : ENDIF 290 : 291 : * CAS A 3 POINTS : 292 : * ----------------- 293 : * famille des 2 vecteurs variations du champ dual 294 : IF (IFO2.NE.0.AND.IFO3.EQ.0) THEN 294 : CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0) 294 : CALL ADCHPO(IFO2,IFO0,IFO20,1D0,-1D0) 295 : * 295 : CALL XTX1(IFO10,F12) 295 : CALL XTX1(IFO20,F22) 296 : * 297 : * Construction d'une famille libre orthogonale 297 : CALL XTY1(IFO10,IFO20,MLMOTF,MLMOTF,F1F2) 297 : CALL XTY1(IFO0 ,IFO10,MLMOTF,MLMOTF,F0F1) 297 : CALL XTY1(IFO0 ,IFO20,MLMOTF,MLMOTF,F0F2) 298 : * 299 : A1 = F12 300 : B1 = F1F2 301 : XJ1 = -F0F1 302 : D1 = F1F2 303 : E1 = F22 304 : XK1 = -F0F2 305 : * 306 : * Calcul des coordonnées de IFO0 dans la famille libre 307 : * et élimination des vecteurs liés 308 : IF (ABS(E1) .LT. xpetit ) GOTO 2200 309 : * 310 : A2 = A1 - (D1 *B1/E1) 311 : XJ2 = XJ1 - (XK1*B1/E1) 312 : * 313 : IF (ABS(A2) .LE. ABS(A1)*1D-6) THEN 314 : * write (6,*) ' act3 a2 a1 ' 315 : GOTO 2200 316 : ENDIF 317 : * 318 : XL = XJ2/A2 319 : * 320 : B3 = B1 321 : XJ3 = XJ1 - (A1*XL) 322 : * 323 : IF (ABS(B3) .LT. xpetit) GOTO 2200 324 : * 325 : XM = XJ3/B3 326 : * 327 : XN = 0 328 : * 329 : * write (6,*) ' act3 xl xm ',xl,xm 330 : GOTO 1210 331 : 332 : 2200 CONTINUE 333 : * le systeme est singulier on essaye avec un vecteur de moins 334 : * 335 : IF (ABS(A1) .LT. xpetit ) GOTO 2300 336 : * 337 : XL = XJ1/A1 338 : * 339 : XM = 0 340 : XN = 0 341 : * 342 : * act3 : reduction a %i1 dimensions 343 : INTERR(1)=1 343 : * CALL ERREUR(859) 344 : GOTO 1210 345 : 346 : 2300 CONTINUE 347 : * act3 : pas d accélération 347 : call erreur(860) 348 : XL = 0 349 : XM = 0 350 : XN = 0 351 : 352 : 1210 CONTINUE 353 : * Vérification qu'on ne retombe pas sur l'un des duals entrés 354 : * Sinon, on réessaye en tenant compte de moins de vecteurs 355 : if (abs(xl-1).lt.1e-1.and.abs(xm).lt.1e-1.and. 356 : > abs(xn).lt.1e-1) goto 2300 357 : if (abs(xm-1).lt.1e-1.and.abs(xn).lt.1e-1.and. 358 : > abs(xl).lt.1e-1) goto 2200 359 : * if (abs(xn-1).lt.1e-1.and.abs(xl).lt.1e-1.and. 360 : * > abs(xm).lt.1e-1) goto 2300 361 : * 361 : CALL ADCHPO(ICHD0,ICHD1,ITEMP5,XL+XM,XM) 362 : * 363 : * ajustement de la nature du champ 364 : * on prend la meme nature que ICHDO 365 : * 366 : MCHPOI = ICHD0 367 : MCHPO1 = ITEMP5 368 : SEGACT MCHPOI,MCHPO1 369 : NAT = JATTRI(/1) 370 : NSOUPO = MCHPO1.IPCHP(/1) 371 : SEGADJ MCHPO1 372 : DO 1220 I=1,NAT 373 : MCHPO1.JATTRI(I)=JATTRI(I) 374 : 1220 CONTINUE 375 : SEGDES MCHPOI,MCHPO1 376 : * 377 : * Affectation 378 : * 378 : CALL ECROBJ('CHPOINT',ITEMP5) 379 : * 380 : * Quelques destructions 381 : * 382 : MLMOTS=MLMOTF 383 : SEGSUP MLMOTS 383 : CALL DTCHPO(IFO10) 383 : CALL DTCHPO(IFO20) 384 : RETURN 385 : ENDIF 386 : 387 : * CAS A 2 POINTS : 388 : * ----------------- 389 : IF (IFO2.EQ.0) THEN 389 : CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0) 390 : * 390 : CALL XTX1(IFO10,F12) 391 : * 392 : * Un seul vecteur = famille libre d'emblée 392 : CALL XTY1(IFO0,IFO10,MLMOTF,MLMOTF,F0F1) 393 : * 394 : A1 = F12 395 : XJ1 = -F0F1 396 : * 397 : IF (ABS(A1) .LT. xpetit ) GOTO 1100 398 : * 399 : XL = XJ1/A1 400 : * 401 : XM = 0 402 : XN = 0 403 : GOTO 1110 404 : 405 : 1100 CONTINUE 406 : * act3 : pas d accélération 406 : call erreur(860) 407 : XL = 0 408 : XM = 0 409 : XN = 0 410 : 411 : 1110 CONTINUE 412 : if (abs(xl-1).lt.1e-1) goto 1100 412 : CALL MUCHPO(ICHD0,XL,ITEMP5,1) 413 : * write (6,*) ' act3 xl vaut ',xl 414 : * 415 : * ajustement de la nature du champ 416 : * on prend la meme nature que ICHDO 417 : * 418 : MCHPOI = ICHD0 419 : MCHPO1 = ITEMP5 420 : SEGACT MCHPOI,MCHPO1 421 : NAT = JATTRI(/1) 422 : NSOUPO = MCHPO1.IPCHP(/1) 423 : SEGADJ MCHPO1 424 : DO 1120 I=1,NAT 425 : MCHPO1.JATTRI(I)=JATTRI(I) 426 : 1120 CONTINUE 427 : SEGDES MCHPOI,MCHPO1 428 : * 429 : * Affectation 430 : * 430 : CALL ECROBJ('CHPOINT',ITEMP5) 431 : * 432 : * QUELQUES DESTRUCTIONS 433 : * 434 : MLMOTS=MLMOTF 435 : SEGSUP MLMOTS 435 : CALL DTCHPO(IFO10) 436 : RETURN 437 : ENDIF 438 : END 439 : 440 : 441 : 442 :
© Cast3M 2003 - All rights reserved.
Disclaimer