[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 | %=======================================================================
|
---|
[1027] | 28 | % Copyright 2008-2018, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
|
---|
[809] | 29 | % http://www.legi.grenoble-inp.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';
|
---|
[515] | 100 | case {'curl','div','strain'}
|
---|
| 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
|
---|
[867] | 107 | %Attr_FF.Role='errorflag';
|
---|
| 108 | %VarAttribute=[VarAttribute {Attr_FF}];
|
---|
[675] | 109 | FieldName(check_remove)=[];
|
---|
[515] | 110 |
|
---|
| 111 | %% loop on subdomains
|
---|
| 112 | for isub=1:NbSubDomain
|
---|
[651] | 113 | nbvec_sub=NbCentre(isub);
|
---|
[515] | 114 | check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
|
---|
[905] | 115 | ind_sel=find(sum(check_range,2)==nb_coord);% select points whose all coordinates are in the prescribed range
|
---|
[515] | 116 | nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
|
---|
| 117 | if check_grid
|
---|
| 118 | EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
|
---|
| 119 | end
|
---|
| 120 | if check_der
|
---|
| 121 | [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
|
---|
| 122 | end
|
---|
[581] | 123 | for ilist=1:length(FieldName)
|
---|
| 124 | switch FieldName{ilist}
|
---|
[515] | 125 | case 'vec(U,V)'
|
---|
| 126 | DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
|
---|
| 127 | DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
| 128 | case 'U'
|
---|
| 129 | DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
|
---|
| 130 | case 'V'
|
---|
| 131 | DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
| 132 | case 'norm(U,V)'
|
---|
| 133 | U=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
|
---|
| 134 | V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
| 135 | DataOut.norm(ind_sel)=sqrt(U.*U+V.*V);
|
---|
| 136 | case 'curl(U,V)'
|
---|
| 137 | DataOut.curl(ind_sel)=DataOut.curl(ind_sel)-EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
| 138 | case 'div(U,V)'
|
---|
| 139 | DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
[938] | 140 | % case 'curl_polar(U,V)'
|
---|
| 141 | % DataOut.curl(ind_sel)=DataOut.curl(ind_sel)+(-EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *(Coord_interp(ind_sel,1).*FieldVar(1:nbvec_sub+3,isub,2)))./Coord_interp(ind_sel,1);
|
---|
[515] | 142 | case 'strain(U,V)'
|
---|
| 143 | DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
| 144 | end
|
---|
| 145 | end
|
---|
| 146 | end
|
---|
[867] | 147 | nbval(nbval==0)=NaN;% to avoid division by zero for averaging
|
---|
[515] | 148 | ListFieldOut=fieldnames(DataOut);
|
---|
| 149 | for ifield=1:numel(ListFieldOut)
|
---|
| 150 | DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval;
|
---|
| 151 | DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx);
|
---|
| 152 | end
|
---|
| 153 |
|
---|
| 154 |
|
---|
| 155 |
|
---|
| 156 |
|
---|