Sources Esope | Cast3M

knrf
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