MATLAB Examples

Interpolation from a regular grid

This is an example of application where the grid where the values are known is a regular distribution of points.

Contents

Creation of the grid where the function is known

% Definition of the geometry
obj_AFP = obj_AFP_V1_00;
obj_AFP.file_geom_definition = 'R_fcts_geom\Files_geom_definition\Cube_3D_centered'; % Shape of the geometry: unit length centered cube
obj_AFP.scale_and_translation_geom_definition = {[10 10 2] [5 0 0]}; % Limits of the geometry: "x": [0 10]; "y": [-5 5]; "z": [-1 1]

% Analysis of the geometry
obj_AFP = Analyse_geometry_definition(obj_AFP);

% Definition of the grid points
Vx = 0:.2:10;
Vy = -5:.2:5;
Vz = -1:.1:1;
[Mx,My,Mz] = meshgrid(Vx,Vy,Vz);
Pts_coord_grid = [Mx(:) My(:) Mz(:)];
obj_AFP.Pts_mean_distance_expr = {[] [.2 .2 .1]}; % This parameter will be used in the interpolation process to determine the size
% of the ponderation functions.
obj_AFP = calc_pts_distr(obj_AFP,'replace',Pts_coord_grid);

Specifications of the points where the values are to be interpolated

% Basic 2D plane
Vix = 1:.1:9;
Viy = -4:.1:4;
Viz = 0;
[Mix,Miy,Miz] = meshgrid(Vix,Viy,Viz);
Pts_coord_interp_basic = [Mix(:) Miy(:) Miz(:)];
Nb_pts_interp = size(Pts_coord_interp_basic,1);
Nb_pts_Vix = length(Vix);
Nb_pts_Viy = length(Viy);

% Rotation of the plane around the coordinate [5 0 0]
Pts_coord_interp_rotated_1 = rotation_pts_by_angle([0 1 0],Pts_coord_interp_basic,-pi/18,[5 0 0]);
Pts_coord_interp_rotated_2 = rotation_pts_by_angle([0 0 1],Pts_coord_interp_rotated_1,pi/20,[5 0 0]);

Specifications of the interpolation parameters

stc_specs.SolOrder = 2;
stc_specs.Kernel_technique = 'FC_variable';
stc_specs.F_Kronecker_enforcement = false;
stc_specs.LargKernel = 1.17*[1 1 1];
stc_specs.Value_points_outside_domain = 'Extrapolation';
stc_specs.L_check_status_inside = false;
stc_specs.L_output_C_NZphih = false;
stc_specs.L_output_C_NZdNh = false;
stc_specs.L_output_C_dNh = true;
stc_specs.L_output_V_interp = false;
stc_specs.L_output_stc_kernel = false;
stc_specs.F_compact_output_cells = false;
stc_specs.L_derivatives_interp_to_compute = true(1,10);
stc_specs.F_use_pts_mean_distance_in_kernel = true;
stc_specs.V_fct_to_interp = [];
stc_specs.F_waitbar = true;

Computation of the interpolation matrix

C_dNh = interp_AFP(obj_AFP,Pts_coord_interp_rotated_2,stc_specs);
M_G1_to_G2_1 = C_dNh{1}; % Direct interpolation matrix
M_G1_to_G2_x = C_dNh{2}; % "x" derivative
M_G1_to_G2_y = C_dNh{3}; % "y" derivative
M_G1_to_G2_z = C_dNh{4}; % "z" derivative
M_G1_to_G2_x2 = C_dNh{5}; % "x^2" derivative
M_G1_to_G2_xy = C_dNh{6}; % "x^y" derivative
M_G1_to_G2_xz = C_dNh{7}; % "x^z" derivative
M_G1_to_G2_y2 = C_dNh{8}; % "y^2" derivative
M_G1_to_G2_yz = C_dNh{9}; % "y^z" derivative
M_G1_to_G2_z2 = C_dNh{10}; % "z^2" derivative

Interpolated values

% Function to interpolate
fct_to_interp = inline('(x-5.8).^3/100 + -2*exp(-(x-6).^2/(2*2.4^2) -(y-2).^2/(2*2.1^2) -(z-.3).^2/(2*0.4^2))','x','y','z');
values_to_interp = fct_to_interp(Pts_coord_grid(:,1),Pts_coord_grid(:,2),Pts_coord_grid(:,3)); % These values represent the
% function known on the grid points.

% Derivative "x" of the function. This is for plotting reference.
fct_to_interp_x = inline('3*(x-5.8).^2/100 + -2*exp(-(x-6).^2/(2*2.4^2) -(y-2).^2/(2*2.1^2) -(z-.3).^2/(2*0.4^2)).*(-2*(x-6)/(2*2.4^2))','x','y','z');
values_to_interp_x = fct_to_interp_x(Pts_coord_grid(:,1),Pts_coord_grid(:,2),Pts_coord_grid(:,3));

