source: trunk/src/calc_field_interp.m @ 681

Last change on this file since 681 was 675, checked in by sommeria, 11 years ago

various bugs corrected

File size: 7.1 KB
Line 
1%'calc_field_interp': calculate fields (velocity, vort, div...) using linear interpolation if requested
2%---------------------------------------------------------------------
3% [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
4%
5% OUTPUT:
6% VarVal: array giving the values of the calculated field
7% ListVarName: corresponding list of variable names
8% VarAttribute: corresponding list of variable attributes, each term #ilist is of the form VarAttribute{ilist}.tag=value
9%
10% INPUT:
11% Coord(nbpoints,2): matrix of x,y coordinates of the input data points
12% Data: inputfield structure
13% FieldName: string representing the field to calculate, or cell array of fields (as displayed in uvmat/FieldName)
14% XI, YI: set of x and y coordinates where the fields need to be linearly interpolated,
15%        if XI, YI are missing, there is no interpolation (case of colors in vector plots)
16%
17function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
18
19%% initialization
20VarVal={};
21ListVarName={};
22VarAttribute={};
23errormsg='';
24InputVarList={};
25if ischar(FieldName),FieldName={FieldName};end
26check_skipped=zeros(size(FieldName));% default, =1 to mark the variables which can be calculated
27check_interp=ones(size(FieldName));% default, =1 to mark the variables which can be interpolated (not ancillary)
28Operator=cell(size(FieldName));
29
30%% analyse the list of input fields: needed variables and requested operations
31for ilist=1:numel(FieldName)
32    Operator{ilist}='';%default empty operator (vec, norm,...)
33    r=regexp(FieldName{ilist},'(?<Operator>(^vec|^norm|^curl|^div|^strain))\((?<UName>.+),(?<VName>.+)\)$','names');% analyse the field name
34    if isempty(r) % the field name is a variable itself
35        ivar=find(strcmp(FieldName{ilist},Data.ListVarName));
36        if isempty(ivar)
37            check_skipped(ilist)=1; %variable not found
38        elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1));
39            if isfield(Data.VarAttribute{ivar},'Role') &&...
40                (strcmp(Data.VarAttribute{ivar}.Role,'ancillary')||strcmp(Data.VarAttribute{ivar}.Role,'warnflag')||strcmp(Data.VarAttribute{ivar}.Role,'errorflag'))
41                check_interp(ilist)=0; % ancillary variable, not interpolated     
42            end
43            InputVarList=[InputVarList FieldName{ilist}];% the variable is added to the list of input variables if it is not already in the list
44        end
45    else
46        if ~isfield(Data,r.UName)||~isfield(Data,r.VName)%needed input variable not found
47            check_skipped(ilist)=1;
48        elseif strcmp(r.Operator,'curl')||strcmp(r.Operator,'div')||strcmp(r.Operator,'strain')
49            Operator{ilist}=r.Operator;
50            switch r.Operator
51                case 'curl'% case of CivX data format
52                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get curl through linear interp: use ProjMode=interp_tps'; return; end
53                    UName{ilist}='vort';
54                    Data.vort=Data.DjUi(:,1,2)-Data.DjUi(:,2,1);
55                case 'div'
56                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get div through linear interp: use ProjMode=interp_tps'; return; end
57                    UName{ilist}='div';
58                    Data.div=Data.DjUi(:,1,1)+Data.DjUi(:,2,2);
59                case 'strain'
60                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get strain through linear interp: use ProjMode=interp_tps'; return; end
61                    UName{ilist}='strain';
62                    Data.strain=Data.DjUi(:,1,2)+Data.DjUi(:,2,1);
63            end
64            InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
65        else % case  'norm' for instance
66            UName{ilist}=r.UName;
67            VName{ilist}=r.VName;
68            if isempty(find(strcmp(r.UName,InputVarList)));
69                InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
70            end
71            if isempty(find(strcmp(r.VName,InputVarList), 1));
72                InputVarList=[InputVarList VName{ilist}]; %the variable is added to the list if it is not already in the list
73            end
74            Operator{ilist}=r.Operator;
75        end
76    end
77end
78
79%% create interpolator for each variable to interpolate
80if exist('XI','var')
81    for ilist=1:numel(InputVarList)
82        F.(InputVarList{ilist})=TriScatteredInterp(Coord,Data.(InputVarList{ilist}),'linear');
83    end
84end
85
86%% perform the linear interpolation for the requested variables
87for ilist=1:numel(FieldName)
88    if ~check_skipped(ilist)
89        nbvar=numel(ListVarName);
90        switch Operator{ilist}
91            case 'vec'
92                if exist('XI','var')
93                    if check_interp(ilist)
94                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
95                    VarVal{nbvar+2}=F.(VName{ilist})(XI,YI);
96                    end
97                else
98                    VarVal{nbvar+1}=Data.(UName{ilist});
99                    VarVal{nbvar+2}=Data.(VName{ilist});
100                end
101                ListVarName{nbvar+1}=UName{ilist};
102                ListVarName{nbvar+2}=VName{ilist};
103                VarAttribute{nbvar+1}.Role='vector_x';
104                VarAttribute{nbvar+2}.Role='vector_y';
105            case 'norm'
106                if exist('XI','var')
107                    if check_interp(ilist)
108                    U2=F.(UName{ilist})(XI,YI).*F.(UName{ilist})(XI,YI);
109                    V2=F.(VName{ilist})(XI,YI).*F.(VName{ilist})(XI,YI);
110                    end
111                else
112                    U2=Data.(UName{ilist}).*Data.(UName{ilist});
113                    V2=Data.(VName{ilist}).*Data.(VName{ilist});
114                end
115                VarVal{nbvar+1}=sqrt(U2+V2);
116                ListVarName{nbvar+1}='norm';
117                VarAttribute{nbvar+1}.Role='scalar';
118            case {'curl','div','strain'}
119                if exist('XI','var')
120                    if check_interp(ilist)
121                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
122                    end
123                else
124                    VarVal{nbvar+1}=Data.(UName{ilist});
125                end
126                ListVarName{nbvar+1}=UName{ilist};
127                VarAttribute{nbvar+1}.Role='scalar';
128            otherwise
129                if ~isempty(FieldName{ilist})
130                    if exist('XI','var')
131                        if check_interp(ilist)
132                        VarVal{nbvar+1}=F.(FieldName{ilist})(XI,YI);
133                        end
134                    else
135                        VarVal{nbvar+1}= Data.(FieldName{ilist});
136                    end
137                    ListVarName{nbvar+1}=FieldName{ilist};
138                    VarAttribute{nbvar+1}.Role='scalar';
139                end
140        end
141    end
142end
143
144%% put an error flag to indicate NaN data
145if exist('XI','var')&&~isempty(VarVal)
146    nbvar=numel(VarVal);
147    ListVarName{nbvar+1}='FF';
148    VarVal{nbvar+1}=isnan(VarVal{nbvar});
149    VarAttribute{nbvar+1}.Role='errorflag';
150end
151
152
153
154
155
Note: See TracBrowser for help on using the repository browser.