Changeset 13008

Timestamp:

08/13/12 10:52:25 (13 years ago)

Author:

Mathieu Morlighem

Message:

CHG: moved stuff around

Location:

issm/trunk-jpl/src/m

Files:

• classes/modellist.m (modified) (1 diff)
• contrib/ecco/MeltingGroundingLines.m (moved) (moved from issm/trunk-jpl/src/m/model/MeltingGroundingLines.m )
• contrib/ecco/PropagateFlagsUntilDistance.m (moved) (moved from issm/trunk-jpl/src/m/model/PropagateFlagsUntilDistance.m ) (3 diffs)
• contrib/hack (added)
• contrib/hack/thicknessevolution.m (moved) (moved from issm/trunk-jpl/src/m/model/thicknessevolution.m )
• contrib/hack/tres.m (moved) (moved from issm/trunk-jpl/src/m/model/tres.m )
• contrib/massbalance/outflow.m (moved) (moved from issm/trunk-jpl/src/m/model/outflow.m )
• geometry/ThicknessCorrection.m (moved) (moved from issm/trunk-jpl/src/m/model/ThicknessCorrection.m )
• geometry/recover_areas.m (moved) (moved from issm/trunk-jpl/src/m/model/recover_areas.m )
• geometry/recover_areas.py (moved) (moved from issm/trunk-jpl/src/m/model/recover_areas.py )
• geometry/slope.m (moved) (moved from issm/trunk-jpl/src/m/model/slope.m )
• interp/ProfileValues.m (moved) (moved from issm/trunk-jpl/src/m/model/ProfileValues.m )
• interp/SectionValues.m (moved) (moved from issm/trunk-jpl/src/m/model/SectionValues.m )
• meca/shear2d.m (moved) (moved from issm/trunk-jpl/src/m/model/shear2d.m )
• model/README (deleted)
• model/inversions/misfit.m (moved) (moved from issm/trunk-jpl/src/m/model/misfit.m )
• model/mesh/bamg.m (modified) (1 diff)
• model/modelsextract.m (deleted)
• model/modelsextractfromdomains.m (deleted)
• model/modis.m (deleted)
• model/plot/radarpower.m (moved) (moved from issm/trunk-jpl/src/m/model/radarpower.m )
• model/processgeometry.m (deleted)
• model/qstat.m (deleted)
• model/regional/regionaltransient2d.m (moved) (moved from issm/trunk-jpl/src/m/model/regionaltransient2d.m )
• model/sia.m (deleted)

issm/trunk-jpl/src/m/classes/modellist.m

@@ -10 +10 @@
         end
         methods
