1 : C KDOM7 SOURCE CHAT 05/01/13 00:54:31 5004 2 : SUBROUTINE KDOM7(IP1,IP2,IP3,IP4,XCEN) 3 : C 4 : C************************************************************************ 5 : C 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : KDOM7 9 : C 9 : C DESCRIPTION : Subroutine called by LEKMA0 10 : C Given a QUA4 (IP1,IP2,IP3,IP4), this subroutine 11 : C compute its center by averaging the centers of 12 : C the two possible configuration it can be 13 : C divided in 14 : C 15 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 16 : C 17 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF 18 : C 19 : C************************************************************************ 20 : C 21 : C 22 : C**** Variables de COOPTIO 23 : C 24 : IMPLICIT INTEGER(I-N) 25 : INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 26 : & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 27 : & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 28 : & ,IECHO, IIMPI, IOSPI 29 : & ,IDIM, IFICLE, IPREFI 30 : C & ,MCOORD 31 : & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 32 : & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 33 : & ,NORINC,NORVAL,NORIND,NORVAD 34 : & ,NUCROU, IPSAUV, IREFOR, ISAFOR 35 : C 36 : INTEGER IP1,IP2,IP3, IP4, I1,I2, IPCEN, IPL, IPR 37 : REAL*8 VOL, XCEN(3),VOLT,DX0,DY0,DZ0,DX1,DY1,DZ1,CEL 38 : & ,RNOX0,RNOY0,RNOZ0,RNOX1,RNOY1,RNOZ1,XCEN1(3),VOLT1,PSCA 39 : & ,RNOT 40 : LOGICAL LOGREG 41 : -INC CCOPTIO 42 : -INC SMCOORD 43 : C 44 : DZ0=0.0D0 45 : DZ1=0.0D0 46 : RNOT=0.0D0 47 : C 48 : C**** We chose as positive direction the one which 49 : C has the maximum modulus IP1IP(1+1) x IP1IP(1+3) 50 : C 51 : DO I1=1,4 52 : IF(I1 .EQ. 1)THEN 53 : C We consider the triangles IP1, IP2, IP4 54 : IPCEN=IP1 55 : IPR=IP2 56 : IPL=IP4 57 : ELSEIF(I1 .EQ. 2)THEN 58 : IPCEN=IP2 59 : IPR=IP3 60 : IPL=IP1 61 : ELSEIF(I1 .EQ. 3)THEN 62 : IPCEN=IP3 63 : IPR=IP4 64 : IPL=IP2 65 : ELSEIF(I1 .EQ. 4)THEN 66 : IPCEN=IP4 67 : IPR=IP1 68 : IPL=IP3 69 : ENDIF 70 : DX0=XCOOR((IPR-1)*(IDIM+1)+1)- 71 : & XCOOR((IPCEN-1)*(IDIM+1)+1) 72 : DY0=XCOOR((IPR-1)*(IDIM+1)+2)- 73 : & XCOOR((IPCEN-1)*(IDIM+1)+2) 74 : IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)- 75 : & XCOOR((IPCEN-1)*(IDIM+1)+3) 76 : DX1=XCOOR((IPL-1)*(IDIM+1)+1)- 77 : & XCOOR((IPCEN-1)*(IDIM+1)+1) 78 : DY1=XCOOR((IPL-1)*(IDIM+1)+2)- 79 : & XCOOR((IPCEN-1)*(IDIM+1)+2) 80 : IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)- 81 : & XCOOR((IPCEN-1)*(IDIM+1)+3) 82 : RNOX1=(DY0*DZ1-DY1*DZ0) 83 : RNOY1=(DZ0*DX1-DZ1*DX0) 84 : RNOZ1=(DX0*DY1-DX1*DY0) 85 : VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0 86 : VOL=VOL/0.5D0 87 : IF(VOL .GT. RNOT)THEN 88 : RNOT=VOL 89 : RNOX0=RNOX1 90 : RNOY0=RNOY1 91 : RNOZ0=RNOZ1 92 : ENDIF 93 : ENDDO 94 : C 95 : C**** First division in triangles 96 : C We consider the triangles IP3, IP1, IP2 and 97 : C IP4, IP1, IP3 98 : C 99 : XCEN1(1)=0.0D0 100 : XCEN1(2)=0.0D0 101 : XCEN1(3)=0.0D0 102 : XCEN(1)=0.0D0 103 : XCEN(2)=0.0D0 104 : XCEN(3)=0.0D0 105 : VOLT=0.0D0 106 : VOLT1=0.0D0 107 : LOGREG=.TRUE. 108 : C 109 : DO I1=1,2 110 : IF(I1 .EQ. 1)THEN 111 : C We consider the triangles IP3, IP1, IP2 112 : IPCEN=IP1 113 : IPR=IP2 114 : IPL=IP3 115 : ELSEIF(I1 .EQ. 2)THEN 116 : IPCEN=IP1 117 : IPR=IP3 118 : IPL=IP4 119 : ENDIF 120 : DX0=XCOOR((IPR-1)*(IDIM+1)+1)- 121 : & XCOOR((IPCEN-1)*(IDIM+1)+1) 122 : DY0=XCOOR((IPR-1)*(IDIM+1)+2)- 123 : & XCOOR((IPCEN-1)*(IDIM+1)+2) 124 : IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)- 125 : & XCOOR((IPCEN-1)*(IDIM+1)+3) 126 : DX1=XCOOR((IPL-1)*(IDIM+1)+1)- 127 : & XCOOR((IPCEN-1)*(IDIM+1)+1) 128 : DY1=XCOOR((IPL-1)*(IDIM+1)+2)- 129 : & XCOOR((IPCEN-1)*(IDIM+1)+2) 130 : IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)- 131 : & XCOOR((IPCEN-1)*(IDIM+1)+3) 132 : RNOX1=(DY0*DZ1-DY1*DZ0) 133 : RNOY1=(DZ0*DX1-DZ1*DX0) 134 : RNOZ1=(DX0*DY1-DX1*DY0) 135 : VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0 136 : VOL=VOL*0.5D0 137 : PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1) 138 : IF((VOL .LE. (RNOT*1.0D-6)) .OR. (PSCA.LE.0))THEN 139 : LOGREG=.FALSE. 140 : ENDIF 141 : VOLT1=VOLT1+VOL 142 : DO I2=1,IDIM,1 143 : CEL=(XCOOR((IPCEN-1)*(IDIM+1)+I2)+ 144 : & XCOOR((IPL-1)*(IDIM+1)+I2)+XCOOR((IPR-1)*(IDIM+1)+I2)) 145 : & /3.0D0 146 : XCEN1(I2)=XCEN1(I2)+(VOL*CEL) 147 : ENDDO 148 : ENDDO 149 : C 150 : C**** If one of the area is negative or almost zero the division is not 151 : C possible (LOGREG=.FALSE.) 152 : C 153 : IF(LOGREG)THEN 154 : VOLT=VOLT1 155 : DO I2=1,IDIM,1 156 : XCEN(I2)=XCEN1(I2) 157 : ENDDO 158 : ENDIF 159 : C 160 : C**** Second division in triangles 161 : C We consider the triangles IP3, IP2, IP4 and 162 : C IP4, IP2, IP1 163 : C 164 : C 165 : VOLT1=0.0D0 166 : LOGREG=.TRUE. 167 : XCEN1(1)=0.0D0 168 : XCEN1(2)=0.0D0 169 : XCEN1(3)=0.0D0 170 : C 171 : DO I1=1,2 172 : IF(I1 .EQ. 1)THEN 173 : C We consider the triangles 174 : IPCEN=IP2 175 : IPR=IP3 176 : IPL=IP4 177 : ELSEIF(I1 .EQ. 2)THEN 178 : IPCEN=IP2 179 : IPR=IP4 180 : IPL=IP1 181 : ENDIF 182 : DX0=XCOOR((IPR-1)*(IDIM+1)+1)- 183 : & XCOOR((IPCEN-1)*(IDIM+1)+1) 184 : DY0=XCOOR((IPR-1)*(IDIM+1)+2)- 185 : & XCOOR((IPCEN-1)*(IDIM+1)+2) 186 : IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)- 187 : & XCOOR((IPCEN-1)*(IDIM+1)+3) 188 : DX1=XCOOR((IPL-1)*(IDIM+1)+1)- 189 : & XCOOR((IPCEN-1)*(IDIM+1)+1) 190 : DY1=XCOOR((IPL-1)*(IDIM+1)+2)- 191 : & XCOOR((IPCEN-1)*(IDIM+1)+2) 192 : IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)- 193 : & XCOOR((IPCEN-1)*(IDIM+1)+3) 194 : RNOX1=(DY0*DZ1-DY1*DZ0) 195 : RNOY1=(DZ0*DX1-DZ1*DX0) 196 : RNOZ1=(DX0*DY1-DX1*DY0) 197 : VOL=((RNOX1**2)+(RNOY1**2)+(RNOZ1**2))**0.5D0 198 : VOL=VOL*0.5D0 199 : PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1) 200 : IF((VOL .LE. (RNOT*1.0D-6)) .OR. (PSCA.LE.0))THEN 201 : LOGREG=.FALSE. 202 : ENDIF 203 : VOLT1=VOLT1+VOL 204 : DO I2=1,IDIM,1 205 : CEL=(XCOOR((IPCEN-1)*(IDIM+1)+I2)+ 206 : & XCOOR((IPL-1)*(IDIM+1)+I2)+XCOOR((IPR-1)*(IDIM+1)+I2)) 207 : & /3.0D0 208 : XCEN1(I2)=XCEN1(I2)+(VOL*CEL) 209 : ENDDO 210 : ENDDO 211 : C 212 : IF(LOGREG)THEN 213 : VOLT=VOLT+VOLT1 214 : DO I2=1,IDIM,1 215 : XCEN(I2)=XCEN(I2)+XCEN1(I2) 216 : ENDDO 217 : ENDIF 218 : C 219 : DO I2=1,IDIM,1 220 : XCEN(I2)=XCEN(I2)/VOLT 221 : ENDDO 222 : C 223 : C**** Finally, we check whether the center of gravity is "inside" 224 : C the domain, i.e. whether the normals CENIPi x CENIP(i+1) 225 : C have the same direction of (RNOX0,RNOY0,RNOZ0) 226 : C 227 : LOGREG=.FALSE. 228 : DO I1=1,4 229 : IF(I1 .EQ. 1)THEN 230 : C We consider the triangles CEN, IP1, IP2 231 : IPR=IP1 232 : IPL=IP2 233 : ELSEIF(I1 .EQ. 2)THEN 234 : IPR=IP2 235 : IPL=IP3 236 : ELSEIF(I1 .EQ. 3)THEN 237 : IPR=IP3 238 : IPL=IP4 239 : ELSEIF(I1 .EQ. 4)THEN 240 : IPR=IP4 241 : IPL=IP1 242 : ENDIF 243 : DX0=XCOOR((IPR-1)*(IDIM+1)+1)- 244 : & XCEN(1) 245 : DY0=XCOOR((IPR-1)*(IDIM+1)+2)- 246 : & XCEN(2) 247 : IF(IDIM. EQ. 3) DZ0=XCOOR((IPR-1)*(IDIM+1)+3)- 248 : & XCEN(3) 249 : DX1=XCOOR((IPL-1)*(IDIM+1)+1)- 250 : & XCEN(1) 251 : DY1=XCOOR((IPL-1)*(IDIM+1)+2)- 252 : & XCEN(2) 253 : IF(IDIM. EQ. 3) DZ1=XCOOR((IPL-1)*(IDIM+1)+3)- 254 : & XCEN(3) 255 : RNOX1=(DY0*DZ1-DY1*DZ0) 256 : RNOY1=(DZ0*DX1-DZ1*DX0) 257 : RNOZ1=(DX0*DY1-DX1*DY0) 258 : PSCA=(RNOX0*RNOX1)+(RNOY0*RNOY1)+(RNOZ0*RNOZ1) 259 : IF(PSCA .LE. 0.0D0)THEN 260 : LOGREG=.TRUE. 261 : ENDIF 262 : ENDDO 263 : C 264 : IF(LOGREG)THEN 265 : C 266 : WRITE(IOIMP,*) 'Subroutine kdom7.eso' 267 : WRITE(IOIMP,*) 'QUA4' 268 : C 269 : C Erreur -107 0 270 : C Liste des noeuds de l'élément : 271 : C 271 : CALL ERREUR(-107) 272 : WRITE(IOIMP,*) IP1, IP2, IP3, IP4 273 : C 274 : C Erreur 845 2 275 : C Maillage donné incorrect ?!!! 276 : C 276 : CALL ERREUR(845) 277 : C 278 : ENDIF 279 : C 280 : RETURN 281 : END 282 : 283 :
© Cast3M 2003 - All rights reserved.
Disclaimer