[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 20082017, LEGI UMR 5519 / CNRS UGA GINP, Grenoble, France


[809]  19  % http://www.legi.grenobleinp.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]  35  function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)


[515]  36 


[644]  37  %% initialization


[579]  38  VarVal={};


[515]  39  ListVarName={};


 40  VarAttribute={};


 41  errormsg='';


[517]  42  InputVarList={};


[579]  43  if ischar(FieldName),FieldName={FieldName};end


 44  check_skipped=zeros(size(FieldName));% default, =1 to mark the variables which can be calculated


[654]  45  check_interp=ones(size(FieldName));% default, =1 to mark the variables which can be interpolated (not ancillary)


[579]  46  Operator=cell(size(FieldName));


[644]  47 


 48  %% analyse the list of input fields: needed variables and requested operations


[579]  49  for 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


 97  end


[644]  98 


 99  %% create interpolator for each variable to interpolate


[517]  100  if exist('XI','var')


 101  for ilist=1:numel(InputVarList)


 102  F.(InputVarList{ilist})=TriScatteredInterp(Coord,Data.(InputVarList{ilist}),'linear');


 103  end


[515]  104  end


[644]  105 


 106  %% perform the linear interpolation for the requested variables


[579]  107  for 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


 162  end


[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 

