# include # include # include "blend.h" int main ( ); void cubic_rs ( double r, double s, int i, double *xi ); void identity_r ( double r, int i, double *xi ); void identity_rs ( double r, double s, int i, double *xi ); void identity_rst ( double r, double s, double t, int i, double *xi ); void quad_rst ( double r, double s, double t, int i, double *xi ); void stretch_r ( double r, int i, double *xi ); void stretch_rs ( double r, double s, int i, double *xi ); void stretch_rst ( double r, double s, double t, int i, double *xi ); void test01 ( ); void test02 ( ); void test03 ( ); void test04 ( ); void test05 ( ); void test06 ( ); void test07 ( ); void test08 ( ); void test09 ( ); void test10 ( ); void test11 ( ); void test12 ( ); void blend_101_test ( ); void blend_102_test ( ); void blend_103_test ( ); void blend_i_0d1_test ( ); void blend_ij_0d1_test ( ); void blend_ij_1d1_test ( ); void blend_ijk_0d1_test ( ); void blend_ijk_1d1_test ( ); void blend_ijk_2d1_test ( ); /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: MAIN is the main program for BLEND_TEST. Discussion: BLEND_TEST tests the BLEND library. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 October 2018 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "BLEND_TEST\n" ); printf ( " C version\n" ); printf ( " Test the BLEND library.\n" ); test01 ( ); test02 ( ); test03 ( ); test04 ( ); test05 ( ); test06 ( ); test07 ( ); test08 ( ); test09 ( ); test10 ( ); test11 ( ); test12 ( ); blend_101_test ( ); blend_102_test ( ); blend_103_test ( ); blend_i_0d1_test ( ); blend_ij_0d1_test ( ); blend_ij_1d1_test ( ); blend_ijk_0d1_test ( ); blend_ijk_1d1_test ( ); blend_ijk_2d1_test ( ); /* Terminate. */ printf ( "\n" ); printf ( "BLEND_TEST:\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void test01 ( ) /******************************************************************************/ /* Purpose: TEST01 tests BLEND_R_0DN on the identity map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double x[1]; printf ( "\n" ); printf ( "TEST01\n" ); printf ( " Identity test on BLEND_R_0DN.\n" ); n = 1; r = 0.0; blend_r_0dn ( r, x, n, identity_r ); printf ( " %8g %8g\n", r, x[0] ); r = 1.0; blend_r_0dn ( r, x, n, identity_r ); printf ( " %8g %8g\n", r, x[0] ); r = 0.5; blend_r_0dn ( r, x, n, identity_r ); printf ( " %8g %8g\n", r, x[0] ); return; } /******************************************************************************/ void test02 ( ) /******************************************************************************/ /* Purpose: TEST02 tests BLEND_RS_0DN on the identity map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double x[2]; printf ( "\n" ); printf ( "TEST02\n" ); printf ( " Identity test on BLEND_RS_0DN.\n" ); n = 2; r = 0.0; s = 0.0; blend_rs_0dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 1.0; s = 0.0; blend_rs_0dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 0.0; s = 1.0; blend_rs_0dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 1.0; s = 1.0; blend_rs_0dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 0.5; s = 0.5; blend_rs_0dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); return; } /******************************************************************************/ void test03 ( ) /******************************************************************************/ /* Purpose: TEST03 tests BLEND_RS_1DN on the identity map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double x[2]; printf ( "\n" ); printf ( "TEST03\n" ); printf ( " Identity test on BLEND_RS_1DN.\n" ); n = 2; r = 0.0; s = 0.0; blend_rs_1dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 1.0; s = 0.0; blend_rs_1dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 0.0; s = 1.0; blend_rs_1dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 1.0; s = 1.0; blend_rs_1dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 0.5; s = 0.5; blend_rs_1dn ( r, s, x, n, identity_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); return; } /******************************************************************************/ void test04 ( ) /******************************************************************************/ /* Purpose: TEST04 tests BLEND_RST_0DN on the identity map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double t; double x[3]; printf ( "\n" ); printf ( "TEST04\n" ); printf ( " Identity test on BLEND_RST_0DN.\n" ); n = 3; r = 0.0; s = 0.0; t = 0.0; blend_rst_0dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 0.0; t = 0.0; blend_rst_0dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 1.0; t = 0.0; blend_rst_0dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 0.0; t = 1.0; blend_rst_0dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 1.0; t = 1.0; blend_rst_0dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.5; s = 0.5; t = 0.5; blend_rst_0dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); return; } /******************************************************************************/ void test05 ( ) /******************************************************************************/ /* Purpose: TEST05 tests BLEND_RST_1DN on the identity map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double t; double x[3]; printf ( "\n" ); printf ( "TEST05\n" ); printf ( " Identity test on BLEND_RST_1DN.\n" ); n = 3; r = 0.0; s = 0.0; t = 0.0; blend_rst_1dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 0.0; t = 0.0; blend_rst_1dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 1.0; t = 0.0; blend_rst_1dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 0.0; t = 1.0; blend_rst_1dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 1.0; t = 1.0; blend_rst_1dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.5; s = 0.5; t = 0.5; blend_rst_1dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); return; } /******************************************************************************/ void test06 ( ) /******************************************************************************/ /* Purpose: TEST06 tests BLEND_RST_2DN on the identity map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double t; double x[3]; printf ( "\n" ); printf ( "TEST06\n" ); printf ( " Identity test on BLEND_RST_2DN.\n" ); n = 3; r = 0.0; s = 0.0; t = 0.0; blend_rst_2dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 0.0; t = 0.0; blend_rst_2dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 1.0; t = 0.0; blend_rst_2dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 0.0; t = 1.0; blend_rst_2dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 1.0; t = 1.0; blend_rst_2dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.5; s = 0.5; t = 0.5; blend_rst_2dn ( r, s, t, x, n, identity_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); return; } /******************************************************************************/ void test07 ( ) /******************************************************************************/ /* Purpose: TEST07 tests BLEND_R_0DN on the stretch map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double x[1]; printf ( "\n" ); printf ( "TEST07\n" ); printf ( " Stretch test on BLEND_R_0DN.\n" ); n = 1; r = 0.0; blend_r_0dn ( r, x, n, stretch_r ); printf ( " %8g %8g\n", r, x[0] ); r = 1.0; blend_r_0dn ( r, x, n, stretch_r ); printf ( " %8g %8g\n", r, x[0] ); r = 0.5; blend_r_0dn ( r, x, n, stretch_r ); printf ( " %8g %8g\n", r, x[0] ); return; } /******************************************************************************/ void test08 ( ) /******************************************************************************/ /* Purpose: TEST08 tests BLEND_RS_0DN on the stretch map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double x[2]; printf ( "\n" ); printf ( "TEST08\n" ); printf ( " Stretch test on BLEND_RS_0DN.\n" ); n = 2; r = 0.0; s = 0.0; blend_rs_0dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 1.0; s = 0.0; blend_rs_0dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 0.0; s = 1.0; blend_rs_0dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 1.0; s = 1.0; blend_rs_0dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 0.5; s = 0.5; blend_rs_0dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); return; } /******************************************************************************/ void test09 ( ) /******************************************************************************/ /* Purpose: TEST09 tests BLEND_RS_1DN on the stretch map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double x[2]; printf ( "\n" ); printf ( "TEST09\n" ); printf ( " Stretch test on BLEND_RS_1DN.\n" ); n = 2; r = 0.0; s = 0.0; blend_rs_1dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 1.0; s = 0.0; blend_rs_1dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 0.0; s = 1.0; blend_rs_1dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 1.0; s = 1.0; blend_rs_1dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); r = 0.5; s = 0.5; blend_rs_1dn ( r, s, x, n, stretch_rs ); printf ( " %8g %8g %8g %8g\n", r, s, x[0], x[1] ); return; } /******************************************************************************/ void test10 ( ) /******************************************************************************/ /* Purpose: TEST10 tests BLEND_RST_0DN on the stretch map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double t; double x[3]; printf ( "\n" ); printf ( "TEST10\n" ); printf ( " Stretch test on BLEND_RST_0DN.\n" ); n = 3; r = 0.0; s = 0.0; t = 0.0; blend_rst_0dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 0.0; t = 0.0; blend_rst_0dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 1.0; t = 0.0; blend_rst_0dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 0.0; t = 1.0; blend_rst_0dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 1.0; t = 1.0; blend_rst_0dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.5; s = 0.5; t = 0.5; blend_rst_0dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); return; } /******************************************************************************/ void test11 ( ) /******************************************************************************/ /* Purpose: TEST11 tests BLEND_RST_1DN on the stretch map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double t; double x[3]; printf ( "\n" ); printf ( "TEST11\n" ); printf ( " Stretch test on BLEND_RST_1DN.\n" ); n = 3; r = 0.0; s = 0.0; t = 0.0; blend_rst_1dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 0.0; t = 0.0; blend_rst_1dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 1.0; t = 0.0; blend_rst_1dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 0.0; t = 1.0; blend_rst_1dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 1.0; t = 1.0; blend_rst_1dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.5; s = 0.5; t = 0.5; blend_rst_1dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); return; } /******************************************************************************/ void test12 ( ) /******************************************************************************/ /* Purpose: TEST12 tests BLEND_RST_2DN on the stretch map. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int n; double r; double s; double t; double x[3]; printf ( "\n" ); printf ( "TEST12\n" ); printf ( " Stretch test on BLEND_RST_2DN.\n" ); n = 3; r = 0.0; s = 0.0; t = 0.0; blend_rst_2dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 0.0; t = 0.0; blend_rst_2dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 1.0; t = 0.0; blend_rst_2dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.0; s = 0.0; t = 1.0; blend_rst_2dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 1.0; s = 1.0; t = 1.0; blend_rst_2dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); r = 0.5; s = 0.5; t = 0.5; blend_rst_2dn ( r, s, t, x, n, stretch_rst ); printf ( " %8g %8g %8g %8g %8g %8g\n", r, s, t, x[0], x[1], x[2] ); return; } /******************************************************************************/ void blend_101_test ( ) /******************************************************************************/ /* Purpose: BLEND_101_TEST tests BLEND_101. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 October 2018 Author: John Burkardt */ { double *d; int i; int n1; double r; n1 = 10; d = r8vec_zeros_new ( n1 ); for ( i = 0; i < n1; i++ ) { if ( i == 0 || i == n1 - 1 ) { d[i] = ( double ) ( i + 1 ); } } printf ( "\n" ); printf ( "BLEND_101_TEST\n" ); printf ( " BLEND_101 blends endpoint values into a list.\n" ); r8vec_print ( n1, d, " Initial data list" ); for ( i = 0; i < n1; i++ ) { if ( i != 0 && i != n1 - 1 ) { r = ( double ) ( i ) / ( double ) ( n1 - 1 ); d[i] = blend_101 ( r, d[0], d[n1-1] ); } } r8vec_print ( n1, d, " Interpolated data list" ); free ( d ); return; } /******************************************************************************/ void blend_102_test ( ) /******************************************************************************/ /* Purpose: BLEND_102_TEST tests BLEND_102. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 October 2008 Author: John Burkardt */ { double *d; int i; int j; int n1; int n2; double r; double s; n1= 5; n2 = 5; d = r8mat_zeros_new ( n1, n2 ); for ( i = 1; i <= n1; i++ ) { for ( j = 1; j <= n2; j++ ) { if ( ( i == 1 || i == n1 ) && ( j == 1 || j == n2 ) ) { d[i-1+(j-1)*n1] = ( double ) ( i + j ); } } } printf ( "\n" ); printf ( "BLEND_102_TEST\n" ); printf ( " BLEND_102 blends corner values into a table.\n" ); r8mat_print ( n1, n2, d, " Initial data array" ); for ( i = 1; i <= n1; i++ ) { r = ( double ) ( i - 1 ) / ( double ) ( n1 - 1 ); for ( j = 1; j <= n2; j++ ) { s = ( double ) ( j - 1 ) / ( double ) ( n2 - 1 ); if ( ( i == 1 || i == n1 ) && ( j == 1 || j == n2 ) ) { continue; } d[i-1+(j-1)*n1] = blend_102 ( r, s, d[0+0*n1], d[0+(n2-1)*n1], d[n1-1+0*n1], d[n1-1+(n2-1)*n1] ); } } r8mat_print ( n1, n2, d, " Interpolated data array" ); free ( d ); return; } /******************************************************************************/ void blend_103_test ( ) /******************************************************************************/ /* Purpose: BLEND_103_TEST tests BLEND_103. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 October 2008 Author: John Burkardt */ { double *d; int i; int j; int k; int n1; int n2; int n3; double r; double s; double t; n1 = 3; n2 = 5; n3 = 4; d = r8block_zeros_new ( n1, n2, n3 ); for ( i = 1; i <= n1; i++ ) { for ( j = 1; j <= n2; j++ ) { for ( k = 1; k <= n3; k++ ) { if ( ( i == 1 || i == n1 ) && ( j == 1 || j == n2 ) && ( k == 1 || k == n3 ) ) { d[i-1+(j-1)*n1+(k-1)*n1*n2] = ( double ) ( i + j + k ); } } } } printf ( "\n" ); printf ( "BLEND_103_TEST\n" ); printf ( " BLEND_103 blends corner values into a table.\n" ); r8block_print ( n1, n2, n3, d, " Initial data array" ); for ( i = 1; i <= n1; i++ ) { r = ( double ) ( i - 1 ) / ( double ) ( n1 - 1 ); for ( j = 1; j <= n2; j++ ) { s = ( double ) ( j - 1 ) / ( double ) ( n2 - 1 ); for ( k = 1; k <= n3; k++ ) { t = ( double ) ( k - 1 ) / ( double ) ( n3 - 1 ); if ( ( i == 1 || i == n1 ) && ( j == 1 || j == n2 ) && ( k == 1 || k == n3 ) ) { continue; } d[i-1+(j-1)*n1+(k-1)*n1*n2] = blend_103 ( r, s, t, d[ 0+ 0 *n1+ 0 *n1*n2], d[ 0+ 0 *n1+(n3-1)*n1*n2], d[ 0+(n2-1)*n1+ 0 *n1*n2], d[ 0+(n2-1)*n1+(n3-1)*n1*n2], d[n1-1+ 0 *n1+ 0 *n1*n2], d[n1-1+ 0 *n1+(n3-1)*n1*n2], d[n1-1+(n2-1)*n1+ 0 *n1*n2], d[n1-1+(n2-1)*n1+(n3-1)*n1*n2] ); } } } r8block_print ( n1, n2, n3, d, " Interpolated data array" ); free ( d ); return; } /******************************************************************************/ void blend_i_0d1_test ( ) /******************************************************************************/ /* Purpose: BLEND_I_0D1 tests BLEND_I_0D1. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int i; int m; double x[5]; m = 5; x[0] = 100.0; x[m-1] = 100.0 + ( m - 1 ) * 5; printf ( "\n" ); printf ( "BLEND_I_0D1\n" ); printf ( " BLEND_I_0D1 interpolates data in a vector.\n" ); printf ( "\n" ); printf ( " X[0] = %g\n", x[0] ); printf ( " X(%d)= %g\n", m - 1, x[m-1] ); printf ( "\n" ); printf ( " Interpolated values:\n" ); printf ( "\n" ); blend_i_0d1 ( x, m ); for ( i = 0; i < m; i++ ) { printf ( " %6d %8g\n", i, x[i] ); } return; } /******************************************************************************/ void blend_ij_0d1_test ( ) /******************************************************************************/ /* Purpose: BLEND_IJ_0D1_TEST tests BLEND_IJ_0D1. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int i; int j; int m1 = 5; int m2 = 4; double r; double s; double temp; double x[20]; printf ( "\n" ); printf ( "BLEND_IJ_0D1_TEST\n" ); printf ( " BLEND_IJ_0D1 interpolates data in a table,\n" ); printf ( " from corner data.\n" ); printf ( "\n" ); printf ( " The table is %d rows by %d columns.\n", m1, m2 ); /* Load data in the corners only. */ i = 0; j = 0; r = ( double ) i / ( double ) ( m1 - 1 ); s = ( double ) j / ( double ) ( m2 - 1 ); cubic_rs ( r, s, 1, &temp ); x[i*m2+j] = temp; i = m1 - 1; j = 0; r = ( double ) i / ( double ) ( m1 - 1 ); s = ( double ) j / ( double ) ( m2 - 1 ); cubic_rs ( r, s, 1, &temp ); x[i*m2+j] = temp; i = 0; j = m2 - 1; r = ( double ) i / ( double ) ( m1 - 1 ); s = ( double ) j / ( double ) ( m2 - 1 ); cubic_rs ( r, s, 1, &temp ); x[i*m2+j] = temp; i = m1 - 1; j = m2 - 1; r = ( double ) i / ( double ) ( m1 - 1 ); s = ( double ) j / ( double ) ( m2 - 1 ); cubic_rs ( r, s, 1, &temp ); x[i*m2+j] = temp; blend_ij_0d1 ( x, m1, m2 ); printf ( "\n" ); printf ( " Values interpolated by BLEND_IJ_0D1:\n" ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g %8g\n", x[i*m2], x[i*m2+1], x[i*m2+2], x[i*m2+3] ); } return; } /******************************************************************************/ void blend_ij_1d1_test ( ) /******************************************************************************/ /* Purpose: BLEND_IJ_1D1_TEST tests BLEND_IJ_1D1. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int i; int j; int m1 = 5; int m2 = 4; double r; double s; double temp; double x[20]; printf ( "\n" ); printf ( "BLEND_IJ_1D1_TEST\n" ); printf ( " BLEND_IJ_1D1 interpolates data in a table,\n" ); printf ( " from edge data.\n" ); printf ( "\n" ); printf ( " The table is %d rows by %d columns.\n", m1, m2 ); /* Load data in the edges only. */ for ( i = 0; i < m1; i++ ) { r = ( double ) i / ( double ) ( m1 - 1 ); j = 0; s = ( double ) j / ( double ) ( m2 - 1 ); cubic_rs ( r, s, 1, &temp ); x[i*m2+j] = temp; j = m2 - 1; s = ( double ) j / ( double ) ( m2 - 1 ); cubic_rs ( r, s, 1, &temp ); x[i*m2+j] = temp; } for ( j = 0; j < m2; j++ ) { s = ( double ) j / ( double ) ( m2 - 1 ); i = 0; r = ( double ) i / ( double ) ( m1 - 1 ); cubic_rs ( r, s, 1, &temp ); x[i*m2+j] = temp; i = m1 - 1; r = ( double ) i / ( double ) ( m1 - 1 ); cubic_rs ( r, s, 1, &temp ); x[i*m2+j] = temp; } blend_ij_1d1 ( x, m1, m2 ); printf ( "\n" ); printf ( " Values interpolated by BLEND_IJ_1D1:\n" ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g %8g\n", x[i*m2], x[i*m2+1], x[i*m2+2], x[i*m2+3] ); } return; } /******************************************************************************/ void blend_ijk_0d1_test ( ) /******************************************************************************/ /* Purpose: BLEND_IJK_0D1_TEST tests BLEND_IJK_0D1. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int i; int j; int k; int m1 = 4; int m2 = 3; int m3 = 3; int num_extreme; double r; double s; double t; double temp; double x[36]; printf ( "\n" ); printf ( "BLEND_IJK_0D1_TEST\n" ); printf ( " BLEND_IJK_0D1 interpolates data in a 3D table,\n" ); printf ( " from corner data.\n" ); printf ( "\n" ); printf ( " The table is %d rows by %d columns by %d layers.\n", m1, m2, m3 ); /* Load data in the faces only. */ for ( i = 0; i < m1; i++ ) { r = ( double ) i / ( double ) ( m1 - 1 ); for ( j = 0; j < m2; j++ ) { s = ( double ) j / ( double ) ( m2 - 1 ); for ( k = 0; k < m3; k++ ) { t = ( double ) k / ( double ) ( m3 - 1 ); num_extreme = 0; if ( i == 0 || i == m1 - 1 ) { num_extreme = num_extreme + 1; } if ( j == 0 || j == m2 - 1 ) { num_extreme = num_extreme + 1; } if ( k == 0 || k == m3 - 1 ) { num_extreme = num_extreme + 1; } if ( num_extreme >= 3 ) { quad_rst ( r, s, t, 1, &temp ); } else { temp = 0.0; } x[(i*m3+j)*m2+k] = temp; } } } printf ( "\n" ); printf ( " Data given to BLEND_IJK_0D1:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } blend_ijk_0d1 ( x, m1, m2, m3 ); printf ( "\n" ); printf ( " Values interpolated by BLEND_IJK_0D1:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } for ( i = 0; i < m1; i++ ) { r = ( double ) i / ( double ) ( m1 - 1 ); for ( j = 0; j < m2; j++ ) { s = ( double ) j / ( double ) ( m2 - 1 ); for ( k = 0; k < m3; k++ ) { t = ( double ) k / ( double ) ( m3 - 1 ); quad_rst ( r, s, t, 1, &temp ); x[(i*m3+j)*m2+k] = temp; } } } printf ( "\n" ); printf ( " Exact data:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } return; } /******************************************************************************/ void blend_ijk_1d1_test ( ) /******************************************************************************/ /* Purpose: BLEND_IJK_1D1_TEST tests BLEND_IJK_1D1. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int i; int j; int k; int m1 = 4; int m2 = 3; int m3 = 3; int num_extreme; double r; double s; double t; double temp; double x[36]; printf ( "\n" ); printf ( "BLEND_IJK_1D1_TEST\n" ); printf ( " BLEND_IJK_1D1 interpolates data in a 3D table,\n" ); printf ( " from edge data.\n" ); printf ( "\n" ); printf ( " The table is %d rows by %d columns by %d layers.\n", m1, m2, m3 ); /* Load data in the faces only. */ for ( i = 0; i < m1; i++ ) { r = ( double ) i / ( double ) ( m1 - 1 ); for ( j = 0; j < m2; j++ ) { s = ( double ) j / ( double ) ( m2 - 1 ); for ( k = 0; k < m3; k++ ) { t = ( double ) k / ( double ) ( m3 - 1 ); num_extreme = 0; if ( i == 0 || i == m1 - 1 ) { num_extreme = num_extreme + 1; } if ( j == 0 || j == m2 - 1 ) { num_extreme = num_extreme + 1; } if ( k == 0 || k == m3 - 1 ) { num_extreme = num_extreme + 1; } if ( num_extreme >= 2 ) { quad_rst ( r, s, t, 1, &temp ); } else { temp = 0.0; } x[(i*m3+j)*m2+k] = temp; } } } printf ( "\n" ); printf ( " Data given to BLEND_IJK_1D1:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } blend_ijk_1d1 ( x, m1, m2, m3 ); printf ( "\n" ); printf ( " Values interpolated by BLEND_IJK_1D1:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } for ( i = 0; i < m1; i++ ) { r = ( double ) i / ( double ) ( m1 - 1 ); for ( j = 0; j < m2; j++ ) { s = ( double ) j / ( double ) ( m2 - 1 ); for ( k = 0; k < m3; k++ ) { t = ( double ) k / ( double ) ( m3 - 1 ); quad_rst ( r, s, t, 1, &temp ); x[(i*m3+j)*m2+k] = temp; } } } printf ( "\n" ); printf ( " Exact data:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } return; } /******************************************************************************/ void blend_ijk_2d1_test ( ) /******************************************************************************/ /* Purpose: BLEND_IJK_2D1_TEST tests BLEND_IJK_2D1. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { int i; int j; int k; int m1 = 4; int m2 = 3; int m3 = 3; double r; double s; double t; double temp; double x[36]; printf ( "\n" ); printf ( "BLEND_IJK_2D1_TEST\n" ); printf ( " BLEND_IJK_2D1 interpolates data in a 3D table,\n" ); printf ( " from face data.\n" ); printf ( "\n" ); printf ( " The table is %d rows by %d columns by %d layers.\n", m1, m2, m3 ); /* Load data in the faces only. */ for ( i = 0; i < m1; i++ ) { r = ( double ) i / ( double ) ( m1 - 1 ); for ( j = 0; j < m2; j++ ) { s = ( double ) j / ( double ) ( m2 - 1 ); for ( k = 0; k < m3; k++ ) { t = ( double ) k / ( double ) ( m3 - 1 ); if ( i == 0 || i == m1 - 1 || j == 0 || j == m2 - 1 || k == 0 || k == m3 - 1 ) { quad_rst ( r, s, t, 1, &temp ); } else { temp = 0.0; } x[(i*m3+j)*m2+k] = temp; } } } printf ( "\n" ); printf ( " Data given to BLEND_IJK_2D1:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } blend_ijk_2d1 ( x, m1, m2, m3 ); printf ( "\n" ); printf ( " Values interpolated by BLEND_IJK_2D1:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } for ( i = 0; i < m1; i++ ) { r = ( double ) i / ( double ) ( m1 - 1 ); for ( j = 0; j < m2; j++ ) { s = ( double ) j / ( double ) ( m2 - 1 ); for ( k = 0; k < m3; k++ ) { t = ( double ) k / ( double ) ( m3 - 1 ); quad_rst ( r, s, t, 1, &temp ); x[(i*m3+j)*m2+k] = temp; } } } printf ( "\n" ); printf ( " Exact data:\n" ); printf ( "\n" ); for ( k = 0; k < m3; k++ ) { printf ( "\n" ); printf ( " Layer K = %d\n", k ); printf ( "\n" ); for ( i = 0; i < m1; i++ ) { printf ( " %8g %8g %8g\n", x[(i*m3+0)*m2+k], x[(i*m3+1)*m2+k], x[(i*m3+2)*m2+k] ); } } return; } /******************************************************************************/ void cubic_rs ( double r, double s, int i, double *xi ) /******************************************************************************/ /* Purpose: CUBIC_RS evaluates a function of R and S used for some tests. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { *xi = 20.0 * ( r * r * s * s * s ); return; } /******************************************************************************/ void quad_rst ( double r, double s, double t, int i, double *xi ) /******************************************************************************/ /* Purpose: QUAD_RST evaluates a function of R, S and T used for some tests. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt */ { *xi = 18.0 * ( r * r + s + t ); return; } /******************************************************************************/ void identity_r ( double r, int i, double *xi ) /******************************************************************************/ /* Purpose: IDENTITY_R returns a data component given (R). Discussion: This is a dummy routine, which simply returns (R). Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt Parameters: Input, double R, the coordinate of a point that lies on the boundary of the cube. Input, int I, the component of the data which is to be returned. Output, double *XI, the I-th component of the data vector X, evaluated at the point (R), which is on an endpoint of the unit line segment. */ { if ( i == 0 ) { *xi = r; } else { printf ( "\n" ); printf ( "IDENTITY_R - Fatal error!\n" ); printf ( " Illegal component index I = %d\n", i ); *xi = 0.0; } return; } /******************************************************************************/ void identity_rs ( double r, double s, int i, double *xi ) /******************************************************************************/ /* Purpose: IDENTITY_RS returns a data component given (R,S). Discussion: This is a dummy routine, which simply returns (R,S). Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt Parameters: Input, double R, S, the coordinates of a point that lies on the boundary of the square. Input, int I, the component of the data which is to be returned. Output, double *XI, the I-th component of the data vector X, evaluated at the point (R,S), which is on a corner, or edge, of the unit square. */ { if ( i == 0 ) { *xi = r; } else if ( i == 1 ) { *xi = s; } else { printf ( "\n" ); printf ( "IDENTITY_RS - Fatal error!\n" ); printf ( " Illegal component index I = %d\n", i ); *xi = 0.0; } return; } /******************************************************************************/ void identity_rst ( double r, double s, double t, int i, double *xi ) /******************************************************************************/ /* Purpose: IDENTITY_RST returns a data component given (R,S,T). Discussion: This is a dummy routine, which simply returns (R,S,T). Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt Parameters: Input, double R, S, T, the coordinates of a point that lies on the boundary of the cube. Input, int I, the component of the data which is to be returned. Output, double *XI, the I-th component of the data vector X, evaluated at the point (R,S), which is on a corner, edge or face of the unit cube. */ { if ( i == 0 ) { *xi = r; } else if ( i == 1 ) { *xi = s; } else if ( i == 2 ) { *xi = t; } else { printf ( "\n" ); printf ( "IDENTITY_RST - Fatal error!\n" ); printf ( " Illegal component index I = %d\n", i ); *xi = 0.0; } return; } /******************************************************************************/ void stretch_r ( double r, int i, double *xi ) /******************************************************************************/ /* Purpose: STRETCH_R returns a data component given (R). Discussion: This routine shifts by 1 and stretches by 2. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt Parameters: Input, double R, the coordinate of a point that lies on the boundary of the cube. Input, int I, the component of the data which is to be returned. Output, double *XI, the I-th component of the data vector X, evaluated at the point (R), which is on an endpoint of the unit line segment. */ { if ( i == 0 ) { *xi = 2.0 * r + 1.0; } else { printf ( "\n" ); printf ( "STRETCH_R - Fatal error\n" ); printf ( " Illegal component index I = %d\n", i ); *xi = 0.0; } return; } /******************************************************************************/ void stretch_rs ( double r, double s, int i, double *xi ) /******************************************************************************/ /* Purpose: STRETCH_RS returns a data component given (R,S). Discussion: This routine shifts by (1,2) and stretches by (3,4). Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt Parameters: Input, double R, S, the coordinates of a point that lies on the boundary of the square. Input, int I, the component of the data which is to be returned. Output, double *XI, the I-th component of the data vector X, evaluated at the point (R,S), which is on a corner, or edge, of the unit square. */ { if ( i == 0 ) { *xi = 3.0 * r + 1.0; } else if ( i == 1 ) { *xi = 4.0 * s + 2.0; } else { printf ( "\n" ); printf ( "STRETCH_RS - Fatal error!\n" ); printf ( " Illegal component index I = %d\n", i ); *xi = 0.0; } return; } /******************************************************************************/ void stretch_rst ( double r, double s, double t, int i, double *xi ) /******************************************************************************/ /* Purpose: STRETCH_RST returns a data component given (R,S,T). Discussion: This routine shifts by (1,2,3) and stretches by (4,5,6) Licensing: This code is distributed under the GNU LGPL license. Modified: 02 December 2013 Author: John Burkardt Parameters: Input, double R, S, T, the coordinates of a point that lies on the boundary of the cube. Input, int I, the component of the data which is to be returned. Output, double *XI, the I-th component of the data vector X, evaluated at the point (R,S), which is on a corner, edge or face of the unit cube. */ { if ( i == 0 ) { *xi = 4.0 * r + 1.0; } else if ( i == 1 ) { *xi = 5.0 * s + 2.0; } else if ( i == 2 ) { *xi = 6.0 * t + 3.0; } else { printf ( "\n" ); printf ( "STRETCH_RST - Fatal error\n" ); printf ( " Illegal component index I = %d\n", i ); *xi = 0.0; } return; }