source: issm/trunk/src/m/model/plot/processdata.m@ 9719

function [data datatype]=processdata(md,data,options);
%PROCESSDATA - process data to be plotted
%
%   datatype = 1 -> elements
%   datatype = 2 -> nodes
%   datatype = 3 -> node quivers
%   datatype = 4 -> patch
%
%   Usage:
%      [data datatype]=processdata(md,data,options);
%
%   See also: PLOTMODEL, PROCESSMESH

%check format
if (iscell(data) | isempty(data) | length(data)==0 | (length(data)==1 & ~isstruct(data) & isnan(data))),
        error('plotmodel error message: data provided is empty');
end

%Process Patch
if isstruct(data) 
        if (isfield(data,'index') & isfield(data,'value')),
                if data.interpolation(1)==P1Enum(),
                        data=data.value;
                        data=data';
                        data=data(:);
                        datatype=4;
                elseif data.interpolation(1)==P0Enum(),
                        data=data.value;
                        datatype=5;
                else
                        error(['interpolation ' data.interpolation(1) ' not supported yet']);
                end
        else
                error('structure other than Patch not supported yet');
        end
else
        %initialize datatype
        datatype=0;
end

%get datatype
datasize=size(data);

%non patch processing
if datatype~=4 & datatype~=5,

        %transpose data if necessary
        if (size(data,2) > size(data,1)),
                data=data';
        end
        datasize=size(data);

        %convert to double if necessary
        if ~isnumeric(data);
                disp('processdata info message: data is not numeric (logical?). Converted to double');
                data=double(data);
        end

        %check length
        if datasize(1)~=md.numberofnodes & datasize(1)~=md.numberofelements & datasize(1)~=md.numberofnodes*6 & (md.mesh.dimension==3 & ~(datasize(1)==md.numberofelements2d | datasize(1)==md.numberofnodes2d))
                error('plotmodel error message: data not supported yet');
        end

        %quiver?
        if datasize(2)>1,
                datatype=3;

                %check number of columns, add zeros if necessary,
                if (md.mesh.dimension==3)
                        if datasize(2)==2,
                                data=[data, zeros(datasize(1),1)];
                        elseif datasize(2)~=3,
                                error('plotmodel error message: data provided should have 2 or 3 columns for quiver plot, and 1 for regular plot');
                        end
                %elseif ((md.mesh.dimension==2) & datasize(2)~=2),
                %       error('plotmodel error message: data provided should have 2 columns for quiver plot, and 1 for regular plot');
                end
        end

        %treat the case datasize(1)=6*nodes
        if datasize(1)==6*md.numberofnodes
                %keep the only norm of data
                data1=data(1:6:md.numberofnodes*6,:);
                data2=data(2:6:md.numberofnodes*6,:);
                data=sqrt(data1.^2+data2.^2);
                datasize(1)=md.numberofnodes;
                %---> go to node data
        end

        %treat the case datasize(1)=nodes2d
        if (md.mesh.dimension==3 & datasize(1)==md.numberofnodes2d),
                data=project3d(md,data,'node');
                datasize(1)=md.numberofnodes;
                %---> go to node data
        end

        %treat the case datasize(1)=nodes2d
        if (md.mesh.dimension==3 & datasize(1)==md.numberofelements2d),
                data=project3d(md,data,'element');
                datasize(1)=md.numberofelements;
                %---> go to node data
        end

        %smoothing?
        if exist(options,'smooth')
                data=averaging(md,data,getfieldvalue(options,'smooth'));
                datasize(1)=md.numberofnodes;
                %---> go to node data
        end
end

%element data
if (datasize(1)==md.numberofelements & datasize(2)==1),

        %Initialize datatype if non patch
        if datatype~=4 & datatype~=5,
                datatype=1;
        end

        %Mask?
        if exist(options,'mask'),
                flags=getfieldvalue(options,'mask');
                pos=find(~flags);
                if length(flags)==md.numberofnodes,
                        [pos2 dummy]=find(ismember(md.elements,pos));
                        data(pos2,:)=NaN;
                elseif length(flags)==md.numberofelements
                        data(pos,:)=NaN;
                else
                        disp('plotmodel warning: mask length not supported yet (supported length are md.numberofnodes and md.numberofelements');
                end
        end

        %log?
        if exist(options,'log'),
                bounds=getfieldvalue(options,'caxis',[min(data(:)) max(data(:))]);
                data(find(data<bounds(1)))=bounds(1);
                if any(data<=0),
                        error('Log option cannot be applied on negative values. Use caxis option (Rignot''s settings: [1.5 max(data)])');
                end
                pos=find(~isnan(data));
                data(pos)=log(data(pos))/log(getfieldvalue(options,'log'));
        end
end

%node data
if (datasize(1)==md.numberofnodes & datasize(2)==1),
        datatype=2;

        %Mask?
        if exist(options,'mask'),
                flags=getfieldvalue(options,'mask');
                pos=find(~flags);
                if length(flags)==md.numberofnodes,
                        data(pos,:)=NaN;
                elseif length(flags)==md.numberofelements
                        data(md.elements(pos,:),:)=NaN;
                else
                        disp('plotmodel warning: mask length not supported yet (supported length are md.numberofnodes and md.numberofelements');
                end
        end

        %log?
        if exist(options,'log'),
                %if any(data<=0),
                %       error('Log option cannot be applied on negative values. Use caxis option (Rignot''s settings: [1.5 max(data)])');
                %end
                data=log(data)/log(getfieldvalue(options,'log'));
        end
end

%layer projection? 
if getfieldvalue(options,'layer',0)>=1,
        data=project2d(md,data,getfieldvalue(options,'layer')); %project onto 2d mesh
end

%control arrow density if quiverplot
if datatype==3 & exist(options,'density')
        databak=data;
        data=NaN*ones(datasize);
        density=getfieldvalue(options,'density');
        data(1:density:end,:)=databak(1:density:end,:);
        clear databak
end

%OK, if datatype=0 error out
if datatype==0,
        error(['data provided not recognized or not supported']);
end