Changeset 24855

Timestamp:

05/14/20 23:11:18 (5 years ago)

Author:

Eric.Larour

Message:

CHG: considerable speed up + 2D connectivity method.

Location:

issm/trunk-jpl/src/m/contrib/larour

Files:

• glacier_inventory.m (modified) (6 diffs)
• isintrianglearray.m (added)
• isintrianglearraytotal.m (added)

issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m

@@ -84 +84 @@
 
                 end % }}}
-                function mesh_connectivity(self,mesh) % {{{
-
-                        plotr=1;
+                function mesh_connectivity(self,mesh,varargin) % {{{
+                        
+                        %retrieve options: 
+                        options=pairoptions(varargin{:}); 
+
+                        subsetregions=getfieldvalue(options,'regions',1:length(self.boxes));
+                        errornotfound=getfieldvalue(options,'errornotfound',1);
 
                         %The way we run this is through the O2 zones defined in boxes. We go through 
@@ -119 +123 @@
 
                         %Go through O2 regions: 
-                        for i=1:length(self.boxes),
+                        for i=subsetregions,
                                 string=self.boxes(i).RGI_CODE; 
                                 disp(['progressing with region ' num2str(i) ' ' string]);
@@ -138 +142 @@
                                         radius=sqrt( (lat0-minlat)^2+(long0-maxlong)^2);
 
-                %some radius adjustment:  {{{
-                switch i,
-                        case 4, radius=40;
-                        case 12, radius=25;
-                        case 33, radius=5;
-                        case 41, radius=75;
-                        case 42, radius=45;
-                        case 68, radius=10;
-                        case 82, radius=30;
-                        otherwise,
-                end % }}}
+                                        %some radius adjustment:  {{{
+                                        switch i,
+                                                case 4, radius=40;
+                                                case 12, radius=25;
+                                                case 33, radius=5;
+                                                case 41, radius=75;
+                                                case 42, radius=45;
+                                                case 68, radius=10;
+                                                case 82, radius=30;
+                                                otherwise,
+                                        end % }}}
 
                                         [lids,rids]=self.gidtolid(glaciers);
@@ -197 +201 @@
                                                         end
                                                 end
-                                                self.glacier_connectivity(glaciers(j))=el;
-                                                if ~found, 
-                                                        error(sprintf('could not find element for glacier %i with lid %i',j,lids(j)));
+                                                if ~found,
+                                                        if errornotfound,
+                                                                error(sprintf('could not find element for glacier %i with lid %i',j,lids(j)));
+                                                        end
                                                 end
+                                                if(found)self.glacier_connectivity(glaciers(j))=el; end;
                                         end
                                 end
@@ -208 +214 @@
                         for j=1:length(self.glacier_connectivity),
                                 el=self.glacier_connectivity(j);
+                                if ~el,continue; end;
+                                count=self.element_connectivity(el,ny);
+                                if count>ny,
+                                        error('need to enlarge connectivity table to at least');
+                                end
+                                self.element_connectivity(el,count+1)=j;
+                                self.element_connectivity(el,ny)=count+1;
+                        end
+
+                        %Reduce the number of columns (we did not initially, so we chose an arbitrary ny:
+                        nmax=max(self.element_connectivity(:,end));
+                        self.element_connectivity=self.element_connectivity(:,[1:nmax,ny]);
+
+
+                end % }}}
+                function mesh_connectivity2d(self,md,proj,varargin) % {{{
+                        
+                        %retrieve options: 
+                        options=pairoptions(varargin{:}); 
+
+                        subsetregions=getfieldvalue(options,'regions',1:self.nregions());
+
+                        %initialize glacier connectivity: 
+                        self.glacier_connectivity=zeros(self.nglaciers,1);
+
+                        %assume maximum width for connectivity table and initialize 
+                        %as sparse matrix: 
+                        ny=self.nglaciers(); self.element_connectivity=sparse(md.mesh.numberofelements,ny);
+                                
+                        disp('Building element and glacier connectivity table: ');
+                        [lat0,long0]=projlatlong(proj);
+                
+                        [mpartition,npartition]=self.partition();
+                        for i=subsetregions,
+                                region=self.regions(i);
+                                disp(sprintf(' progress for region: %s',region.name));
+
+                                %project lat,long: 
+                                [xlid,ylid]=gdaltransform(region.CenLon,region.CenLat,'EPSG:4326',proj);
+
+                                %go through lids: 
+                                x0=xlid; y0=ylid;
+                                x1=md.mesh.x(md.mesh.elements(:,1)); y1=md.mesh.y(md.mesh.elements(:,1));
+                                x2=md.mesh.x(md.mesh.elements(:,2)); y2=md.mesh.y(md.mesh.elements(:,2));
+                                x3=md.mesh.x(md.mesh.elements(:,3)); y3=md.mesh.y(md.mesh.elements(:,3));
+                                in=isintrianglearraytotal(x0,x1,x2,x3,y0,y1,y2,y3);
+                                [els,glacs]=find(in);
+                                self.glacier_connectivity(mpartition(i)-1+glacs)=els;
+                        end
+
+                        %build element connectivity table: 
+                        for j=1:length(self.glacier_connectivity),
+                                el=self.glacier_connectivity(j);
+                                if ~el,continue; end;
                                 count=self.element_connectivity(el,ny);
                                 if count>ny,
@@ -332 +392 @@
                 end
                 %}}}
+                function listboxes(self) % {{{
+
+                        for i=1:length(self.boxes),
+                                b=self.boxes(i);
+                                disp(sprintf('Region #%i: %s (%s)',i,b.FULL_NAME,b.RGI_CODE));
+                        end
+                        
+                end
+                %}}}
         end
 end