我试图找出立方体的每个面的表面积和相应的向外单位法线。此操作在有限元网格上完成，因此我使用形状（基础）函数将立方体的每个表面转换为等参数形式，然后尝试提取区域和法线。
这是代码：
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
每个表面的面积应该是1个单位，对于这种情况，但是这个代码不是这样返回的。法线也是不正确的。输出是
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。这可能永远不是一个立方体，它可以是任何不规则的六面体，因此这些区域总是不相等，因此我们需要计算它们中的每一个。湾面可能沿不同方向定向，因此等参数变换是必要的。
这是解决此问题的正确方法吗？如果有人能帮我解决这个问题会很高兴。我也尝试使用对角线的矢量乘积来计算面积和单位法线，但是当结构不规则时它们不起作用。这是一个不规则的六面体1的示例图片。一个普通的立方体大致立方体看起来像这样。 2
正常：取三个顶点并计算它们形成的向量的叉积。
区域：在XY，YZ和ZX上应用鞋带公式，然后采用三个结果的欧几里德范数。