source: trunk/src/calc_field_tps.m @ 849

Last change on this file since 849 was 809, checked in by g7moreau, 10 years ago
  • Add license
File size: 8.0 KB
RevLine 
[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%=======================================================================
28% Copyright 2008-2014, LEGI UMR 5519 / CNRS UJF 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
[651]45function [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
55errormsg='';
56   
57%% nbre of subdomains
58if 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);
64else
65    nb_coord=size(Coord_interp,2);
66    nb_sites=size(Coord_interp,1);
67end
68NbSubDomain=size(Coord_tps,3);
69nbval=zeros(nb_sites,1);
70
71%% list of operations
72check_grid=0;
73check_der=0;
[675]74check_vec=0;
75check_remove=false(size(FieldName));
76VarAttribute={};
[581]77for 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]106end
107Attr_FF.Role='errorflag';
108VarAttribute=[VarAttribute {Attr_FF}];
[675]109FieldName(check_remove)=[];
[515]110
111%% loop on subdomains
112for 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)');
115    ind_sel=find(sum(check_range,2)==nb_coord);
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
[675]123%     ListVar={};
[581]124    for ilist=1:length(FieldName)
[675]125%         var_count=numel(ListVar);
[581]126        switch FieldName{ilist}
[515]127            case 'vec(U,V)'
[675]128%                 ListVar=[ListVar {'U', 'V'}];
129%                 VarAttribute{var_count+1}.Role='vector_x';
130%                 VarAttribute{var_count+2}.Role='vector_y';
[515]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'
[675]134%                 ListVar=[ListVar {'U'}];
135%                 VarAttribute{var_count+1}.Role='scalar';
[515]136                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
137            case 'V'
[675]138%                 ListVar=[ListVar {'V'}];
139%                 VarAttribute{var_count+1}.Role='scalar';
[515]140                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
141            case 'norm(U,V)'
[675]142%                 ListVar=[ListVar {'norm'}];
143%                 VarAttribute{var_count+1}.Role='scalar';
[515]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)'
[675]148%                 ListVar=[ListVar {'curl'}];
149%                 VarAttribute{var_count+1}.Role='scalar';
[515]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)'
[675]152%                 ListVar=[ListVar {'div'}];
153%                 VarAttribute{var_count+1}.Role='scalar';
[515]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)'
[675]156%                 ListVar=[ListVar {'strain'}];
157%                 VarAttribute{var_count+1}.Role='scalar';
[515]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
161end
162DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
163nbval(nbval==0)=1;% to avoid division by zero for averaging
164ListFieldOut=fieldnames(DataOut);
165for ifield=1:numel(ListFieldOut)
166    DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval;
167    DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx);
168end
169
170
171
172
Note: See TracBrowser for help on using the repository browser.