check.for

边界元前处理程序,检查边界网格及其法向方向

      PROGRAM MAIN

	PARAMETER (N=1120,NE=1120)
      REAL*4   COOR(N,3),q1,q2,q3
	INTEGER  NORD(NE,8),KODT(NE),k1,k2,k3,k4,k5


	OPEN(1,FILE='lilinchock0.DAT',STATUS='OLD')
	OPEN(3,FILE='drawroll.SCR',STATUS='old')
	
 


	DO I=1,N
      READ(1,*)J,COOR(J,:)
      ENDDO
   	
	     DO I=1,NE
      READ(1,*) J,NORD(I,1:4)    !,KODT(J)
      ENDDO
 


   

 	write(*,*)"======="
	
	CLOSE(1)
	CALL DRAW(N,NE,COOR,NORD,KODT)

	END
       

C----------- 输出CAD脚本文件，用来检查网格划分是否正确 ---------------
      
      SUBROUTINE DRAW(N,NE,COOR,NORD,KODT)   
	COMMON/S6/ FJACOB,COSBX,COSBY,COSBZ
	REAL*4  COOR(N,3),XK(8),YK(8),ZK(8),FN(8)
	INTEGER NORD(NE,8),NODAL(8),KODT(NE) 

	DO KK=1,NE
  	WRITE(3,1002) ( COOR(NORD(KK,I),1),
	1                COOR(NORD(KK,I),2),
	1                COOR(NORD(KK,I),3), I=1,4 )            
      ENDDO

1002  FORMAT(7H3DFACE ,SP,E12.6,1H,,SP,E12.6,1H,,SP,E12.6,1X,SP,
	1 E12.6,',',SP,E12.6,',',SP,E12.6,1X,SP,E12.6,',',SP,E12.6,',',  
	1 SP,E12.6,1X,SP,E12.6,',',SP,E12.6,',',SP,E12.6,' ')

C--------------- 生成单元外法矢延长线 ---------------------------
	
	PI=ATAN(1.0)*4.

	DO 260 I=1,ne
      S=20.
	IF(KODT(I).EQ.3.OR.KODT(I).EQ.2) S=50.
	DO 261 J=1,8
	NODAL(J)=NORD(I,J)
	IF(NODAL(J).EQ.0) GOTO 262
	FN(J)=0.5
	XK(J)=COOR(NODAL(J),1)
	YK(J)=COOR(NODAL(J),2)
	ZK(J)=COOR(NODAL(J),3)
	GOTO 261
262   FN(J)=0.
	XK(J)=0.
	YK(J)=0.
	ZK(J)=0.
261   CONTINUE
	FN(1)=0.25-0.5*(FN(5)+FN(8))
	FN(2)=0.25-0.5*(FN(5)+FN(6))
	FN(3)=0.25-0.5*(FN(6)+FN(7))
	FN(4)=0.25-0.5*(FN(7)+FN(8))
 	
      XM=FN(1)*XK(1)+FN(2)*XK(2)+FN(3)*XK(3)+FN(4)*XK(4)+
 	1  FN(5)*XK(5)+FN(6)*XK(6)+FN(7)*XK(7)+FN(8)*XK(8)
      YM=FN(1)*YK(1)+FN(2)*YK(2)+FN(3)*YK(3)+FN(4)*YK(4)+
  	1  FN(5)*YK(5)+FN(6)*YK(6)+FN(7)*YK(7)+FN(8)*YK(8)
      ZM=FN(1)*ZK(1)+FN(2)*ZK(2)+FN(3)*ZK(3)+FN(4)*ZK(4)+
 	1  FN(5)*ZK(5)+FN(6)*ZK(6)+FN(7)*ZK(7)+FN(8)*ZK(8)
	CALL JACOB(NODAL,0.0,0.0,XK,YK,ZK)

	ANGLX=XM+S*COSBX
	ANGLY=YM+S*COSBY
	ANGLZ=ZM+S*COSBZ

      WRITE(3,1003) XM,YM,ZM,ANGLX,ANGLY,ANGLZ
260   CONTINUE
1003  FORMAT(5HLINE ,SP,E12.6,1H,,SP,E12.6,1H,,SP,E12.6,1X,
     1               SP,E12.6,',',SP,E12.6,',',SP,E12.6,' ')    
	CLOSE(3)

	RETURN 
	END

