1 : C DIAGO SOURCE FANDEUR 11/07/19 21:15:57 7039 2 : C---------------------------------------------------------------------| 3 : C | 4 : LOGICAL FUNCTION DIAGO(IP,JP,XCRIT) 5 : C | 6 : C CETTE FONCTION LOGIQUE VERIFIE QUE LE SEGMENT >IP,JP! | 7 : C CREE EST VALIDE. | 8 : C -1- CE SEGMENT NE DOIT PAS ETRE LA DIAGONALE D'UNE | 9 : C FACETTE CARREE DEJA EXISTANTE. | 10 : C -2- CE SEGMENT NE DOIT PAS ETRE TROP LONG | 11 : C | 12 : C---------------------------------------------------------------------| 13 : C 14 : IMPLICIT INTEGER(I-N) 15 : IMPLICIT REAL*8(A-H,O-Z) 16 : -INC CCOPTIO 17 : -INC TDEMAIT 18 : LOGICAL REPONS,IN,in2,vervol 19 : C 20 : C WRITE(6,1000) 21 : C1000 FORMAT(' *** DIAGO ***') 22 : C 23 : REPONS=.FALSE. 24 : DIAGO=.FALSE. 25 : IFDIAG=0 26 : C 27 : C LE SEGMENT APPARTIENT-IL A UNE FACETTE CARREE ? 28 : C QUESTION SUBSIDIAIRE EXISTE-T-IL DEJA ??? | 29 : C ----------------------------------------------- 30 : DO 100 I=1,40 31 : JF=NPF(I,IP) 32 : IF (JF.EQ.0) GOTO 200 33 : DO 110 J=1,40 34 : IF (NPF(J,JP).EQ.0) GOTO 100 35 : IF (JF.NE.NPF(J,JP)) GOTO 110 36 : * LE SEGMENT EST UNE ARETE D'UNE FACETTE TRIANGULAIRE ==> OK 37 : IF (NFC(4,JF).EQ.0) RETURN 38 : IF (IP.EQ.NFC(1,JF).AND.JP.EQ.NFC(3,JF)) REPONS=.TRUE. 39 : IF (IP.EQ.NFC(3,JF).AND.JP.EQ.NFC(1,JF)) REPONS=.TRUE. 40 : IF (IP.EQ.NFC(2,JF).AND.JP.EQ.NFC(4,JF)) REPONS=.TRUE. 41 : IF (IP.EQ.NFC(4,JF).AND.JP.EQ.NFC(2,JF)) REPONS=.TRUE. 42 : IF (REPONS) GOTO 200 43 : * LE SEGMENT EST UNE ARETE D'UNE FACETTE QUADRANGULAIRE ==> OK 44 : RETURN 45 : 110 CONTINUE 46 : 100 CONTINUE 47 : 200 CONTINUE 48 : C 49 : DIAGO=REPONS 50 : C 51 : IF (.NOT.REPONS) GOTO 300 52 : C 53 : WRITE(6,1010)IP,JP,JF,nfc(1,jf),nfc(2,jf),nfc(3,jf) 54 : 1010 FORMAT(' LE SEGMENT >',I3,',',I3,'! APPARTIENT A LA FACETTE ', 55 : # I3,3i5) 56 : C 57 : RETURN 58 : * 59 : * TEST QUE L'ARETE N'EST PAS INSCRITE DANS UNE FACETTE 60 : * 61 : 300 CONTINUE 62 : * test pour voir si l'arete est dans le volume a mailler 63 : * premier test rapide par les volumes aux facettes issues du sommet 64 : * en cas d'echec test supplementaire 65 : DO 900 I=1,40 66 : JF=NPF(I,IP) 67 : IF (JF.EQ.0) GOTO 901 68 : if (vol(jp,nfc(1,jf),nfc(2,jf),nfc(3,jf)).gt.0.d0) goto 950 69 : if (nfc(4,jf).eq.0) goto 900 70 : if (vol(jp,nfc(1,jf),nfc(3,jf),nfc(4,jf)).gt.0.d0) goto 950 71 : 900 continue 72 : 901 continue 73 : DO 910 I=1,40 74 : JF=NPF(I,JP) 75 : IF (JF.EQ.0) GOTO 911 76 : if (vol(ip,nfc(1,jf),nfc(2,jf),nfc(3,jf)).gt.0.d0) goto 950 77 : if (nfc(4,jf).eq.0) goto 910 78 : if (vol(ip,nfc(1,jf),nfc(3,jf),nfc(4,jf)).gt.0.d0) goto 950 79 : 910 continue 80 : 911 continue 81 : goto 960 82 : 950 continue 83 : * write (6,*) ' tests complementaires ' 84 : do 3 i=1,4 85 : xyz(i,nptmax+1)=0.99d0*xyz(i,ip)+0.01d0*xyz(i,jp) 86 : xyz(i,nptmax+2)=0.01d0*xyz(i,ip)+0.99d0*xyz(i,jp) 87 : 3 continue 88 : diago=.not.vervol(nptmax+1,ip,ip,-1,-1,-1,-1,-1,-1) 89 : if (diago) return 90 : diago=.not.vervol(nptmax+1,jp,jp,-1,-1,-1,-1,-1,-1) 91 : if (diago) return 92 : 960 continue 93 : XIP=XYZ(1,IP) 94 : YIP=XYZ(2,IP) 95 : ZIP=XYZ(3,IP) 96 : XJP=XYZ(1,JP) 97 : YJP=XYZ(2,JP) 98 : ZJP=XYZ(3,JP) 99 : DO 310 I=1,40 100 : IF=NPF(I,IP) 101 : IF (IF.EQ.0) GOTO 315 102 : II=ISUCC(IF,IP) 103 : IF (II.EQ.JP) GOTO 310 104 : JJ=IPRED(IF,IP) 105 : IF (JJ.EQ.JP) GOTO 310 106 : XII=XYZ(1,II)-XIP 107 : YII=XYZ(2,II)-YIP 108 : ZII=XYZ(3,II)-ZIP 109 : VII=SQRT(XII**2+YII**2+ZII**2) 110 : XII=XII/VII 111 : YII=YII/VII 112 : ZII=ZII/VII 113 : XJJ=XYZ(1,JJ)-XIP 114 : YJJ=XYZ(2,JJ)-YIP 115 : ZJJ=XYZ(3,JJ)-ZIP 116 : VJJ=SQRT(XJJ**2+YJJ**2+ZJJ**2) 117 : XJJ=XJJ/VJJ 118 : YJJ=YJJ/VJJ 119 : ZJJ=ZJJ/VJJ 120 : * (IP,JP)^(IP,II) 121 : XN1=(YJP-YIP)* ZII -(ZJP-ZIP)* YII 122 : YN1=(ZJP-ZIP)* XII -(XJP-XIP)* ZII 123 : ZN1=(XJP-XIP)* YII -(YJP-YIP)* XII 124 : VN1=SQRT(XN1**2+YN1**2+ZN1**2) 125 : * (IP,JP)^(IP,JJ) 126 : XN2=(YJP-YIP)* ZJJ -(ZJP-ZIP)* YJJ 127 : YN2=(ZJP-ZIP)* XJJ -(XJP-XIP)* ZJJ 128 : ZN2=(XJP-XIP)* YJJ -(YJP-YIP)* XJJ 129 : VN2=SQRT(XN2**2+YN2**2+ZN2**2) 130 : IF((XJP-XIP)*(XII+XJJ)+(YJP-YIP)*(YII+YJJ)+ 131 : # (ZJP-ZIP)*(ZII+ZJJ).LT.0.d0) GOTO 310 132 : SCAL=XN1*XN2+YN1*YN2+ZN1*ZN2 133 : IF (SCAL.LT.-xcrit*VN1*VN2) THEN 134 : if (iimpi.ne.0) 135 : > WRITE (6,*) ' DIAGO-1 ',IP,JP,' EST DANS ',IF, 136 : > nfc(1,if),nfc(2,if),nfc(3,if) 137 : DIAGO=.TRUE. 138 : ifdiag=if 139 : RETURN 140 : ENDIF 141 : 310 CONTINUE 142 : 315 CONTINUE 143 : DO 320 I=1,40 144 : IF=NPF(I,JP) 145 : IF (IF.EQ.0) GOTO 325 146 : II=ISUCC(IF,JP) 147 : IF (II.EQ.IP) GOTO 320 148 : JJ=IPRED(IF,JP) 149 : IF (JJ.EQ.IP) GOTO 320 150 : XII=XYZ(1,II)-XJP 151 : YII=XYZ(2,II)-YJP 152 : ZII=XYZ(3,II)-ZJP 153 : VII=SQRT(XII**2+YII**2+ZII**2) 154 : XII=XII/VII 155 : YII=YII/VII 156 : ZII=ZII/VII 157 : XJJ=XYZ(1,JJ)-XJP 158 : YJJ=XYZ(2,JJ)-YJP 159 : ZJJ=XYZ(3,JJ)-ZJP 160 : VJJ=SQRT(XJJ**2+YJJ**2+ZJJ**2) 161 : XJJ=XJJ/VJJ 162 : YJJ=YJJ/VJJ 163 : ZJJ=ZJJ/VJJ 164 : * (JP,IP)^(JP,II) 165 : XN1=(YIP-YJP)* ZII -(ZIP-ZJP)* YII 166 : YN1=(ZIP-ZJP)* XII -(XIP-XJP)* ZII 167 : ZN1=(XIP-XJP)* YII -(YIP-YJP)* XII 168 : VN1=SQRT(XN1**2+YN1**2+ZN1**2) 169 : * (JP,IP)^(JP,JJ) 170 : XN2=(YIP-YJP)* ZJJ -(ZIP-ZJP)* YJJ 171 : YN2=(ZIP-ZJP)* XJJ -(XIP-XJP)* ZJJ 172 : ZN2=(XIP-XJP)* YJJ -(YIP-YJP)* XJJ 173 : VN2=SQRT(XN2**2+YN2**2+ZN2**2) 174 : IF((XIP-XJP)*(XII+XJJ)+(YIP-YJP)*(YII+YJJ)+ 175 : # (ZIP-ZJP)*(ZII+ZJJ).LT.0d0) GOTO 320 176 : SCAL=XN1*XN2+YN1*YN2+ZN1*ZN2 177 : IF (SCAL.LT.-xcrit*VN1*VN2) THEN 178 : if (iimpi.ne.0) 179 : > WRITE (6,*) ' DIAGO-2 ',JP,IP,' EST DANS ',IF, 180 : > nfc(1,if),nfc(2,if),nfc(3,if) 181 : ifdiag=if 182 : DIAGO=.TRUE. 183 : RETURN 184 : ENDIF 185 : 320 CONTINUE 186 : 325 CONTINUE 187 : * VERIFICATION DES DISTANCES AVEC LES AUTRES ARETES RENTRANTES 188 : * RECHERCHE SI ARETE EXISTE DEJA 189 : DO 500 I=1,40 190 : IF=NPF(I,IP) 191 : IF (IF.EQ.0) GOTO 501 192 : DO 502 J=1,40 193 : JF=NPF(J,JP) 194 : IF (JF.EQ.0) GOTO 500 195 : IF (IF.EQ.JF) GOTO 800 196 : 502 CONTINUE 197 : 500 CONTINUE 198 : 501 CONTINUE 199 : XP1=XYZ(1,IP) 200 : YP1=XYZ(2,IP) 201 : ZP1=XYZ(3,IP) 202 : XP2=XYZ(1,JP) 203 : YP2=XYZ(2,JP) 204 : ZP2=XYZ(3,JP) 205 : DP=(XP1-XP2)**2+(YP1-YP2)**2+(ZP1-ZP2)**2 206 : XP3=(XP1+XP2)/2.d0 207 : YP3=(YP1+YP2)/2.d0 208 : ZP3=(ZP1+ZP2)/2.d0 209 : * TEST AVEC LES AUTRES ARETES POSSIBLES 210 : goto 800 211 : do 520 ipt=ip,jp,jp-ip 212 : DO 510 I=1,40 213 : IF=NPF(I,IPT) 214 : if (IF.EQ.0) goto 520 215 : kp=isucc(if,ipt) 216 : IF (IP.EQ.KP) GOTO 510 217 : IF (JP.EQ.KP) GOTO 510 218 : do 511 J=1,40 219 : jf=npf(j,kp) 220 : if (jF.EQ.0) goto 510 221 : lp=isucc(jf,kp) 222 : IF (IP.EQ.LP) GOTO 511 223 : IF (JP.EQ.LP) GOTO 511 224 : XQ1=XYZ(1,KP) 225 : YQ1=XYZ(2,KP) 226 : ZQ1=XYZ(3,KP) 227 : XQ2=XYZ(1,LP) 228 : YQ2=XYZ(2,LP) 229 : ZQ2=XYZ(3,LP) 230 : DQ=(XQ1-XQ2)**2+(YQ1-YQ2)**2+(ZQ1-ZQ2)**2 231 : XQ3=(XQ1+XQ2)/2.d0 232 : YQ3=(YQ1+YQ2)/2.d0 233 : ZQ3=(ZQ1+ZQ2)/2.d0 234 : dmin=min(dp,dq)/64d0 235 : * DT=(XP1-XQ1)**2+(YP1-YQ1)**2+(ZP1-ZQ1)**2 236 : * IF (DT.LE.dmin) GOTO 810 237 : * DT=(XP1-XQ2)**2+(YP1-YQ2)**2+(ZP1-ZQ2)**2 238 : * IF (DT.LE.dmin) GOTO 810 239 : * DT=(XP1-XQ3)**2+(YP1-YQ3)**2+(ZP1-ZQ3)**2 240 : * IF (DT.LE.dmin) GOTO 810 241 : * DT=(XP2-XQ1)**2+(YP2-YQ1)**2+(ZP2-ZQ1)**2 242 : * IF (DT.LE.dmin) GOTO 810 243 : * DT=(XP2-XQ2)**2+(YP2-YQ2)**2+(ZP2-ZQ2)**2 244 : * IF (DT.LE.dmin) GOTO 810 245 : * DT=(XP2-XQ3)**2+(YP2-YQ3)**2+(ZP2-ZQ3)**2 246 : * IF (DT.LE.dmin) GOTO 810 247 : DT=(XP3-XQ1)**2+(YP3-YQ1)**2+(ZP3-ZQ1)**2 248 : IF (DT.LE.dmin) GOTO 810 249 : DT=(XP3-XQ2)**2+(YP3-YQ2)**2+(ZP3-ZQ2)**2 250 : IF (DT.LE.dmin) GOTO 810 251 : DT=(XP3-XQ3)**2+(YP3-YQ3)**2+(ZP3-ZQ3)**2 252 : IF (DT.LE.dmin) GOTO 810 253 : 511 continue 254 : 510 CONTINUE 255 : 520 CONTINUE 256 : 800 CONTINUE 257 : DIAGO=.FALSE. 258 : RETURN 259 : 810 CONTINUE 260 : DIAGO=.TRUE. 261 : IFDIAG=JF 262 : if (iimpi.ne.0) 263 : >WRITE (6,*) ' ARETE TROP PROCHES ',IP,JP,KP,LP,dt,dmin 264 : * WRITE (6,*) 'xp1,yp1,zp1 ',xp1,yp1,zp1 265 : * WRITE (6,*) 'xp2,yp2,zp2 ',xp2,yp2,zp2 266 : * WRITE (6,*) 'xp3,yp3,zp3 ',xp3,yp3,zp3 267 : * WRITE (6,*) 'xq1,yq1,zq1 ',xq1,yq1,zq1 268 : * WRITE (6,*) 'xq2,yq2,zq2 ',xq2,yq2,zq2 269 : * WRITE (6,*) 'xq3,yq3,zq3 ',xq3,yq3,zq3 270 : RETURN 271 : END 272 : 273 :
© Cast3M 2003 - All rights reserved.
Disclaimer