1 : C INCLU4 SOURCE PASCAL 09/01/22 21:15:35 6263 2 : SUBROUTINE INCLU4(IPT1,IPT2,IPEX,XCRIT) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : -INC CCREEL 6 : -INC CCOPTIO 7 : -INC SMELEME 8 : -INC SMCOORD 9 : * 10 : * IL EST SUPPOSER QUE LE IPT2 ne contiennent quer des TET4 et que ipt1 11 : * soit des poi1, 12 : * 13 : * IPEX(I)=1 veut dire que le noeud I est interne à un element de ipt2 14 : * 15 : * POUR DECIDER SI UN POINT EST A L'INTERIEUR D'UN ELEMENT ON CALCULE 16 : * LES COORDONNEES BARYCENTRIQUES DU POINT ET IL FAUT QU'ELLES SOIENT 17 : * TOUTES POSITIVES OU QUE CELLES QUI SOIENT NEGATIVES SOIENT D'UN 18 : * ORDRE DE GRANDEUR TRES INFERIEUR AUX AUTRES ( 0.0001 FOIS ). 19 : * 20 : SEGMENT ISEG1 21 : REAL*8 XLIM(2,NBEL),YLIM(2,NBEL),ZLIM(2,NBEL) 22 : ENDSEGMENT 23 : SEGMENT ISEG3 24 : INTEGER NIZO(NZO+1) 25 : ENDSEGMENT 26 : SEGMENT ISEG4 27 : INTEGER NUMZO(NZO) 28 : ENDSEGMENT 29 : SEGMENT ISEG5 30 : INTEGER NNMEL(ILON),IDEJ(NZO) 31 : ENDSEGMENT 32 : SEGMENT ISEG6 33 : REAL*8 AM(4,4),P(4),AL(4),A(4),B(4),C(4),D(4) 34 : REAL*8 XA(3,4),XPU(3) 35 : ENDSEGMENT 36 : SEGMENT IPEX(XCOOR(/1)/(IDIM+1)) 37 : * 38 : * ON CALCULE LA TAILLE MAXI D'UN ELEMENT DANS TOUTES LES DIRECTIONS 39 : * AFIN DE CREER UN ZONAGE DE L'ESPACE. EN MEME TEMPS ON CALCULE 40 : * LA DIMENSION HORS TOUT DU MAILLAGE ET ON COMMENCE A CREER 41 : * UNE NUMEROTATION LOCALE DES POINTS DU MAILLAGE IPT 42 : * 43 : IDIM1=IDIM+1 44 : XPREC=-XCRIT 45 : MELEME= IPT2 46 : SEGACT MELEME 47 : IF(ITYPEL.NE.23) THEN 47 : CALL ERREUR(5) 48 : RETURN 49 : ENDIF 50 : NBEL = NUM(/2) 51 : NBNN=NUM(/1) 52 : * WRITE(6,FMT='('' NBEL NBNN '',2I6)') NBEL,NBNN 53 : SEGINI ISEG1 54 : ILOC=0 55 : XZO=0.D0 56 : YZO=0.D0 57 : ZZO=0.D0 58 : XZA=XGRAND 59 : YZA=XGRAND 60 : ZZA=XGRAND 61 : XTOMI=XGRAND 62 : XTOMA=-XGRAND 63 : YTOMI=XGRAND 64 : YTOMA=-XGRAND 65 : ZTOMI=XGRAND 66 : ZTOMA=-XGRAND 67 : DO 1 I1=1,NBEL 68 : XMI=XGRAND 69 : YMI=XGRAND 70 : ZMI=XGRAND 71 : YMA=-XGRAND 72 : XMA=-XGRAND 73 : ZMA=-XGRAND 74 : DO 2 I2 = 1,NBNN 75 : IB=NUM(I2,I1) 76 : IA=(IB-1)*IDIM1 77 : IF(XCOOR(IA+1).LT.XMI) XMI=XCOOR(IA+1) 78 : IF(XCOOR(IA+1).GT.XMA) XMA=XCOOR(IA+1) 79 : IF(XCOOR(IA+2).LT.YMI) YMI=XCOOR(IA+2) 80 : IF(XCOOR(IA+2).GT.YMA) YMA=XCOOR(IA+2) 81 : * IF( IDIM.EQ.3 ) THEN 82 : IF(XCOOR(IA+3).LT.ZMI) ZMI=XCOOR(IA+3) 83 : IF(XCOOR(IA+3).GT.ZMA) ZMA=XCOOR(IA+3) 84 : * ENDIF 85 : 2 CONTINUE 86 : XLIM(1,I1)=XMI 87 : XLIM(2,I1)=XMA 88 : YLIM(1,I1)=YMI 89 : YLIM(2,I1)=YMA 90 : XZO=MAX (XZO,XMA-XMI) 91 : YZO=MAX (YZO,YMA-YMI) 92 : XZA=MIN(XZA,XMA-XMI) 93 : YZA=MIN(YZA,YMA-YMI) 94 : IF(XMI.LT.XTOMI) XTOMI=XMI 95 : IF(XMA.GT.XTOMA) XTOMA=XMA 96 : IF(YMI.LT.YTOMI) YTOMI=YMI 97 : IF(YMA.GT.YTOMA) YTOMA=YMA 98 : * IF(IDIM.EQ.3) THEN 99 : ZLIM(1,I1)=ZMI 100 : ZLIM(2,I1)=ZMA 101 : ZZO=MAX(ZZO,ZMA-ZMI) 102 : ZZA=MIN(ZZA,ZMA-ZMI) 103 : IF(ZMI.LT.ZTOMI) ZTOMI=ZMI 104 : IF(ZMA.GT.ZTOMA) ZTOMA=ZMA 105 : * ENDIF 106 : 1 CONTINUE 107 : * WRITE(6,FMT='(''XZO YZO '',4E12.5)') XZO,YZO 108 : * WRITE(6,FMT='(''XZA YZA '',4E12.5)') XZA,YZA 109 : * WRITE(6,FMT='(''XTOMI XTOMA '',4E12.5)') XTOMI,XTOMA 110 : * WRITE(6,FMT='(''YTOMI YTOMA '',4E12.5)') YTOMI,YTOMA 111 : XPR=MIN(XZO*1.D-2,(XTOMA-XTOMI)/2.D+4) 112 : YPR=MIN(YZO*1.D-2,(YTOMA-YTOMI)/2.D+4) 113 : C XZO=XZO*1.1 114 : C YZO=YZO*1.1 115 : XZA=XZA*0.97 116 : YZA=YZA*0.97 117 : XTOMI= XTOMI - (XTOMA-XTOMI)/1.D+4 118 : XTOMA= XTOMA + (XTOMA-XTOMI)/1.D+4 119 : YTOMI= YTOMI - (YTOMA-YTOMI)/1.D+4 120 : YTOMA= YTOMA + (YTOMA-YTOMI)/1.D+4 121 : C XZO=MIN ( XZO, XTOMA-XTOMI) 122 : C YZO=MIN ( YZO, YTOMA-YTOMI) 123 : XZA=MIN ( XZA, XTOMA-XTOMI) 124 : YZA=MIN ( YZA, YTOMA-YTOMI) 125 : NXZO=INT((XTOMA-XTOMI)/XZA) + 1 126 : NYZO=INT((YTOMA-YTOMI)/YZA) + 1 127 : XZO=XZA 128 : YZO=YZA 129 : NZZO=1 130 : * WRITE(6,FMT='('' NXZO NYZO'',2I7)') NXZO,NYZO 131 : * IF(IDIM.EQ.3) THEN 132 : ZPR=MIN(ZZO*1.D-2,(ZTOMA-ZTOMI)/2.D+4) 133 : C ZZO=ZZO*1.1 134 : ZZA=ZZA*0.97 135 : ZTOMI= ZTOMI - (ZTOMA-ZTOMI)/1.D+4 136 : ZTOMA= ZTOMA + (ZTOMA-ZTOMI)/1.D+4 137 : C ZZO=MIN ( ZZO, ZTOMA-ZTOMI) 138 : ZZA=MIN ( ZZA, ZTOMA-ZTOMI) 139 : NZZO=INT((ZTOMA-ZTOMI)/ZZA)+ 1 140 : ZZO=ZZA 141 : * WRITE(6,FMT='('' zz0,zzA,ztomi,ztoma'',4e12.5)') 142 : * $ xzo,xza,ztomi,ztoma 143 : * ENDIF 144 : * WRITE(6,FMT='('' XTOMI XTOMA YTOMI YTOMA '',4E12.5 )') 145 : * $ XTOMI, XTOMA, YTOMI ,YTOMA 146 : NXDEP=MIN(NXZO,10) 147 : NYDEP=MIN(NYZO,10) 148 : * IF(IDIM.EQ.2) THEN 149 : * IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN 150 : * XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/90 151 : * NXZO=MAX(INT(NXZO/XY),NXDEP) 152 : * NYZO=MAX(INT(NYZO/XY),NYDEP) 153 : * IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN 154 : * XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/60 155 : * NXZO=MAX(INT(NXZO/XY),NXDEP) 156 : * NYZO=MAX(INT(NYZO/XY),NYDEP) 157 : * ENDIF 158 : * XZO=(XTOMA-XTOMI)/NXZO 159 : * YZO=(YTOMA-YTOMI)/NYZO 160 : * NXZO=(XTOMA-XTOMI)/XZO +1 161 : * NYZO=(YTOMA-YTOMI)/YZO +1 162 : * ENDIF 163 : C WRITE(6,FMT='('' XZO NXZO YZO NYZO '' , E12.5,I5,E12.5,I5)') 164 : C $ XZO ,NXZO, YZO, NYZO 165 : * ELSE 166 : NZDEP=MIN(NZZO,10) 167 : * WRITE(6,FMT='('' XZO NXZO YZO NYZO ZZO NZZO'' , E12.5,I7,/, 168 : * $ E12.5,I7,E12.5,I7)') 169 : C $ XZO ,NXZO, YZO, NYZO,ZZO,NZZO 170 : IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NXZO NYZO NZZO '' 171 : $,4I7) ') NXZO,NYZO,NZZO 172 : IF(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO).GT.25000.) THEN 173 : XYZ =(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO))**0.3333/25. 174 : NXZO=MAX(INT(FLOAT(NXZO)/XYZ),NXDEP) 175 : NYZO=MAX(INT(FLOAT(NYZO)/XYZ),NYDEP) 176 : NZZO=MAX(INT(FLOAT(NZZO)/XYZ),NZDEP) 177 : IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NXZO NYZO NZZO '' 178 : $,4I7) ') NXZO,NYZO,NZZO 179 : IF(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO).GT.20000.) THEN 180 : XYZ =(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO))**0.3333/25. 181 : NXZO=MAX(INT(FLOAT(NXZO)/XYZ),NXDEP) 182 : NYZO=MAX(INT(FLOAT(NYZO)/XYZ),NYDEP) 183 : NZZO=MAX(INT(FLOAT(NZZO)/XYZ),NZDEP) 184 : IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NXZO NYZO NZZO '' 185 : $,4I7) ') NXZO,NYZO,NZZO 186 : IF(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO).GT.20000.) THEN 187 : XYZ =(FLOAT(NXZO)*FLOAT(NYZO)*FLOAT(NZZO))**0.3333/25. 188 : NXZO=MAX(INT(FLOAT(NXZO)/XYZ),NXDEP) 189 : NYZO=MAX(INT(FLOAT(NYZO)/XYZ),NYDEP) 190 : NZZO=MAX(INT(FLOAT(NZZO)/XYZ),NZDEP) 191 : ENDIF 192 : ENDIF 193 : XZO=(XTOMA-XTOMI)/FLOAT(NXZO) 194 : YZO=(YTOMA-YTOMI)/FLOAT(NYZO) 195 : ZZO=(ZTOMA-ZTOMI)/FLOAT(NZZO) 196 : NXZO=INT((XTOMA-XTOMI)/XZO)+1 197 : NYZO=INT((YTOMA-YTOMI)/YZO)+1 198 : NZZO=INT((ZTOMA-ZTOMI)/ZZO)+1 199 : ENDIF 200 : * ENDIF 201 : * 202 : * ON VEUT CONSTRUIRE LA LISTE DES ELEMENTS TOUCHANT UNE ZONE 203 : * POUR CELA ON COMMENCE PAR COMPTER COMBIEN D'ELEMENT TOUCHENT 204 : * CHAQUE ZONE ET EN MEME TEMPS ON STOCKE LES ZONES TOUCHEES 205 : * PAR CHAQUE ELEMENT ET LEUR NOMBRE 206 : * 207 : 208 : NZO=NXZO*NYZO*NZZO 209 : IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NZO NXZO NYZO NZZO '' 210 : $,4I7) ') NZO,NXZO,NYZO,NZZO 211 : NXYZO=NXZO*NYZO 212 : * IDI=4 213 : * IF(IDIM.EQ.3) THEN 214 : * IDI=8 215 : * ENDIF 216 : SEGINI ISEG3 217 : SEGINI ISEG4 218 : DO 3 I1=1,NBEL 219 : NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1 220 : NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1 221 : NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1 222 : NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1 223 : * IF(IDIM.EQ.3) THEN 224 : NIZ1Z=INT((ZLIM(1,I1)-ZTOMI-ZPR)/ZZO) +1 225 : NIZ2Z=INT((ZLIM(2,I1)-ZTOMI+ZPR)/ZZO) +1 226 : DO 200 L3=NIZ1Z,NIZ2Z 227 : DO 200 L1=NIZ1Y,NIZ2Y 228 : DO 200 L2=NIZ1X,NIZ2X 229 : NIZA = L2 + ( L1-1) * NXZO + ( L3-1)*NXYZO 230 : NUMZO(NIZA) = NUMZO(NIZA) +1 231 : 200 CONTINUE 232 : * ELSE 233 : * DO 201 L1=NIZ1Y,NIZ2Y 234 : * DO 201 L2=NIZ1X,NIZ2X 235 : * NIZA = L2 + ( L1-1) * NXZO 236 : * NUMZO(NIZA) = NUMZO(NIZA) +1 237 : * 201 CONTINUE 238 : * ENDIF 239 : 3 CONTINUE 240 : * 241 : * CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES 242 : * ELEMENTS CONCERNEES PAR UNE ZONE 243 : * 244 : ILON=0 245 : NIZO(1)=1 246 : DO 202 L1=1,NZO 247 : NIZO(L1+1)=NIZO(L1)+NUMZO(L1) 248 : ILON=ILON+ NUMZO(L1) 249 : 202 CONTINUE 250 : * WRITE(6,FMT='('' ILON '',I5)') ILON 251 : * WRITE(6,109) (KKK,NUMZO(KKK),(NELZO(KI,KKK),KI=1,4),KKK=1,NBEL) 252 : * 109 FORMAT(I6,I5,4I5) 253 : * WRITE(6,110)( NIZO(KI),KI=1,NZO+1) 254 : 110 FORMAT(16I5) 255 : SEGINI ISEG5 256 : DO 5 I1=1,NBEL 257 : NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1 258 : NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1 259 : NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1 260 : NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1 261 : * IF(IDIM.EQ.3) THEN 262 : NIZ1Z=INT((ZLIM(1,I1)-ZTOMI-ZPR)/ZZO) +1 263 : NIZ2Z=INT((ZLIM(2,I1)-ZTOMI+ZPR)/ZZO) +1 264 : DO 205 L3=NIZ1Z,NIZ2Z 265 : DO 205 L1=NIZ1Y,NIZ2Y 266 : DO 205 L2=NIZ1X,NIZ2X 267 : NIZA = L2 + ( L1-1) * NXZO + ( L3-1)*NXYZO 268 : IAD=NIZO(NIZA)+IDEJ(NIZA) 269 : NNMEL(IAD)=I1 270 : IDEJ(NIZA)=IDEJ(NIZA)+1 271 : 205 CONTINUE 272 : * ELSE 273 : * DO 203 L1=NIZ1Y,NIZ2Y 274 : * DO 203 L2=NIZ1X,NIZ2X 275 : * NIZA = L2 + ( L1-1) * NXZO 276 : * IAD=NIZO(NIZA)+IDEJ(NIZA) 277 : * NNMEL(IAD)=I1 278 : * IDEJ(NIZA)=IDEJ(NIZA)+1 279 : * 203 CONTINUE 280 : * ENDIF 281 : 5 CONTINUE 282 : * 283 : * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT POUR CHAQUE 284 : * POINT DE L'OBJETR MAILLAGE IPT1, ON COMMENCE PAR LE METTRE SOUS 285 : * FORME D'ELEMENTS DE TYPE POI1 286 : * 287 : SEGSUP ISEG1,ISEG4 288 : SEGACT IPT1 289 : IF(IPT1.ITYPEL.NE.1) THEN 289 : CALL CHANGE(IPT1,1) 290 : ENDIF 291 : SEGINI IPEX 292 : SEGINI ISEG6 293 : C WRITE(6,FMT='('' AVANT BOUCLE 10'')') 294 : DO 10 I=1,IPT1.NUM(/2) 295 : IP=IPT1.NUM(1,I) 296 : XPU(1)=XCOOR((IP-1)*IDIM1+1) 297 : XPU(2)=XCOOR((IP-1)*IDIM1+2) 298 : XPU(3)=XCOOR((IP-1)*IDIM1+3) 299 : * write(6,fmt='('' point x y '',i5,2e12.5)') iP,XPU(1),xpu(2) 300 : IF(XPU(1).LT.XTOMI.OR.XPU(1).GT.XTOMA) GO TO 10 301 : IF(XPU(2).LT.YTOMI.OR.XPU(2).GT.YTOMA) GO TO 10 302 : * IF(IDIM.EQ.3) THEN 303 : IF(XPU(3).LT.ZTOMI.OR.XPU(3).GT.ZTOMA) GO TO 10 304 : * ENDIF 305 : INDZO=INT((XPU(1)-XTOMI)/XZO)+ 1 +INT((XPU(2)-YTOMI)/YZO)*NXZO 306 : # +INT((XPU(3)-ZTOMI)/ZZO)*NXZO*NYZO 307 : * IF(IDIM.EQ.3) INDZO=INDZO+INT((XPU(3)-ZTOMI)/ZZO)*NXZO*NYZO 308 : IDEB=NIZO(INDZO) 309 : IFIN=NIZO(INDZO+1)-1 310 : * write(6,fmt='('' ideb ifin'',2i5)') ideb,ifin 311 : IF(IDEB.GT.IFIN) GO TO 10 312 : IEL=0 313 : DO 11 KK=IDEB,IFIN 314 : * 315 : * ON CALCULE LES COORDONNEES BARYCENTRIQUES ( AU PLUS 4 ) 316 : * 317 : K=NNMEL(KK) 318 : * write(6,fmt='('' k '',i5)') k 319 : J1=NUM(1,K) 320 : J2=NUM(2,K) 321 : J3=NUM(3,K) 322 : J1IDIM=(J1-1)*IDIM1 323 : J2IDIM=(J2-1)*IDIM1 324 : J3IDIM=(J3-1)*IDIM1 325 : * IF(IDIM.EQ.3) THEN 326 : J4=NUM(4,K) 327 : J4IDIM=(J4-1)*IDIM1 328 : * ENDIF 329 : P(1)=1.D0 330 : DO 12 K1=1,NBNN 331 : AM(1,K1)=1.D0 332 : 12 CONTINUE 333 : DO 13 K1=1,IDIM 334 : P(K1+1)=XPU(K1) 335 : AM(K1+1,1)=XCOOR(J1IDIM+K1) 336 : AM(K1+1,2)=XCOOR(J2IDIM+K1) 337 : AM(K1+1,3)=XCOOR(J3IDIM+K1) 338 : * IF(IDIM.EQ.3) THEN 339 : AM(K1+1,4)=XCOOR(J4IDIM+K1) 340 : * ENDIF 341 : 13 CONTINUE 342 : * IF(IDIM.EQ.2) THEN 343 : * X1=AM(2,1) 344 : * X2=AM(2,2) 345 : * X3=AM(2,3) 346 : * Y1=AM(3,1) 347 : * Y2=AM(3,2) 348 : * Y3=AM(3,3) 349 : * X=P(2) 350 : * Y=P(3) 351 : * DETAM=X1*Y2+X2*Y3+X3*Y1-Y1*X2-Y2*X3-Y3*X1 352 : * A(1)=X2*Y3-X3*Y2 353 : * A(2)=X3*Y1-X1*Y3 354 : * A(3)=X1*Y2-X2*Y1 355 : * B(1)=Y2-Y3 356 : * B(2)=Y3-Y1 357 : * B(3)=Y1-Y2 358 : * C(1)=X3-X2 359 : * C(2)=X1-X3 360 : * C(3)=X2-X1 361 : * DO 14 IK=1,NBNN 362 : * AL(IK)=(A(IK)+B(IK)*X+C(IK)*Y)/DETAM 363 : * 14 CONTINUE 364 : * AL(4)=1.D0 365 : * ELSE 366 : X1=AM(2,1) 367 : X2=AM(2,2) 368 : X3=AM(2,3) 369 : X4=AM(2,4) 370 : Y1=AM(3,1) 371 : Y2=AM(3,2) 372 : Y3=AM(3,3) 373 : Y4=AM(3,4) 374 : Z1=AM(4,1) 375 : Z2=AM(4,2) 376 : Z3=AM(4,3) 377 : Z4=AM(4,4) 378 : X=P(2) 379 : Y=P(3) 380 : Z=P(4) 381 : DETAM=X2*Y3*Z4+X3*Y4*Z2+X4*Y2*Z3-X4*Y3*Z2-X2*Y4*Z3-X3*Y2*Z4-X1*Y3* 382 : 1Z4-X3*Y4*Z1-X4*Y1*Z3+X4*Y3*Z1+X3*Y1*Z4+X1*Y4*Z3+X1*Y2*Z4+X4*Y1*Z2+ 383 : 2X2*Y4*Z1-X4*Y2*Z1-X2*Y1*Z4-X1*Y4*Z2-X1*Y2*Z3-X3*Y1*Z2-X2*Y3*Z1+X3* 384 : 3Y2*Z1+X2*Y1*Z3+X1*Y3*Z2 385 : A(1)=X2*Y3*Z4+X3*Y4*Z2+X4*Y2*Z3-X4*Y3*Z2-X2*Y4*Z3-X3*Y2*Z4 386 : A(2)=X4*Y3*Z1+X3*Y1*Z4+X1*Y4*Z3-X1*Y3*Z4-X3*Y4*Z1-X4*Y1*Z3 387 : A(3)=X1*Y2*Z4+X4*Y1*Z2+X2*Y4*Z1-X4*Y2*Z1-X2*Y1*Z4-X1*Y4*Z2 388 : A(4)=X3*Y2*Z1+X2*Y1*Z3+X1*Y3*Z2-X1*Y2*Z3-X3*Y1*Z2-X2*Y3*Z1 389 : B(1)=Y4*Z3-Y3*Z4+Y2*Z4-Y4*Z2+Y3*Z2-Y2*Z3 390 : B(2)=Y3*Z4-Y4*Z3+Y4*Z1-Y1*Z4+Y1*Z3-Y3*Z1 391 : B(3)=Y4*Z2-Y2*Z4+Y1*Z4-Y4*Z1+Y2*Z1-Y1*Z2 392 : B(4)=Y2*Z3-Y3*Z2+Y3*Z1-Y1*Z3+Y1*Z2-Y2*Z1 393 : C(1)=X3*Z4-X4*Z3+X4*Z2-X2*Z4+X2*Z3-X3*Z2 394 : C(2)=X4*Z3-X3*Z4+X1*Z4-X4*Z1+X3*Z1-X1*Z3 395 : C(3)=X2*Z4-X4*Z2+X4*Z1-X1*Z4+X1*Z2-X2*Z1 396 : C(4)=X3*Z2-X2*Z3+X1*Z3-X3*Z1+X2*Z1-X1*Z2 397 : D(1)=X4*Y3-X3*Y4+X2*Y4-X4*Y2+X3*Y2-X2*Y3 398 : D(2)=X3*Y4-X4*Y3+X4*Y1-X1*Y4+X1*Y3-X3*Y1 399 : D(3)=X4*Y2-X2*Y4+X1*Y4-X4*Y1+X2*Y1-X1*Y2 400 : D(4)=X2*Y3-X3*Y2+X3*Y1-X1*Y3+X1*Y2-X2*Y1 401 : DO 15 II=1,NBNN 402 : AL(II)=(A(II)+B(II)*X+C(II)*Y+D(II)*Z)/DETAM 403 : 15 CONTINUE 404 : * ENDIF 405 : * write(6,fmt='(''K al '',I5,4e12.3)')K, 406 : * $ al(1),al(2),al(3),al(4) 407 : IF( AL(1).GT.XPREC.AND.AL(2).GT.XPREC.AND.AL(3).GT.XPREC. 408 : $ AND.AL(4).GT.XPREC) THEN 409 : * 410 : * LE POINT EST INTERNE A L'ELEMENT 411 : * 412 : * IEL=IEL+1 413 : IPEX(IP)=1 414 : * WRITE(6,*) IP 415 : GO TO 10 416 : ENDIF 417 : 11 CONTINUE 418 : 10 CONTINUE 419 : SEGSUP ISEG5,ISEG6,ISEG3 420 : SEGDES IPT1,MELEME 421 : RETURN 422 : END 423 : 424 : 425 : 426 : 427 : 428 :
© Cast3M 2003 - All rights reserved.
Disclaimer