C MATRAN SOURCE GOUNAND 25/11/12 21:15:40 12399 SUBROUTINE MATRAN(MTABD,MCHPOI) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C------------------------------------------------------------------ C PROJET : CASTEM 2000 C C NOM : MATRAN C C DESCRIPTION : Cette subroutine cree un chpoint 'FACE' sur MTABD C contenant la matrice de passage du repaire global C dans le repaire de chaque face (voir CALJQB) C C LANGAGE : FORTRAN 77 + ESOPE 2000 C C AUTEUR : R. MOREL, DRN/DMT/SEMT/TTMF C C----------------------------------------------------------------- C C APPELES (E/S) : LEKTAB C C APPELES (Calcul): CALJBR C C----------------------------------------------------------------- C C ENTREES : C MTABD : Table de sous-type domaine C C SORTIES : C MCHPOI : Pointeur du CHPOIN contenant les C Matrices decrite ci-dessus. C C----------------------------------------------------------------- C C HISTORIQUE (Anomalies et modifications eventuelles) C C HISTORIQUE : 23.10.98, Creation C C----------------------------------------------------------------- C -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC SMTABLE POINTEUR MTABD.MTABLE -INC SMELEME POINTEUR MELEMC.MELEME,MELEMD.MELEME POINTEUR MELEMF.MELEME -INC SMCOORD -INC SMCHPOI -INC SMLENTI -INC SIZFFB PARAMETER (NBO=5) REAL*8 CFT,NORMX,NORMY,NORMZ,A2J(2,2,1),T,TX,TY,TZ,AIRE REAL*8 X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3,SCAL,SGN,A3J(3,3,1) INTEGER NP,NPG,NES CHARACTER*8 TYPE,TYPC,NOM0,LIST1(NBO),LIST2(NBO) C Les elements complets ne sont pas tous implementes C On se ramene aux elements facep sans le centre de la face DATA LIST1/'SEG3 ','TRI4 ','TRI7 ','QUA5 ', & 'QUA9 '/ DATA LIST2/'SEG2 ','TRI3 ','TRI6 ','QUA4 ', & 'QUA8 '/ C CALL LEKTAB(MTABD,'FACEL',MELEMD) SEGACT MELEMD IPT2=MELEMD CALL LEKTAB(MTABD,'FACE',MELEMF) CALL KRIPAD(MELEMF,MLENTI) SEGACT MELEMF NC=IDIM*IDIM TYPE='FACE' CALL CRCHPT(TYPE,MELEMF,NC,1,MCHPOI) CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) CALL LEKTAB(MTABD,'FACEP',MELEME) SEGACT MELEME C SEGACT MCHPOI*MOD MSOUPO=MCHPOI.IPCHP(1) SEGACT MSOUPO*MOD NBSOUS=LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 IF (IDIM.EQ.2) THEN C NOCOMP(1)='UX ' NOCOMP(2)='UY ' NOCOMP(3)='RX ' NOCOMP(4)='RY ' DO 1 L=1,NBSOUS IPT1=MELEME IF(NBSOUS.NE.1)IPT1=LISOUS(L) SEGACT IPT1 NP=IPT1.NUM(/1)-1 NEL=IPT1.NUM(/2) C Toutes les faces sont des SEG3 donc NP=2 IF(NP.NE.2)RETURN C C C On considere la face sans son centre comme un elt fini SEG2 C TYPE = NOMS(IPT1.ITYPEL)//' ' CALL OPTLI(IP,LIST1,TYPE,NBO) IF (IP .EQ. 0) CALL ERREUR(5) TYPE = LIST2(IP) CALL KALPBG(TYPE,'FONFORM0',IZFFM) SEGACT IZFFM*MOD IZHR=KZHR(1) SEGACT IZHR*MOD NES=1 NPG=1 C Boucle sur toutes les faces pour le paquet L DO 10 K=1,NEL NF=LECT(IPT1.NUM(NP+1,K)) C REMPLISSAGE DE XYZ C ------------------ C DO 12 I=1,NP J=IPT1.NUM(I,K) DO 121 N=1,IDIM XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 121 CONTINUE 12 CONTINUE C C CALJBR calcul le jacobien du passage de l'elt de ref. C a l'elt. reel C A2J=Jacobien AIRE=detA2J CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP, &NPG,IAXI,AIRE,A2J,SGEN) C NORMX=A2J(1,2,1) NORMY=A2J(2,2,1) C On verifie que n est dans le meme sens que vec(13) de FACEL C Calcul de vec(13) C FACEL est repere par IPT2 J1=IPT2.NUM(1,NF) J2=IPT2.NUM(2,NF) JJ2=LECT(J2) J3=IPT2.NUM(3,NF) C IF (J1.eq.J3) THEN X1=XCOOR((J1-1)*(IDIM+1)+1) Y1=XCOOR((J1-1)*(IDIM+1)+2) X2=XCOOR((J2-1)*(IDIM+1)+1) Y2=XCOOR((J2-1)*(IDIM+1)+2) SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY SGN=SIGN(1.D0,SCAL) C ELSE X1=XCOOR((J1-1)*(IDIM+1)+1) Y1=XCOOR((J1-1)*(IDIM+1)+2) X3=XCOOR((J3-1)*(IDIM+1)+1) Y3=XCOOR((J3-1)*(IDIM+1)+2) SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY SGN=SIGN(1.D0,SCAL) ENDIF MPOVAL.VPOCHA(JJ2,1)=A2J(1,1,1) MPOVAL.VPOCHA(JJ2,2)=A2J(2,1,1) MPOVAL.VPOCHA(JJ2,3)=SGN*NORMX MPOVAL.VPOCHA(JJ2,4)=SGN*NORMY C C 10 CONTINUE C SEGDES IPT1 1 CONTINUE C C ELSE C NOCOMP(1)='UX ' NOCOMP(2)='UY ' NOCOMP(3)='UZ ' NOCOMP(4)='RX ' NOCOMP(5)='RY ' NOCOMP(6)='RZ ' NOCOMP(7)='MX ' NOCOMP(8)='MY ' NOCOMP(9)='MZ ' DO 21 L=1,NBSOUS IPT1=MELEME IF(NBSOUS.NE.1)IPT1=LISOUS(L) SEGACT IPT1 NP=IPT1.NUM(/1)-1 NEL=IPT1.NUM(/2) C Les elts complets ne sont pas implementes, on utilise C l'elt. face sans le centre C TYPE = NOMS(IPT1.ITYPEL)//' ' CALL OPTLI(IP,LIST1,TYPE,NBO) IF (IP .EQ. 0) CALL ERREUR(5) TYPE = LIST2(IP) C Calcul des fonctions de forme sur l'elt. de reference C CALL KALPBG(TYPE,'FONFORM0',IZFFM) SEGACT IZFFM*MOD IZHR=KZHR(1) SEGACT IZHR*MOD NES=2 NPG=1 C Boucle sur les faces du paquet L DO 210 K=1,NEL NF=LECT(IPT1.NUM(NP+1,K)) C C REMPLISSAGE DE XYZ C ------------------ C C write(6,*)'num',(ipt1.num(ii,k),ii=1,np+1),' NF=',nf DO 212 I=1,NP J=IPT1.NUM(I,K) DO 2121 N=1,IDIM XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) 2121 CONTINUE 212 CONTINUE C C C calcul du jacobien idem dim 2 C CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP, & NPG,IAXI,AIRE,A3J,SGEN) NORMX=A3J(1,3,1) NORMY=A3J(2,3,1) NORMZ=A3J(3,3,1) C C On verifie que n est dans le meme sens que vec(13) de FACEL C Calcul de vec(13) C FACEL est repere par IPT2 J1=IPT2.NUM(1,NF) J2=IPT2.NUM(2,NF) JJ2=LECT(J2) J3=IPT2.NUM(3,NF) C IF (J1.eq.J3) THEN X1=XCOOR((J1-1)*(IDIM+1)+1) Y1=XCOOR((J1-1)*(IDIM+1)+2) Z1=XCOOR((J1-1)*(IDIM+1)+3) X2=XCOOR((J2-1)*(IDIM+1)+1) Y2=XCOOR((J2-1)*(IDIM+1)+2) Z2=XCOOR((J2-1)*(IDIM+1)+3) C SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY+(Z2-Z1)*NORMZ SGN=SIGN(1.D0,SCAL) C ELSE X1=XCOOR((J1-1)*(IDIM+1)+1) Y1=XCOOR((J1-1)*(IDIM+1)+2) Z1=XCOOR((J1-1)*(IDIM+1)+3) X3=XCOOR((J3-1)*(IDIM+1)+1) Y3=XCOOR((J3-1)*(IDIM+1)+2) Z3=XCOOR((J3-1)*(IDIM+1)+3) SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY+(Z3-Z1)*NORMZ SGN=SIGN(1.D0,SCAL) ENDIF C TX=A3J(1,1,1) TY=A3J(2,1,1) TZ=A3J(3,1,1) T=SQRT(TX*TX+TY*TY+TZ*TZ) TX=TX / T TY=TY / T TZ=TZ / T MPOVAL.VPOCHA(JJ2,1)=TX MPOVAL.VPOCHA(JJ2,2)=TY MPOVAL.VPOCHA(JJ2,3)=TZ MPOVAL.VPOCHA(JJ2,4)=NORMY*TZ-NORMZ*TY MPOVAL.VPOCHA(JJ2,5)=NORMZ*TX-NORMX*TZ MPOVAL.VPOCHA(JJ2,6)=NORMX*TY-NORMY*TX MPOVAL.VPOCHA(JJ2,7)=SGN*NORMX MPOVAL.VPOCHA(JJ2,8)=SGN*NORMY MPOVAL.VPOCHA(JJ2,9)=SGN*NORMZ C 210 CONTINUE C SEGDES IPT1 21 CONTINUE ENDIF CALL ECMO(MTABD,'MATROT','CHPOINT ',MCHPOI) SEGSUP MLENTI SEGDES MELEMD SEGDES MELEME SEGDES MELEMF SEGDES MCHPOI SEGDES MSOUPO SEGDES MPOVAL RETURN END