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-2020, 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','DUDX','DUDY','DVDX','DVDY'}
|
---|
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 | FieldName(check_remove)=[];
|
---|
108 |
|
---|
109 | %% loop on subdomains
|
---|
110 | for isub=1:NbSubDomain
|
---|
111 | nbvec_sub=NbCentre(isub);
|
---|
112 | check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
|
---|
113 | ind_sel=find(sum(check_range,2)==nb_coord);% select points whose all coordinates are in the prescribed range
|
---|
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
|
---|
121 | for ilist=1:length(FieldName)
|
---|
122 | switch FieldName{ilist}
|
---|
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);
|
---|
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);
|
---|
148 | end
|
---|
149 | end
|
---|
150 | end
|
---|
151 | nbval(nbval==0)=NaN;% to avoid division by zero for averaging
|
---|
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 |
|
---|