1 |
|
---|
2 | %'calc_field_tps': defines fields (velocity, vort, div...) from civ data and calculate them with tps interpolation
|
---|
3 | %---------------------------------------------------------------------
|
---|
4 | % [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbCentre,SubRange,FieldVar,FieldName,Coord_interp)
|
---|
5 | %
|
---|
6 | % OUTPUT:
|
---|
7 | % DataOut: structure representing the output fields
|
---|
8 | % VarAttribute: cell array of structures coontaining the variable attributes
|
---|
9 | % errormsg: error msg , = '' by default
|
---|
10 | %
|
---|
11 | % INPUT:
|
---|
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]
|
---|
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)
|
---|
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]
|
---|
26 |
|
---|
27 | %=======================================================================
|
---|
28 | % Copyright 2008-2016, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
|
---|
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 |
|
---|
45 | function [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbCentre,SubRange,FieldVar,FieldName,Coord_interp)
|
---|
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;
|
---|
74 | check_vec=0;
|
---|
75 | check_remove=false(size(FieldName));
|
---|
76 | VarAttribute={};
|
---|
77 | for ilist=1:length(FieldName)
|
---|
78 | FieldNameType=regexprep(FieldName{ilist},'(.+','');% detect the char string before the parenthesis
|
---|
79 | VarAttributeNew={};
|
---|
80 | switch FieldNameType
|
---|
81 | case 'vec'
|
---|
82 | check_grid=1;
|
---|
83 | DataOut.U=zeros(nb_sites,1);
|
---|
84 | DataOut.V=zeros(nb_sites,1);
|
---|
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
|
---|
92 | check_grid=1;
|
---|
93 | DataOut.(FieldNameType)=zeros(nb_sites,1);
|
---|
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';
|
---|
100 | case {'curl','div','strain'}
|
---|
101 | check_der=1;
|
---|
102 | DataOut.(FieldNameType)=zeros(nb_sites,1);
|
---|
103 | VarAttributeNew{1}.Role='scalar';
|
---|
104 | end
|
---|
105 | VarAttribute=[VarAttribute VarAttributeNew];
|
---|
106 | end
|
---|
107 | %Attr_FF.Role='errorflag';
|
---|
108 | %VarAttribute=[VarAttribute {Attr_FF}];
|
---|
109 | FieldName(check_remove)=[];
|
---|
110 |
|
---|
111 | %% loop on subdomains
|
---|
112 | for isub=1:NbSubDomain
|
---|
113 | nbvec_sub=NbCentre(isub);
|
---|
114 | check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
|
---|
115 | ind_sel=find(sum(check_range,2)==nb_coord);% select points whose all coordinates are in the prescribed range
|
---|
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
|
---|
123 | for ilist=1:length(FieldName)
|
---|
124 | %Operator{ilist}='';%default empty operator (vec, norm,...)
|
---|
125 | %r=regexp(FieldName{ilist},'(?<Operator>(^vec|^norm|^curl|^div|^strain))\((?<UName>.+),(?<VName>.+)\)$','names');% TODO, replace U, V
|
---|
126 | switch FieldName{ilist}
|
---|
127 | case 'vec(U,V)'
|
---|
128 | % ListVar=[ListVar {'U', 'V'}];
|
---|
129 | % VarAttribute{var_count+1}.Role='vector_x';
|
---|
130 | % VarAttribute{var_count+2}.Role='vector_y';
|
---|
131 | DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
|
---|
132 | DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
133 | case 'U'
|
---|
134 | % ListVar=[ListVar {'U'}];
|
---|
135 | % VarAttribute{var_count+1}.Role='scalar';
|
---|
136 | DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
|
---|
137 | case 'V'
|
---|
138 | % ListVar=[ListVar {'V'}];
|
---|
139 | % VarAttribute{var_count+1}.Role='scalar';
|
---|
140 | DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
141 | case 'norm(U,V)'
|
---|
142 | % ListVar=[ListVar {'norm'}];
|
---|
143 | % VarAttribute{var_count+1}.Role='scalar';
|
---|
144 | U=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
|
---|
145 | V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
146 | DataOut.norm(ind_sel)=sqrt(U.*U+V.*V);
|
---|
147 | case 'curl(U,V)'
|
---|
148 | % ListVar=[ListVar {'curl'}];
|
---|
149 | % VarAttribute{var_count+1}.Role='scalar';
|
---|
150 | DataOut.curl(ind_sel)=DataOut.curl(ind_sel)-EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
151 | case 'div(U,V)'
|
---|
152 | % ListVar=[ListVar {'div'}];
|
---|
153 | % VarAttribute{var_count+1}.Role='scalar';
|
---|
154 | DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
155 | case 'strain(U,V)'
|
---|
156 | % ListVar=[ListVar {'strain'}];
|
---|
157 | % VarAttribute{var_count+1}.Role='scalar';
|
---|
158 | DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
|
---|
159 | end
|
---|
160 | end
|
---|
161 | end
|
---|
162 | %DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
|
---|
163 | nbval(nbval==0)=NaN;% to avoid division by zero for averaging
|
---|
164 | ListFieldOut=fieldnames(DataOut);
|
---|
165 | for ifield=1:numel(ListFieldOut)
|
---|
166 | DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval;
|
---|
167 | DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx);
|
---|
168 | end
|
---|
169 |
|
---|
170 |
|
---|
171 |
|
---|
172 |
|
---|