source: trunk/src/calc_field_interp.m @ 1002

Last change on this file since 1002 was 1002, checked in by sommeria, 8 years ago

correlation_time and correlation_x corrected

File size: 8.2 KB
RevLine 
[654]1%'calc_field_interp': calculate fields (velocity, vort, div...) using linear interpolation if requested
[515]2%---------------------------------------------------------------------
[579]3% [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
[515]4%
5% OUTPUT:
[575]6% VarVal: array giving the values of the calculated field
[579]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
[515]9%
10% INPUT:
[654]11% Coord(nbpoints,2): matrix of x,y coordinates of the input data points
[579]12% Data: inputfield structure
13% FieldName: string representing the field to calculate, or cell array of fields (as displayed in uvmat/FieldName)
[654]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)
[809]16
17%=======================================================================
[977]18% Copyright 2008-2017, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
[809]19%   http://www.legi.grenoble-inp.fr
20%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
[654]21%
[809]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
[579]35function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
[515]36
[644]37%% initialization
[579]38VarVal={};
[515]39ListVarName={};
40VarAttribute={};
41errormsg='';
[517]42InputVarList={};
[579]43if ischar(FieldName),FieldName={FieldName};end
44check_skipped=zeros(size(FieldName));% default, =1 to mark the variables which can be calculated
[654]45check_interp=ones(size(FieldName));% default, =1 to mark the variables which can be interpolated (not ancillary)
[579]46Operator=cell(size(FieldName));
[644]47
48%% analyse the list of input fields: needed variables and requested operations
[579]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
[863]52    if isempty(r) % no operator: the field name is a variable itself
[650]53        ivar=find(strcmp(FieldName{ilist},Data.ListVarName));
[863]54        if isempty(ivar)% the requested variable does not exist
[579]55            check_skipped(ilist)=1; %variable not found
[863]56        elseif isempty(find(strcmp(FieldName{ilist},InputVarList), 1));% the variable exists and has not been already selected
[864]57            if exist('XI','var')&& isfield(Data.VarAttribute{ivar},'Role') &&...
[863]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                check_skipped(ilist)=1; %variable not used
61            else
62                InputVarList=[InputVarList FieldName{ilist}];% the variable is added to the list of input variables
[521]63            end
[517]64        end
65    else
[579]66        if ~isfield(Data,r.UName)||~isfield(Data,r.VName)%needed input variable not found
[521]67            check_skipped(ilist)=1;
[579]68        elseif strcmp(r.Operator,'curl')||strcmp(r.Operator,'div')||strcmp(r.Operator,'strain')
69            Operator{ilist}=r.Operator;
70            switch r.Operator
71                case 'curl'% case of CivX data format
[675]72                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get curl through linear interp: use ProjMode=interp_tps'; return; end
[579]73                    UName{ilist}='vort';
74                    Data.vort=Data.DjUi(:,1,2)-Data.DjUi(:,2,1);
75                case 'div'
[675]76                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get div through linear interp: use ProjMode=interp_tps'; return; end
[579]77                    UName{ilist}='div';
78                    Data.div=Data.DjUi(:,1,1)+Data.DjUi(:,2,2);
79                case 'strain'
[675]80                    if ~isfield(Data,'DjUi'), errormsg='field DjUi needed to get strain through linear interp: use ProjMode=interp_tps'; return; end
[579]81                    UName{ilist}='strain';
82                    Data.strain=Data.DjUi(:,1,2)+Data.DjUi(:,2,1);
83            end
[1002]84            InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if iTriScatteredInterpt is not already in the list
[988]85        else % case 'norm' for instance
[521]86            UName{ilist}=r.UName;
87            VName{ilist}=r.VName;
88            if isempty(find(strcmp(r.UName,InputVarList)));
89                InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
90            end
[579]91            if isempty(find(strcmp(r.VName,InputVarList), 1));
[521]92                InputVarList=[InputVarList VName{ilist}]; %the variable is added to the list if it is not already in the list
93            end
94            Operator{ilist}=r.Operator;
[517]95        end
[515]96    end
97end
[644]98
99%% create interpolator for each variable to interpolate
[517]100if exist('XI','var')
101    for ilist=1:numel(InputVarList)
102        F.(InputVarList{ilist})=TriScatteredInterp(Coord,Data.(InputVarList{ilist}),'linear');
103    end
[515]104end
[644]105
106%% perform the linear interpolation for the requested variables
[579]107for ilist=1:numel(FieldName)
[521]108    if ~check_skipped(ilist)
[533]109        nbvar=numel(ListVarName);
110        switch Operator{ilist}
111            case 'vec'
[517]112                if exist('XI','var')
[667]113                    if check_interp(ilist)
[533]114                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
115                    VarVal{nbvar+2}=F.(VName{ilist})(XI,YI);
[654]116                    end
[517]117                else
[533]118                    VarVal{nbvar+1}=Data.(UName{ilist});
119                    VarVal{nbvar+2}=Data.(VName{ilist});
[517]120                end
[533]121                ListVarName{nbvar+1}=UName{ilist};
122                ListVarName{nbvar+2}=VName{ilist};
123                VarAttribute{nbvar+1}.Role='vector_x';
124                VarAttribute{nbvar+2}.Role='vector_y';
125            case 'norm'
126                if exist('XI','var')
[667]127                    if check_interp(ilist)
[533]128                    U2=F.(UName{ilist})(XI,YI).*F.(UName{ilist})(XI,YI);
129                    V2=F.(VName{ilist})(XI,YI).*F.(VName{ilist})(XI,YI);
[654]130                    end
[533]131                else
132                    U2=Data.(UName{ilist}).*Data.(UName{ilist});
133                    V2=Data.(VName{ilist}).*Data.(VName{ilist});
134                end
135                VarVal{nbvar+1}=sqrt(U2+V2);
136                ListVarName{nbvar+1}='norm';
[517]137                VarAttribute{nbvar+1}.Role='scalar';
[579]138            case {'curl','div','strain'}
139                if exist('XI','var')
[667]140                    if check_interp(ilist)
[579]141                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
[654]142                    end
[579]143                else
144                    VarVal{nbvar+1}=Data.(UName{ilist});
145                end
146                ListVarName{nbvar+1}=UName{ilist};
147                VarAttribute{nbvar+1}.Role='scalar';
[533]148            otherwise
[579]149                if ~isempty(FieldName{ilist})
[533]150                    if exist('XI','var')
[667]151                        if check_interp(ilist)
[579]152                        VarVal{nbvar+1}=F.(FieldName{ilist})(XI,YI);
[654]153                        end
[533]154                    else
[579]155                        VarVal{nbvar+1}= Data.(FieldName{ilist});
[533]156                    end
[579]157                    ListVarName{nbvar+1}=FieldName{ilist};
[533]158                    VarAttribute{nbvar+1}.Role='scalar';
159                end
160        end
[515]161    end
162end
[644]163
164%% put an error flag to indicate NaN data
[863]165% if exist('XI','var')&&~isempty(VarVal)
166%     nbvar=numel(VarVal);
167%     ListVarName{nbvar+1}='FF';
168%     VarVal{nbvar+1}=isnan(VarVal{nbvar});
169%     VarAttribute{nbvar+1}.Role='errorflag';
170% end
[515]171
172
173
174
175
Note: See TracBrowser for help on using the repository browser.