1 : C SENS1 SOURCE CHAT 05/01/13 03:13:33 5004 2 : SUBROUTINE SENS1 3 : c 4 : c sous routine castem 2000 5 : c 6 : c utilisée pour l'opérateur SENS 7 : c determine dans quel sens tourne un contour en 2 D 8 : c 9 : c synthaxe de l'opérateur TAB2 = SENS TAB1 ; 10 : c 11 : c tab1 est issue de CCON et contient des maillages indicé par 12 : c des entiers. ces maillages doivent etre des contours fermes 13 : c orientes. 14 : c 15 : c tab2 contient des entiers (+/-1) indice par les meme entiers 16 : c +1 indique que le contour est dans le sens trigonometrique 17 : c -1 indique que le contour est dans le sens inverse 18 : c 19 : c langage Fortran 77 + esope 20 : c auteur A de Gayffier 21 : c 22 : IMPLICIT INTEGER(I-N) 23 : IMPLICIT REAL*8 (A-H,O-Z) 24 : CHARACTER CHARI,CHARR 25 : LOGICAL LOGII,LOGIR 26 : c 27 : -INC CCREEL 28 : -INC SMTABLE 29 : -INC SMELEME 30 : -INC CCOPTIO 31 : -INC SMCOORD 32 : c 33 : c declaration du segment INCR 34 : c 35 : SEGMENT INCR 36 : INTEGER NUM0(NBELEM) 37 : ENDSEGMENT 38 : c 39 : c operateur disponible seulement den dimension 2 40 : IF (IDIM .NE. 2 ) THEN 40 : CALL ERREUR(19) 41 : RETURN 42 : ENDIF 43 : c 43 : CALL LIROBJ ('TABLE',MTABLE,1,IRET) 44 : IF (IRET .EQ. 0) RETURN 45 : SEGACT ,MTABLE 46 : c 47 : c compte combien d'indice entier 48 : NBENT = 0 49 : DO 5 I=1,MTABTI(/2) 50 : IF (MTABTI(I) .EQ. 'ENTIER' ) NBENT = NBENT +1 51 : 5 CONTINUE 52 : c 53 : M = 0 54 : SEGINI ,MTAB1 55 : c 56 : c boucle sur tous les contours 57 : c 58 : DO 10 INDTAB=1,NBENT 59 : c 59 : CALL ACCTAB(MTABLE,'ENTIER',INDTAB,XVALI,CHARI,LOGII,0, 60 : & 'MAILLAGE',IVALR,XVALR,CHARR,LOGIR,IRETR) 61 : IF (IERR .NE. 0 ) RETURN 62 : c 63 : c 64 : MELEME = IRETR 65 : SEGACT MELEME 66 : c 67 : c quelques controles sur le maillage 68 : c 69 : IF (LISOUS(/1) .NE. 0 .OR. (ITYPEL .NE. 2 70 : & .AND. ITYPEL .NE. 3)) THEN 70 : CALL ERREUR(26) 71 : SEGDES MELEME 72 : SEGSUP MTAB1 73 : RETURN 74 : ENDIF 75 : c 76 : c initialisation de la boucle sur les segments2 77 : c 78 : NBELEM = NUM(/2) 79 : SEGINI INCR 80 : I = 1 81 : IF (ITYPEL .EQ. 2 ) THEN 82 : I1 = NUM(1,1) 83 : I2 = NUM(2,1) 84 : ELSE IF ( ITYPEL .EQ. 3 ) THEN 85 : I1 = NUM(1,1) 86 : I2 = NUM(3,1) 87 : ENDIF 88 : STHETA = 0.D0 89 : ISEG2 = 0 90 : c 91 : c Initialisation du tableau NUM0 92 : c 93 : DO 9 J=1,NBELEM 94 : NUM0(J) = NUM(1,J) 95 : 9 CONTINUE 96 : c 97 : c debut de la boucle sur les segments seg2 du maillage 98 : c 99 : 15 CONTINUE 100 : c 101 : DO 20 J=1,NBELEM 102 : IF ( NUM0(J) .EQ. I2 ) THEN 103 : c on a trouve le segment suivant 104 : ISEG2 = J 105 : NUM0(J) = 0 ; 106 : IF (ITYPEL .EQ. 2 ) THEN 107 : I3 = NUM(2,J) 108 : ELSE IF ( ITYPEL .EQ. 3 ) THEN 109 : I3 = NUM(3,J) 110 : ENDIF 111 : GOTO 30 112 : ENDIF 113 : 20 CONTINUE 114 : c 115 : c on n'est pas parvenu a boucler le contour 116 : c 117 : c 118 : c dans le cas axisymmetrique le contour n'est pas necessairement fermé 119 : c 120 : IF (IFOUR .EQ. 0) THEN 121 : c 122 : c on verifie que le point est sur l'axe OZ 123 : c 124 : XI1 = ABS(XCOOR((I2-1)*(IDIM+1) +1 )) 125 : DENS = MAX(XCOOR((I2-1)*(IDIM+1)),1.D-10) 126 : c print *,'i2=',i2,'r=',xi1,'dens=',dens 127 : IF ( XI1 .GT. (DENS/100.d0)) GOTO 25 128 : c 129 : c on cherche l'autre point sur l'axe 130 : c 131 : DO 23 J=1,NBELEM 132 : c 133 : XI2 =ABS(XCOOR((NUM(1,J)-1)*(IDIM+1) +1 )) 134 : DENS =MAX(XCOOR((NUM(1,J)-1)*(IDIM+1)),1.D-10) 135 : c 136 : IF (XI2.LT.(DENS/100.d0)) THEN 137 : I3 = NUM(1,J) 138 : * print *, 'on a trouvé le second point',i3,'r=',xi2 139 : GOTO 30 140 : ENDIF 141 : 23 CONTINUE 142 : * print *, 'pas trouvé le second point sur l axe' 143 : c 144 : ENDIF 145 : C 146 : 25 CONTINUE 147 : c 148 : c gestion de l'erreur contour non ferme 149 : c 149 : CALL ERREUR(28) 150 : SEGSUP MTAB1 151 : SEGDES MELEME,MTABLE 152 : RETURN 153 : c 154 : 30 CONTINUE 155 : c 156 : c I sert a compter le nombre de segments pris en compte 157 : c 158 : I = I + 1 159 : c 160 : c calcul de l'angle entre les deux segments 161 : c 162 : XI1 = XCOOR((I1-1)*(IDIM+1) +1 ) 163 : YI1 = XCOOR((I1-1)*(IDIM+1) +2 ) 164 : XI2 = XCOOR((I2-1)*(IDIM+1) +1 ) 165 : YI2 = XCOOR((I2-1)*(IDIM+1) +2 ) 166 : XI3 = XCOOR((I3-1)*(IDIM+1) +1 ) 167 : YI3 = XCOOR((I3-1)*(IDIM+1) +2 ) 168 : X12 = XI2-XI1 169 : Y12 = YI2-YI1 170 : X23 = XI3-XI2 171 : Y23 = YI3-YI2 172 : D12 = SQRT( X12*X12+Y12*Y12) 173 : D23 = SQRT ( X23*X23+Y23*Y23) 174 : SITHET = ( X12*Y23 - X23*Y12 ) / D12 / D23 175 : COTHET = ( X12*X23 + Y12*Y23 ) / D12 / D23 176 : IF ( (ABS(SITHET)) .LT. 1.D-10 ) THEN 177 : THETA = 0.D0 178 : ELSE 179 : IF ( ABS(SITHET) .GT. 1.D0 ) THEN 180 : THETA = SIGN((XPI/2.D0),SITHET) 181 : ELSE 182 : THETA = ASIN( SITHET ) 183 : ENDIF 184 : ENDIF 185 : IF (COTHET.LT.0.D0 .AND. SITHET.GT.0.D0) THETA=XPI-THETA 186 : IF (COTHET.LT.0.D0 .AND. SITHET.LT.0.D0) THETA=-XPI-THETA 187 : STHETA = STHETA + THETA 188 : I1 = I2 189 : I2 = I3 190 : c 191 : c recherche du point suivant si on n'est pas parvenu 192 : c au point initial 193 : c 194 : IF ( ISEG2 .NE. 1 ) GOTO 15 195 : c 196 : c on controle que tous les segments ont ete examinés 197 : c 198 : IF (I.NE.(NBELEM+1) .AND. I.NE.(NBELEM+2)) THEN 199 : c nbelem+2 correspond au cas axis avec un contour non fermé 199 : CALL ERREUR(28) 200 : SEGSUP MTAB1 201 : SEGDES MELEME,MTABLE 202 : RETURN 203 : ENDIF 204 : c 205 : c determination du sens du contour la somme vaut +/- XPI 206 : c 207 : IF ( ABS(STHETA - XPI*2.D0 ) .LT. 1.D-3 ) THEN 208 : IVALR = 1 209 : ELSE IF ( ABS(STHETA + XPI*2.D0 ) .LT. 1.D-3 ) THEN 210 : IVALR = -1 211 : ELSE 212 : IVALR = 0 213 : c il y un problème: c'est sur! 214 : c print *, 'somme des angles en radian',stheta 215 : INTERR(1)=INDTAB 215 : CALL ERREUR(718) 216 : ENDIF 217 : c 218 : c ecriture dans la table du resultat 219 : c 219 : CALL ECCTAB(MTAB1,'ENTIER',INDTAB,XVALI,CHARI,LOGII,IOBIN, 220 : & 'ENTIER',IVALR,XVALR,CHARR,LOGIR,IRETR) 221 : SEGDES MELEME 222 : SEGSUP INCR 223 : c 224 : 10 CONTINUE 225 : c 226 : c il n'y a plus rien dans la table 226 : CALL ECROBJ('TABLE',MTAB1) 227 : SEGDES MTABLE,MTAB1 228 : RETURN 229 : c 230 : END 231 : 232 : 233 : 234 : 235 : 236 : 237 : 238 : 239 :
© Cast3M 2003 - All rights reserved.
Disclaimer