Changeset 579

Timestamp:

Mar 6, 2013, 9:20:51 PM (12 years ago)

Author:

sommeria

Message:

bugs corrected for reading vort div

Location:

trunk/src

Files:

• calc_field_interp.m (modified) (5 diffs)
• read_civxdata.m (modified) (2 diffs)
• uvmat.m (modified) (1 diff)

trunk/src/calc_field_interp.m

@@ -2 +2 @@
 % for projection with linear interpolation
 %---------------------------------------------------------------------
-% [DataOut,VarAttribute,errormsg]=calc_field_interp(Coord_tps,NbSites,SubRange,FieldVar,Operation,Coord_interp)
+% [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
 %
 % OUTPUT:
 % VarVal: array giving the values of the calculated field
+% ListVarName: corresponding list of variable names
+% VarAttribute: corresponding list of variable attributes, each term #ilist is of the form VarAttribute{ilist}.tag=value
 %
 % INPUT:
-% Coord_tps:
-% NbSites
-% SubRange
-% FieldVar
-% Operation: cell array representing the list of operations (eg div, rot..)
-% Coord_interp: coordiantes of sites on which the fields need to be calculated
+% Coord(nbpoints,2): matrix of x,y coordinates of theinput data points
+% Data: inputfield structure
+% FieldName: string representing the field to calculate, or cell array of fields (as displayed in uvmat/FieldName)
+% XI, YI: set of x and y coordiantes where the fields need to be linearly interpolated
 
-function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,Operation,XI,YI)
+function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
 
-VarVal=[];
+VarVal={};
 ListVarName={};
 VarAttribute={};
 errormsg='';
 InputVarList={};
-if ischar(Operation),Operation={Operation};end
-check_skipped=zeros(size(Operation));
-Operator=cell(size(Operation));
-for ilist=1:numel(Operation)
-    r=regexp(Operation{ilist},'(?<Operator>(^vec|^norm))\((?<UName>.+),(?<VName>.+)\)$','names');
-    if isempty(r) % the operation is the variable
-        if ~isfield(Data,Operation{ilist})
-            check_skipped(ilist)=1;
+if ischar(FieldName),FieldName={FieldName};end
+check_skipped=zeros(size(FieldName));% default, =1 to mark the variables which can be calculated
+Operator=cell(size(FieldName));
+for ilist=1:numel(FieldName)
+    Operator{ilist}='';%default empty operator (vec, norm,...)
+    r=regexp(FieldName{ilist},'(?<Operator>(^vec|^norm|^curl|^div|^strain))\((?<UName>.+),(?<VName>.+)\)$','names');% analyse the field name
+    if isempty(r) % the field name is a variable itself
+        if ~isfield(Data,FieldName{ilist})
+            check_skipped(ilist)=1; %variable not found
         else
-            if isempty(find(strcmp(Operation{ilist},InputVarList)));
-                InputVarList=[InputVarList Operation{ilist}];% the variable is added to the list if it is not already in the list
+            if isempty(find(strcmp(FieldName{ilist},InputVarList), 1));
+                InputVarList=[InputVarList FieldName{ilist}];% the variable is added to the list of input variables if it is not already in the list
             end
-            Operator{ilist}='';
         end
     else
-        if ~isfield(Data,r.UName)||~isfield(Data,r.VName)
+        if ~isfield(Data,r.UName)||~isfield(Data,r.VName)%needed input variable not found
             check_skipped(ilist)=1;
+        elseif strcmp(r.Operator,'curl')||strcmp(r.Operator,'div')||strcmp(r.Operator,'strain')
+            Operator{ilist}=r.Operator;
+            switch r.Operator
+                case 'curl'% case of CivX data format
+                    UName{ilist}='vort';
+                    Data.vort=Data.DjUi(:,1,2)-Data.DjUi(:,2,1);
+                case 'div'
+                    UName{ilist}='div';
+                    Data.div=Data.DjUi(:,1,1)+Data.DjUi(:,2,2);
+                case 'strain'
+                    UName{ilist}='strain';
+                    Data.strain=Data.DjUi(:,1,2)+Data.DjUi(:,2,1);
+            end
+            InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
         else
             UName{ilist}=r.UName;
