program polyhedron
IMPLICIT NONE
real(8) coord(3,8)
INTEGER face, INPT, I, ino,k, ii
REAL(8) XI(3), dNdxi(8,3), ZERO, ONE, MONE, EIGHT
REAL(8) VJACOB(3,3),XII(8,3),norm(3), TWO, THREE, FOUR, HALF, SIX
PARAMETER(ZERO=0.D0,ONE=1.D0,MONE=-1.D0,EIGHT=8.D0)
PARAMETER(TWO=2.D0,THREE=3.D0,FOUR=4.D0,HALF=0.5D0,SIX=6.D0)
REAL(8) dXdXi,dXdEta,dXdZeta,dYdXi,dYdEta,dYdZeta,dZdXi,dZdEta,area_isd
REAL(8) dZdZeta, dA, mag, normal(3,1),xLocal(4),yLocal(4),zLocal(4)
!COORDINATES OF THE CUBE
coord(1,1)=1.00
coord(2,1)=1.00
coord(3,1)=1.00
coord(1,2)=1.00
coord(2,2)=0.00
coord(3,2)=1.00
coord(1,3)=1.00
coord(2,3)=1.00
coord(3,3)=0.00
coord(1,4)=1.00
coord(2,4)=0.00
coord(3,4)=0.00
coord(1,5)=0.00
coord(2,5)=1.00
coord(3,5)=1.00
coord(1,6)=0.00
coord(2,6)=0.00
coord(3,6)=1.00
coord(1,7)=0.00
coord(2,7)=1.00
coord(3,7)=0.00
coord(1,8)=0.00
coord(2,8)=0.00
coord(3,8)=0.00
do face=1,6 !Loop over the faces
area_isd=0.0
call xintSurf3D4pt(face,xLocal,yLocal,zLocal) !get local points
do ii=1,4
call computeSurf3D(xLocal(ii),yLocal(ii),zLocal(ii),face,coord,dA,norm) !compute area and normal
area_isd=area_isd+dA
end do
write(*,*) 'face', face, 'area', area_isd
write(*,*) 'norm', norm
end do
end program polyhedron
计算局部雅可比矩阵和法线的子程序为:
subroutine computeSurf3D(xLocal,yLocal,zLocal,face,coords,dA,normal)
IMPLICIT NONE
integer face,stat,i,j,k
real(8) xLocal,yLocal,zLocal,dA,dshxi(8,3),zero,dsh(8,3),one
real(8) coords(3,8),two,eighth,mapJ(3,3),mag,normal(3)
real(8) dXdXi,dXdEta,dXdZeta,dYdXi,dYdEta,dYdZeta,dZdXi,dZdEta
real(8) dZdZeta
parameter(one=1.d0,two=2.d0,eighth=1.d0/8.d0,zero=0.d0)
!Hex shape function derivatives
dshxi(1,1) = -eighth*(one - yLocal)*(one - zLocal)
dshxi(1,2) = -eighth*(one - xLocal)*(one - zLocal)
dshxi(1,3) = -eighth*(one - xLocal)*(one - yLocal)
dshxi(2,1) = eighth*(one - yLocal)*(one - zLocal)
dshxi(2,2) = -eighth*(one + xLocal)*(one - zLocal)
dshxi(2,3) = -eighth*(one + xLocal)*(one - yLocal)
dshxi(3,1) = eighth*(one + yLocal)*(one - zLocal)
dshxi(3,2) = eighth*(one + xLocal)*(one - zLocal)
dshxi(3,3) = -eighth*(one + xLocal)*(one + yLocal)
dshxi(4,1) = -eighth*(one + yLocal)*(one - zLocal)
dshxi(4,2) = eighth*(one - xLocal)*(one - zLocal)
dshxi(4,3) = -eighth*(one - xLocal)*(one + yLocal)
dshxi(5,1) = -eighth*(one - yLocal)*(one + zLocal)
dshxi(5,2) = -eighth*(one - xLocal)*(one + zLocal)
dshxi(5,3) = eighth*(one - xLocal)*(one - yLocal)
dshxi(6,1) = eighth*(one - yLocal)*(one + zLocal)
dshxi(6,2) = -eighth*(one + xLocal)*(one + zLocal)
dshxi(6,3) = eighth*(one + xLocal)*(one - yLocal)
dshxi(7,1) = eighth*(one + yLocal)*(one + zLocal)
dshxi(7,2) = eighth*(one + xLocal)*(one + zLocal)
dshxi(7,3) = eighth*(one + xLocal)*(one + yLocal)
dshxi(8,1) = -eighth*(one + yLocal)*(one + zLocal)
dshxi(8,2) = eighth*(one - xLocal)*(one + zLocal)
dshxi(8,3) = eighth*(one - xLocal)*(one + yLocal)
dXdXi = zero
dXdEta = zero
dXdZeta = zero
dYdXi = zero
dYdEta = zero
dYdZeta = zero
dZdXi = zero
dZdEta = zero
dZdZeta = zero
do k=1,8
dXdXi = dXdXi + dshxi(k,1)*coords(1,k)
dXdEta = dXdEta + dshxi(k,2)*coords(1,k)
dXdZeta = dXdZeta + dshxi(k,3)*coords(1,k)
dYdXi = dYdXi + dshxi(k,1)*coords(2,k)
dYdEta = dYdEta + dshxi(k,2)*coords(2,k)
dYdZeta = dYdZeta + dshxi(k,3)*coords(2,k)
dZdXi = dZdXi + dshxi(k,1)*coords(3,k)
dZdEta = dZdEta + dshxi(k,2)*coords(3,k)
dZdZeta = dZdZeta + dshxi(k,3)*coords(3,k)
enddo
! Jacobian of the mapping
!
if((face.eq.1).or.(face.eq.2)) then
! zeta = constant on this face
dA = dsqrt((dYdXi*dZdEta - dYdEta*dZdXi)**2+(dXdXi*dZdEta - dXdEta*dZdXi)**2+(dXdXi*dYdEta - dXdEta*dYdXi)**2)
elseif((face.eq.3).or.(face.eq.4)) then
! eta = constant on this face
dA = dsqrt((dYdXi*dZdZeta - dYdZeta*dZdXi)**2+(dXdXi*dZdZeta - dXdZeta*dZdXi)**2+(dXdXi*dYdZeta - dXdZeta*dYdXi)**2)
elseif((face.eq.5).or.(face.eq.6)) then
! xi = constant on this face
dA = dsqrt((dYdEta*dZdZeta - dYdZeta*dZdEta)**2+(dXdEta*dZdZeta - dXdZeta*dZdEta)**2+(dXdEta*dYdZeta - dXdZeta*dYdEta)**2)
endif
!
if((face.eq.1).or.(face.eq.2)) then
! zeta = constant on this face
normal(1) = dYdXi*dZdEta - dYdEta*dZdXi
normal(2) = dXdXi*dZdEta - dXdEta*dZdXi
normal(3) = dXdXi*dYdEta - dXdEta*dYdXi
if(face.eq.1) normal = -normal
elseif((face.eq.3).or.(face.eq.4)) then
! eta = constant on this face
normal(1) = dYdXi*dZdZeta - dYdZeta*dZdXi
normal(2) = dXdXi*dZdZeta - dXdZeta*dZdXi
normal(3) = dXdXi*dYdZeta - dXdZeta*dYdXi
if(face.eq.3) normal = -normal
elseif((face.eq.5).or.(face.eq.6)) then
! xi = constant on this face
normal(1) = dYdEta*dZdZeta - dYdZeta*dZdEta
normal(2) = dXdEta*dZdZeta - dXdZeta*dZdEta
normal(3) = dXdEta*dYdZeta - dXdZeta*dYdEta
if(face.eq.5) normal = -normal
endif
mag = dsqrt(normal(1)**two+normal(2)**two+normal(3)**two)
normal(1) = normal(1)/mag
normal(2) = normal(2)/mag
normal(3) = normal(3)/mag
end subroutine computeSurf3D
局部高斯点可从该子程序获得:
subroutine xintSurf3D4pt(face,xLocal,yLocal,zLocal)
implicit none
integer face
real(8) xLocal(4),yLocal(4),zLocal(4),w(4),one,three
parameter(one=1.d0,three=3.d0)
! Gauss pt locations in master element
!
if(face.eq.1) then
xLocal(1) = -dsqrt(one/three)
yLocal(1) = -dsqrt(one/three)
zLocal(1) = -one
xLocal(2) = dsqrt(one/three)
yLocal(2) = -dsqrt(one/three)
zLocal(2) = -one
xLocal(3) = dsqrt(one/three)
yLocal(3) = dsqrt(one/three)
zLocal(3) = -one
xLocal(4) = -dsqrt(one/three)
yLocal(4) = dsqrt(one/three)
zLocal(4) = -one
elseif(face.eq.2) then
xLocal(1) = -dsqrt(one/three)
yLocal(1) = -dsqrt(one/three)
zLocal(1) = one
xLocal(2) = dsqrt(one/three)
yLocal(2) = -dsqrt(one/three)
zLocal(2) = one
xLocal(3) = dsqrt(one/three)
yLocal(3) = dsqrt(one/three)
zLocal(3) = one
xLocal(4) = -dsqrt(one/three)
yLocal(4) = dsqrt(one/three)
zLocal(4) = one
elseif(face.eq.3) then
xLocal(1) = -dsqrt(one/three)
yLocal(1) = -one
zLocal(1) = -dsqrt(one/three)
xLocal(2) = dsqrt(one/three)
yLocal(2) = -one
zLocal(2) = -dsqrt(one/three)
xLocal(3) = dsqrt(one/three)
yLocal(3) = -one
zLocal(3) = dsqrt(one/three)
xLocal(4) = -dsqrt(one/three)
yLocal(4) = -one
zLocal(4) = dsqrt(one/three)
elseif(face.eq.4) then
xLocal(1) = one
yLocal(1) = -dsqrt(one/three)
zLocal(1) = -dsqrt(one/three)
xLocal(2) = one
yLocal(2) = dsqrt(one/three)
zLocal(2) = -dsqrt(one/three)
xLocal(3) = one
yLocal(3) = dsqrt(one/three)
zLocal(3) = dsqrt(one/three)
xLocal(4) = one
yLocal(4) = -dsqrt(one/three)
zLocal(4) = dsqrt(one/three)
elseif(face.eq.5) then
xLocal(1) = -dsqrt(one/three)
yLocal(1) = one
zLocal(1) = -dsqrt(one/three)
xLocal(2) = dsqrt(one/three)
yLocal(2) = one
zLocal(2) = -dsqrt(one/three)
xLocal(3) = dsqrt(one/three)
yLocal(3) = one
zLocal(3) = dsqrt(one/three)
xLocal(4) = -dsqrt(one/three)
yLocal(4) = one
zLocal(4) = dsqrt(one/three)
elseif(face.eq.6) then
xLocal(1) = -one
yLocal(1) = -dsqrt(one/three)
zLocal(1) = -dsqrt(one/three)
xLocal(2) = -one
yLocal(2) = dsqrt(one/three)
zLocal(2) = -dsqrt(one/three)
xLocal(3) = -one
yLocal(3) = dsqrt(one/three)
zLocal(3) = dsqrt(one/three)
xLocal(4) = -one
yLocal(4) = -dsqrt(one/three)
zLocal(4) = dsqrt(one/three)
endif
end subroutine xintSurf3D4pt
face 1 area 0.57735026918962573
norm 1.0000000000000000 0.0000000000000000E+000 0.0000000000000000E+000
face 2 area 0.57735026918962573
norm -1.0000000000000000 -0.0000000000000000E+000 -0.0000000000000000E+000
face 3 area 1.0000000000000000
norm 0.0000000000000000E+000 -0.0000000000000000E+000 1.0000000000000000
face 4 area 0.57735026918962573
norm -0.0000000000000000E+000 0.0000000000000000E+000 -1.0000000000000000
face 5 area 1.1547005383792515
norm -0.0000000000000000E+000 0.86602540378443871 0.50000000000000000
face 6 area 1.4142135623730951
norm 0.0000000000000000E+000 -0.70710678118654746 -0.70710678118654746
注:a.这可能永远不是一个立方体,它可以是任何不规则的六面体,因此面积不会总是相等的,因此我们需要计算每一个。
b、 这些面可能沿不同的方向定向,因此等参变换是必要的。
1
. 一个规则的立方体看起来像这样。
2
