[515]  1 


[575]  2  %'calc_field_tps': defines fields (velocity, vort, div...) from civ data and calculate them with tps interpolation


[515]  3  %


[651]  4  % [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbCentre,SubRange,FieldVar,FieldName,Coord_interp)


[515]  5  %


 6  % OUTPUT:


 7  % DataOut: structure representing the output fields


[675]  8  % VarAttribute: cell array of structures coontaining the variable attributes


 9  % errormsg: error msg , = '' by default


[515]  10  %


 11  % INPUT:


[651]  12  % Coord_tps: coordinates of the centres, of dimensions [nb_point,nb_coord,nb_subdomain], where


 13  % nb_point is the max number of data point in a subdomain,


 14  % nb_coord the space dimension,


 15  % nb_subdomain the nbre of subdomains used for tps


 16  % NbCentre: nbre of tps centres for each subdomain, of dimension nb_subdomain


 17  % SubRange: coordinate range for each subdomain, of dimensions [nb_coord,2,nb_subdomain]


[675]  18  % FieldVar: array representing the input fields as tps weights with dimension (nbvec_sub+3,NbSubDomain,nb_dim)


 19  % nbvec_sub= max nbre of vectors in a subdomain


 20  % NbSubDomain =nbre of subdomains


 21  % nb_dim: nbre of dimensions for vector components (x> 1, y>2)


[651]  22  % FieldName: cell array representing the list of operations (eg div(U,V), rot(U,V))


 23  % Coord_interp: coordinates of sites on which the fields need to be calculated of dimensions


 24  % [nb_site,nb_coord] for an array of interpolation sites


 25  % [nb_site_y,nb_site_x,nb_coord] for interpolation on a plane grid of size [nb_site_y,nb_site_x]


[515]  26 


[809]  27  %=======================================================================


[1071]  28  % Copyright 20082020, LEGI UMR 5519 / CNRS UGA GINP, Grenoble, France


[809]  29  % http://www.legi.grenobleinp.fr


 30  % Joel.Sommeria  Joel.Sommeria (A) legi.cnrs.fr


 31  %


 32  % This file is part of the toolbox UVMAT.


 33  %


 34  % UVMAT is free software; you can redistribute it and/or modify


 35  % it under the terms of the GNU General Public License as published


 36  % by the Free Software Foundation; either version 2 of the license,


 37  % or (at your option) any later version.


 38  %


 39  % UVMAT is distributed in the hope that it will be useful,


 40  % but WITHOUT ANY WARRANTY; without even the implied warranty of


 41  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


 42  % GNU General Public License (see LICENSE.txt) for more details.


 43  %=======================================================================


 44 


[651]  45  function [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbCentre,SubRange,FieldVar,FieldName,Coord_interp)


[515]  46 


 47  %list of defined scalars to display in menus (in addition to 'ima_cor').


 48  % a type is associated to each scalar:


 49  % 'discrete': related to the individual velocity vectors, not interpolated by patch


 50  % 'vel': calculated from velocity components, continuous field (interpolated with velocity)


 51  % 'der': needs spatial derivatives


 52  % 'var': the scalar name corresponds to a field name in the netcdf files


 53  % a specific variable name for civ1 and civ2 fields are also associated, if


 54  % the scalar is calculated from other fields, as explicited below


 55  errormsg='';


 56 


 57  %% nbre of subdomains


 58  if ndims(Coord_interp)==3


 59  nb_coord=size(Coord_interp,3);


 60  npx=size(Coord_interp,2);


 61  npy=size(Coord_interp,1);


 62  nb_sites=npx*npy;


 63  Coord_interp=reshape(Coord_interp,nb_sites,nb_coord);


 64  else


 65  nb_coord=size(Coord_interp,2);


 66  nb_sites=size(Coord_interp,1);


 67  end


 68  NbSubDomain=size(Coord_tps,3);


 69  nbval=zeros(nb_sites,1);


 70 


 71  %% list of operations


 72  check_grid=0;


 73  check_der=0;


[675]  74  check_vec=0;


 75  check_remove=false(size(FieldName));


 76  VarAttribute={};


[581]  77  for ilist=1:length(FieldName)


[675]  78  FieldNameType=regexprep(FieldName{ilist},'(.+','');% detect the char string before the parenthesis


 79  VarAttributeNew={};


[581]  80  switch FieldNameType


[515]  81  case 'vec'


[521]  82  check_grid=1;


[515]  83  DataOut.U=zeros(nb_sites,1);


 84  DataOut.V=zeros(nb_sites,1);


[675]  85  VarAttributeNew{1}.Role='vector_x';


 86  VarAttributeNew{2}.Role='vector_y';


 87  check_vec=1;


 88  case {'U','V'}


 89  if check_vec% no new data needed


 90  check_remove(ilist)=1;


 91  else


[515]  92  check_grid=1;


[581]  93  DataOut.(FieldNameType)=zeros(nb_sites,1);


[675]  94  VarAttributeNew{1}.Role='scalar';


 95  end


 96  case 'norm'


 97  check_grid=1;


 98  DataOut.(FieldNameType)=zeros(nb_sites,1);


 99  VarAttributeNew{1}.Role='scalar';


[1090]  100  case {'curl','div','strain','DUDX','DUDY','DVDX','DVDY'}


[515]  101  check_der=1;


[581]  102  DataOut.(FieldNameType)=zeros(nb_sites,1);


[675]  103  VarAttributeNew{1}.Role='scalar';


[515]  104  end


[675]  105  VarAttribute=[VarAttribute VarAttributeNew];


[515]  106  end


[675]  107  FieldName(check_remove)=[];


[515]  108 


 109  %% loop on subdomains


 110  for isub=1:NbSubDomain


[651]  111  nbvec_sub=NbCentre(isub);


[515]  112  check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');


[905]  113  ind_sel=find(sum(check_range,2)==nb_coord);% select points whose all coordinates are in the prescribed range


[515]  114  nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)


 115  if check_grid


 116  EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'


 117  end


 118  if check_der


 119  [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'


 120  end


[581]  121  for ilist=1:length(FieldName)


 122  switch FieldName{ilist}


[515]  123  case 'vec(U,V)'


 124  DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);


 125  DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);


 126  case 'U'


 127  DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);


 128  case 'V'


 129  DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);


 130  case 'norm(U,V)'


 131  U=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);


 132  V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);


 133  DataOut.norm(ind_sel)=sqrt(U.*U+V.*V);


 134  case 'curl(U,V)'


 135  DataOut.curl(ind_sel)=DataOut.curl(ind_sel)EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);


 136  case 'div(U,V)'


 137  DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);


 138  case 'strain(U,V)'


 139  DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);


[1090]  140  case 'DUDX(U,V)'


 141  DataOut.DUDX(ind_sel)=DataOut.DUDX(ind_sel)+EMDX *FieldVar(1:nbvec_sub+3,isub,1);


 142  case 'DUDY(U,V)'


 143  DataOut.DUDY(ind_sel)=DataOut.DUDY(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1);


 144  case 'DVDX(U,V)'


 145  DataOut.DVDX(ind_sel)=DataOut.DVDX(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,2);


 146  case 'DVDY(U,V)'


 147  DataOut.DVDY(ind_sel)=DataOut.DVDY(ind_sel)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);


[515]  148  end


 149  end


 150  end


[867]  151  nbval(nbval==0)=NaN;% to avoid division by zero for averaging


[515]  152  ListFieldOut=fieldnames(DataOut);


 153  for ifield=1:numel(ListFieldOut)


 154  DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval;


 155  DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx);


 156  end


 157 


 158 


 159 


 160 

