C KDOM11 SOURCE PV 21/10/31 21:15:01 11158 SUBROUTINE KDOM11(MELF,MELFL,MELFP,MCHPSU,MCHPNO,MCHPMR) C C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : KDOM11 C C DESCRIPTION : Subroutine called by KDOM10 in the case of EULER C model. C We create the CHPOINTS C MCHPSU = surfaces (faces dimension) C MCHPNO = normals (oriented from the first C to the second point of FACEL) C MCHPMR = rotation matrix C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF C C************************************************************************ C C INPUTS : C C MELF : meleme 'FACE' C MELFL : meleme 'FACEL' C MELFP : meleme 'FACEP' C C OUTPUTS : C C MCHPSU : mchpoi 'XXSURFAC' C MCHPNO : mchpoi 'XXNORMAF' C MCHPMR : mchpoi 'MATROT' C C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) C INTEGER IGEOM, MCHPSU, MCHPNO, MCHPMR & ,JGN, JGM, NBSOUS, ISOUS, NP, NEL, IEL, IPO, NLCF & , NF, IFAC REAL*8 XF,YF,X1,Y1, DSURFX, DSURFY, XSURF, YSURF, SURF, A & ,ORIENT & ,ZF, DXP, DYP, DZP, DXPM1, DYPM1, DZPM1, DSURFZ, ZSURF & ,RMX, RMY, RMZ, RRX, RRY, RRZ CHARACTER*8 TYPE -INC PPARAM -INC CCOPTIO -INC SMELEME POINTEUR MELFL.MELEME,MELFP.MELEME,MELF.MELEME -INC SMCHPOI POINTEUR MPOVSU.MPOVAL, MPOVNO.MPOVAL, MPOVMR.MPOVAL -INC SMLENTI -INC SMLMOTS -INC SMCOORD -INC CCREEL C C**** Corresp. FACE C CALL KRIPAD(MELF,MLENTI) C SEGINI MLENTI C C**** Champoint surfaces C JGN=4 JGM=1 SEGINI MLMOTS MLMOTS.MOTS(1)='SCAL' TYPE='FACE ' CALL KRCHP1(TYPE,MELF,MCHPSU,MLMOTS) IF(IERR.NE.0)GOTO 9999 CALL LICHT(MCHPSU,MPOVSU,TYPE,IGEOM) IF(IERR.NE.0)GOTO 9999 C SEGACT MPOVSU SEGSUP MLMOTS C C**** Champoint normales C JGN=4 JGM=IDIM SEGINI MLMOTS MLMOTS.MOTS(1)='UX' MLMOTS.MOTS(2)='UY' IF(IDIM .EQ. 3) MLMOTS.MOTS(3)='UZ' TYPE='FACE ' CALL KRCHP1(TYPE,MELF,MCHPNO,MLMOTS) IF(IERR.NE.0)GOTO 9999 CALL LICHT(MCHPNO,MPOVNO,TYPE,IGEOM) IF(IERR.NE.0)GOTO 9999 C SEGACT MPOVNO SEGSUP MLMOTS C C**** Champoint matrice de rotation C JGN=4 JGM=IDIM*IDIM SEGINI MLMOTS IF(IDIM.EQ.2)THEN MLMOTS.MOTS(1)='RX' MLMOTS.MOTS(2)='RY' MLMOTS.MOTS(3)='MX' MLMOTS.MOTS(4)='MY' * Normale en M * vect(M,U) = z ELSE MLMOTS.MOTS(1)='UX' MLMOTS.MOTS(2)='UY' MLMOTS.MOTS(3)='UZ' MLMOTS.MOTS(4)='RX' MLMOTS.MOTS(5)='RY' MLMOTS.MOTS(6)='RZ' MLMOTS.MOTS(7)='MX' MLMOTS.MOTS(8)='MY' MLMOTS.MOTS(9)='MZ' * Normale en M * vect(M,R) = U ENDIF CALL KRCHP1(TYPE,MELF,MCHPMR,MLMOTS) IF(IERR.NE.0)GOTO 9999 CALL LICHT(MCHPMR,MPOVMR,TYPE,IGEOM) IF(IERR.NE.0)GOTO 9999 C SEGACT MPOVMR C SEGACT MELFP NBSOUS=MELFP.LISOUS(/1) IF(NBSOUS.EQ.0)NBSOUS=1 SEGACT MELFL C********************** C CAS DE LA DIMENSION 2 C********************** surf=0.d0 xsurf=0.d0 ysurf=0.d0 zsurf=0.d0 IF (IDIM.EQ.2) THEN C DO ISOUS=1,NBSOUS,1 IPT1 = MELFP IF(NBSOUS .NE. 1)THEN IPT1 = MELFP.LISOUS(ISOUS) SEGACT IPT1 ENDIF C NP=IPT1.NUM(/1)-1 NEL=IPT1.NUM(/2) IF(NP .NE. 2)THEN WRITE(IOIMP,*) 'Subroutine knrf.eso' CALL ERREUR(5) ENDIF C DO IEL=1,NEL,1 NF=IPT1.NUM(NP+1,IEL) XF = XCOOR((NF-1)*(IDIM+1)+1) YF = XCOOR((NF-1)*(IDIM+1)+2) * IPO=IPT1.NUM(1,IEL) X1=XCOOR((IPO-1)*(IDIM+1)+1) Y1=XCOOR((IPO-1)*(IDIM+1)+2) DSURFX=Y1-YF DSURFY=XF-X1 XSURF=DSURFX YSURF=DSURFY SURF=((DSURFX*DSURFX)+(DSURFY*DSURFY))**0.5D0 * IPO=IPT1.NUM(2,IEL) X1=XCOOR((IPO-1)*(IDIM+1)+1) Y1=XCOOR((IPO-1)*(IDIM+1)+2) DSURFX=YF-Y1 DSURFY=X1-XF XSURF=XSURF+DSURFX YSURF=YSURF+DSURFY SURF=SURF+(((DSURFX*DSURFX)+(DSURFY*DSURFY))**0.5D0) C NLCF=MLENTI.LECT(NF) MPOVSU.VPOCHA(NLCF,1)=SURF C C************* Orientation selon FACEL C IPO=MELFL.NUM(1,NLCF) X1=XCOOR((IPO-1)*(IDIM+1)+1) Y1=XCOOR((IPO-1)*(IDIM+1)+2) DSURFX=XF-X1 DSURFY=YF-Y1 ORIENT=SIGN(1.0D0,((DSURFX*XSURF)+(DSURFY*YSURF))) C MPOVNO.VPOCHA(NLCF,1)=XSURF/SURF*ORIENT MPOVNO.VPOCHA(NLCF,2)=YSURF/SURF*ORIENT C MPOVMR.VPOCHA(NLCF,3)=XSURF/SURF*ORIENT MPOVMR.VPOCHA(NLCF,4)=YSURF/SURF*ORIENT MPOVMR.VPOCHA(NLCF,1)=-1*YSURF/SURF*ORIENT MPOVMR.VPOCHA(NLCF,2)=XSURF/SURF*ORIENT C ENDDO IF(NBSOUS .NE. 1) SEGDES IPT1 ENDDO C********************** C CAS DE LA DIMENSION 3 C********************** ELSE C DO ISOUS=1,NBSOUS,1 IPT1 = MELFP IF(NBSOUS .NE. 1)THEN IPT1 = MELFP.LISOUS(ISOUS) SEGACT IPT1 ENDIF C NP=IPT1.NUM(/1)-1 NEL=IPT1.NUM(/2) C DO IEL=1,NEL,1 NF=IPT1.NUM(NP+1,IEL) XF = XCOOR((NF-1)*(IDIM+1)+1) YF = XCOOR((NF-1)*(IDIM+1)+2) ZF = XCOOR((NF-1)*(IDIM+1)+3) * IPO=IPT1.NUM(NP,IEL) DXP = XCOOR((IPO-1)*(IDIM+1)+1) - XF DYP = XCOOR((IPO-1)*(IDIM+1)+2) - YF DZP = XCOOR((IPO-1)*(IDIM+1)+3) - ZF XSURF=0.0D0 YSURF=0.0D0 ZSURF=0.0D0 DO IFAC=1,NP,1 DXPM1 = DXP DYPM1 = DYP DZPM1 = DZP IPO=IPT1.NUM(IFAC,IEL) DXP = XCOOR((IPO-1)*(IDIM+1)+1) - XF DYP = XCOOR((IPO-1)*(IDIM+1)+2) - YF DZP = XCOOR((IPO-1)*(IDIM+1)+3) - ZF DSURFX = 0.5D0 * ((DYPM1 * DZP) - (DZPM1 * DYP)) DSURFY = 0.5D0 * ((DZPM1 * DXP) - (DXPM1 * DZP)) DSURFZ = 0.5D0 * ((DXPM1 * DYP) - (DYPM1 * DXP)) XSURF=XSURF+DSURFX YSURF=YSURF+DSURFY ZSURF=ZSURF+DSURFZ ENDDO * SURF=(XSURF*XSURF)+(YSURF*YSURF)+(ZSURF*ZSURF) surf=max(surf,xpetit) SURF=sqrt(SURF) NLCF=MLENTI.LECT(NF) MPOVSU.VPOCHA(NLCF,1)=SURF C C************* Orientation selon FACEL C IPO=MELFL.NUM(1,NLCF) DSURFX=XF-XCOOR((IPO-1)*(IDIM+1)+1) DSURFY=YF-XCOOR((IPO-1)*(IDIM+1)+2) DSURFZ=ZF-XCOOR((IPO-1)*(IDIM+1)+3) ORIENT=SIGN(1.0D0,((DSURFX*XSURF)+(DSURFY*YSURF)+(DSURFZ & *ZSURF))) C MPOVNO.VPOCHA(NLCF,1)=XSURF/SURF*ORIENT MPOVNO.VPOCHA(NLCF,2)=YSURF/SURF*ORIENT MPOVNO.VPOCHA(NLCF,3)=ZSURF/SURF*ORIENT C C************* MATROT C C Normal C MPOVMR.VPOCHA(NLCF,7)=XSURF/SURF*ORIENT MPOVMR.VPOCHA(NLCF,8)=YSURF/SURF*ORIENT MPOVMR.VPOCHA(NLCF,9)=ZSURF/SURF*ORIENT C First direction (RX,RY,RZ) is normal to the C normal and FP (P = first point of FACEP) C IPO=IPT1.NUM(1,IEL) DXP = XCOOR((IPO-1)*(IDIM+1)+1) - XF DYP = XCOOR((IPO-1)*(IDIM+1)+2) - YF DZP = XCOOR((IPO-1)*(IDIM+1)+3) - ZF DSURFX = (ZSURF * DYP) - (YSURF * DZP) DSURFY = (XSURF * DZP) - (ZSURF * DXP) DSURFZ = (YSURF * DXP) - (XSURF * DYP) C C DZP=modulus of the RX,RY,RZ C DZP=(((DSURFX*DSURFX)+(DSURFY*DSURFY)+ & (DSURFZ*DSURFZ))**0.5D0) C MPOVMR.VPOCHA(NLCF,4)=DSURFX/DZP MPOVMR.VPOCHA(NLCF,5)=DSURFY/DZP MPOVMR.VPOCHA(NLCF,6)=DSURFZ/DZP C C (UX,UY,UZ,RX,RY,RZ,MX,MY,MZ)=(1,2,3,4,5,6,7,8,9) C M,R,U is a right-hand normal frame C RMX=MPOVMR.VPOCHA(NLCF,7) RMY=MPOVMR.VPOCHA(NLCF,8) RMZ=MPOVMR.VPOCHA(NLCF,9) RRX=MPOVMR.VPOCHA(NLCF,4) RRY=MPOVMR.VPOCHA(NLCF,5) RRZ=MPOVMR.VPOCHA(NLCF,6) C MPOVMR.VPOCHA(NLCF,1)=(RMY*RRZ) - (RMZ*RRY) MPOVMR.VPOCHA(NLCF,2)=(RMZ*RRX) - (RMX*RRZ) MPOVMR.VPOCHA(NLCF,3)=(RMX*RRY) - (RMY*RRX) ENDDO IF(NBSOUS .NE. 1) SEGDES IPT1 ENDDO ENDIF C SEGDES MPOVSU SEGDES MPOVNO SEGDES MPOVMR SEGDES MELFP C SEGDES MELFL SEGSUP MLENTI C 9999 RETURN END