source: trunk/src/calc_field_interp.m @ 648

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

get_field updated, several bugs corrected,open_uvmat suppressd

File size: 6.1 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
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
27Operator=cell(size(FieldName));
28
29%% analyse the list of input fields: needed variables and requested operations
30for ilist=1:numel(FieldName)
31    Operator{ilist}='';%default empty operator (vec, norm,...)
32    r=regexp(FieldName{ilist},'(?<Operator>(^vec|^norm|^curl|^div|^strain))\((?<UName>.+),(?<VName>.+)\)$','names');% analyse the field name
33    if isempty(r) % the field name is a variable itself
34        ivar=strcmp(FieldName{ilist},Data.ListVarName);
35        if isempty(ivar)
36            check_skipped(ilist)=1; %variable not found
37        else
38            if isfield(Data.VarAttribute{ivar},'Role')&&strcmp(Data.VarAttribute{ivar}.Role,'ancillary')
39                check_skipped(ilist)=1; % ancillary variable, not found interpolated     
40            elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1));
41                InputVarList=[InputVarList FieldName{ilist}];% the variable is added to the list of input variables if it is not already in the list
42            end
43        end
44    else
45        if ~isfield(Data,r.UName)||~isfield(Data,r.VName)%needed input variable not found
46            check_skipped(ilist)=1;
47        elseif strcmp(r.Operator,'curl')||strcmp(r.Operator,'div')||strcmp(r.Operator,'strain')
48            Operator{ilist}=r.Operator;
49            switch r.Operator
50                case 'curl'% case of CivX data format
51                    UName{ilist}='vort';
52                    Data.vort=Data.DjUi(:,1,2)-Data.DjUi(:,2,1);
53                case 'div'
54                    UName{ilist}='div';
55                    Data.div=Data.DjUi(:,1,1)+Data.DjUi(:,2,2);
56                case 'strain'
57                    UName{ilist}='strain';
58                    Data.strain=Data.DjUi(:,1,2)+Data.DjUi(:,2,1);
59            end
60            InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
61        else
62            UName{ilist}=r.UName;
63            VName{ilist}=r.VName;
64            if isempty(find(strcmp(r.UName,InputVarList)));
65                InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
66            end
67            if isempty(find(strcmp(r.VName,InputVarList), 1));
68                InputVarList=[InputVarList VName{ilist}]; %the variable is added to the list if it is not already in the list
69            end
70            Operator{ilist}=r.Operator;
71        end
72    end
73end
74
75%% create interpolator for each variable to interpolate
76if exist('XI','var')
77    for ilist=1:numel(InputVarList)
78        F.(InputVarList{ilist})=TriScatteredInterp(Coord,Data.(InputVarList{ilist}),'linear');
79    end
80end
81
82%% perform the linear interpolation for the requested variables
83for ilist=1:numel(FieldName)
84    if ~check_skipped(ilist)
85        nbvar=numel(ListVarName);
86        switch Operator{ilist}
87            case 'vec'
88                if exist('XI','var')
89                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
90                    VarVal{nbvar+2}=F.(VName{ilist})(XI,YI);
91                else
92                    VarVal{nbvar+1}=Data.(UName{ilist});
93                    VarVal{nbvar+2}=Data.(VName{ilist});
94                end
95                ListVarName{nbvar+1}=UName{ilist};
96                ListVarName{nbvar+2}=VName{ilist};
97                VarAttribute{nbvar+1}.Role='vector_x';
98                VarAttribute{nbvar+2}.Role='vector_y';
99            case 'norm'
100                if exist('XI','var')
101                    U2=F.(UName{ilist})(XI,YI).*F.(UName{ilist})(XI,YI);
102                    V2=F.(VName{ilist})(XI,YI).*F.(VName{ilist})(XI,YI);
103                else
104                    U2=Data.(UName{ilist}).*Data.(UName{ilist});
105                    V2=Data.(VName{ilist}).*Data.(VName{ilist});
106                end
107                VarVal{nbvar+1}=sqrt(U2+V2);
108                ListVarName{nbvar+1}='norm';
109                VarAttribute{nbvar+1}.Role='scalar';
110            case {'curl','div','strain'}
111                if exist('XI','var')
112                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
113                else
114                    VarVal{nbvar+1}=Data.(UName{ilist});
115                end
116                ListVarName{nbvar+1}=UName{ilist};
117                VarAttribute{nbvar+1}.Role='scalar';
118            otherwise
119                if ~isempty(FieldName{ilist})
120                    if exist('XI','var')
121                        VarVal{nbvar+1}=F.(FieldName{ilist})(XI,YI);
122                    else
123                        VarVal{nbvar+1}= Data.(FieldName{ilist});
124                    end
125                    ListVarName{nbvar+1}=FieldName{ilist};
126                    VarAttribute{nbvar+1}.Role='scalar';
127                end
128        end
129    end
130end
131
132%% put an error flag to indicate NaN data
133if exist('XI','var')&&~isempty(VarVal)
134    nbvar=numel(ListVarName);
135    ListVarName{nbvar+1}='FF';
136    VarVal{nbvar+1}=isnan(VarVal{nbvar});
137    VarAttribute{nbvar+1}.Role='errorflag';
138end
139
140
141
142
143
Note: See TracBrowser for help on using the repository browser.