📄 wmlintpakimauniform3.cpp
字号:
// Magic Software, Inc.
// http://www.magic-software.com
// http://www.wild-magic.com
// Copyright (c) 2003. All Rights Reserved
//
// The Wild Magic Library (WML) source code is supplied under the terms of
// the license agreement http://www.magic-software.com/License/WildMagic.pdf
// and may not be copied or disclosed except in accordance with the terms of
// that agreement.
#include "WmlIntpAkimaUniform3.h"
#include "WmlMath.h"
using namespace Wml;
//----------------------------------------------------------------------------
template <class Real>
IntpAkimaUniform3<Real>::IntpAkimaUniform3 (int iXBound, int iYBound,
int iZBound, Real fXMin, Real fXSpacing, Real fYMin, Real fYSpacing,
Real fZMin, Real fZSpacing, Real*** aaafF)
{
// At least a 3x3x3 block of data points are needed to construct the
// estimates of the boundary derivatives.
assert( iXBound >= 3 && iYBound >= 3 && iZBound >= 3 && aaafF );
assert( fXSpacing > (Real)0.0 && fYSpacing > (Real)0.0
&& fZSpacing > (Real)0.0 );
m_iXBound = iXBound;
m_iYBound = iYBound;
m_iZBound = iZBound;
m_fXMin = fXMin;
m_fXSpacing = fXSpacing;
m_fYMin = fYMin;
m_fYSpacing = fYSpacing;
m_fZMin = fZMin;
m_fZSpacing = fZSpacing;
m_aaafF = aaafF;
int iXBoundM1 = iXBound - 1;
int iYBoundM1 = iYBound - 1;
int iZBoundM1 = iZBound - 1;
m_iQuantity = iXBound*iYBound*iZBound;
m_fXMax = fXMin + fXSpacing*iXBoundM1;
m_fYMax = fYMin + fYSpacing*iYBoundM1;
m_fZMax = fZMin + fZSpacing*iZBoundM1;
// compute slopes
Real fInvDX = ((Real)1.0)/fXSpacing;
Real fInvDY = ((Real)1.0)/fYSpacing;
Real fInvDZ = ((Real)1.0)/fZSpacing;
Real fInvDXDY = fInvDX*fInvDY;
Real fInvDXDZ = fInvDX*fInvDZ;
Real fInvDYDZ = fInvDY*fInvDZ;
Real fInvDXDYDZ = fInvDX*fInvDYDZ;
// xslope[z][y][x]
Real*** aaafXSlope;
Allocate3D(iXBound+3,iYBound,iZBound,aaafXSlope);
// yslope[z][x][y]
Real*** aaafYSlope;
Allocate3D(iYBound+3,iXBound,iZBound,aaafYSlope);
// zslope[y][x][z]
Real*** aaafZSlope;
Allocate3D(iZBound+3,iXBound,iYBound,aaafZSlope);
int iX, iY, iZ;
for (iZ = 0; iZ < iZBound; iZ++)
{
for (iY = 0; iY < iYBound; iY++)
{
for (iX = 0; iX < iXBoundM1; iX++)
{
aaafXSlope[iZ][iY][iX+2] = (aaafF[iZ][iY][iX+1] -
aaafF[iZ][iY][iX])*fInvDX;
}
aaafXSlope[iZ][iY][1] = ((Real)2.0)*aaafXSlope[iZ][iY][2] -
aaafXSlope[iZ][iY][3];
aaafXSlope[iZ][iY][0] = ((Real)2.0)*aaafXSlope[iZ][iY][1] -
aaafXSlope[iZ][iY][2];
aaafXSlope[iZ][iY][iXBound+1] =
((Real)2.0)*aaafXSlope[iZ][iY][iXBound] -
aaafXSlope[iZ][iY][iXBound-1];
aaafXSlope[iZ][iY][iXBound+2] =
((Real)2.0)*aaafXSlope[iZ][iY][iXBound+1] -
aaafXSlope[iZ][iY][iXBound];
}
}
for (iZ = 0; iZ < iZBound; iZ++)
{
for (iX = 0; iX < iXBound; iX++)
{
for (iY = 0; iY < iYBoundM1; iY++)
{
aaafYSlope[iZ][iX][iY+2] = (aaafF[iZ][iY+1][iX] -
aaafF[iZ][iY][iX])*fInvDY;
}
aaafYSlope[iZ][iX][1] = ((Real)2.0)*aaafYSlope[iZ][iX][2] -
aaafYSlope[iZ][iX][3];
aaafYSlope[iZ][iX][0] = ((Real)2.0)*aaafYSlope[iZ][iX][1] -
aaafYSlope[iZ][iX][2];
aaafYSlope[iZ][iX][iYBound+1] =
((Real)2.0)*aaafYSlope[iZ][iX][iYBound] -
aaafYSlope[iZ][iX][iYBound-1];
aaafYSlope[iZ][iX][iYBound+2] =
((Real)2.0)*aaafYSlope[iZ][iX][iYBound+1] -
aaafYSlope[iZ][iX][iYBound];
}
}
for (iY = 0; iY < iYBound; iY++)
{
for (iX = 0; iX < iXBound; iX++)
{
for (iZ = 0; iZ < iZBoundM1; iZ++)
{
aaafZSlope[iY][iX][iZ+2] = (aaafF[iZ+1][iY][iX] -
aaafF[iZ][iY][iX])*fInvDZ;
}
aaafZSlope[iY][iX][1] = ((Real)2.0)*aaafZSlope[iY][iX][2] -
aaafZSlope[iY][iX][3];
aaafZSlope[iY][iX][0] = ((Real)2.0)*aaafZSlope[iY][iX][1] -
aaafZSlope[iY][iX][2];
aaafZSlope[iY][iX][iZBound+1] =
((Real)2.0)*aaafZSlope[iY][iX][iZBound] -
aaafZSlope[iY][iX][iZBound-1];
aaafZSlope[iY][iX][iZBound+2] =
((Real)2.0)*aaafZSlope[iY][iX][iZBound+1] -
aaafZSlope[iY][iX][iZBound];
}
}
// construct first-order derivatives
Real*** aaafFX;
Allocate3D(iXBound,iYBound,iZBound,aaafFX);
Real*** aaafFY;
Allocate3D(iXBound,iYBound,iZBound,aaafFY);
Real*** aaafFZ;
Allocate3D(iXBound,iYBound,iZBound,aaafFZ);
for (iZ = 0; iZ < iZBound; iZ++)
{
for (iY = 0; iY < iYBound; iY++)
{
for (iX = 0; iX < iXBound; iX++)
{
aaafFX[iZ][iY][iX] = ComputeDerivative(
aaafXSlope[iZ][iY]+iX);
}
}
}
for (iZ = 0; iZ < iZBound; iZ++)
{
for (iX = 0; iX < iXBound; iX++)
{
for (iY = 0; iY < iYBound; iY++)
{
aaafFY[iZ][iY][iX] = ComputeDerivative(
aaafYSlope[iZ][iX]+iY);
}
}
}
for (iY = 0; iY < iYBound; iY++)
{
for (iX = 0; iX < iXBound; iX++)
{
for (iZ = 0; iZ < iZBound; iZ++)
{
aaafFZ[iZ][iY][iX] = ComputeDerivative(
aaafZSlope[iY][iX]+iZ);
}
}
}
// construct second-order derivatives
Real*** aaafFXY;
Allocate3D(iXBound,iYBound,iZBound,aaafFXY);
Real*** aaafFXZ;
Allocate3D(iXBound,iYBound,iZBound,aaafFXZ);
Real*** aaafFYZ;
Allocate3D(iXBound,iYBound,iZBound,aaafFYZ);
int iX0 = iXBoundM1, iX1 = iX0-1, iX2 = iX1-1;
int iY0 = iYBoundM1, iY1 = iY0-1, iY2 = iY1-1;
int iZ0 = iZBoundM1, iZ1 = iZ0-1, iZ2 = iZ1-1;
for (iZ = 0; iZ < iZBound; iZ++)
{
// corners of z-slice
aaafFXY[iZ][0][0] = ((Real)0.25)*fInvDXDY*(
((Real)9.0)*aaafF[iZ][0][0]
-((Real)12.0)*aaafF[iZ][0][1]
+ ((Real)3.0)*aaafF[iZ][0][2]
-((Real)12.0)*aaafF[iZ][1][0]
+((Real)16.0)*aaafF[iZ][1][1]
- ((Real)4.0)*aaafF[iZ][1][2]
+ ((Real)3.0)*aaafF[iZ][2][0]
- ((Real)4.0)*aaafF[iZ][2][1]
+ aaafF[iZ][2][2]);
aaafFXY[iZ][0][iXBoundM1] = ((Real)0.25)*fInvDXDY*(
((Real)9.0)*aaafF[iZ][0][iX0]
-((Real)12.0)*aaafF[iZ][0][iX1]
+ ((Real)3.0)*aaafF[iZ][0][iX2]
-((Real)12.0)*aaafF[iZ][1][iX0]
+((Real)16.0)*aaafF[iZ][1][iX1]
- ((Real)4.0)*aaafF[iZ][1][iX2]
+ ((Real)3.0)*aaafF[iZ][2][iX0]
- ((Real)4.0)*aaafF[iZ][2][iX1]
+ aaafF[iZ][2][iX2]);
aaafFXY[iZ][iYBoundM1][0] = ((Real)0.25)*fInvDXDY*(
((Real)9.0)*aaafF[iZ][iY0][0]
-((Real)12.0)*aaafF[iZ][iY0][1]
+ ((Real)3.0)*aaafF[iZ][iY0][2]
-((Real)12.0)*aaafF[iZ][iY1][0]
+((Real)16.0)*aaafF[iZ][iY1][1]
- ((Real)4.0)*aaafF[iZ][iY1][2]
+ ((Real)3.0)*aaafF[iZ][iY2][0]
- ((Real)4.0)*aaafF[iZ][iY2][1]
+ aaafF[iZ][iY2][2]);
aaafFXY[iZ][iYBoundM1][iXBoundM1] = ((Real)0.25)*fInvDXDY*(
((Real)9.0)*aaafF[iZ][iY0][iX0]
-((Real)12.0)*aaafF[iZ][iY0][iX1]
+ ((Real)3.0)*aaafF[iZ][iY0][iX2]
-((Real)12.0)*aaafF[iZ][iY1][iX0]
+((Real)16.0)*aaafF[iZ][iY1][iX1]
- ((Real)4.0)*aaafF[iZ][iY1][iX2]
+ ((Real)3.0)*aaafF[iZ][iY2][iX0]
- ((Real)4.0)*aaafF[iZ][iY2][iX1]
+ aaafF[iZ][iY2][iX2]);
// x-edges of z-slice
for (iX = 1; iX < iXBoundM1; iX++)
{
aaafFXY[iZ][0][iX] = ((Real)0.25)*fInvDXDY*(
((Real)3.0)*(aaafF[iZ][0][iX-1] - aaafF[iZ][0][iX+1]) -
((Real)4.0)*(aaafF[iZ][1][iX-1] - aaafF[iZ][1][iX+1]) +
(aaafF[iZ][2][iX-1] - aaafF[iZ][2][iX+1]));
aaafFXY[iZ][iYBoundM1][iX] = ((Real)0.25)*fInvDXDY*(
((Real)3.0)*(aaafF[iZ][iY0][iX-1] - aaafF[iZ][iY0][iX+1]) -
((Real)4.0)*(aaafF[iZ][iY1][iX-1] - aaafF[iZ][iY1][iX+1]) +
(aaafF[iZ][iY2][iX-1] - aaafF[iZ][iY2][iX+1]));
}
// y-edges of z-slice
for (iY = 1; iY < iYBoundM1; iY++)
{
aaafFXY[iZ][iY][0] = ((Real)0.25)*fInvDXDY*(
((Real)3.0)*(aaafF[iZ][iY-1][0] - aaafF[iZ][iY+1][0]) -
((Real)4.0)*(aaafF[iZ][iY-1][1] - aaafF[iZ][iY+1][1]) +
(aaafF[iZ][iY-1][2] - aaafF[iZ][iY+1][2]));
aaafFXY[iZ][iY][iXBoundM1] = ((Real)0.25)*fInvDXDY*(
((Real)3.0)*(aaafF[iZ][iY-1][iX0] - aaafF[iZ][iY+1][iX0]) -
((Real)4.0)*(aaafF[iZ][iY-1][iX1] - aaafF[iZ][iY+1][iX1]) +
(aaafF[iZ][iY-1][iX2] - aaafF[iZ][iY+1][iX2]));
}
// interior of z-slice
for (iY = 1; iY < iYBoundM1; iY++)
{
for (iX = 1; iX < iXBoundM1; iX++)
{
aaafFXY[iZ][iY][iX] = ((Real)0.25)*fInvDXDY*(
aaafF[iZ][iY-1][iX-1] - aaafF[iZ][iY-1][iX+1] -
aaafF[iZ][iY+1][iX-1] + aaafF[iZ][iY+1][iX+1]);
}
}
}
for (iY = 0; iY < iYBound; iY++)
{
// corners of z-slice
aaafFXZ[0][iY][0] = ((Real)0.25)*fInvDXDZ*(
((Real)9.0)*aaafF[0][iY][0]
-((Real)12.0)*aaafF[0][iY][1]
+ ((Real)3.0)*aaafF[0][iY][2]
-((Real)12.0)*aaafF[1][iY][0]
+((Real)16.0)*aaafF[1][iY][1]
- ((Real)4.0)*aaafF[1][iY][2]
+ ((Real)3.0)*aaafF[2][iY][0]
- ((Real)4.0)*aaafF[2][iY][1]
+ aaafF[2][iY][2]);
aaafFXZ[0][iY][iXBoundM1] = ((Real)0.25)*fInvDXDZ*(
((Real)9.0)*aaafF[0][iY][iX0]
-((Real)12.0)*aaafF[0][iY][iX1]
+ ((Real)3.0)*aaafF[0][iY][iX2]
-((Real)12.0)*aaafF[1][iY][iX0]
+((Real)16.0)*aaafF[1][iY][iX1]
- ((Real)4.0)*aaafF[1][iY][iX2]
+ ((Real)3.0)*aaafF[2][iY][iX0]
- ((Real)4.0)*aaafF[2][iY][iX1]
+ aaafF[2][iY][iX2]);
aaafFXZ[iZBoundM1][iY][0] = ((Real)0.25)*fInvDXDZ*(
((Real)9.0)*aaafF[iZ0][iY][0]
-((Real)12.0)*aaafF[iZ0][iY][1]
+ ((Real)3.0)*aaafF[iZ0][iY][2]
-((Real)12.0)*aaafF[iZ1][iY][0]
+((Real)16.0)*aaafF[iZ1][iY][1]
- ((Real)4.0)*aaafF[iZ1][iY][2]
+ ((Real)3.0)*aaafF[iZ2][iY][0]
- ((Real)4.0)*aaafF[iZ2][iY][1]
+ aaafF[iZ2][iY][2]);
aaafFXZ[iZBoundM1][iY][iXBoundM1] = ((Real)0.25)*fInvDXDZ*(
((Real)9.0)*aaafF[iZ0][iY][iX0]
-((Real)12.0)*aaafF[iZ0][iY][iX1]
+ ((Real)3.0)*aaafF[iZ0][iY][iX2]
-((Real)12.0)*aaafF[iZ1][iY][iX0]
+((Real)16.0)*aaafF[iZ1][iY][iX1]
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -