source:trunk/src/calc_field_interp.m@579

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

bugs corrected for reading vort div

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