BS3IG
This function evaluates the integral of a tensor-product spline in three dimensions over a three-dimensional rectangle, given its tensor-product B-spline representation.
Function Return Value
BS3IG — Integral of the spline over the three-dimensional rectangle (A, B) by (C, D) by (E, F). (Output)
Required Arguments
A — Lower limit of the X-variable. (Input)
B — Upper limit of the X-variable. (Input)
C — Lower limit of the Y-variable. (Input)
D — Upper limit of the Y-variable. (Input)
E — Lower limit of the Z-variable. (Input)
F — Upper limit of the Z-variable. (Input)
KXORD — Order of the spline in the X-direction. (Input)
KYORD — Order of the spline in the Y-direction. (Input)
KZORD — Order of the spline in the Z-direction. (Input)
XKNOT — Array of length NXCOEF + KXORD containing the knot sequence in the X-direction. (Input)
XKNOT must be nondecreasing.
YKNOT — Array of length NYCOEF + KYORD containing the knot sequence in the Y-direction. (Input)
YKNOT must be nondecreasing.
ZKNOT — Array of length NZCOEF + KZORD containing the knot sequence in the Z-direction. (Input)
ZKNOT must be nondecreasing.
NXCOEF — Number of B-spline coefficients in the X-direction. (Input)
NYCOEF — Number of B-spline coefficients in the Y-direction. (Input)
NZCOEF — Number of B-spline coefficients in the Z-direction. (Input)
BSCOEF — Array of length NXCOEF * NYCOEF * NZCOEF containing the tensor-product
B-spline coefficients. (Input)
BSCOEF is treated internally as a matrix of size NXCOEF by NYCOEF by NZCOEF.
FORTRAN 90 Interface
Generic: BS3IG (A, B, C, D, E, F, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, NXCOEF, NYCOEF, NZCOEF, BSCOEF)
Specific: The specific interface names are S_BS3IG and D_BS3IG.
FORTRAN 77 Interface
Single: BS3IG (A, B, C, D, E, F, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, NXCOEF, NYCOEF, NZCOEF, BSCOEF)
Double: The double precision function name is DBS3IG.
Description
The routine BS3IG computes the integral of a tensor-product three-dimensional spline, given its
B-spline representation. Specifically, given the knot sequence tx = XKNOT, ty = YKNOT, tz = ZKNOT, the order kx = KXORD, ky = KYORD, kz = KZORD, the coefficients β = BSCOEF, the number of coefficients nx = NXCOEF, ny = NYCOEF, nz = NZCOEF, and a three-dimensional rectangle [a, b] by [c, d] by [e, f], BS3IG returns the value
where
This routine uses the identity (22) on page 151 of de Boor (1978). It assumes (for all knot sequences) that the first and last k knots are stacked, that is, t1 = … = tk and tn + 1 = … = tn + k, where k is the order of the spline in the x, y, or z direction.
Comments
1. Workspace may be explicitly provided, if desired, by use of B23IG/DB23IG. The reference is:
CALL B23IG(A, B, C, D, E, F, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, NXCOEF, NYCOEF, NZCOEF, BSCOEF, WK)
The additional argument is:
WK — Work array of length 4 * (MAX(KXORD, KYORD, KZORD) + 1) + NYCOEF + NZCOEF.
2. Informational errors
|
Type |
Code |
Description |
|
3 |
1 |
The lower limit of the X-integration is less than XKNOT(KXORD). |
|
3 |
2 |
The upper limit of the X-integration is greater than XKNOT(NXCOEF + 1). |
|
3 |
3 |
The lower limit of the Y-integration is less than YKNOT(KYORD). |
|
3 |
4 |
The upper limit of the Y-integration is greater than YKNOT(NYCOEF + 1). |
|
3 |
5 |
The lower limit of the Z- integration is less than ZKNOT(KZORD). |
|
3 |
6 |
The upper limit of the Z-integration is greater than ZKNOT(NZCOEF + 1). |
|
4 |
13 |
Multiplicity of the knots cannot exceed the order of the spline. |
|
4 |
14 |
The knots must be nondecreasing. |
Example
We integrate the three-dimensional tensor-product quartic (kx = 5) by linear (ky = 2) by quadratic (kz = 3) spline which interpolates x3 + xyz at the points
{(i/10, j/5, m/7) : = -10, …, 10, j = 0, …, 5, and m = 0,…, 7}
over the rectangle [0, 1] × [.5, 1] × [0, .5]. The exact answer is 11/128.
USE BS3IG_INT
USE BS3IN_INT
USE BSNAK_INT
USE UMACH_INT
IMPLICIT NONE
! SPECIFICATIONS FOR PARAMETERS
INTEGER KXORD, KYORD, KZORD, LDF, MDF, NXDATA, NXKNOT,&
NYDATA, NYKNOT, NZDATA, NZKNOT
PARAMETER (KXORD=5, KYORD=2, KZORD=3, NXDATA=21, NYDATA=6,&
NZDATA=8, LDF=NXDATA, MDF=NYDATA,&
NXKNOT=NXDATA+KXORD, NYKNOT=NYDATA+KYORD,&
NZKNOT=NZDATA+KZORD)
!
INTEGER I, J, K, NOUT, NXCOEF, NYCOEF, NZCOEF
REAL A, B, BSCOEF(NXDATA,NYDATA,NZDATA), C , D, E,&
F, FDATA(LDF,MDF,NZDATA), FF, FIG, FLOAT, G, H, RI,&
RJ, VAL, X, XDATA(NXDATA), XKNOT(NXKNOT), Y,&
YDATA(NYDATA), YKNOT(NYKNOT), Z, ZDATA(NZDATA),&
ZKNOT(NZKNOT)
INTRINSIC FLOAT
! Define function
F(X,Y,Z) = X*X*X + X*Y*Z
! Set up interpolation points
DO 10 I=1, NXDATA
XDATA(I) = FLOAT(I-11)/10.0
10 CONTINUE
! Generate knot sequence
CALL BSNAK (NXDATA, XDATA, KXORD, XKNOT)
! Set up interpolation points
DO 20 I=1, NYDATA
YDATA(I) = FLOAT(I-1)/FLOAT(NYDATA-1)
20 CONTINUE
! Generate knot sequence
CALL BSNAK (NYDATA, YDATA, KYORD, YKNOT)
! Set up interpolation points
DO 30 I=1, NZDATA
ZDATA(I) = FLOAT(I-1)/FLOAT(NZDATA-1)
30 CONTINUE
! Generate knot sequence
CALL BSNAK (NZDATA, ZDATA, KZORD, ZKNOT)
! Generate FDATA
DO 50 K=1, NZDATA
DO 40 I=1, NYDATA
DO 40 J=1, NXDATA
FDATA(J,I,K) = F(XDATA(J),YDATA(I),ZDATA(K))
40 CONTINUE
50 CONTINUE
! Get output unit number
CALL UMACH (2, NOUT)
! Interpolate
CALL BS3IN (XDATA, YDATA, ZDATA, FDATA, KXORD, KYORD, KZORD, XKNOT, &
YKNOT, ZKNOT, BSCOEF)
!
NXCOEF = NXDATA
NYCOEF = NYDATA
NZCOEF = NZDATA
A = 0.0
B = 1.0
C = 0.5
D = 1.0
E = 0.0
FF = 0.5
! Integrate
VAL = BS3IG(A,B,C ,D,E,FF,KXORD,KYORD,KZORD,XKNOT,YKNOT,ZKNOT,&
NXCOEF,NYCOEF,NZCOEF,BSCOEF)
! Calculate integral directly
G = .5*(B**4-A**4)
H = (B-A)*(B+A)
RI = G*(D-C )
RJ = .5*H*(D-C )*(D+C )
FIG = .5*(RI*(FF-E)+.5*RJ*(FF-E)*(FF+E))
! Print results
WRITE (NOUT,99999) VAL, FIG, FIG - VAL
99999 FORMAT (' Computed Integral = ', F10.5, /, ' Exact Integral '&
, '= ', F10.5,/, ' Error '&
, '= ', F10.6, /)
END
Computed Integral = 0.08594
Exact Integral = 0.08594
Error = 0.000000
