1 : C LIGNE SOURCE CHAT 07/10/30 21:16:34 5942 2 : 3 : C FABRIQUE DES LIGNES CONSTITUES D ELEM A 2 OU 3 NOEUDS 4 : C 5 : C ILIGGI indique quel est le type de la ligne support : 6 : C 1 = SEGMENT DE DROITE (operateur DROI) 7 : C 2 = 8 : C 3 = ARC DE CERCLE (operateur CERC) 9 : C 4 = 10 : C 5 = <-- indisponible 11 : C 6 = <-- indisponible 12 : C 7 = ARC DE PARABOLE (operateur PARA) 13 : C 8 = <-- indisponible : ELLIPSE (operateur ELLI) 14 : C 9 = PORTION DE CUBIQUE (operateur CUBP) 15 : C 10 = PORTION DE CUBIQUE (operateur CUBT) 16 : C 11 = ARC DE CERCLE (operateur CERC3) 17 : C 12 = 18 : C 19 : C 09/2003 : 20 : C Modifications suite a la mise en place du cas IDIM=1. Ajout de 21 : C tests de coherence. Indisponibilite des operateurs PARA,CERC,CER3, 22 : C CUBP,CUBT dans le cas IDIM=1. Ajout de commentaires. 23 : 24 : SUBROUTINE LIGNE(ILIGGI,IFONC,DEN1OX,DEN2OX,INBR) 25 : 26 : IMPLICIT INTEGER(I-N) 27 : IMPLICIT REAL*8 (A-H,O-Z) 28 : 29 : -INC CCREEL 30 : -INC CCOPTIO 31 : -INC CCGEOME 32 : -INC SMELEME 33 : -INC SMCOORD 34 : 35 : PARAMETER (NITER=5) 36 : logical ltelq 37 : SEGMENT TABPAR(NBELEM) 38 : SEGMENT / STRAV / (DIST2(500),XPOIN(500),YPOIN(500), 39 : . ZPOIN(500),DIST(500),U(NSTRAV,500)) 40 : 41 : C :::::: AJOUT PIFFARD : 16/01/86 POUR CUBIQUES ET CER3 :::::::: 42 : DIMENSION NDP(12) 43 : DIMENSION TU(4),CU(4,4),DLONGT(NITER+1),DMOY(NITER) 44 : CHARACTER*(4) ITTEMP*8 45 : CHARACTER*4 MCLE(2) 46 : CHARACTER*4 MUNIF(1) 47 : DATA NDP / 2,2,3,3,0,0,3,0,4,4,3,3 / 48 : DATA MCLE / 'DINI','DFIN' / 49 : DATA MUNIF / 'UNIF' / 50 : C :::::: AJOUT PIFFARD : 16/01/86 POUR CUBIQUES :::::::: 51 : DATA CU / 1.D0, 0.D0, -5.5D0, -1.0D0, 52 : . 0.D0, 0.D0, 9.0D0, 4.5D0, 53 : . 0.D0, 0.D0, -4.5D0, -9.0D0, 54 : . 0.D0, 1.D0, 1.0D0, 5.5D0 / 55 : 56 : ILIG=ILIGGI 57 : C En DIMEnsion 1, seul l'operateur DROIte peut etre utilise 58 : IF ((IDIM.EQ.1).AND.(ILIG.NE.1)) THEN 59 : INTERR(1)=IDIM 59 : CALL ERREUR(709) 60 : RETURN 61 : ENDIF 62 : C Quelques verifications 62 : IF (ILCOUR.EQ.0) CALL ERREUR(16) 63 : INOMB=NDP(ILIG) 63 : IF (INOMB.EQ.0) CALL ERREUR(19) 64 : ITYPELT=KDEGRE(ILCOUR) 65 : C Pas de verif. sur ITYPELT (2 ou 3 uniquement) 65 : IF (ITYPELT.EQ.0) CALL ERREUR(16) 66 : C Arc de CUBIQUE (CUBP et CUBT) uniquement en SEG2 67 : IF (((ILIG.EQ.9).OR.(ILIG.EQ.10)).AND.(ITYPELT.EQ.3)) 67 : . CALL ERREUR(16) 67 : IF ((ILIG.EQ.12).AND.(ITYPELT.NE.3)) CALL ERREUR(16) 68 : NBNN=NBNNE(ITYPELT) 68 : IF (NBNN.EQ.0) CALL ERREUR(16) 69 : IF (IERR.NE.0) RETURN 70 : 71 : C Initialisations dans le cas de CUBP et CUBT 72 : IF (ILIG.EQ.9.OR.ILIG.EQ.10) THEN 73 : DO i=1,NITER 74 : DLONGT(i)=0.D0 75 : DMOY(i)=0.D0 76 : ENDDO 77 : DLONGT(NITER+1)=0.D0 78 : ENDIF 79 : 80 : DEN1=DEN1OX 81 : DEN2=DEN2OX 82 : 83 : C Lecture des donnees facultatives (decoupage, densites...) 84 : IF (IFONC.EQ.0) GOTO 400 85 : C Y a-t-il un facteur de decoupage ? 86 : IF (ILIG.NE.10) THEN 87 : INBR=0 87 : CALL MESLIR(-172) 87 : CALL LIRENT(INBR,0,IRETOU) 88 : ELSE 89 : C Si CUBIque a vecteurs TGTS, il FAUT un facteur de decoupage. 89 : CALL MESLIR(-172) 89 : CALL LIRENT(INBR,1,IRETOU) 90 : IF (IERR.NE.0) RETURN 91 : ENDIF 92 : C Y a t-il des densites imposees ? 93 : IMPOI=0 94 : IMPOF=0 94 : 413 CALL MESLIR(-171) 94 : CALL LIRMOT(MCLE,2,IRETOU,0) 95 : IF (IRETOU.EQ.1) THEN 95 : CALL MESLIR(-170) 95 : CALL LIRREE(DEN1,1,IRETOU) 96 : IF (IERR.NE.0) RETURN 97 : IMPOI=1 98 : GOTO 413 99 : ELSEIF (IRETOU.EQ.2) THEN 99 : CALL MESLIR(-169) 99 : CALL LIRREE(DEN2,1,IRETOU) 100 : IF (IERR.NE.0) RETURN 101 : IMPOF=1 102 : GOTO 413 103 : ENDIF 104 : 400 CONTINUE 105 : C Y-a-t-il le mot UNIF pour les cubiques ? 106 : IF (ILIG.EQ.9) THEN 107 : IUNIF=0 107 : CALL LIRMOT(MUNIF,1,IUNIF,0) 108 : ENDIF 109 : C INOMB : nombre de "points" a lire pour definir la ligne 110 : C Lecture du premier POINT (ou LIGNE) 111 : C Quel est le type de la premiere donnee (POINT OU LIGNE) ? 112 : IPOIN3=0 113 : ITP1=0 114 : ITP2=0 115 : ITTEMP=' ' 115 : CALL MESLIR(-168) 115 : CALL QUETYP(ITTEMP,1,IRETOU) 116 : IF (IERR.NE.0) RETURN 116 : CALL MESLIR(-168) 116 : CALL LIROBJ(ITTEMP,IPT1,1,IRETOU) 117 : IF (ITTEMP.EQ.'MAILLAGE') THEN 118 : SEGACT IPT1 118 : CALL EXTRPO(IPT1,1,IPOIN1) 118 : CALL EXTRPO(IPT1,2,IPOIN3) 119 : ITP1=1 120 : ELSE 121 : IF (ITTEMP.NE.'POINT ') THEN 122 : MOTERR(1:8)=ITTEMP 122 : CALL ERREUR(39) 123 : RETURN 124 : ENDIF 125 : IPOIN1=IPT1 126 : IPOIN3=IPT1 127 : ENDIF 128 : C Lecture des troisieme et quatrieme POINTs si necessaire 129 : IF (INOMB.EQ.2) GOTO 503 129 : CALL MESLIR(-167) 129 : CALL LIROBJ('POINT ',IPCEN,1,IRETOU) 130 : IF (IERR.NE.0) RETURN 131 : IF (INOMB.EQ.3) GOTO 503 131 : CALL LIROBJ('POINT ',IPCEN2,1,IRETOU) 132 : IF (IERR.NE.0) RETURN 132 : 503 CALL MESLIR(-166) 133 : C Lecture eventuelle du second POINT ou LIGNE 133 : CALL LIROBJ('MAILLAGE',IP2,0,IRETOU) 134 : IF (IRETOU.NE.0) THEN 135 : IPT6=IP2 136 : SEGACT IPT6 136 : CALL EXTRPO(IP2,2,IPOIN2) 137 : ITP2=1 138 : ELSE 138 : CALL MESLIR(-166) 138 : CALL LIROBJ('POINT ',IPOIN2,0,IRETOU) 139 : IF (IRETOU.EQ.0) THEN 140 : IPOIN2=IPOIN3 140 : IF (IPOIN1.EQ.IPOIN2) CALL ERREUR(20) 141 : ENDIF 141 : IF (IPOIN1.EQ.IPOIN2) CALL ERREUR(432) 142 : ENDIF 143 : 144 : C Affectation des densites suivant les donnees 145 : IBP1=(IDIM+1)*IPOIN1 146 : IBP2=(IDIM+1)*IPOIN2 147 : IF (IFONC.NE.0) THEN 148 : IF (IMPOI.EQ.0) DEN1=XCOOR(IBP1) 149 : IF (IMPOF.EQ.0) DEN2=XCOOR(IBP2) 150 : ENDIF 150 : IF ((INBR.LE.0).AND.((DEN1*DEN2).LE.0.D0)) CALL ERREUR(17) 151 : IF (IERR.NE.0) RETURN 152 : DENI=DEN1 153 : DECA=DEN2-DEN1 154 : C ----- AJOUT PIFFARD : 23/01/86 POUR CALCUL D'UNE TOL FONCT. DES DENS 155 : C ----- AJOUT GOUNAND : 30/10/97 : ELLE NE SERT NULLE PART ! 156 : C TOL=(DEN1+DEN2)/20.D0 157 : 158 : C Branchement suivant le type de LIGNE 159 : C ------------- 160 : GOTO (11,11,13,13,11,11,13,17,18,18,19,13),ILIG 161 : 162 : C Segment de DROITE (DROI) 163 : C -------------------------- 164 : 11 XDIS=XCOOR(IBP2-IDIM)-XCOOR(IBP1-IDIM) 165 : IF (IDIM.GE.2) THEN 166 : YDIS=XCOOR(IBP2-IDIM+1)-XCOOR(IBP1-IDIM+1) 167 : ZDIS=0.D0 168 : IF (IDIM.GE.3) ZDIS=XCOOR(IBP2-IDIM+2)-XCOOR(IBP1-IDIM+2) 169 : DIS=SQRT(XDIS*XDIS+YDIS*YDIS+ZDIS*ZDIS) 170 : ELSE 171 : DIS=ABS(XDIS) 172 : ENDIF 173 : C Modif. pour permettre de faire des montagnes 174 : C Si decoupage impose, on accepte points confondus 175 : IF ((INBR.NE.0).AND.(DIS.EQ.0.D0)) DIS=1.D0 175 : IF (DIS.EQ.0.D0) CALL ERREUR(21) 176 : IF (IERR.NE.0) RETURN 177 : DLONG=DIS 178 : DEN1=DEN1/DIS 179 : DEN2=DEN2/DIS 180 : XIN=XCOOR(IBP1-IDIM) 181 : IF (IDIM.EQ.1) THEN 182 : YIN=0.D0 183 : ZIN=0.D0 184 : ELSE IF (IDIM.EQ.2) THEN 185 : YIN=XCOOR(IBP1-IDIM+1) 186 : ZIN=0. 187 : ELSE IF (IDIM.EQ.3) THEN 188 : YIN=XCOOR(IBP1-IDIM+1) 189 : ZIN=XCOOR(IBP1-IDIM+2) 190 : ENDIF 191 : GOTO 200 192 : C Arc de CERCLE (CERC) ou de PARABOLE (PARA) 193 : C -------------------------------------------- 194 : 13 IBPC=(IDIM+1)*IPCEN 195 : XCEN=XCOOR(IBPC-IDIM) 196 : YCEN=XCOOR(IBPC-IDIM+1) 197 : ZCEN=0.D0 198 : IF (IDIM.GE.3) ZCEN=XCOOR(IBPC-IDIM+2) 199 : IF (ILIG.NE.7) THEN 200 : C Cas de l'arc de cercle 201 : C On ne recalcule pas le centre des paraboles 202 : C ----- 23/01/86 AJOUT PIFFARD : RECALCUL DU CENTRE ----- 203 : XVP=XCOOR(IBP2-IDIM)-XCOOR(IBP1-IDIM) 204 : XMP=(XCOOR(IBP2-IDIM)+XCOOR(IBP1-IDIM))*0.5D0 205 : YVP=XCOOR(IBP2-IDIM+1)-XCOOR(IBP1-IDIM+1) 206 : YMP=(XCOOR(IBP2-IDIM+1)+XCOOR(IBP1-IDIM+1))*0.5D0 207 : ZVP=0.D0 208 : ZMP=0.D0 209 : IF (IDIM.GE.3) THEN 210 : ZVP=XCOOR(IBP2-IDIM+2)-XCOOR(IBP1-IDIM+2) 211 : ZMP=(XCOOR(IBP2-IDIM+2)+XCOOR(IBP1-IDIM+2))*0.5D0 212 : ENDIF 213 : SNOR=SQRT(XVP*XVP+YVP*YVP+ZVP*ZVP) 214 : XVP=XVP/SNOR 215 : YVP=YVP/SNOR 216 : ZVP=ZVP/SNOR 217 : RO=(XMP-XCEN)*XVP+(YMP-YCEN)*YVP+(ZMP-ZCEN)*ZVP 218 : XCEN=XCEN+RO*XVP 219 : YCEN=YCEN+RO*YVP 220 : ZCEN=ZCEN+RO*ZVP 221 : C ----- FIN DE L'AJOUT --------------------------------- 222 : ENDIF 223 : XV1=XCOOR(IBP1-IDIM)-XCEN 224 : XV2=XCOOR(IBP2-IDIM)-XCEN 225 : YV1=XCOOR(IBP1-IDIM+1)-YCEN 226 : YV2=XCOOR(IBP2-IDIM+1)-YCEN 227 : ZV1=0.D0 228 : ZV2=0.D0 229 : IF (IDIM.GE.3) THEN 230 : ZV1=XCOOR(IBP1-IDIM+2)-ZCEN 231 : ZV2=XCOOR(IBP2-IDIM+2)-ZCEN 232 : ENDIF 233 : IF (ILIG.NE.7) THEN 234 : C Arc de CERCLE 235 : RAY1=SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) 236 : RAY2=SQRT(XV2*XV2+YV2*YV2+ZV2*ZV2) 236 : IF (RAY1/RAY2.LT.0.99D0.OR.RAY1/RAY2.GT.1.01D0) CALL ERREUR(21) 237 : IF (IERR.NE.0) RETURN 238 : RAY=SQRT(RAY1*RAY2) 239 : SCAL=XV1*XV2+YV1*YV2+ZV1*ZV2 240 : XVECT=YV1*ZV2-YV2*ZV1 241 : YVECT=XV2*ZV1-XV1*ZV2 242 : ZVECT=XV1*YV2-XV2*YV1 243 : VECT=SQRT(XVECT*XVECT+YVECT*YVECT+ZVECT*ZVECT) 244 : ANG=ATAN2(VECT,SCAL) 245 : IF (IIMPI.EQ.1) WRITE (IOIMP,9143) ANG 246 : 9143 FORMAT(' ANGLE DE L ARC EN RADIAN ',G12.5) 247 : IF (ABS(ANG).GT.XPI) THEN 247 : CALL ERREUR(313) 248 : RETURN 249 : ENDIF 250 : SINA=SIN(ANG) 251 : DLONG=ANG*RAY 252 : DEN1=DEN1/DLONG 253 : DEN2=DEN2/DLONG 254 : XV3=(YVECT*ZV1-YV1*ZVECT)/(RAY*RAY*SINA) 255 : YV3=(XV1*ZVECT-XVECT*ZV1)/(RAY*RAY*SINA) 256 : ZV3=(XVECT*YV1-XV1*YVECT)/(RAY*RAY*SINA) 257 : ELSE 258 : C Arc de PARABOLE 259 : XIN=XCOOR(IBP1-IDIM) 260 : YIN=XCOOR(IBP1-IDIM+1) 261 : XDIS=0.5D0*(XCOOR(IBP2-IDIM)-XIN) 262 : YDIS=0.5D0*(XCOOR(IBP2-IDIM+1)-YIN) 263 : IF (IDIM.GE.3) THEN 264 : ZIN=XCOOR(IBP1-IDIM+2) 265 : ZDIS=0.5D0*(XCOOR(IBP2-IDIM+2)-ZIN) 266 : ELSE 267 : ZIN=0.D0 268 : ZDIS=0.D0 269 : ENDIF 270 : DIS=2*SQRT(XDIS*XDIS+YDIS*YDIS+ZDIS*ZDIS) 271 : XIN=XIN+XDIS 272 : YIN=YIN+YDIS 273 : ZIN=ZIN+ZDIS 274 : DLONG=0.5D0*(DIS+SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1) 275 : . +SQRT(XV2*XV2+YV2*YV2+ZV2*ZV2)) 276 : DEN1=DEN1/DLONG 277 : DEN2=DEN2/DLONG 278 : XV=-0.25D0*(XV1+XV2) 279 : YV=-0.25D0*(YV1+YV2) 280 : ZV=-0.25D0*(ZV1+ZV2) 281 : ENDIF 282 : GOTO 200 283 : C NON disponible 284 : C ---------------- 285 : 17 GOTO 200 286 : C Arc de CUBIQUE (CUBP,CUBT) 287 : C ---------------- 288 : C ::::::::: AJOUT PIFFARD 16/01/86 : POUR CUBIQUES ::::: 289 : C ...... RECUPERATION DES 4 PTS DE BASE SI CUBP ..... 290 : C ...... DES 2 PTS DE BASE + TGTES SI CUBT ..... 291 : 18 XINI=XCOOR(IBP1-IDIM) 292 : YINI=XCOOR(IBP1-IDIM+1) 293 : ZINI=0.D0 294 : IF (IDIM.GE.3) ZINI=XCOOR(IBP1-IDIM+2) 295 : IBP3=(IDIM+1)*IPCEN 296 : XINT1=XCOOR(IBP3-IDIM) 297 : YINT1=XCOOR(IBP3-IDIM+1) 298 : ZINT1=0.D0 299 : IF (IDIM.GE.3) ZINT1=XCOOR(IBP3-IDIM+2) 300 : IBP4=(IDIM+1)*IPCEN2 301 : XINT2=XCOOR(IBP4-IDIM) 302 : YINT2=XCOOR(IBP4-IDIM+1) 303 : ZINT2=0.D0 304 : IF (IDIM.GE.3) ZINT2=XCOOR(IBP4-IDIM+2) 305 : XFIN=XCOOR(IBP2-IDIM) 306 : YFIN=XCOOR(IBP2-IDIM+1) 307 : ZFIN=0.D0 308 : IF (IDIM.GE.3) ZFIN=XCOOR(IBP2-IDIM+2) 309 : C Calcul de la distance entre les 2 ou 4 pts pour determiner 310 : C le nombre de segments a creer si il n'est pas donne 311 : ZDIS=0.D0 312 : IF (ILIG.EQ.9) THEN 313 : XDIS=XFIN-XINT2-XINT1-XINI 314 : YDIS=YFIN-YINT2-YINT1-YINI 315 : IF (IDIM.GE.3) ZDIS=ZFIN-ZINT2-ZINT1-ZINI 316 : ELSE 317 : XDIS=XFIN-XINI 318 : YDIS=YFIN-YINI 319 : IF (IDIM.GE.3) ZDIS=ZFIN-ZINI 320 : ENDIF 321 : DIS=SQRT(XDIS*XDIS+YDIS*YDIS+ZDIS*ZDIS) 321 : IF (DIS.EQ.0) CALL ERREUR(21) 322 : IF (IERR.NE.0) RETURN 323 : DLONG=DIS 324 : DEN1=DEN1/DIS 325 : DEN2=DEN2/DIS 326 : GOTO 200 327 : C Arc de CERCLE (CERC3) 328 : C ----------------------- 329 : C Ajout PIFFARD : JANV 86 POUR CERCLE PASSANT PAR 3 PTS 330 : C Remanie par VERPEAUX on se branche sur le cercle normal 331 : 19 IBPC=(IDIM+1)*IPCEN 332 : XOM1=XCOOR(IBP1-IDIM) 333 : YOM1=XCOOR(IBP1-IDIM+1) 334 : ZOM1=0.D0 335 : IF (IDIM.GE.3) ZOM1=XCOOR(IBP1-IDIM+2) 336 : XOM2=XCOOR(IBPC-IDIM) 337 : YOM2=XCOOR(IBPC-IDIM+1) 338 : ZOM2=0.D0 339 : IF (IDIM.GE.3) ZOM2=XCOOR(IBPC-IDIM+2) 340 : XOM3=XCOOR(IBP2-IDIM) 341 : YOM3=XCOOR(IBP2-IDIM+1) 342 : ZOM3=0.D0 343 : IF (IDIM.GE.3) ZOM3=XCOOR(IBP2-IDIM+2) 344 : XM1=XOM2-XOM1 345 : YM1=YOM2-YOM1 346 : ZM1=ZOM2-ZOM1 347 : SNOR=SQRT(XM1*XM1+YM1*YM1+ZM1*ZM1) 348 : XM1=XM1/SNOR 349 : YM1=YM1/SNOR 350 : ZM1=ZM1/SNOR 351 : XM2=XOM3-XOM2 352 : YM2=YOM3-YOM2 353 : ZM2=ZOM3-ZOM2 354 : SNOR=SQRT(XM2*XM2+YM2*YM2+ZM2*ZM2) 355 : XM2=XM2/SNOR 356 : YM2=YM2/SNOR 357 : ZM2=ZM2/SNOR 358 : XVECT=YM1*ZM2-YM2*ZM1 359 : YVECT=ZM1*XM2-ZM2*XM1 360 : ZVECT=XM1*YM2-XM2*YM1 361 : VECT=SQRT(XVECT*XVECT+YVECT*YVECT+ZVECT*ZVECT) 362 : XVECT=XVECT/VECT 363 : YVECT=YVECT/VECT 364 : ZVECT=ZVECT/VECT 365 : XV1=YM1*ZVECT-YVECT*ZM1 366 : YV1=ZM1*XVECT-ZVECT*XM1 367 : ZV1=XM1*YVECT-XVECT*YM1 368 : XM3=XOM3-XOM1 369 : YM3=YOM3-YOM1 370 : ZM3=ZOM3-ZOM1 371 : RO=(XM2*XM3+YM2*YM3+ZM2*ZM3)/(XM2*XV1+YM2*YV1+ZM2*ZV1) 372 : XCEN=0.5D0*(XOM1+XOM2+RO*XV1) 373 : YCEN=0.5D0*(YOM1+YOM2+RO*YV1) 374 : ZCEN=0.5D0*(ZOM1+ZOM2+RO*ZV1) 375 : IF (IIMPI.EQ.1) WRITE(IOIMP,6032) XCEN,YCEN,ZCEN 376 : 6032 FORMAT(2X,'CENTRE DU CERCLE =',3G13.5) 377 : RAY=SQRT((XOM1-XCEN)**2+(YOM1-YCEN)**2+(ZOM1-ZCEN)**2) 378 : IF (IIMPI.EQ.1) WRITE(IOIMP,6033) RAY 379 : 6033 FORMAT(2X,'RAYON DU CERCLE =',G13.5) 380 : XV1=XCOOR(IBP1-IDIM)-XCEN 381 : XV2=XCOOR(IBP2-IDIM)-XCEN 382 : YV1=XCOOR(IBP1-IDIM+1)-YCEN 383 : YV2=XCOOR(IBP2-IDIM+1)-YCEN 384 : ZV1=0.D0 385 : ZV2=0.D0 386 : IF (IDIM.GE.3) THEN 387 : ZV1=XCOOR(IBP1-IDIM+2)-ZCEN 388 : ZV2=XCOOR(IBP2-IDIM+2)-ZCEN 389 : ENDIF 390 : SCAL=XV1*XV2+YV1*YV2+ZV1*ZV2 391 : XVE=YV1*ZV2-YV2*ZV1 392 : YVE=XV2*ZV1-XV1*ZV2 393 : ZVE=XV1*YV2-XV2*YV1 394 : VE=XVE*XVECT+YVE*YVECT+ZVE*ZVECT 395 : ANG=ATAN2(VE,SCAL) 396 : IF (ANG.LT.0.D0) ANG=ANG+(2.D0*XPI) 397 : IF (IIMPI.EQ.1) WRITE(IOIMP,9143) ANG 398 : SINA=SIN(ANG) 399 : DLONG=ANG*RAY 400 : DEN1=DEN1/DLONG 401 : DEN2=DEN2/DLONG 402 : XV3=YVECT*ZV1-YV1*ZVECT 403 : YV3=XV1*ZVECT-XVECT*ZV1 404 : ZV3=XVECT*YV1-XV1*YVECT 405 : ILIG=3 406 : GOTO 200 407 : 408 : C Decoupage de la ligne (base sur l'abscisse curviligne) 409 : C ----------------------- 409 : 200 CALL DECOUP(INBR,DEN1,DEN2,APROG,NBELEM,DENI,DECA,DLONG) 410 : IF (IIMPI.EQ.1) WRITE(IOIMP,1000) NBELEM,APROG 411 : 1000 FORMAT (/,' ELEMENT ',I6,' RAISON ',G12.5) 411 : IF ((ILIG.EQ.9).AND.(NBELEM.LT.3)) CALL ERREUR(21) 412 : IF (IERR.NE.0) RETURN 413 : 414 : NBSOUS=0 415 : NBREF=0 416 : SEGINI MELEME 417 : ITYPEL=ITYPELT 418 : NX=NBELEM-1 419 : SEGINI TABPAR 420 : C Initialisation segment de travail (arc de CUBIQUE) 421 : IF ((ILIG.EQ.9).OR.(ILIG.EQ.10)) THEN 422 : NSTRAV=NITER+1 423 : SEGINI STRAV 424 : ENDIF 425 : 426 : C Creation des elements (objet MAILLAGE MELEME) 427 : C Calcul des densites de chacun des points 428 : DIN=DEN1 429 : IPOO=XCOOR(/1)/(IDIM+1) 430 : NUM(1,1)=IPOIN1 431 : IF (NX.NE.0) THEN 432 : DO i=1,NX 433 : DIN=DIN*APROG 434 : TABPAR(i)=DIN 435 : IPOO=IPOO+1 436 : NUM(2,i)=IPOO 437 : IF (ITYPEL.EQ.3) THEN 438 : IPOO=IPOO+1 439 : NUM(3,i)=IPOO 440 : C NUM(3,i)=IPOO+1 441 : ENDIF 442 : NUM(1,i+1)=IPOO 443 : ENDDO 444 : ENDIF 445 : DIN=DIN*APROG 446 : IF (ITYPEL.EQ.3) THEN 447 : TABPAR(NBELEM)=DIN 448 : IPOO=IPOO+1 449 : NUM(2,NBELEM)=IPOO 450 : ENDIF 451 : NUM(ITYPEL,NBELEM)=IPOIN2 452 : 453 : C Creation des POINTs 454 : C --------------------- 455 : C :::::: AJOUT PIFFARD POUR CUBIQUES :::::: 456 : IF ((ILIG.EQ.9).OR.(ILIG.EQ.10)) THEN 457 : NX=NX+1 458 : IDEP=2 459 : XPOIN(1)=XINI 460 : YPOIN(1)=YINI 461 : ZPOIN(1)=ZINI 462 : DIST(1)=0.D0 463 : ELSE 464 : IDEP=1 465 : ENDIF 466 : DPAR=0 467 : IADR=XCOOR(/1)/(IDIM+1) 468 : IF (NX.EQ.0) GOTO 300 469 : NBPTS=IADR+(NX-IDEP+1)*(ITYPEL-1) 470 : SEGADJ MCOORD 471 : DO 301 i=IDEP,NX 472 : DIN=TABPAR(i) 473 : DPAR=DPAR+DIN 474 : C Creation des POINTS milieux (generation de SEG3) 475 : IF (ITYPEL.EQ.2) GOTO 302 476 : DPAR1=DPAR-DIN*0.5D0 477 : GOTO (1001,1001,1003,1003,1001,1001,1007,1008,1009,1009, 478 : . 1010,1003),ILIG 479 : C Segment de DROITE 480 : 1001 XCOOR(IADR*(IDIM+1)+1)=XIN+DPAR1*XDIS 481 : IF (IDIM.GE.2) THEN 482 : XCOOR(IADR*(IDIM+1)+2)=YIN+DPAR1*YDIS 483 : IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZIN+DPAR1*ZDIS 484 : ENDIF 485 : GOTO 1020 486 : C Arc de CERCLE 487 : 1003 PARA=ANG*DPAR1 488 : COSA=COS(PARA) 489 : SINA=SIN(PARA) 490 : IF (ILIG.EQ.12) THEN 491 : PAR2=ANG*DIN/2.D0 492 : XDD=2.D0-COS(PAR2)-(TAN(PAR2)/RAY/(2.D0*RAY/SIN(PAR2)-1.D0)) 493 : ELSE 494 : XDD=1.D0 495 : ENDIF 496 : XCOOR(IADR*(IDIM+1)+1)=XCEN+(COSA*XV1+SINA*XV3)*XDD 497 : XCOOR(IADR*(IDIM+1)+2)=YCEN+(COSA*YV1+SINA*YV3)*XDD 498 : IF (IDIM.GE.3) THEN 499 : XCOOR(IADR*(IDIM+1)+3)=ZCEN+(COSA*ZV1+SINA*ZV3)*XDD 500 : ENDIF 501 : GOTO 1020 502 : C Arc de PARABOLE 503 : 1007 PARA=2*(DPAR1-0.5D0) 504 : XCOOR(IADR*(IDIM+1)+1)=XIN+PARA*XDIS+(1-PARA**2)*XV 505 : XCOOR(IADR*(IDIM+1)+2)=YIN+PARA*YDIS+(1-PARA**2)*YV 506 : IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZIN+PARA*ZDIS+ 507 : . (1-PARA*PARA)*ZV 508 : GOTO 1020 509 : C NON disponible 510 : 1008 GOTO 1020 511 : C Arc de CUBIQUE (CUBT et CUBP) 511 : 1009 CALL ERREUR(16) 512 : RETURN 513 : C Erreur anormale car ILG=11 --> ILIG=3 513 : 1010 CALL ERREUR(5) 514 : RETURN 515 : C Densite du point milieu cree 516 : 1020 XCOOR((IADR+1)*(IDIM+1))=DENI+DECA*DPAR1 517 : IADR=IADR+1 518 : 519 : C Creation des points finaux de chaque SEG2 ou SEG3 520 : C :::::: AJOUT PIFFARD POUR CUBIQUES (GENERATION DE SEG2 UNIQUEMENT) ::: 521 : 302 GOTO (1101,1101,1103,1103,1101,1101,1107,1108,1109,1109, 522 : . 1110,1103),ILIG 523 : C Segment de DROITE 524 : 1101 XCOOR(IADR*(IDIM+1)+1)=XIN+DPAR*XDIS 525 : IF (IDIM.GE.2) THEN 526 : XCOOR(IADR*(IDIM+1)+2)=YIN+DPAR*YDIS 527 : IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZIN+DPAR*ZDIS 528 : ENDIF 529 : GOTO 1120 530 : C Arc de CERCLE 531 : 1103 PARA=ANG*DPAR 532 : COSA=COS(PARA) 533 : SINA=SIN(PARA) 534 : XCOOR(IADR*(IDIM+1)+1)=XCEN+COSA*XV1+SINA*XV3 535 : XCOOR(IADR*(IDIM+1)+2)=YCEN+COSA*YV1+SINA*YV3 536 : IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZCEN+COSA*ZV1+SINA*ZV3 537 : GOTO 1120 538 : C Arc de PARABOLE 539 : 1107 PARA=2*(DPAR-0.5D0) 540 : XCOOR(IADR*(IDIM+1)+1)=XIN+PARA*XDIS+(1-PARA*PARA)*XV 541 : XCOOR(IADR*(IDIM+1)+2)=YIN+PARA*YDIS+(1-PARA*PARA)*YV 542 : IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZIN+PARA*ZDIS+ 543 : . (1-PARA*PARA)*ZV 544 : GOTO 1120 545 : 1108 GOTO 1120 546 : C Arc de CUBIQUE 547 : C ::::: AJOUT PIFFARD : CUBIQUES . 16/01/86 :::::: 548 : 1109 UP=(i-1.D0)/NX 549 : U(1,i)=UP 550 : F1=2*UP**3-3*UP**2+1 551 : F2=-2*UP**3+3*UP**2 552 : F3=UP**3-2*UP**2+UP 553 : F4=UP**3-UP**2 554 : IF (ILIG.EQ.9) THEN 555 : C ...... CUBIQUES PASSANT PAR 4 POINTS ....... 556 : DO J=1,4 557 : TU(J)=+CU(J,1)*XINI+CU(J,2)*XINT1+CU(J,3)*XINT2+CU(J,4) 558 : . *XFIN 559 : ENDDO 560 : XPOIN(i)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4) 561 : DO J=1,4 562 : TU(J)=+CU(J,1)*YINI+CU(J,2)*YINT1+CU(J,3)*YINT2+CU(J,4) 563 : . *YFIN 564 : ENDDO 565 : YPOIN(i)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4) 566 : IF (IDIM.GE.3) THEN 567 : DO J=1,4 568 : TU(J)=+CU(J,1)*ZINI+CU(J,2)*ZINT1+CU(J,3)*ZINT2+CU(J,4) 569 : . *ZFIN 570 : ENDDO 571 : ZPOIN(i)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4) 572 : ENDIF 573 : ELSE 574 : C ...... CUBIQUES PASSANT PAR 2 POINTS + TGTES ....... 575 : XPOIN(i)=F1*XINI+F2*XFIN+F3*XINT1+F4*XINT2 576 : YPOIN(i)=F1*YINI+F2*YFIN+F3*YINT1+F4*YINT2 577 : IF (IDIM.GE.2) ZPOIN(i)=F1*ZINI+F2*ZFIN+F3*ZINT1+F4*ZINT2 578 : ENDIF 579 : C ... SI IUNIF = 0 PAS DE CORRECTION SUR LES PARAM DE LA CUBIQUE .... 580 : IF (IUNIF.EQ.1) THEN 581 : DIST(I)=DIST(I-1)+SQRT((XPOIN(I)-XPOIN(I-1))**2 582 : . +(YPOIN(I)-YPOIN(I-1))**2+(ZPOIN(I)-ZPOIN(I-1))**2) 583 : IF (IIMPI.EQ.1) WRITE(IOIMP,6012) DIST(I) 584 : 6012 FORMAT(2X,'DISTANCE CUMULEE ENTRE PTS =',G13.5) 585 : DLONGT(1)=DIST(I) 586 : GOTO 301 587 : ELSE 588 : XCOOR(IADR*(IDIM+1)+1)=XPOIN(I) 589 : XCOOR(IADR*(IDIM+1)+2)=YPOIN(I) 590 : IF(IDIM.GT.2) XCOOR(IADR*(IDIM+1)+3)=ZPOIN(I) 591 : ENDIF 592 : GOTO 1120 593 : C Erreur anormale car ILIG=11 --> ILIG=3 593 : 1110 CALL ERREUR(5) 594 : RETURN 595 : C Densite du point cree 596 : 1120 XCOOR((IADR+1)*(IDIM+1))=DENI+DECA*DPAR 597 : IADR=IADR+1 598 : 301 CONTINUE 599 : C Fin de la boucle de creation des points 600 : NBPTS=IADR 601 : SEGADJ MCOORD 602 : 603 : 300 DIN=TABPAR(NBELEM) 604 : DPAR=DPAR+DIN 605 : C Rajout du point milieu du dernier SEG3 (non effectue dans la boucle) 606 : IF (ITYPEL.EQ.3) THEN 607 : NBPTS=IADR+1 608 : SEGADJ MCOORD 609 : DPAR1=DPAR-DIN*0.5D0 610 : GOTO (1201,1201,1203,1203,1201,1201,1207,1208, 611 : . 1201,1201,1201,1203),ILIG 612 : 1201 XCOOR(IADR*(IDIM+1)+1)=XIN+DPAR1*XDIS 613 : IF (IDIM.GE.2) XCOOR(IADR*(IDIM+1)+2)=YIN+DPAR1*YDIS 614 : IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZIN+DPAR1*ZDIS 615 : GOTO 1220 616 : 1203 PARA=ANG*DPAR1 617 : COSA=COS(PARA) 618 : SINA=SIN(PARA) 619 : IF (ILIG.EQ.12) THEN 620 : PAR2=ANG*DIN/2.D0 621 : XDD=2.D0-COS(PAR2)-(TAN(PAR2)/RAY/(2.D0*RAY/SIN(PAR2)-1.D0)) 622 : ELSE 623 : XDD=1.D0 624 : ENDIF 625 : XCOOR(IADR*(IDIM+1)+1)=XCEN+(XV1*COSA+XV3*SINA)*XDD 626 : XCOOR(IADR*(IDIM+1)+2)=YCEN+(YV1*COSA+YV3*SINA)*XDD 627 : IF (IDIM.GE.3) THEN 628 : XCOOR(IADR*(IDIM+1)+3)=ZCEN+(ZV1*COSA+ZV3*SINA)*XDD 629 : ENDIF 630 : GOTO 1220 631 : 1207 PARA=2*(DPAR1-0.5D0) 632 : XCOOR(IADR*(IDIM+1)+1)=XIN+PARA*XDIS+(1-PARA*PARA)*XV 633 : XCOOR(IADR*(IDIM+1)+2)=YIN+PARA*YDIS+(1-PARA*PARA)*YV 634 : IF (IDIM.GE.3) 635 : . XCOOR(IADR*(IDIM+1)+3)=ZIN+PARA*ZDIS+(1-PARA*PARA)*ZV 636 : GOTO 1220 637 : 1208 GOTO 1220 638 : 1220 XCOOR((IADR+1)*(IDIM+1))=DENI+DECA*DPAR1 639 : IADR=IADR+1 640 : ENDIF 641 : 642 : C Cas particulier de l'arc de CUBIQUE (CUBP) avec option UNIF 643 : C :::::::::::::::: AJOUT PIFFARD POUR CUBIQUE PAR 4 POINTS :::::::::::: 644 : C .... SI CUBIQUE ET SI IUNIF = 1 CORRECTION POUR PTS UNIFORMES. 645 : IF ((ILIG.EQ.9).AND.(IUNIF.EQ.1)) THEN 646 : C .... ON RECALCULE LES POINTS AVEC LA CORRECTION SUR LES PARAM U . 647 : C .. ON FAIT 5 ITERATIONS MAXI (ON PEUT SORTIR AVT PAR DMOY) 648 : DLONGT(1)=DLONGT(1)+SQRT((XFIN-XPOIN(NX))**2 649 : & +(YFIN-YPOIN(NX))**2 650 : & +(ZFIN-ZPOIN(NX))**2) 651 : DMOY(1)=DLONGT(1)/NX 652 : IF (IIMPI.EQ.1) WRITE(IOIMP,6003) DMOY(1) 653 : 6003 FORMAT(2X,'DISTANCE MOYENNE =',G13.5) 654 : DO KT=1,NITER 655 : IF (IIMPI.EQ.1) WRITE(IOIMP,6000) KT,NITER 656 : 6000 FORMAT(2X,'ITERATION',I2,' SUR ',I2,' DEMANDEES') 657 : DO I=IDEP,NX 658 : U(KT+1,I)=U(KT,I)*(I-1.D0)/NX*DLONGT(KT)/DIST(I) 659 : UP=U(KT+1,I) 660 : F1=2*UP**3-3*UP**2+1 661 : F2=-2*UP**3+3*UP**2 662 : F3=UP**3-2*UP**2+UP 663 : F4=UP**3-UP**2 664 : DO J=1,4 665 : TU(J)=+CU(J,1)*XINI+CU(J,2)*XINT1 666 : & +CU(J,3)*XINT2+CU(J,4)*XFIN 667 : ENDDO 668 : XPOIN(I)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4) 669 : DO J=1,4 670 : TU(J)=+CU(J,1)*YINI+CU(J,2)*YINT1 671 : & +CU(J,3)*YINT2+CU(J,4)*YFIN 672 : ENDDO 673 : YPOIN(I)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4) 674 : ZPOIN(I)=0.D0 675 : IF (IDIM.GE.3) THEN 676 : DO J=1,4 677 : TU(J)=+CU(J,1)*ZINI+CU(J,2)*ZINT1 678 : & +CU(J,3)*ZINT2+CU(J,4)*ZFIN 679 : ENDDO 680 : ZPOIN(I)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4) 681 : ENDIF 682 : DIST2(I)=SQRT((XPOIN(I)-XPOIN(I-1))**2 683 : & +(YPOIN(I)-YPOIN(I-1))**2+(ZPOIN(I)-ZPOIN(I-1))**2) 684 : IF (IIMPI.EQ.1) WRITE(IOIMP,6001) DIST2(I) 685 : 6001 FORMAT(2X,'DISTANCE ENTRE 2 PTS =',G13.5) 686 : DIST(I)=DIST(I-1)+DIST2(I) 687 : DLONGT(KT+1)=DLONGT(KT+1)+DIST2(I) 688 : ENDDO 689 : DLONGT(KT+1)=DLONGT(KT+1)+SQRT((XFIN-XPOIN(NX))**2 690 : & +(YFIN-YPOIN(NX))**2+(ZFIN-ZPOIN(NX))**2) 691 : IF (IIMPI.EQ.1) WRITE(IOIMP,6002) DLONGT(KT+1) 692 : 6002 FORMAT(2X,'DISTANCE TOTALE ENTRE PTS =',G13.5) 693 : DMOY(KT)=DLONGT(KT+1)/NX 694 : IF (IIMPI.EQ.1) WRITE(IOIMP,6003) DMOY(KT) 695 : C6003 FORMAT(2X,'DISTANCE MOYENNE =',G13.5) 696 : CRIT=0.05D0*DMOY(KT) 697 : ICHANG=0 698 : DO IL=IDEP,NX 699 : IF (ABS(DIST2(IL)-DMOY(KT)).GT.CRIT) ICHANG=ICHANG+1 700 : ENDDO 701 : IF (ICHANG.EQ.0) GOTO 307 702 : ENDDO 703 : C .... J'AI FINI LA CORRECTION SI : 704 : C .. - JE SUIS ARRIVE AU BOUT DES ITERATIONS 705 : C . - J'AI TROUVE LA BONNE REPARTITION 706 : C ---- STOCKAGE DES POINTS DS LA BASE GIBI ----- 707 : 307 DPAR=0.D0 708 : NBPTS=IADR+NX-IDEP+1 709 : SEGADJ MCOORD 710 : DO i=IDEP,NX 711 : DPAR=DPAR+TABPAR(i) 712 : XCOOR(IADR*(IDIM+1)+1)=XPOIN(i) 713 : XCOOR(IADR*(IDIM+1)+2)=YPOIN(i) 714 : IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZPOIN(i) 715 : XCOOR((IADR+1)*(IDIM+1))=DENI+DECA+DPAR 716 : IADR=IADR+1 717 : ENDDO 718 : ENDIF 719 : 720 : C Un peu de menage ! 721 : SEGSUP TABPAR 722 : IF ((ILIG.EQ.9).OR.(ILIG.EQ.10)) SEGSUP STRAV 723 : 724 : IPT2=MELEME 725 : DO i=1,NUM(/2) 726 : ICOLOR(i)=IDCOUL 727 : ENDDO 728 : IF (ITP1.NE.0) THEN 729 : ltelq=.false. 729 : CALL FUSE(IPT1,MELEME,IPT2,ltelq) 730 : SEGDES IPT1,MELEME 731 : ENDIF 732 : MELEME=IPT2 733 : IF (ITP2.NE.0) THEN 734 : ltelq=.false. 734 : CALL FUSE(MELEME,IP2,IPT2,ltelq) 735 : SEGDES IPT6,MELEME 736 : ENDIF 737 : MELEME=IPT2 738 : SEGDES MELEME 739 : IF (IERR.NE.0) SEGSUP MELEME 739 : CALL ECROBJ('MAILLAGE',MELEME) 740 : 741 : RETURN 742 : END 743 : 744 : 745 : 746 : 747 :
© Cast3M 2003 - All rights reserved.
Disclaimer