Skip to content

Commit f13e799

Browse files
perazzperazz
committed
Add scattered-point eval to fitpack_grid_surface
- Add eval(x,y) and eval(x(:),y(:)) for arbitrary scattered points using bispeu (no sorting requirement, O(1) workspace) - Rename previous gridded eval to eval_ongrid for consistency with fitpack_surface, which already had this split - Both surface types now have matching interfaces: eval (scattered, bispeu) and eval_ongrid (grid, bispev) - Add test_grid_surface_scattered_eval verifying scattered, single- point, and gridded evaluation on z = x^2 + y^2 - Update existing tests to use eval_ongrid where grid output is needed
1 parent 2631abf commit f13e799

File tree

3 files changed

+127
-7
lines changed

3 files changed

+127
-7
lines changed

‎src/fitpack_grid_surfaces.f90‎

Lines changed: 49 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@
1919
! **************************************************************************************************
2020
module fitpack_grid_surfaces
2121
use fitpack_core, only: FITPACK_SUCCESS,FP_REAL,FP_SIZE,FP_FLAG,FP_COMM,zero,IOPT_NEW_SMOOTHING,IOPT_OLD_FIT, &
22-
IOPT_NEW_LEASTSQUARES,bispev,fitpack_error_handling,get_smoothing,regrid, &
22+
IOPT_NEW_LEASTSQUARES,bispev,bispeu,fitpack_error_handling,get_smoothing,regrid, &
2323
parder,pardeu,FITPACK_INPUT_ERROR, &
2424
dblint,profil,pardtc, &
2525
FP_COMM_SIZE,FP_COMM_PACK,FP_COMM_EXPAND
@@ -70,10 +70,15 @@ module fitpack_grid_surfaces
7070
procedure :: least_squares => surface_fit_least_squares
7171
procedure :: interpolate => surface_fit_interpolating
7272

73-
!> Evaluate gridded domain at given x,y coordinates
73+
!> Evaluate at scattered (x,y) points
74+
procedure, private :: gridsurf_eval_one
75+
procedure, private :: gridsurf_eval_many
76+
generic :: eval => gridsurf_eval_one,gridsurf_eval_many
77+
78+
!> Evaluate on a grid domain
7479
procedure, private :: gridded_eval_one
7580
procedure, private :: gridded_eval_many
76-
generic :: eval => gridded_eval_one,gridded_eval_many
81+
generic :: eval_ongrid => gridded_eval_one,gridded_eval_many
7782

7883
!> Evaluate derivatives at given coordinates
7984
procedure, private :: gridded_derivatives_gridded
@@ -290,6 +295,47 @@ integer function surf_new_fit(this,x,y,z,smoothing,order)
290295

291296
end function surf_new_fit
292297

298+
!> Evaluate surface at a list of scattered (x(i),y(i)) points using bispeu
299+
function gridsurf_eval_many(this,x,y,ierr) result(f)
300+
class(fitpack_grid_surface), intent(inout) :: this
301+
real(FP_REAL), intent(in) :: x(:),y(size(x))
302+
real(FP_REAL) :: f(size(x))
303+
integer(FP_FLAG), optional, intent(out) :: ierr
304+
305+
integer(FP_FLAG) :: ier
306+
integer(FP_SIZE) :: min_lwrk
307+
real(FP_REAL), allocatable :: wrk(:)
308+
309+
! bispeu workspace: lwrk >= kx+ky+2
310+
min_lwrk = sum(this%order) + 2
311+
allocate(wrk(min_lwrk))
312+
313+
call bispeu(tx=this%t(:,1),nx=this%knots(1), &
314+
ty=this%t(:,2),ny=this%knots(2), &
315+
c=this%c, &
316+
kx=this%order(1),ky=this%order(2), &
317+
x=x,y=y,z=f,m=size(x), &
318+
wrk=wrk,lwrk=min_lwrk, &
319+
ier=ier)
320+
321+
call fitpack_error_handling(ier,ierr,'evaluate grid surface at scattered points')
322+
323+
end function gridsurf_eval_many
324+
325+
!> Evaluate surface at a single (x,y) point using bispeu
326+
real(FP_REAL) function gridsurf_eval_one(this,x,y,ierr) result(f)
327+
class(fitpack_grid_surface), intent(inout) :: this
328+
real(FP_REAL), intent(in) :: x,y
329+
integer(FP_FLAG), optional, intent(out) :: ierr
330+
331+
real(FP_REAL) :: z1(1)
332+
333+
z1 = gridsurf_eval_many(this,[x],[y],ierr)
334+
f = z1(1)
335+
336+
end function gridsurf_eval_one
337+
338+
!> Evaluate surface on a grid of x(:) x y(:) points using bispev
293339
function gridded_eval_many(this,x,y,ierr) result(f)
294340
class(fitpack_grid_surface), intent(inout) :: this
295341
real(FP_REAL), intent(in) :: x(:),y(:) ! Evaluation points

‎test/fitpack_curve_tests.f90‎

Lines changed: 77 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@ module fitpack_curve_tests
5252
public :: test_surface_integral
5353
public :: test_cross_section
5454
public :: test_derivative_spline
55+
public :: test_grid_surface_scattered_eval
5556

5657

5758
contains
@@ -1004,8 +1005,8 @@ logical function test_gridded_fit(iunit) result(success)
10041005
exit
10051006
end if
10061007

1007-
! Evaluation of the spline approximation
1008-
fit_z = surf%eval(daregr_x,daregr_y,ierr)
1008+
! Evaluation of the spline approximation on grid
1009+
fit_z = surf%eval_ongrid(daregr_x,daregr_y,ierr)
10091010

