立方体の各面の表面積と、対応する外側の単位法線を調べようとしています。この操作は有限要素メッシュで行われるため、立方体の各面を shape(basis) 関数を使用してアイソパラメトリック形式に変換してから、領域と法線を抽出しようとしました。
コードは次のとおりです。
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
注: これは常に立方体であるとは限らず、不規則な六面体である可能性があるため、面積は常に等しくなるとは限らないため、それぞれを計算する必要があります。b. 面はさまざまな方向に向いている可能性があるため、アイソパラメトリック変換が必要です。
これは、この問題にアプローチする正しい方法ですか? 誰かがこれを理解するのを手伝ってくれたらうれしいです。また、各面の対角線のベクトル積を使用して面積と単位法線を計算しようとしましたが、構造が不規則な場合は機能しません。これは、不規則な六面体の写真の例です1。通常のルーブは大体立方体でこんな感じ。2
