SectionValues

PURPOSE ^

SECTIONVALUES - compute the value of a field on a section

SYNOPSIS ^

function [index,X,Y,Z,S,data_interp]=SectionValues(md,data,infile,resolution)

DESCRIPTION ^

SECTIONVALUES - compute the value of a field on a section

   This routine gets the value of a given field of the model on points
   given by filname (Argus type file)

   Usage:
      [elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [index,X,Y,Z,S,data_interp]=SectionValues(md,data,infile,resolution)
0002 %SECTIONVALUES - compute the value of a field on a section
0003 %
0004 %   This routine gets the value of a given field of the model on points
0005 %   given by filname (Argus type file)
0006 %
0007 %   Usage:
0008 %      [elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution)
0009 
0010 %read infile:
0011 contempt=expread(infile,1);
0012 nods=contempt.nods;
0013 x=contempt.x;
0014 y=contempt.y;
0015 
0016 %get the specified resolution
0017 if isnumeric(resolution(1))
0018     res_h=resolution(1);
0019 else
0020     error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]')
0021 end
0022 if strcmpi(md.type,'3d')
0023     if (length(resolution)==2 & isnumeric(resolution(2)))
0024         res_v=resolution(2);
0025     else
0026         error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]')
0027     end
0028 end
0029 
0030 %initialization
0031 X=[]; %X-coordinate
0032 Y=[]; %Y-coordinate
0033 S=0;  %curvilinear coordinate
0034 
0035 for i=1:nods-1
0036 
0037     x_start=x(i);
0038     x_end=x(i+1);
0039     y_start=y(i);
0040     y_end=y(i+1);
0041     s_start=S(end);
0042 
0043     length_segment=sqrt((x_end-x_start)^2+(y_end-y_start)^2);
0044     portion=ceil(length_segment/res_h);
0045 
0046     x_segment=zeros(portion,1);
0047     y_segment=zeros(portion,1);
0048     s_segment=zeros(portion,1);
0049 
0050     for j=1:portion
0051         x_segment(j)=x_start+(j-1)*(x_end-x_start)/portion;
0052         y_segment(j)=y_start+(j-1)*(y_end-y_start)/portion;
0053         s_segment(j)=s_start+j*length_segment/portion;
0054     end
0055 
0056     %plug into X and Y
0057     X=[X;x_segment];
0058     Y=[Y;y_segment];
0059     S=[S;s_segment];
0060 end
0061 X(end+1)=x(nods);
0062 Y(end+1)=y(nods);
0063 
0064 %Number of grids:
0065 numberofgrids=size(X,1);
0066 
0067 %Compute Z
0068 Z=zeros(numberofgrids,1);
0069 
0070 %New mesh and Data interpolation
0071 if strcmpi(md.type,'2d')
0072 
0073     %Interpolation of data on specified points
0074     data_interp=griddata_mesh_to_mesh(md.elements,md.x,md.y,data,X,Y);
0075 
0076     %Compute index
0077     index=[1:1:(numberofgrids-1);2:1:numberofgrids]';
0078 
0079 else
0080 
0081     %vertically extrude mesh
0082 
0083     %Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
0084     offset=10^-10;
0085     bed=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,md.bed,X,Y,Z,'node')+offset;
0086     surface=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,md.surface,X,Y,Z,'node')-offset;
0087 
0088     %Some useful parameters
0089     layers=ceil(res_v/md.numlayers);
0090     gridsperlayer=numberofgrids;
0091     gridstot=gridsperlayer*layers;
0092     elementsperlayer=gridsperlayer-1;
0093     elementstot=(gridsperlayer-1)*(layers-1);
0094 
0095     %initialization
0096     X3=zeros(gridsperlayer*layers,1); Y3=zeros(gridsperlayer*layers,1); Z3=zeros(gridsperlayer*layers,1); S3=zeros(gridsperlayer*layers,1); index3=zeros(elementstot,4);
0097 
0098     %Get new coordinates in 3d
0099     for i=1:layers
0100         X3(i:layers:end)=X;
0101         Y3(i:layers:end)=Y;
0102         Z3(i:layers:end)=bed+(i-1)*(surface-bed)/(layers-1);
0103         S3(i:layers:end)=S;
0104 
0105         if i<layers %Build index3 with quads
0106             index3((i-1)*elementsperlayer+1:i*elementsperlayer,:)=[i:layers:gridstot-layers; i+1:layers:gridstot-layers; i+layers+1:layers:gridstot; i+layers:layers:gridstot]';
0107         end
0108     end
0109 
0110     %Interpolation of data on specified points
0111     data_interp=griddata_mesh_to_mesh_3d(md.elements,md.x,md.y,md.z,data,X3,Y3,Z3);
0112 
0113     %build outputs
0114     X=X3; Y=Y3; Z=Z3;  S=S3; index=index3;
0115 end

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