10101011
if (.not.FITPACK_SUCCESS(ierr)) then
10111012
success = .false.
@@ -2354,7 +2355,7 @@ logical function test_derivative_spline() result(success)
23542355
return
23552356
end if
23562357

2357-
fgrid = dsurf%eval(xeval, yeval, ierr)
2358+
fgrid = dsurf%eval_ongrid(xeval, yeval, ierr)
23582359
if (.not.FITPACK_SUCCESS(ierr)) then
23592360
print *, '[test_derivative_spline] gridded ds/dx eval failed: ', FITPACK_MESSAGE(ierr)
23602361
return
@@ -2402,7 +2403,7 @@ logical function test_derivative_spline() result(success)
24022403
return
24032404
end if
24042405

2405-
fgrid = dsurf%eval(xeval, yeval, ierr)
2406+
fgrid = dsurf%eval_ongrid(xeval, yeval, ierr)
24062407
if (.not.FITPACK_SUCCESS(ierr)) then
24072408
print *, '[test_derivative_spline] gridded ds/dy eval failed: ', FITPACK_MESSAGE(ierr)
24082409
return
@@ -2464,6 +2465,78 @@ logical function test_derivative_spline() result(success)
24642465

24652466
end function test_derivative_spline
24662467

2468+
!> Test scattered-point evaluation on fitpack_grid_surface
2469+
logical function test_grid_surface_scattered_eval() result(success)
2470+
2471+
type(fitpack_grid_surface) :: gsurf
2472+
integer(FP_FLAG) :: ierr
2473+
real(FP_REAL), parameter :: TOL = 1.0e-10_FP_REAL
2474+
2475+
integer, parameter :: NX = 15, NY = 15
2476+
real(FP_REAL) :: xg(NX), yg(NY), zg(NY,NX)
2477+
real(FP_REAL) :: xp(5), yp(5), fp(5), exact(5)
2478+
real(FP_REAL) :: fgrid(3,3), xeval(3), yeval(3)
2479+
integer :: j
2480+
2481+
success = .false.
2482+
2483+
! Grid data for z = x^2 + y^2 on [0,1] x [0,1]
2484+
xg = linspace(zero, one, NX)
2485+
yg = linspace(zero, one, NY)
2486+
do j = 1, NX
2487+
zg(:,j) = xg(j)**2 + yg**2
2488+
end do
2489+
2490+
ierr = gsurf%new_fit(xg, yg, zg, smoothing=zero)
2491+
if (.not.FITPACK_SUCCESS(ierr)) then
2492+
print *, '[test_grid_surface_scattered_eval] fit failed: ', FITPACK_MESSAGE(ierr)
2493+
return
2494+
end if
2495+
2496+
! --- Scattered eval: arbitrary (x,y) pairs ---
2497+
xp = [0.1_FP_REAL, 0.25_FP_REAL, 0.5_FP_REAL, 0.75_FP_REAL, 0.9_FP_REAL]
2498+
yp = [0.9_FP_REAL, 0.3_FP_REAL, 0.5_FP_REAL, 0.1_FP_REAL, 0.6_FP_REAL]
2499+
exact = xp**2 + yp**2
2500+
2501+
fp = gsurf%eval(xp, yp, ierr)
2502+
if (.not.FITPACK_SUCCESS(ierr)) then
2503+
print *, '[test_grid_surface_scattered_eval] scattered eval failed: ', FITPACK_MESSAGE(ierr)
2504+
return
2505+
end if
2506+
2507+
if (maxval(abs(fp - exact)) > TOL) then
2508+
print *, '[test_grid_surface_scattered_eval] scattered eval error:', maxval(abs(fp - exact))
2509+
return
2510+
end if
2511+
2512+
! --- Single-point eval ---
2513+
if (abs(gsurf%eval(0.3_FP_REAL, 0.7_FP_REAL) - (0.09_FP_REAL + 0.49_FP_REAL)) > TOL) then
2514+
print *, '[test_grid_surface_scattered_eval] single-point eval error'
2515+
return
2516+
end if
2517+
2518+
! --- Gridded eval (eval_ongrid) still works ---
2519+
xeval = linspace(0.2_FP_REAL, 0.8_FP_REAL, 3)
2520+
yeval = linspace(0.2_FP_REAL, 0.8_FP_REAL, 3)
2521+
2522+
fgrid = gsurf%eval_ongrid(xeval, yeval, ierr)
2523+
if (.not.FITPACK_SUCCESS(ierr)) then
2524+
print *, '[test_grid_surface_scattered_eval] eval_ongrid failed: ', FITPACK_MESSAGE(ierr)
2525+
return
2526+
end if
2527+
2528+
! fgrid(j,i) = xeval(i)^2 + yeval(j)^2
2529+
do j = 1, 3
2530+
if (maxval(abs(fgrid(j,:) - (xeval**2 + yeval(j)**2))) > TOL) then
2531+
print *, '[test_grid_surface_scattered_eval] eval_ongrid mismatch at j=', j
2532+
return
2533+
end if
2534+
end do
2535+
2536+
success = .true.
2537+
2538+
end function test_grid_surface_scattered_eval
2539+
24672540
! ODE-style reciprocal error weight
24682541
elemental real(FP_REAL) function rewt(RTOL,ATOL,x)
24692542
real(FP_REAL), intent(in) :: RTOL,ATOL,x

‎test/test.f90‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,7 @@ subroutine run_interface_tests()
6767
call add_test(test_surface_integral())
6868
call add_test(test_cross_section())
6969
call add_test(test_derivative_spline())
70+
call add_test(test_grid_surface_scattered_eval())
7071

7172
end subroutine run_interface_tests
7273

0 commit comments

Comments
 (0)