C********************** 计算雅可比 **************************
	
	SUBROUTINE JACOB(NODAL,XW,YW,XK,YK,ZK)
	COMMON/S6/ FJACOB,COSBX,COSBY,COSBZ
	INTEGER NODAL(8)
      REAL*4  XK(8),YK(8),ZK(8),DFXW(8),DFYW(8)
      REAL    XD1,YD1,ZD1,XD2,YD2,ZD2

 	DO 263 J=5,8
 	DFXW(J)=0.
263   DFYW(J)=0.
	IF(NODAL(5).NE.0) DFXW(5)=-XW*(1.-YW)
	IF(NODAL(6).NE.0) DFYW(6)=-YW*(1.+XW)
	IF(NODAL(7).NE.0) DFXW(7)=-XW*(1.+YW)
	IF(NODAL(8).NE.0) DFYW(8)=-YW*(1.-XW)
	IF(NODAL(5).NE.0) DFYW(5)=-0.5*(1.-XW*XW)
	IF(NODAL(6).NE.0) DFXW(6)= 0.5*(1.-YW*YW)
	IF(NODAL(7).NE.0) DFYW(7)= 0.5*(1.-XW*XW)
	IF(NODAL(8).NE.0) DFXW(8)=-0.5*(1.-YW*YW)
	DFXW(1)=-0.25*(1.-YW)-0.5*(DFXW(5)+DFXW(8))
	DFXW(2)= 0.25*(1.-YW)-0.5*(DFXW(5)+DFXW(6))
	DFXW(3)= 0.25*(1.+YW)-0.5*(DFXW(6)+DFXW(7))
	DFXW(4)=-0.25*(1.+YW)-0.5*(DFXW(7)+DFXW(8))
	DFYW(1)=-0.25*(1.-XW)-0.5*(DFYW(5)+DFYW(8))
	DFYW(2)=-0.25*(1.+XW)-0.5*(DFYW(5)+DFYW(6))
	DFYW(3)= 0.25*(1.+XW)-0.5*(DFYW(6)+DFYW(7))
	DFYW(4)= 0.25*(1.-XW)-0.5*(DFYW(7)+DFYW(8))
	XD1=DFXW(1)*XK(1)+DFXW(2)*XK(2)+DFXW(3)*XK(3)+DFXW(4)*XK(4)
     1   +DFXW(5)*XK(5)+DFXW(6)*XK(6)+DFXW(7)*XK(7)+DFXW(8)*XK(8)
	YD1=DFXW(1)*YK(1)+DFXW(2)*YK(2)+DFXW(3)*YK(3)+DFXW(4)*YK(4)
     1   +DFXW(5)*YK(5)+DFXW(6)*YK(6)+DFXW(7)*YK(7)+DFXW(8)*YK(8)
	ZD1=DFXW(1)*ZK(1)+DFXW(2)*ZK(2)+DFXW(3)*ZK(3)+DFXW(4)*ZK(4)
     1   +DFXW(5)*ZK(5)+DFXW(6)*ZK(6)+DFXW(7)*ZK(7)+DFXW(8)*ZK(8)
	XD2=DFYW(1)*XK(1)+DFYW(2)*XK(2)+DFYW(3)*XK(3)+DFYW(4)*XK(4)
     1   +DFYW(5)*XK(5)+DFYW(6)*XK(6)+DFYW(7)*XK(7)+DFYW(8)*XK(8)
	YD2=DFYW(1)*YK(1)+DFYW(2)*YK(2)+DFYW(3)*YK(3)+DFYW(4)*YK(4)
     1   +DFYW(5)*YK(5)+DFYW(6)*YK(6)+DFYW(7)*YK(7)+DFYW(8)*YK(8)
	ZD2=DFYW(1)*ZK(1)+DFYW(2)*ZK(2)+DFYW(3)*ZK(3)+DFYW(4)*ZK(4)
     1   +DFYW(5)*ZK(5)+DFYW(6)*ZK(6)+DFYW(7)*ZK(7)+DFYW(8)*ZK(8)

	G1=YD1*ZD2-ZD1*YD2
	G2=ZD1*XD2-XD1*ZD2
	G3=XD1*YD2-YD1*XD2

	FJACOB=SQRT(G1*G1+G2*G2+G3*G3)
	COSBX=G1/FJACOB
	COSBY=G2/FJACOB
	COSBZ=G3/FJACOB

	RETURN
	END