C KNRF SOURCE GOUNAND 25/11/12 21:15:31 12399 SUBROUTINE KNRF(MTABD,MCHELM,MCHPOI,MCHPO2) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C C************************************************************************* C C OBJET : Cree deux CHAMPOINT faces et un mchaml C CHPC1 contenant la longueur des faces C ( faux dans le cas 3d ) C CHPC2 contenant la normale aux faces. C MCHAML orientation (normale entrante ou sortante) C SYNTAXE : CHPC1 CHPC2 MCHAML = KNRF OBJDOM ; C C OBJDOM : TABLE de SOUSTYPE DOMAINE C C REMARQUE : C On verifie que la normale est dans le meme sens que pt1->pt3 C des elements FACEL C Si pt1=pt3 on oriente la normale vers l'exterieur cad pt1->pt2 C C Non teste si axisymetrie C************************************************************************* C MODIFICATION du 20/11/98 C C Le calcul de la normale s'effectue a l'aide du jacobien calcule C par CALJBR (au lieu de CALJQB) C C************************************************************************* -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC SMTABLE POINTEUR MTABD.MTABLE -INC SMELEME POINTEUR MELEMC.MELEME,MELEMD.MELEME,MELEMT.MELEME POINTEUR MELEMF.MELEME -INC SMCOORD -INC SMCHPOI -INC SMCHAML -INC SMLENTI -INC SIZFFB -INC CCREEL C PARAMETER (NBO=5) REAL*8 A2J(2,2,1),A3J(3,3,1),CFT,NORMX,NORMY,NORMZ REAL*8 LGR,LGR1,LGR2 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 C Acquisition des connectivites FACEL C write(6,*)' KNRF ' TYPE=' ' CALL ACMO(MTABD,'FACEL',TYPE,MELEMD) C write(6,*)' FACEL type=',type IF(TYPE.NE.'MAILLAGE')GO TO 90 SEGACT MELEMD IPT2=MELEMD C C Acquisition des connectivites FACE C TYPE=' ' CALL ACMO(MTABD,'FACE',TYPE,MELEMF) IF(TYPE.NE.'MAILLAGE')GO TO 90 CALL KRIPAD(MELEMF,MLENTI) SEGACT MELEMF C C Creation du CHPOIN a 3 comp. pour la normale aux faces C NC=IDIM TYPE='FACE' CALL CRCHPT(TYPE,MELEMF,NC,1,MCHPOI) CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM) C C Creation du CHPOIN pour la longueur C NC=1 TYPE='FACE' CALL CRCHPT(TYPE,MELEMF,NC,2,MCHPO2) CALL LICHT(MCHPO2,MPOVA2,TYPC,IGEOM) SEGDES MELEMF C IAXI = 0 IF(IFOMOD.EQ.0)IAXI=2 C C Acquisition des connectivites FACEP C TYPE=' ' CALL ACMO(MTABD,'FACEP',TYPE,MELEME) IF(TYPE.NE.'MAILLAGE')GO TO 90 SEGACT MELEME C NBSOUS=LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 C********************** C CAS DE LA DIMENSION 2 C********************** IF (IDIM.EQ.2) THEN C 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 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=GR(/1) NPG=GR(/3) 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) 7 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)=SGN*NORMX MPOVAL.VPOCHA(JJ2,2)=SGN*NORMY C Calcul de la longueur C CFT=1.D0 IF(IAXI.NE.0)THEN TAMP1=(XYZ(1,2)+XYZ(1,1))*0.5D0 TAMP2=(XYZ(2,2)+XYZ(2,1))*0.5D0 ENDIF C Axisymetrie / ox (non teste) IF (IAXI.EQ.1) CFT=2.D0*XPI*TAMP2 C Axisymetrie / oy (non teste non plus) IF (IAXI.EQ.2) CFT=2.D0*XPI*TAMP1 C MPOVA2.VPOCHA(JJ2,1)=AIRE*CFT C C 10 CONTINUE C SEGDES IPT1 1 CONTINUE C********************** C CAS DE LA DIMENSION 3 C********************** ELSE C 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=GR(/1) NPG=GR(/3) 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 DO 212 I=1,NP J=IPT1.NUM(I,K) DO 2121 N=1,IDIM XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N) C PRINT*,XY3(N,I) 2121 CONTINUE 212 CONTINUE 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) C SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY+(Z3-Z1)*NORMZ SGN=SIGN(1.D0,SCAL) ENDIF C MPOVAL.VPOCHA(JJ2,1)=SGN*NORMX MPOVAL.VPOCHA(JJ2,2)=SGN*NORMY MPOVAL.VPOCHA(JJ2,3)=SGN*NORMZ C C Calcul de la surface de la face C MPOVA2.VPOCHA(JJ2,1)=AIRE 210 CONTINUE C SEGDES IPT1 21 CONTINUE C ENDIF C C Orientation des normales (partie non modifiee) C SEGDES MELEMD,MELEME SEGDES MCHPOI,MPOVAL,MCHPO2,MPOVA2 C TYPE=' ' CALL ACMO(MTABD,'ELTFA',TYPE,MELEMT) IF(MELEMT.EQ.0)GO TO 90 TYPE=' ' CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC) IF(MELEMC.EQ.0)GO TO 90 SEGACT MELEMC TYPE=' ' CALL ACMO(MTABD,'FACEL',TYPE,MELEMF) IF(MELEMF.EQ.0)GO TO 90 SEGACT MELEMF CALL CRCHPE(MELEMT,1,MCHELM) SEGACT MCHELM NBSOUS=IMACHE(/1) C K0=0 DO 41 L=1,NBSOUS MCHAML=ICHAML(L) SEGACT MCHAML MELVAL=IELVAL(1) SEGACT MELVAL*MOD IPT1=IMACHE(L) SEGACT IPT1 NF=IPT1.NUM(/1) NEL=IPT1.NUM(/2) DO 42 K=1,NEL K0=K0+1 NE=MELEMC.NUM(1,K0) DO 421 I=1,NF KF=IPT1.NUM(I,K) KF=LECT(KF) N1=MELEMF.NUM(1,KF) N3=MELEMF.NUM(3,KF) C IF(NE.EQ.N1)THEN C La normale est sortante A=1.D0 ELSEIF(NE.EQ.N3)THEN C La normale est entrante A=-1.D0 ELSE WRITE(6,*)'Problemes dans KNRF NE,N1,N3=',NE,N1,N3 A=0. ENDIF C VELCHE(I,K)=A C 421 CONTINUE 42 CONTINUE SEGDES MELVAL,MCHAML SEGDES IPT1 41 CONTINUE SEGSUP MLENTI SEGDES MCHELM SEGDES MELEMC,MELEMF C Normales aux faces par element C write(6,*)' faces elt Normales aire',MCHELM,MCHPOI,MCHPO2 C CALL ECROBJ('MCHAML ',MCHELM) C Normales aux faces C CALL ECROBJ('CHPOINT ',MCHPOI) C aire des faces C CALL ECROBJ('CHPOINT ',MCHPO2) RETURN C 90 CONTINUE WRITE(6,*)' Interruption anormale de KNRF' RETURN END