source: trunk/src/calc_field_tps.m @ 675

Last change on this file since 675 was 675, checked in by sommeria, 11 years ago

various bugs corrected

File size: 7.1 KB
Line 
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
27function [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbCentre,SubRange,FieldVar,FieldName,Coord_interp)
28
29%list of defined scalars to display in menus (in addition to 'ima_cor').
30% a type is associated to each scalar:
31%              'discrete': related to the individual velocity vectors, not interpolated by patch
32%              'vel': calculated from velocity components, continuous field (interpolated with velocity)
33%              'der': needs spatial derivatives
34%              'var': the scalar name corresponds to a field name in the netcdf files
35% a specific variable name for civ1 and civ2 fields are also associated, if
36% the scalar is calculated from other fields, as explicited below
37errormsg='';
38   
39%% nbre of subdomains
40if ndims(Coord_interp)==3
41    nb_coord=size(Coord_interp,3);
42    npx=size(Coord_interp,2);
43    npy=size(Coord_interp,1);
44    nb_sites=npx*npy;
45    Coord_interp=reshape(Coord_interp,nb_sites,nb_coord);
46else
47    nb_coord=size(Coord_interp,2);
48    nb_sites=size(Coord_interp,1);
49end
50NbSubDomain=size(Coord_tps,3);
51nbval=zeros(nb_sites,1);
52
53%% list of operations
54check_grid=0;
55check_der=0;
56check_vec=0;
57check_remove=false(size(FieldName));
58VarAttribute={};
59for ilist=1:length(FieldName)
60    FieldNameType=regexprep(FieldName{ilist},'(.+','');% detect the char string before the parenthesis
61    VarAttributeNew={};
62    switch FieldNameType
63        case 'vec'
64            check_grid=1;
65            DataOut.U=zeros(nb_sites,1);
66            DataOut.V=zeros(nb_sites,1);
67            VarAttributeNew{1}.Role='vector_x';
68            VarAttributeNew{2}.Role='vector_y';
69            check_vec=1;
70        case {'U','V'}
71            if check_vec% no new data needed
72                check_remove(ilist)=1;
73            else
74            check_grid=1;
75            DataOut.(FieldNameType)=zeros(nb_sites,1);
76            VarAttributeNew{1}.Role='scalar';
77            end
78        case 'norm'
79            check_grid=1;
80            DataOut.(FieldNameType)=zeros(nb_sites,1);
81            VarAttributeNew{1}.Role='scalar';
82        case {'curl','div','strain'}
83            check_der=1;
84            DataOut.(FieldNameType)=zeros(nb_sites,1);
85            VarAttributeNew{1}.Role='scalar';
86    end
87    VarAttribute=[VarAttribute VarAttributeNew];
88end
89Attr_FF.Role='errorflag';
90VarAttribute=[VarAttribute {Attr_FF}];
91FieldName(check_remove)=[];
92
93%% loop on subdomains
94for isub=1:NbSubDomain
95    nbvec_sub=NbCentre(isub);
96    check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
97    ind_sel=find(sum(check_range,2)==nb_coord);
98    nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
99    if check_grid
100        EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
101    end
102    if check_der
103        [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
104    end
105%     ListVar={};
106    for ilist=1:length(FieldName)
107%         var_count=numel(ListVar);
108        switch FieldName{ilist}
109            case 'vec(U,V)'
110%                 ListVar=[ListVar {'U', 'V'}];
111%                 VarAttribute{var_count+1}.Role='vector_x';
112%                 VarAttribute{var_count+2}.Role='vector_y';
113                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
114                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
115            case 'U'
116%                 ListVar=[ListVar {'U'}];
117%                 VarAttribute{var_count+1}.Role='scalar';
118                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
119            case 'V'
120%                 ListVar=[ListVar {'V'}];
121%                 VarAttribute{var_count+1}.Role='scalar';
122                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
123            case 'norm(U,V)'
124%                 ListVar=[ListVar {'norm'}];
125%                 VarAttribute{var_count+1}.Role='scalar';
126                U=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
127                V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
128                DataOut.norm(ind_sel)=sqrt(U.*U+V.*V);
129            case 'curl(U,V)'
130%                 ListVar=[ListVar {'curl'}];
131%                 VarAttribute{var_count+1}.Role='scalar';
132                DataOut.curl(ind_sel)=DataOut.curl(ind_sel)-EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
133            case 'div(U,V)'
134%                 ListVar=[ListVar {'div'}];
135%                 VarAttribute{var_count+1}.Role='scalar';
136                DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);
137            case 'strain(U,V)'
138%                 ListVar=[ListVar {'strain'}];
139%                 VarAttribute{var_count+1}.Role='scalar';
140                DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
141        end
142    end
143end
144DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
145nbval(nbval==0)=1;% to avoid division by zero for averaging
146ListFieldOut=fieldnames(DataOut);
147for ifield=1:numel(ListFieldOut)
148    DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval;
149    DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx);
150end
151
152
153
154
Note: See TracBrowser for help on using the repository browser.