+                function md_list=modelsextract(md,flags,minel,varargin) % {{{
+                        %modelsextract - extract several self contained models according to a list of element flags.
+                        %
+                        %   The difference between this routine and the modelextract.m routine (without an 's') is that 
+                        %   as many models are extracted as there are closed contours defined in area. 
+                        %   This routine is needed for example when doing data assimilation of ice shelves in Antarctica. 
+                        %   Many independent ice shelves are present, and we don't want data assimilation on one ice shelf 
+                        %   to be hindered by another totally independent ice shelf.
+                        %
+                        %   Usage:
+                        %      md_list=modelsextract(md,elementfalgs,minel);
+                        %
+                        %   Examples:
+                        %      md_list=modelsextract(md,md.mask.elementonfloatingice,1000);
+                        %
+                        %   See also: EXTRUDE, COLLAPSE, MODELEXTRACT
+
+                        disp('selecting pools of elements');
+                        %go through flags and build as many independent element flags as there are groups of connected 1s
+                        %in flags.
+
+                        %2D or 3D?
+                        if md.mesh.dimension==3,
+                                numberofelements=md.mesh.numberofelements2d; %this will be forgotten when we get out.
+                                flags=project2d(md,flags,1);
+                        else
+                                numberofelements=md.mesh.numberofelements;
+                        end
+
+                        %recover extra arguments: 
+                        distance=0;
+                        if nargin==4,
+                                distance=varargin{1};
+                        end
+
+                        flag_list=cell(0,1);
+
+                        for i=1:size(flags,1),
+
+                                if (flags(i)),
+
+                                        %ok, we are sure element i is part of a new pool.
+                                        pool=zeros(numberofelements,1);
+                                        pool=PropagateFlagsFromConnectivity(md.mesh.elementconnectivity,pool,i,flags);
+                                        flag_list{end+1,1}=pool;
+
+                                        %speed up rest of computation by taking pool out of flags: 
+                                        pos=find(pool);flags(pos)=0;
+
+                                end
+                        end
+
+                        %go through flag_list and discard any pool of less than minel elements: 
+                        ex_pos=[];
+                        for i=1:length(flag_list),
+                                if length(find(flag_list{i}))<minel,
+                                        ex_pos=[ex_pos; i];
+                                end
+                        end
+                        flag_list(ex_pos)=[];
+
+                        %now, if distance was specified, expand the flag_list by distance km: 
+                        if distance,
+                                for i=1:length(flag_list),
+                                        flag_list{i}=PropagateFlagsUntilDistance(md,flag_list{i},distance);
+                                end
+                        end
+
+                        %now, go use the pools of flags to extract models: 
+                        disp(['extracting ' num2str(size(flag_list,1)) ' models']);
+                        models=cell(0,1);
+
+                        for i=1:size(flag_list,1),
+                                disp(['   ' num2str(i) '/' num2str(size(flag_list,1))]);
+                                if md.mesh.dimension==3,
+                                        flags2d=flag_list{i};
+                                        realflags=project3d(md,flags2d,'element');
+                                else
+                                        realflags=flag_list{i};
+                                end
+                                models{end+1,1}=modelextract(md,realflags);
+                        end
+
+                        %return model list
+                        md_list=modellist(models);
+
+                end %end of this function }}}
+                function md_list=modelsextractfromdomains(md,directory) % {{{
+                        %modelsextractfromdomains- extract several self contained models according to a list of domains
+                        %
+                        %   Usage:
+                        %      md_list=modelsextractfromdomains(md,'Basins/');
+                        %
+                        %   Examples:
+                        %      md_list=modelsextract(md,'Basins/');
+                        %
+                        %   See also: MODELSEXTRACTS, MODELEXTRACT
+
+                        %go into directory and get list of files.
+                        cd(directory);
+                        basins=listfiles;
+                        cd ..
+
+                        models=cell(0,1);
+                        for i=1:length(basins),
+                                models{end+1,1}=modelextract(md,[directory '/' basins{i}]);
+                        end
+
+                        %return model list: 
+                        md_list=modellist(models);
+
+                end % }}}
                 function obj = modellist(varargin) % {{{
 

issm/trunk-jpl/src/m/contrib/ecco/PropagateFlagsUntilDistance.m

@@ -7 +7 @@
 %
 
-
-        
 new_flags=flags;
 
@@ -24 +22 @@
 sum_conn=sum(conn,2);
 border_elements=flag_elements(find(sum_conn>=1));
-
 
 %average x and y over elements: 
@@ -49 +46 @@
         end
 
-
         %check which of these new elements are more than distance away from the border elements
         for i=1:length(new_elements),

issm/trunk-jpl/src/m/model/mesh/bamg.m

@@ -339 +339 @@
         error('Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)');
 end
+end 
+
+function geom=processgeometry(geom,tol,outline); % {{{
+
+%Deal with edges
+disp('Checking Edge crossing...');
+i=0;
+while (i<size(geom.Edges,1)),
+
+        %edge counter
+        i=i+1;
+
+        %Get coordinates
+        x1=geom.Vertices(geom.Edges(i,1),1);
+        y1=geom.Vertices(geom.Edges(i,1),2);
+        x2=geom.Vertices(geom.Edges(i,2),1);
+        y2=geom.Vertices(geom.Edges(i,2),2);
+        color1=geom.Edges(i,3);
+
+        j=i; %test edges located AFTER i only
+        while (j<size(geom.Edges,1)),
+
+                %edge counter
+                j=j+1;
+
+                %Skip if the two edges already have a vertex in common
+                if any(ismember(geom.Edges(i,1:2),geom.Edges(j,1:2))),
+                        continue
+                end
+
+                %Get coordinates
+                x3=geom.Vertices(geom.Edges(j,1),1);
+                y3=geom.Vertices(geom.Edges(j,1),2);
+                x4=geom.Vertices(geom.Edges(j,2),1);
+                y4=geom.Vertices(geom.Edges(j,2),2);
+                color2=geom.Edges(j,3);
+
+                %Check if the two edges are crossing one another
+                if SegIntersect([x1 y1; x2 y2],[x3 y3; x4 y4]),
+
+                        %Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
+                        x=det([det([x1 y1; x2 y2])  x1-x2;det([x3 y3; x4 y4])  x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+                        y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+
+                        %Add vertex to the list of vertices
+                        geom.Vertices(end+1,:)=[x y min(color1,color2)];
+                        id=size(geom.Vertices,1);
+
+                        %Update edges i and j
+                        edgei=geom.Edges(i,:);
+                        edgej=geom.Edges(j,:);
+                        geom.Edges(i,:)    =[edgei(1) id       edgei(3)];
+                        geom.Edges(end+1,:)=[id       edgei(2) edgei(3)];
+                        geom.Edges(j,:)    =[edgej(1) id       edgej(3)];
+                        geom.Edges(end+1,:)=[id       edgej(2) edgej(3)];
+
+                        %update current edge second tip
+                        x2=x; y2=y;
+                end
+        end
+
+end
+
+%Check point outside
+disp('Checking for points outside the domain...');
+i=0;
+num=0;
+while (i<size(geom.Vertices,1)),
+
+        %vertex counter
+        i=i+1;
+
+        %Get coordinates
+        x=geom.Vertices(i,1);
+        y=geom.Vertices(i,2);
+        color=geom.Vertices(i,3);
+
+        %Check that the point is inside the domain
+        if (color~=1 & ~ContourToNodes(x,y,outline(1),1)),
+
+                %Remove points from list of Vertices
+                num=num+1;
+                geom.Vertices(i,:)=[];
+
+                %update edges
+                [posedges dummy]=find(geom.Edges==i);
+                geom.Edges(posedges,:)=[];
+                posedges=find(geom.Edges>i);
+                geom.Edges(posedges)=geom.Edges(posedges)-1;
+
+                %update counter
+                i=i-1;
+        end
+end
+if num,
+        disp(['WARNING: ' num2str(num) ' points outside the domain outline have been removed']);
+end
+
+%Check point spacing
+if ~isnan(tol),
+        disp('Checking point spacing...');
+        i=0;
+        while (i<size(geom.Vertices,1)),
+
+                %vertex counter
+                i=i+1;
+
+                %Get coordinates
+                x1=geom.Vertices(i,1);
+                y1=geom.Vertices(i,2);
+
+                j=i; %test edges located AFTER i only
+                while (j<size(geom.Vertices,1)),
+
+                        %vertex counter
+                        j=j+1;
+
+                        %Get coordinates
+                        x2=geom.Vertices(j,1);
+                        y2=geom.Vertices(j,2);
+
+                        %Check whether the two vertices are too close
+                        if ((x2-x1)^2+(y2-y1)^2<tol^2)
+
+                                %Remove points from list of Vertices
+                                geom.Vertices(j,:)=[];
+
+                                %update edges
+                                posedges=find(ismember(geom.Edges,j));
+                                geom.Edges(posedges)=i;
+                                posedges=find(geom.Edges>j);
+                                geom.Edges(posedges)=geom.Edges(posedges)-1;
+
+                                %update counter
+                                j=j-1;
+
+                        end
+                end
+        end
+end
+%remove empty edges
+geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[];
+end % }}}