1 : C KNRF SOURCE MAGN 10/05/18 21:15:54 6675 2 : SUBROUTINE KNRF(MTABD,MCHELM,MCHPOI,MCHPO2) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C 6 : C************************************************************************* 7 : C 8 : C OBJET : Cree deux CHAMPOINT faces et un mchaml 9 : C CHPC1 contenant la longueur des faces 10 : C ( faux dans le cas 3d ) 11 : C CHPC2 contenant la normale aux faces. 12 : C MCHAML orientation (normale entrante ou sortante) 13 : C SYNTAXE : CHPC1 CHPC2 MCHAML = KNRF OBJDOM ; 14 : C 15 : C OBJDOM : TABLE de SOUSTYPE DOMAINE 16 : C 17 : C REMARQUE : 18 : C On verifie que la normale est dans le meme sens que pt1->pt3 19 : C des elements FACEL 20 : C Si pt1=pt3 on oriente la normale vers l'exterieur cad pt1->pt2 21 : C 22 : C Non teste si axisymetrie 23 : C************************************************************************* 24 : C MODIFICATION du 20/11/98 25 : C 26 : C Le calcul de la normale s'effectue a l'aide du jacobien calcule 27 : C par CALJBR (au lieu de CALJQB) 28 : C 29 : C************************************************************************* 30 : -INC CCOPTIO 31 : -INC CCGEOME 32 : -INC SMTABLE 33 : POINTEUR MTABD.MTABLE 34 : -INC SMELEME 35 : POINTEUR MELEMC.MELEME,MELEMD.MELEME,MELEMT.MELEME 36 : POINTEUR MELEMF.MELEME 37 : -INC SMCOORD 38 : -INC SMCHPOI 39 : -INC SMCHAML 40 : -INC SMLENTI 41 : -INC SIZFFB 42 : -INC CCREEL 43 : C 44 : PARAMETER (NBO=5) 45 : REAL*8 A2J(2,2,1),A3J(3,3,1),CFT,NORMX,NORMY,NORMZ 46 : REAL*8 LGR,LGR1,LGR2 47 : CHARACTER*8 TYPE,TYPC,NOM0,LIST1(NBO),LIST2(NBO) 48 : C Les elements complets ne sont pas tous implementes 49 : C On se ramene aux elements facep sans le centre de la face 50 : DATA LIST1/'SEG3 ','TRI4 ','TRI7 ','QUA5 ', 51 : & 'QUA9 '/ 52 : DATA LIST2/'SEG2 ','TRI3 ','TRI6 ','QUA4 ', 53 : & 'QUA8 '/ 54 : C 55 : C Acquisition des connectivites FACEL 56 : C write(6,*)' KNRF ' 57 : TYPE=' ' 57 : CALL ACMO(MTABD,'FACEL',TYPE,MELEMD) 58 : C write(6,*)' FACEL type=',type 59 : IF(TYPE.NE.'MAILLAGE')GO TO 90 60 : SEGACT MELEMD 61 : IPT2=MELEMD 62 : C 63 : C Acquisition des connectivites FACE 64 : C 65 : TYPE=' ' 65 : CALL ACMO(MTABD,'FACE',TYPE,MELEMF) 66 : IF(TYPE.NE.'MAILLAGE')GO TO 90 66 : CALL KRIPAD(MELEMF,MLENTI) 67 : SEGACT MELEMF 68 : C 69 : C Creation du CHPOIN a 3 comp. pour la normale aux faces 70 : C 71 : NC=IDIM 72 : TYPE='FACE' 72 : CALL CRCHPT(TYPE,MELEMF,NC,MCHPOI) 72 : CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) 73 : C 74 : C Creation du CHPOIN pour la longueur 75 : C 76 : NC=1 77 : TYPE='FACE' 77 : CALL CRCHPT(TYPE,MELEMF,NC,MCHPO2) 77 : CALL LICHT(MCHPO2,MPOVA2,TYPC,IGEOM) 78 : SEGDES MELEMF 79 : C 80 : IAXI = 0 81 : IF(IFOMOD.EQ.0)IAXI=2 82 : C 83 : C Acquisition des connectivites FACEP 84 : C 85 : TYPE=' ' 85 : CALL ACMO(MTABD,'FACEP',TYPE,MELEME) 86 : IF(TYPE.NE.'MAILLAGE')GO TO 90 87 : SEGACT MELEME 88 : C 89 : NBSOUS=LISOUS(/1) 90 : IF(NBSOUS.EQ.0)NBSOUS=1 91 : C********************** 92 : C CAS DE LA DIMENSION 2 93 : C********************** 94 : IF (IDIM.EQ.2) THEN 95 : C 96 : DO 1 L=1,NBSOUS 97 : IPT1=MELEME 98 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 99 : SEGACT IPT1 100 : NP=IPT1.NUM(/1)-1 101 : NEL=IPT1.NUM(/2) 102 : C Toutes les faces sont des SEG3 donc NP=2 103 : IF(NP.NE.2)RETURN 104 : C 105 : C On considere la face sans son centre comme un elt fini SEG2 106 : C 107 : TYPE = NOMS(IPT1.ITYPEL)//' ' 107 : CALL OPTLI(IP,LIST1,TYPE,NBO) 107 : IF (IP .EQ. 0) CALL ERREUR(5) 108 : TYPE = LIST2(IP) 108 : CALL KALPBG(TYPE,'FONFORM0',IZFFM) 109 : SEGACT IZFFM*MOD 110 : IZHR=KZHR(1) 111 : SEGACT IZHR*MOD 112 : NES=GR(/1) 113 : NPG=GR(/3) 114 : C Boucle sur toutes les faces pour le paquet L 115 : DO 10 K=1,NEL 116 : NF=LECT(IPT1.NUM(NP+1,K)) 117 : C REMPLISSAGE DE XYZ 118 : C ------------------ 119 : C 120 : DO 12 I=1,NP 121 : J=IPT1.NUM(I,K) 122 : DO 12 N=1,IDIM 123 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 124 : 12 CONTINUE 125 : C 126 : C CALJBR calcul le jacobien du passage de l'elt de ref. 127 : C a l'elt. reel 128 : C A2J=Jacobien AIRE=detA2J 128 : CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP, 129 : &NPG,IAXI,AIRE,A2J,SGEN) 130 : C 131 : NORMX=A2J(1,2,1) 132 : NORMY=A2J(2,2,1) 133 : C On verifie que n est dans le meme sens que vec(13) de FACEL 134 : C Calcul de vec(13) 135 : C FACEL est repere par IPT2 136 : J1=IPT2.NUM(1,NF) 137 : J2=IPT2.NUM(2,NF) 138 : JJ2=LECT(J2) 139 : 7 J3=IPT2.NUM(3,NF) 140 : C 141 : IF (J1.eq.J3) THEN 142 : X1=XCOOR((J1-1)*(IDIM+1)+1) 143 : Y1=XCOOR((J1-1)*(IDIM+1)+2) 144 : X2=XCOOR((J2-1)*(IDIM+1)+1) 145 : Y2=XCOOR((J2-1)*(IDIM+1)+2) 146 : SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY 147 : SGN=SIGN(1.D0,SCAL) 148 : C 149 : ELSE 150 : X1=XCOOR((J1-1)*(IDIM+1)+1) 151 : Y1=XCOOR((J1-1)*(IDIM+1)+2) 152 : X3=XCOOR((J3-1)*(IDIM+1)+1) 153 : Y3=XCOOR((J3-1)*(IDIM+1)+2) 154 : SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY 155 : SGN=SIGN(1.D0,SCAL) 156 : ENDIF 157 : MPOVAL.VPOCHA(JJ2,1)=SGN*NORMX 158 : MPOVAL.VPOCHA(JJ2,2)=SGN*NORMY 159 : C Calcul de la longueur 160 : C 161 : CFT=1.D0 162 : IF(IAXI.NE.0)THEN 163 : TAMP1=(XYZ(1,2)+XYZ(1,1))*0.5D0 164 : TAMP2=(XYZ(2,2)+XYZ(2,1))*0.5D0 165 : ENDIF 166 : C Axisymetrie / ox (non teste) 167 : IF (IAXI.EQ.1) CFT=2.D0*XPI*TAMP2 168 : C Axisymetrie / oy (non teste non plus) 169 : IF (IAXI.EQ.2) CFT=2.D0*XPI*TAMP1 170 : C 171 : MPOVA2.VPOCHA(JJ2,1)=AIRE*CFT 172 : C 173 : C 174 : 10 CONTINUE 175 : C 176 : SEGDES IPT1 177 : 1 CONTINUE 178 : C********************** 179 : C CAS DE LA DIMENSION 3 180 : C********************** 181 : ELSE 182 : C 183 : DO 21 L=1,NBSOUS 184 : IPT1=MELEME 185 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 186 : SEGACT IPT1 187 : NP=IPT1.NUM(/1)-1 188 : NEL=IPT1.NUM(/2) 189 : C Les elts complets ne sont pas implementes, on utilise 190 : C l'elt. face sans le centre 191 : C 192 : TYPE = NOMS(IPT1.ITYPEL)//' ' 192 : CALL OPTLI(IP,LIST1,TYPE,NBO) 192 : IF (IP .EQ. 0) CALL ERREUR(5) 193 : TYPE = LIST2(IP) 194 : C Calcul des fonctions de forme sur l'elt. de reference 195 : C 195 : CALL KALPBG(TYPE,'FONFORM0',IZFFM) 196 : SEGACT IZFFM*MOD 197 : IZHR=KZHR(1) 198 : SEGACT IZHR*MOD 199 : NES=GR(/1) 200 : NPG=GR(/3) 201 : C Boucle sur les faces du paquet L 202 : DO 210 K=1,NEL 203 : NF=LECT(IPT1.NUM(NP+1,K)) 204 : C 205 : C REMPLISSAGE DE XYZ 206 : C ------------------ 207 : C 208 : DO 212 I=1,NP 209 : J=IPT1.NUM(I,K) 210 : DO 212 N=1,IDIM 211 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 212 : C PRINT*,XY3(N,I) 213 : 212 CONTINUE 214 : C 215 : C calcul du jacobien idem dim 2 216 : C 216 : CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP, 217 : & NPG,IAXI,AIRE,A3J,SGEN) 218 : NORMX=A3J(1,3,1) 219 : NORMY=A3J(2,3,1) 220 : NORMZ=A3J(3,3,1) 221 : C 222 : C On verifie que n est dans le meme sens que vec(13) de FACEL 223 : C Calcul de vec(13) 224 : C FACEL est repere par IPT2 225 : J1=IPT2.NUM(1,NF) 226 : J2=IPT2.NUM(2,NF) 227 : JJ2=LECT(J2) 228 : J3=IPT2.NUM(3,NF) 229 : C 230 : IF (J1.eq.J3) THEN 231 : X1=XCOOR((J1-1)*(IDIM+1)+1) 232 : Y1=XCOOR((J1-1)*(IDIM+1)+2) 233 : Z1=XCOOR((J1-1)*(IDIM+1)+3) 234 : X2=XCOOR((J2-1)*(IDIM+1)+1) 235 : Y2=XCOOR((J2-1)*(IDIM+1)+2) 236 : Z2=XCOOR((J2-1)*(IDIM+1)+3) 237 : C 238 : SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY+(Z2-Z1)*NORMZ 239 : SGN=SIGN(1.D0,SCAL) 240 : C 241 : ELSE 242 : X1=XCOOR((J1-1)*(IDIM+1)+1) 243 : Y1=XCOOR((J1-1)*(IDIM+1)+2) 244 : Z1=XCOOR((J1-1)*(IDIM+1)+3) 245 : X3=XCOOR((J3-1)*(IDIM+1)+1) 246 : Y3=XCOOR((J3-1)*(IDIM+1)+2) 247 : Z3=XCOOR((J3-1)*(IDIM+1)+3) 248 : C 249 : SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY+(Z3-Z1)*NORMZ 250 : SGN=SIGN(1.D0,SCAL) 251 : ENDIF 252 : C 253 : MPOVAL.VPOCHA(JJ2,1)=SGN*NORMX 254 : MPOVAL.VPOCHA(JJ2,2)=SGN*NORMY 255 : MPOVAL.VPOCHA(JJ2,3)=SGN*NORMZ 256 : C 257 : C Calcul de la surface de la face 258 : C 259 : MPOVA2.VPOCHA(JJ2,1)=AIRE 260 : 210 CONTINUE 261 : C 262 : SEGDES IPT1 263 : 21 CONTINUE 264 : C 265 : ENDIF 266 : C 267 : C Orientation des normales (partie non modifiee) 268 : C 269 : SEGDES MELEMD,MELEME 270 : SEGDES MCHPOI,MPOVAL,MCHPO2,MPOVA2 271 : C 272 : TYPE=' ' 272 : CALL ACMO(MTABD,'ELTFA',TYPE,MELEMT) 273 : IF(MELEMT.EQ.0)GO TO 90 274 : 275 : TYPE=' ' 275 : CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC) 276 : IF(MELEMC.EQ.0)GO TO 90 277 : SEGACT MELEMC 278 : 279 : TYPE=' ' 279 : CALL ACMO(MTABD,'FACEL',TYPE,MELEMF) 280 : IF(MELEMF.EQ.0)GO TO 90 281 : SEGACT MELEMF 281 : CALL CRCHPE(MELEMT,1,MCHELM) 282 : SEGACT MCHELM 283 : NBSOUS=IMACHE(/1) 284 : C 285 : K0=0 286 : DO 41 L=1,NBSOUS 287 : MCHAML=ICHAML(L) 288 : SEGACT MCHAML 289 : MELVAL=IELVAL(1) 290 : SEGACT MELVAL*MOD 291 : IPT1=IMACHE(L) 292 : SEGACT IPT1 293 : NF=IPT1.NUM(/1) 294 : NEL=IPT1.NUM(/2) 295 : DO 42 K=1,NEL 296 : K0=K0+1 297 : NE=MELEMC.NUM(1,K0) 298 : DO 42 I=1,NF 299 : KF=IPT1.NUM(I,K) 300 : KF=LECT(KF) 301 : N1=MELEMF.NUM(1,KF) 302 : N3=MELEMF.NUM(3,KF) 303 : C 304 : IF(NE.EQ.N1)THEN 305 : C La normale est sortante 306 : A=1.D0 307 : ELSEIF(NE.EQ.N3)THEN 308 : C La normale est entrante 309 : A=-1.D0 310 : ELSE 311 : WRITE(6,*)'Problemes dans KNRF NE,N1,N3=',NE,N1,N3 312 : A=0. 313 : ENDIF 314 : C 315 : VELCHE(I,K)=A 316 : C 317 : 42 CONTINUE 318 : SEGDES MELVAL,MCHAML 319 : SEGDES IPT1 320 : 41 CONTINUE 321 : SEGSUP MLENTI 322 : SEGDES MCHELM 323 : SEGDES MELEMC,MELEMF 324 : C Normales aux faces par element 325 : C write(6,*)' faces elt Normales aire',MCHELM,MCHPOI,MCHPO2 325 : C CALL ECROBJ('MCHAML ',MCHELM) 326 : C Normales aux faces 326 : C CALL ECROBJ('CHPOINT ',MCHPOI) 327 : C aire des faces 327 : C CALL ECROBJ('CHPOINT ',MCHPO2) 328 : RETURN 329 : C 330 : 90 CONTINUE 331 : WRITE(6,*)' Interruption anormale de KNRF' 332 : RETURN 333 : END 334 : 335 : 336 : 337 : 338 : 339 : 340 : 341 : 342 : 343 : 344 : 345 : 346 : 347 : 348 : 349 : 350 : 351 :
© Cast3M 2003 - All rights reserved.
Disclaimer