griddata_mesh_to_mesh

PURPOSE ^

GRIDDATA_MESH_TO_MESH - interpolates data from a 2d mesh to a 2d mesh

SYNOPSIS ^

function data_prime=griddata_mesh_to_mesh(index,x,y,data,x_prime,y_prime,varargin);

DESCRIPTION ^

GRIDDATA_MESH_TO_MESH - interpolates data from a  2d mesh to a 2d mesh

   griddata_mesh_to_mesh: interpolate data (from mesh with triangulation index,x,y)
   into mesh with triangulation (x_prime,y_prime).

   Usage:
      data_prime=griddata_mesh_to_mesh(index,x,y,data,x_prime,y_prime,varargin)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function data_prime=griddata_mesh_to_mesh(index,x,y,data,x_prime,y_prime,varargin);
0002 %GRIDDATA_MESH_TO_MESH - interpolates data from a  2d mesh to a 2d mesh
0003 %
0004 %   griddata_mesh_to_mesh: interpolate data (from mesh with triangulation index,x,y)
0005 %   into mesh with triangulation (x_prime,y_prime).
0006 %
0007 %   Usage:
0008 %      data_prime=griddata_mesh_to_mesh(index,x,y,data,x_prime,y_prime,varargin)
0009 
0010     global client_server_mode IMDataCounter uset
0011 
0012     %restrict mesh_prime grids to the ones which are inside the mesh hull
0013     in=inhull([x_prime y_prime],[x y]);
0014     pos_hull=find(in);
0015 
0016     %keep only mesh_prime grids inside the hull
0017     x_prime_hull=x_prime(pos_hull);
0018     y_prime_hull=y_prime(pos_hull);
0019 
0020     %call griddata_mesh_to_mesh_hulled on this new sub mesh from mesh 2
0021     if strcmpi(client_server_mode,'yes'),
0022         r_data_prime_hull=GriddataMeshToMeshHulled_client(index,x,y,data,x_prime_hull,y_prime_hull);
0023         data_prime_hull=IMdb('select matrix from r_data_prime_hull');
0024         IMdb('drop r_data_prime_hull');
0025     else
0026         data_prime_hull=griddata_mesh_to_mesh_hulled(index,x,y,data,x_prime_hull,y_prime_hull,varargin);
0027     end
0028 
0029     %return output
0030     data_prime=zeros(length(x_prime),1); data_prime(:)=NaN;
0031     data_prime(pos_hull)=data_prime_hull;
0032 end
0033 
0034 function data_prime=griddata_mesh_to_mesh_hulled(index,x,y,data,x_prime,y_prime,varargin);
0035 
0036     %Analyse data: is it given per node or per element?
0037     if isempty(varargin{1})
0038         element_data=(length(data)==size(index,1));%data is given for a given element
0039         grid_data=(length(data)==length(x));       %data is given for a given grid
0040         if element_data==0 & grid_data==0, error('griddata_mesh_to_mesh_hulled error message: unknown size of data, specify element or node'), end
0041     else
0042         if strcmpi(varargin{1},'node')
0043             element_data=0;
0044             grid_data=1;
0045         elseif strcmpi(varargin{1},'element')
0046             element_data=1;
0047             grid_data=0;
0048         else
0049             error('griddata_mesh_to_mesh_hulled error message: wrong input argument. Must be ''node'' or ''element''')
0050         end
0051     end
0052 
0053     %Some parameters and initialization of data_prime
0054     nel=length(index);
0055     data_prime=zeros(length(x_prime),1); data_prime(:)=NaN;
0056 
0057     %compute some quantities that will speed up the process
0058     x3x1=x(index(:,1))-x(index(:,3));
0059     y3y1=y(index(:,1))-y(index(:,3));
0060     x3x2=x(index(:,2))-x(index(:,3));
0061     y3y2=y(index(:,2))-y(index(:,3));
0062     x3=x(index(:,3));
0063     y3=y(index(:,3));
0064 
0065     %triangles jacobian determinant (2*area)
0066     delta=x(index(:,2)).*y(index(:,3))-y(index(:,2)).*x(index(:,3))-x(index(:,1)).*y(index(:,3))+y(index(:,1)).*x(index(:,3))+x(index(:,1)).*y(index(:,2))-y(index(:,1)).*x(index(:,2));
0067 
0068     %Get area coordinates:
0069     for n=1:nel,
0070         if mod(n,10000)==0,
0071             disp(sprintf('\r%s%.2f%s','      interpolation progress:  ',n/nel*100,' %'));
0072         end
0073         %first area coordinate
0074         area_1=(y3y2(n)*(x_prime-x3(n))-x3x2(n)*(y_prime-y3(n)))/delta(n);
0075         %second area coordinate
0076         area_2=(x3x1(n)*(y_prime-y3(n))-y3y1(n)*(x_prime-x3(n)))/delta(n);
0077         %third area coordinate
0078         area_3=1-area_1-area_2;
0079         %get grids in mesh 2 that have all area coordinates positive (meaning they belong to mesh 1 triangle n:
0080         pos=find((area_1>=0) & (area_2>=0) & (area_3>=0));
0081         if ~isempty(pos),
0082             if grid_data
0083                 %mesh 2 grids in pos belong to mesh 1 element n. interpolate value linearly.
0084                 data_prime(pos)=area_1(pos)*data(index(n,1))+area_2(pos)*data(index(n,2))+area_3(pos)*data(index(n,3));
0085             elseif element_data
0086                 data_prime(pos)=data(n);
0087             end
0088         end
0089     end
0090     if nel>1000,
0091         disp(sprintf('\r%s%.2f%s','      interpolation progress:  ',100,' %'));
0092     end
0093 end

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003