griddata_mesh_to_mesh_3d

PURPOSE ^

GRIDDATA_MESH_TO_MESH_3D - interpolates data from a 3d mesh to a 3d mesh

SYNOPSIS ^

function data_prime=griddata_mesh_to_mesh_3d(index,x,y,z,data,x_prime,y_prime,z_prime,varargin);

DESCRIPTION ^

GRIDDATA_MESH_TO_MESH_3D - interpolates data from a 3d mesh to a 3d mesh

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

   Usage:
      data_prime=griddata_mesh_to_mesh_3d(index,x,y,z,data,x_prime,y_prime,z_prime,varargin)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function data_prime=griddata_mesh_to_mesh_3d(index,x,y,z,data,x_prime,y_prime,z_prime,varargin);
0002 %GRIDDATA_MESH_TO_MESH_3D - interpolates data from a 3d mesh to a 3d mesh
0003 %
0004 %   griddata_mesh_to_mesh: interpolate data (from mesh with triangulation index,x,y,z)
0005 %   into mesh with triangulation (x_prime,y_prime,z_prime).
0006 %
0007 %   Usage:
0008 %      data_prime=griddata_mesh_to_mesh_3d(index,x,y,z,data,x_prime,y_prime,z_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 z_prime],[x y z]);
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     z_prime_hull=z_prime(pos_hull);
0020 
0021     %call griddata_mesh_to_mesh_hulled on this new sub mesh from mesh 2
0022     if strcmpi(client_server_mode,'yes'),
0023         r_data_prime_hull=GriddataMeshToMeshHulled_client(index,x,y,data,x_prime_hull,y_prime_hull);
0024         data_prime_hull=IMdb('select matrix from r_data_prime_hull');
0025         IMdb('drop r_data_prime_hull');
0026     else
0027         data_prime_hull=griddata_mesh_to_mesh_hulled(index,x,y,z,data,x_prime_hull,y_prime_hull,z_prime_hull,varargin);
0028     end
0029 
0030     %return output
0031     data_prime=zeros(length(x_prime),1); data_prime(:)=NaN;
0032     data_prime(pos_hull)=data_prime_hull;
0033 end
0034 
0035 function data_prime=griddata_mesh_to_mesh_hulled(index,x,y,z,data,x_prime,y_prime,z_prime,varargin);
0036 
0037     %Analyse data: is it given per node or per element?
0038     if isempty(varargin{1})
0039         element_data=(length(data)==size(index,1));%data is given for a given element
0040         grid_data=(length(data)==length(x));       %data is given for a given grid
0041         if element_data==0 & grid_data==0, error('griddata_mesh_to_mesh_hulled error message: unknown size of data, specify element or node'), end
0042     else
0043         if strcmpi(varargin{1},'node')
0044             element_data=0;
0045             grid_data=1;
0046         elseif strcmpi(varargin{1},'element')
0047             element_data=1;
0048             grid_data=0;
0049         else
0050             error('griddata_mesh_to_mesh_hulled error message: wrong input argument. Must be ''node'' or ''element''')
0051         end
0052     end
0053 
0054     %Some parameters and initialization of data_prime
0055     nel=length(index);
0056     data_prime=zeros(length(x_prime),1); data_prime(:)=NaN;
0057 
0058     %compute some quantities that will speed up the process
0059     x3x1=x(index(:,1))-x(index(:,3));
0060     y3y1=y(index(:,1))-y(index(:,3));
0061     x3x2=x(index(:,2))-x(index(:,3));
0062     y3y2=y(index(:,2))-y(index(:,3));
0063     x3=x(index(:,3));
0064     y3=y(index(:,3));
0065 
0066     %triangles jacobian determinant (2*area)
0067     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));
0068 
0069     %Get area coordinates:
0070     for n=1:nel,
0071         if mod(n,10000)==0,
0072             disp(sprintf('\r%s%.2f%s','      interpolation progress:  ',n/nel*100,' %'));
0073         end
0074         %first area coordinate
0075         area_1=(y3y2(n)*(x_prime-x3(n))-x3x2(n)*(y_prime-y3(n)))/delta(n);
0076         %second area coordinate
0077         area_2=(x3x1(n)*(y_prime-y3(n))-y3y1(n)*(x_prime-x3(n)))/delta(n);
0078         %third area coordinate
0079         area_3=1-area_1-area_2;
0080 
0081         %get grids in mesh prime that have all area coordinates positive (meaning they are on the same extruded column):
0082         pos_horiz=find((area_1>=0) & (area_2>=0) & (area_3>=0));
0083         if ~isempty(pos_horiz),
0084 
0085             %Figure out in which extruded penta it is located
0086             lower=area_1(pos_horiz)*z(index(n,1))+area_2(pos_horiz)*z(index(n,2))+area_3(pos_horiz)*z(index(n,3));
0087             upper=area_1(pos_horiz)*z(index(n,4))+area_2(pos_horiz)*z(index(n,5))+area_3(pos_horiz)*z(index(n,6));
0088             xi=2*(z_prime(pos_horiz)-lower)./(upper-lower)-1;
0089             pos=find(z_prime(pos_horiz)>=lower & z_prime(pos_horiz)<=upper);
0090             pos_vert=pos_horiz(pos);
0091             xi=xi(pos);
0092 
0093             if ~isempty(pos_vert)
0094                 if grid_data
0095                     %mesh 2 grids in pos belong to mesh 1 element n. interpolate value linearly.
0096                     data_prime(pos_vert)=area_1(pos_vert).*(1-xi)/2*data(index(n,1))+area_2(pos_vert).*(1-xi)/2*data(index(n,2))+area_3(pos_vert).*(1-xi)/2*data(index(n,3))...
0097                             +area_1(pos_vert).*(1+xi)/2*data(index(n,4))+area_2(pos_vert).*(1+xi)/2*data(index(n,5))+area_3(pos_vert).*(1+xi)/2*data(index(n,6));
0098                 elseif element_data
0099                     data_prime(pos_vert)=data(n);
0100                 end
0101             end
0102         end
0103     end
0104     if nel>1000,
0105         disp(sprintf('\r%s%.2f%s','      interpolation progress:  ',100,' %'));
0106     end
0107 end

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