1 : C MATRAN SOURCE MAGN 10/05/18 21:16:19 6675 2 : SUBROUTINE MATRAN(MTABD,MCHPOI) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : C------------------------------------------------------------------ 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : MATRAN 9 : C 10 : C DESCRIPTION : Cette subroutine cree un chpoint 'FACE' sur MTABD 11 : C contenant la matrice de passage du repaire global 12 : C dans le repaire de chaque face (voir CALJQB) 13 : C 14 : C LANGAGE : FORTRAN 77 + ESOPE 2000 15 : C 16 : C AUTEUR : R. MOREL, DRN/DMT/SEMT/TTMF 17 : C 18 : C----------------------------------------------------------------- 19 : C 20 : C APPELES (E/S) : LEKTAB 21 : C 22 : C APPELES (Calcul): CALJBR 23 : C 24 : C----------------------------------------------------------------- 25 : C 26 : C ENTREES : 27 : C MTABD : Table de sous-type domaine 28 : C 29 : C SORTIES : 30 : C MCHPOI : Pointeur du CHPOIN contenant les 31 : C Matrices decrite ci-dessus. 32 : C 33 : C----------------------------------------------------------------- 34 : C 35 : C HISTORIQUE (Anomalies et modifications eventuelles) 36 : C 37 : C HISTORIQUE : 23.10.98, Creation 38 : C 39 : C----------------------------------------------------------------- 40 : C 41 : -INC CCOPTIO 42 : -INC CCGEOME 43 : -INC SMTABLE 44 : POINTEUR MTABD.MTABLE 45 : -INC SMELEME 46 : POINTEUR MELEMC.MELEME,MELEMD.MELEME 47 : POINTEUR MELEMF.MELEME 48 : -INC SMCOORD 49 : -INC SMCHPOI 50 : -INC SMLENTI 51 : -INC SIZFFB 52 : PARAMETER (NBO=5) 53 : REAL*8 CFT,NORMX,NORMY,NORMZ,A2J(2,2,1),T,TX,TY,TZ,AIRE 54 : REAL*8 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,SCAL,SGN,A3J(3,3,1) 55 : INTEGER NP,NPG,NES 56 : CHARACTER*8 TYPE,TYPC,NOM0,LIST1(NBO),LIST2(NBO) 57 : C Les elements complets ne sont pas tous implementes 58 : C On se ramene aux elements facep sans le centre de la face 59 : DATA LIST1/'SEG3 ','TRI4 ','TRI7 ','QUA5 ', 60 : & 'QUA9 '/ 61 : DATA LIST2/'SEG2 ','TRI3 ','TRI6 ','QUA4 ', 62 : & 'QUA8 '/ 63 : C 63 : CALL LEKTAB(MTABD,'FACEL',MELEMD) 64 : SEGACT MELEMD 65 : IPT2=MELEMD 65 : CALL LEKTAB(MTABD,'FACE',MELEMF) 65 : CALL KRIPAD(MELEMF,MLENTI) 66 : SEGACT MELEMF 67 : NC=IDIM*IDIM 68 : TYPE='FACE' 68 : CALL CRCHPT(TYPE,MELEMF,NC,MCHPOI) 68 : CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) 68 : CALL LEKTAB(MTABD,'FACEP',MELEME) 69 : SEGACT MELEME 70 : C 71 : SEGACT MCHPOI*MOD 72 : MSOUPO=MCHPOI.IPCHP(1) 73 : SEGACT MSOUPO*MOD 74 : NBSOUS=LISOUS(/1) 75 : IF(NBSOUS.EQ.0)NBSOUS=1 76 : IF (IDIM.EQ.2) THEN 77 : C 78 : NOCOMP(1)='UX ' 79 : NOCOMP(2)='UY ' 80 : NOCOMP(3)='RX ' 81 : NOCOMP(4)='RY ' 82 : DO 1 L=1,NBSOUS 83 : IPT1=MELEME 84 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 85 : SEGACT IPT1 86 : NP=IPT1.NUM(/1)-1 87 : NEL=IPT1.NUM(/2) 88 : C Toutes les faces sont des SEG3 donc NP=2 89 : IF(NP.NE.2)RETURN 90 : C 91 : C 92 : C On considere la face sans son centre comme un elt fini SEG2 93 : C 94 : TYPE = NOMS(IPT1.ITYPEL)//' ' 94 : CALL OPTLI(IP,LIST1,TYPE,NBO) 94 : IF (IP .EQ. 0) CALL ERREUR(5) 95 : TYPE = LIST2(IP) 95 : CALL KALPBG(TYPE,'FONFORM0',IZFFM) 96 : SEGACT IZFFM*MOD 97 : IZHR=KZHR(1) 98 : SEGACT IZHR*MOD 99 : NES=1 100 : NPG=1 101 : C Boucle sur toutes les faces pour le paquet L 102 : DO 10 K=1,NEL 103 : NF=LECT(IPT1.NUM(NP+1,K)) 104 : C REMPLISSAGE DE XYZ 105 : C ------------------ 106 : C 107 : DO 12 I=1,NP 108 : J=IPT1.NUM(I,K) 109 : DO 12 N=1,IDIM 110 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 111 : 12 CONTINUE 112 : C 113 : C CALJBR calcul le jacobien du passage de l'elt de ref. 114 : C a l'elt. reel 115 : C A2J=Jacobien AIRE=detA2J 115 : CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP, 116 : &NPG,IAXI,AIRE,A2J,SGEN) 117 : C 118 : NORMX=A2J(1,2,1) 119 : NORMY=A2J(2,2,1) 120 : C On verifie que n est dans le meme sens que vec(13) de FACEL 121 : C Calcul de vec(13) 122 : C FACEL est repere par IPT2 123 : J1=IPT2.NUM(1,NF) 124 : J2=IPT2.NUM(2,NF) 125 : JJ2=LECT(J2) 126 : J3=IPT2.NUM(3,NF) 127 : C 128 : IF (J1.eq.J3) THEN 129 : X1=XCOOR((J1-1)*(IDIM+1)+1) 130 : Y1=XCOOR((J1-1)*(IDIM+1)+2) 131 : X2=XCOOR((J2-1)*(IDIM+1)+1) 132 : Y2=XCOOR((J2-1)*(IDIM+1)+2) 133 : SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY 134 : SGN=SIGN(1.D0,SCAL) 135 : C 136 : ELSE 137 : X1=XCOOR((J1-1)*(IDIM+1)+1) 138 : Y1=XCOOR((J1-1)*(IDIM+1)+2) 139 : X3=XCOOR((J3-1)*(IDIM+1)+1) 140 : Y3=XCOOR((J3-1)*(IDIM+1)+2) 141 : SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY 142 : SGN=SIGN(1.D0,SCAL) 143 : ENDIF 144 : MPOVAL.VPOCHA(JJ2,1)=A2J(1,1,1) 145 : MPOVAL.VPOCHA(JJ2,2)=A2J(2,1,1) 146 : MPOVAL.VPOCHA(JJ2,3)=SGN*NORMX 147 : MPOVAL.VPOCHA(JJ2,4)=SGN*NORMY 148 : C 149 : C 150 : 10 CONTINUE 151 : C 152 : SEGDES IPT1 153 : 1 CONTINUE 154 : C 155 : C 156 : ELSE 157 : C 158 : NOCOMP(1)='UX ' 159 : NOCOMP(2)='UY ' 160 : NOCOMP(3)='UZ ' 161 : NOCOMP(4)='RX ' 162 : NOCOMP(5)='RY ' 163 : NOCOMP(6)='RZ ' 164 : NOCOMP(7)='MX ' 165 : NOCOMP(8)='MY ' 166 : NOCOMP(9)='MZ ' 167 : DO 21 L=1,NBSOUS 168 : IPT1=MELEME 169 : IF(NBSOUS.NE.1)IPT1=LISOUS(L) 170 : SEGACT IPT1 171 : NP=IPT1.NUM(/1)-1 172 : NEL=IPT1.NUM(/2) 173 : C Les elts complets ne sont pas implementes, on utilise 174 : C l'elt. face sans le centre 175 : C 176 : TYPE = NOMS(IPT1.ITYPEL)//' ' 176 : CALL OPTLI(IP,LIST1,TYPE,NBO) 176 : IF (IP .EQ. 0) CALL ERREUR(5) 177 : TYPE = LIST2(IP) 178 : C Calcul des fonctions de forme sur l'elt. de reference 179 : C 179 : CALL KALPBG(TYPE,'FONFORM0',IZFFM) 180 : SEGACT IZFFM*MOD 181 : IZHR=KZHR(1) 182 : SEGACT IZHR*MOD 183 : NES=2 184 : NPG=1 185 : C Boucle sur les faces du paquet L 186 : DO 210 K=1,NEL 187 : NF=LECT(IPT1.NUM(NP+1,K)) 188 : C 189 : C REMPLISSAGE DE XYZ 190 : C ------------------ 191 : C 192 : C write(6,*)'num',(ipt1.num(ii,k),ii=1,np+1),' NF=',nf 193 : DO 212 I=1,NP 194 : J=IPT1.NUM(I,K) 195 : DO 212 N=1,IDIM 196 : XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 197 : 212 CONTINUE 198 : C 199 : C 200 : C calcul du jacobien idem dim 2 201 : C 201 : CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP, 202 : & NPG,IAXI,AIRE,A3J,SGEN) 203 : NORMX=A3J(1,3,1) 204 : NORMY=A3J(2,3,1) 205 : NORMZ=A3J(3,3,1) 206 : C 207 : C On verifie que n est dans le meme sens que vec(13) de FACEL 208 : C Calcul de vec(13) 209 : C FACEL est repere par IPT2 210 : J1=IPT2.NUM(1,NF) 211 : J2=IPT2.NUM(2,NF) 212 : JJ2=LECT(J2) 213 : J3=IPT2.NUM(3,NF) 214 : C 215 : IF (J1.eq.J3) THEN 216 : X1=XCOOR((J1-1)*(IDIM+1)+1) 217 : Y1=XCOOR((J1-1)*(IDIM+1)+2) 218 : Z1=XCOOR((J1-1)*(IDIM+1)+3) 219 : X2=XCOOR((J2-1)*(IDIM+1)+1) 220 : Y2=XCOOR((J2-1)*(IDIM+1)+2) 221 : Z2=XCOOR((J2-1)*(IDIM+1)+3) 222 : C 223 : SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY+(Z2-Z1)*NORMZ 224 : SGN=SIGN(1.D0,SCAL) 225 : C 226 : ELSE 227 : X1=XCOOR((J1-1)*(IDIM+1)+1) 228 : Y1=XCOOR((J1-1)*(IDIM+1)+2) 229 : Z1=XCOOR((J1-1)*(IDIM+1)+3) 230 : X3=XCOOR((J3-1)*(IDIM+1)+1) 231 : Y3=XCOOR((J3-1)*(IDIM+1)+2) 232 : Z3=XCOOR((J3-1)*(IDIM+1)+3) 233 : SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY+(Z3-Z1)*NORMZ 234 : SGN=SIGN(1.D0,SCAL) 235 : ENDIF 236 : C 237 : TX=A3J(1,1,1) 238 : TY=A3J(2,1,1) 239 : TZ=A3J(3,1,1) 240 : T=SQRT(TX*TX+TY*TY+TZ*TZ) 241 : TX=TX / T 242 : TY=TY / T 243 : TZ=TZ / T 244 : MPOVAL.VPOCHA(JJ2,1)=TX 245 : MPOVAL.VPOCHA(JJ2,2)=TY 246 : MPOVAL.VPOCHA(JJ2,3)=TZ 247 : MPOVAL.VPOCHA(JJ2,4)=NORMY*TZ-NORMZ*TY 248 : MPOVAL.VPOCHA(JJ2,5)=NORMZ*TX-NORMX*TZ 249 : MPOVAL.VPOCHA(JJ2,6)=NORMX*TY-NORMY*TX 250 : MPOVAL.VPOCHA(JJ2,7)=SGN*NORMX 251 : MPOVAL.VPOCHA(JJ2,8)=SGN*NORMY 252 : MPOVAL.VPOCHA(JJ2,9)=SGN*NORMZ 253 : C 254 : 210 CONTINUE 255 : C 256 : SEGDES IPT1 257 : 21 CONTINUE 258 : ENDIF 258 : CALL ECMO(MTABD,'MATROT','CHPOINT ',MCHPOI) 259 : SEGSUP MLENTI 260 : SEGDES MELEMD 261 : SEGDES MELEME 262 : SEGDES MELEMF 263 : SEGDES MCHPOI 264 : SEGDES MSOUPO 265 : SEGDES MPOVAL 266 : RETURN 267 : END 268 : 269 : 270 : 271 : 272 : 273 : 274 : 275 :
© Cast3M 2003 - All rights reserved.
Disclaimer