source: trunk/src/calc_field_interp.m @ 809

Last change on this file since 809 was 809, checked in by g7moreau, 7 years ago
  • Add license
File size: 8.0 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
17%=======================================================================
18% Copyright 2008-2014, LEGI UMR 5519 / CNRS UJF G-INP, Grenoble, France
19%   http://www.legi.grenoble-inp.fr
20%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
21%
22%     This file is part of the toolbox UVMAT.
23%
24%     UVMAT is free software; you can redistribute it and/or modify
25%     it under the terms of the GNU General Public License as published
26%     by the Free Software Foundation; either version 2 of the license,
27%     or (at your option) any later version.
28%
29%     UVMAT is distributed in the hope that it will be useful,
30%     but WITHOUT ANY WARRANTY; without even the implied warranty of
31%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
32%     GNU General Public License (see LICENSE.txt) for more details.
33%=======================================================================
34
35function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
36
37%% initialization
38VarVal={};
39ListVarName={};
40VarAttribute={};
41errormsg='';
42InputVarList={};
43if ischar(FieldName),FieldName={FieldName};end
44check_skipped=zeros(size(FieldName));% default, =1 to mark the variables which can be calculated
45check_interp=ones(size(FieldName));% default, =1 to mark the variables which can be interpolated (not ancillary)
46Operator=cell(size(FieldName));
47
48%% analyse the list of input fields: needed variables and requested operations
49for ilist=1:numel(FieldName)
50    Operator{ilist}='';%default empty operator (vec, norm,...)
51    r=regexp(FieldName{ilist},'(?<Operator>(^vec|^norm|^curl|^div|^strain))\((?<UName>.+),(?<VName>.+)\)$','names');% analyse the field name
52    if isempty(r) % the field name is a variable itself
53        ivar=find(strcmp(FieldName{ilist},Data.ListVarName));
54        if isempty(ivar)
55            check_skipped(ilist)=1; %variable not found
56        elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1));
57            if isfield(Data.VarAttribute{ivar},'Role') &&...
58                (strcmp(Data.VarAttribute{ivar}.Role,'ancillary')||strcmp(Data.VarAttribute{ivar}.Role,'warnflag')||strcmp(Data.VarAttribute{ivar}.Role,'errorflag'))
59                check_interp(ilist)=0; % ancillary variable, not interpolated     
60            end
61            InputVarList=[InputVarList FieldName{ilist}];% the variable is added to the list of input variables if it is not already in the list
62        end
63    else
64        if ~isfield(Data,r.UName)||~isfield(Data,r.VName)%needed input variable not found
65            check_skipped(ilist)=1;
66        elseif strcmp(r.Operator,'curl')||strcmp(r.Operator,'div')||strcmp(r.Operator,'strain')
67            Operator{ilist}=r.Operator;
68            switch r.Operator
69                case 'curl'% case of CivX data format
70                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get curl through linear interp: use ProjMode=interp_tps'; return; end
71                    UName{ilist}='vort';
72                    Data.vort=Data.DjUi(:,1,2)-Data.DjUi(:,2,1);
73                case 'div'
74                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get div through linear interp: use ProjMode=interp_tps'; return; end
75                    UName{ilist}='div';
76                    Data.div=Data.DjUi(:,1,1)+Data.DjUi(:,2,2);
77                case 'strain'
78                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get strain through linear interp: use ProjMode=interp_tps'; return; end
79                    UName{ilist}='strain';
80                    Data.strain=Data.DjUi(:,1,2)+Data.DjUi(:,2,1);
81            end
82            InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
83        else % case  'norm' for instance
84            UName{ilist}=r.UName;
85            VName{ilist}=r.VName;
86            if isempty(find(strcmp(r.UName,InputVarList)));
87                InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
88            end
89            if isempty(find(strcmp(r.VName,InputVarList), 1));
90                InputVarList=[InputVarList VName{ilist}]; %the variable is added to the list if it is not already in the list
91            end
92            Operator{ilist}=r.Operator;
93        end
94    end
95end
96
97%% create interpolator for each variable to interpolate
98if exist('XI','var')
99    for ilist=1:numel(InputVarList)
100        F.(InputVarList{ilist})=TriScatteredInterp(Coord,Data.(InputVarList{ilist}),'linear');
101    end
102end
103
104%% perform the linear interpolation for the requested variables
105for ilist=1:numel(FieldName)
106    if ~check_skipped(ilist)
107        nbvar=numel(ListVarName);
108        switch Operator{ilist}
109            case 'vec'
110                if exist('XI','var')
111                    if check_interp(ilist)
112                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
113                    VarVal{nbvar+2}=F.(VName{ilist})(XI,YI);
114                    end
115                else
116                    VarVal{nbvar+1}=Data.(UName{ilist});
117                    VarVal{nbvar+2}=Data.(VName{ilist});
118                end
119                ListVarName{nbvar+1}=UName{ilist};
120                ListVarName{nbvar+2}=VName{ilist};
121                VarAttribute{nbvar+1}.Role='vector_x';
122                VarAttribute{nbvar+2}.Role='vector_y';
123            case 'norm'
124                if exist('XI','var')
125                    if check_interp(ilist)
126                    U2=F.(UName{ilist})(XI,YI).*F.(UName{ilist})(XI,YI);
127                    V2=F.(VName{ilist})(XI,YI).*F.(VName{ilist})(XI,YI);
128                    end
129                else
130                    U2=Data.(UName{ilist}).*Data.(UName{ilist});
131                    V2=Data.(VName{ilist}).*Data.(VName{ilist});
132                end
133                VarVal{nbvar+1}=sqrt(U2+V2);
134                ListVarName{nbvar+1}='norm';
135                VarAttribute{nbvar+1}.Role='scalar';
136            case {'curl','div','strain'}
137                if exist('XI','var')
138                    if check_interp(ilist)
139                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
140                    end
141                else
142                    VarVal{nbvar+1}=Data.(UName{ilist});
143                end
144                ListVarName{nbvar+1}=UName{ilist};
145                VarAttribute{nbvar+1}.Role='scalar';
146            otherwise
147                if ~isempty(FieldName{ilist})
148                    if exist('XI','var')
149                        if check_interp(ilist)
150                        VarVal{nbvar+1}=F.(FieldName{ilist})(XI,YI);
151                        end
152                    else
153                        VarVal{nbvar+1}= Data.(FieldName{ilist});
154                    end
155                    ListVarName{nbvar+1}=FieldName{ilist};
156                    VarAttribute{nbvar+1}.Role='scalar';
157                end
158        end
159    end
160end
161
162%% put an error flag to indicate NaN data
163if exist('XI','var')&&~isempty(VarVal)
164    nbvar=numel(VarVal);
165    ListVarName{nbvar+1}='FF';
166    VarVal{nbvar+1}=isnan(VarVal{nbvar});
167    VarAttribute{nbvar+1}.Role='errorflag';
168end
169
170
171
172
173
Note: See TracBrowser for help on using the repository browser.