Changeset 12605

Timestamp:

07/03/12 13:44:06 (13 years ago)

Author:

seroussi

Message:

added possibility to average on 1 layer for 3d meshes

File:

• issm/trunk-jpl/src/m/model/averaging.m (modified) (2 diffs)

issm/trunk-jpl/src/m/model/averaging.m

@@ -1 @@
-function average=averaging(md,data,iterations)
+function average=averaging(md,data,iterations,varargin)
 %AVERAGING - smooths the input over the mesh
 %
@@ -8 +8 @@
 %   by taking the average of the element around a node weighted by the
 %   elements volume
+%   For 3d mesh, a last argument can be added to specify the layer to be averaged on.
 %
 %   Usage:
 %      smoothdata=averaging(md,data,iterations)
+%      smoothdata=averaging(md,data,iterations,layer)
 %
 %   Examples:
 %      velsmoothed=averaging(md,md.initialization.vel,4);
 %      pressure=averaging(md,md.initialization.pressure,0);
+%      temperature=averaging(md,md.initialization.temperature,1,1);
 
-if length(data)~=md.mesh.numberofelements & length(data)~=md.mesh.numberofvertices
+if ((nargin~=4) & (nargin~=3)),
+        error('averaging error message');
+end
+if (length(data)~=md.mesh.numberofelements & length(data)~=md.mesh.numberofvertices),
         error('averaging error message: data not supported yet');
+end
+if md.mesh.dimension==3 & nargin==4,
+        if varargin{1}<=0 | varargin{1}>md.mesh.numberoflayers,
+                error('layer should be between 1 and md.mesh.numberoflayers');
+        end
+        layer=varargin{1};
+else
+        layer=0;
 end
 
 %initialization
-weights=zeros(md.mesh.numberofvertices,1);
-data=data(:);
+if layer==0,
+        weights=zeros(md.mesh.numberofvertices,1);
+        data=data(:);
+else 
+        weights=zeros(md.mesh.numberofvertices2d,1);
+        data=data((layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:);
+end
 
-%load some variables (it is much faster if the variab;es are loaded from md once for all)
-index=md.mesh.elements;
-numberofnodes=md.mesh.numberofvertices;
-numberofelements=md.mesh.numberofelements;
+%load some variables (it is much faster if the variabes are loaded from md once for all)
+if layer==0,
+        index=md.mesh.elements;
+        numberofnodes=md.mesh.numberofvertices;
+        numberofelements=md.mesh.numberofelements;
+else
+        index=md.mesh.elements2d;
+        numberofnodes=md.mesh.numberofvertices2d;
+        numberofelements=md.mesh.numberofelements2d;
+end
 
 %build some variables
 line=index(:);
-if md.mesh.dimension==3
+if md.mesh.dimension==3 & layer==0,
         rep=6;
         areas=GetAreas(index,md.mesh.x,md.mesh.y,md.mesh.z);
+elseif md.mesh.dimension==2,
+        rep=3;
+        areas=GetAreas(index,md.mesh.x,md.mesh.y);
 else
         rep=3;
-        areas=GetAreas(index,md.mesh.x,md.mesh.y);
+        areas=GetAreas(index,md.mesh.x2d,md.mesh.y2d);
 end
 summation=1/rep*ones(rep,1);