@@ -45 +59 @@
                 InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
             end
-            if isempty(find(strcmp(r.VName,InputVarList)));
+            if isempty(find(strcmp(r.VName,InputVarList), 1));
                 InputVarList=[InputVarList VName{ilist}]; %the variable is added to the list if it is not already in the list
             end
@@ -58 +72 @@
     end
 end
-for ilist=1:numel(Operation)
+for ilist=1:numel(FieldName)
     if ~check_skipped(ilist)
         nbvar=numel(ListVarName);
@@ -85 +99 @@
                 ListVarName{nbvar+1}='norm';
                 VarAttribute{nbvar+1}.Role='scalar';
+            case {'curl','div','strain'}
+                if exist('XI','var')
+                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
+                else
+                    VarVal{nbvar+1}=Data.(UName{ilist});
+                end
+                ListVarName{nbvar+1}=UName{ilist};
+                VarAttribute{nbvar+1}.Role='scalar';
             otherwise
-                if ~isempty(Operation{ilist})
+                if ~isempty(FieldName{ilist})
                     if exist('XI','var')
-                        VarVal{nbvar+1}=F.(Operation{ilist})(XI,YI);
+                        VarVal{nbvar+1}=F.(FieldName{ilist})(XI,YI);
                     else
-                        VarVal{nbvar+1}= Data.(Operation{ilist});
+                        VarVal{nbvar+1}= Data.(FieldName{ilist});
                     end
-                    ListVarName{nbvar+1}=Operation{ilist};
+                    ListVarName{nbvar+1}=FieldName{ilist};
                     VarAttribute{nbvar+1}.Role='scalar';
                 end
@@ -99 +121 @@
 end
 % put an error flag to indicate NaN data
-if exist('XI','var')
+if exist('XI','var')&&~isempty(VarVal)
     nbvar=numel(ListVarName);
     ListVarName{nbvar+1}='FF';

trunk/src/read_civxdata.m

@@ -93 +93 @@
     Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
     if strcmp(role{var_ind(ivar)},'vector_x')
-        Field.VarAttribute{ivar}.Operation=FieldNames;
+        Field.VarAttribute{ivar}.FieldName=FieldNames;
         if testinterp
             Field.VarAttribute{ivar}.FieldRequest='interp_lin';
@@ -213 +213 @@
 testder=0;
 for ilist=1:length(FieldNames)
-        testder=~isempty(regexp(FieldNames{ilist},'(^curl|^div|strain)', 'once'));%test need for derivatives
+        testder=~isempty(regexp(FieldNames{ilist},'(^curl|^div|^strain)', 'once'));%test need for derivatives
         if testder, break;end
 end      

trunk/src/uvmat.m

@@ -2464 +2464 @@
         if strcmp(CellInfo{imax(ind)}.CoordType,'tps')
             CoordName=UvData.Field.ListVarName{CellInfo{imax(ind)}.CoordIndex};% X,Y coordinates in a single variable
-            CoordMax(ind,NbDim)=max(UvData.Field.(CoordName)(1:end-3,1,:));
-            CoordMax(ind,NbDim-1)=max(UvData.Field.(CoordName)(1:end-3,2,:));
-            CoordMin(ind,NbDim)=min(UvData.Field.(CoordName)(1:end-3,1,:));
-            CoordMin(ind,NbDim-1)=min(UvData.Field.(CoordName)(1:end-3,2,:));
+            CoordMax(ind,NbDim)=max(max(UvData.Field.(CoordName)(1:end-3,1,:),[],1),[],3);% max of x component (2D case)
+            CoordMax(ind,NbDim-1)=max(max(UvData.Field.(CoordName)(1:end-3,2,:),[],1),[],3);% max of y component (2D case)
+            CoordMin(ind,NbDim)=min(min(UvData.Field.(CoordName)(1:end-3,1,:),[],1),[],3);
+            CoordMin(ind,NbDim-1)=min(min(UvData.Field.(CoordName)(1:end-3,2,:),[],1),[],3);% min of y component (2D case)
         else
             XName=UvData.Field.ListVarName{CellInfo{imax(ind)}.CoordIndex(end)};