values_on_G2_1 = M_G1_to_G2_1*values_to_interp; % Direct interpolation on G2
values_on_G2_x = M_G1_to_G2_x*values_to_interp; % Interpolation on G2 of df/dx
values_on_G2_y = M_G1_to_G2_y*values_to_interp; % Interpolation on G2 of df/dy
values_on_G2_z = M_G1_to_G2_z*values_to_interp; % Interpolation on G2 of df/dz
values_on_G2_x2 = M_G1_to_G2_x2*values_to_interp; % Interpolation on G2 of d2f/dx2
values_on_G2_xy = M_G1_to_G2_xy*values_to_interp; % Interpolation on G2 of d2f/dxdy
values_on_G2_xz = M_G1_to_G2_xz*values_to_interp; % Interpolation on G2 of d2f/dxdz
values_on_G2_y2 = M_G1_to_G2_y2*values_to_interp; % Interpolation on G2 of d2f/dy2
values_on_G2_yz = M_G1_to_G2_yz*values_to_interp; % Interpolation on G2 of d2f/dydz
values_on_G2_z2 = M_G1_to_G2_z2*values_to_interp; % Interpolation on G2 of d2f/dz2

Display of the function

% Color display of the function to interpolate
stc_specs_plotc.Display_colorbar = true;
plotc(obj_AFP,values_to_interp,stc_specs_plotc);
xlabel('x');
ylabel('y');
zlabel('z');
title('Function to interpolate')
% Interpolation on the slice
h_fig = figure;
h_surf = surf(reshape(Pts_coord_interp_rotated_2(:,1),Nb_pts_Viy,Nb_pts_Vix,1),...
              reshape(Pts_coord_interp_rotated_2(:,2),Nb_pts_Viy,Nb_pts_Vix,1),...
              reshape(Pts_coord_interp_rotated_2(:,3),Nb_pts_Viy,Nb_pts_Vix,1),...
              reshape(values_on_G2_1,Nb_pts_Viy,Nb_pts_Vix,1));

% Preparation of isosurface
obj_AFP = isosurface(obj_AFP,[],{0}); % Long computation

% Plotting of the isosurface
stc_specs_isosurfaces.h_fig = h_fig;
stc_specs_isosurfaces.F_display_colorbar = false;
[obj_AFP,V_handles_isosurfaces,~] = isosurface(obj_AFP,values_to_interp,{5},stc_specs_isosurfaces);
set(V_handles_isosurfaces,'FaceAlpha',.7);

colorbar; % We add the colorbar here to keep all the colors of the distribution (not only the colors of the isosurfaces)
set(gca,'CLim',[min(values_to_interp) max(values_to_interp)]); % The colormap of the "surf" distribution should be extended
% to the values on which the isosurfaces are computed to have a correspondence of these colors.

xlabel('x');
ylabel('y');
zlabel('z');
title('Surf: interpolated function')

Display of the "x" derivative

% Color display of the function to interpolate
stc_specs_plotc.Display_colorbar = true;
plotc(obj_AFP,values_to_interp_x,stc_specs_plotc);
xlabel('x');
ylabel('y');
zlabel('z');
title('"x" derivative of the function to interpolate')
% Interpolation on the slice
h_fig = figure;
h_surf = surf(reshape(Pts_coord_interp_rotated_2(:,1),Nb_pts_Viy,Nb_pts_Vix,1),...
              reshape(Pts_coord_interp_rotated_2(:,2),Nb_pts_Viy,Nb_pts_Vix,1),...
              reshape(Pts_coord_interp_rotated_2(:,3),Nb_pts_Viy,Nb_pts_Vix,1),...
              reshape(values_on_G2_x,Nb_pts_Viy,Nb_pts_Vix,1));
% It should be reminded that these values on the slice are obtained by the only knowledge of the basic function, not his
% derivative.

% The determination of neighbors for isosurface plotting was done previously.

% Plotting of the isosurface
stc_specs_isosurfaces.h_fig = h_fig;
stc_specs_isosurfaces.F_display_colorbar = false;
[obj_AFP,V_handles_isosurfaces,~] = isosurface(obj_AFP,values_to_interp_x,{5},stc_specs_isosurfaces);
set(V_handles_isosurfaces,'FaceAlpha',.5);

colorbar; % We add the colorbar here to keep all the colors of the distribution (not only the colors of the isosurfaces)
set(gca,'CLim',[min(values_to_interp_x) max(values_to_interp_x)]); % The colormap of the "surf" distribution should be extended
% to the values on which the isosurfaces are computed to have a correspondence of these colors.

xlabel('x');
ylabel('y');
zlabel('z');
title('Surf: interpolation for the "x" derivative')

Copyright 2013 Mathieu Gendron