Changeset 3260

Timestamp:

03/11/10 13:33:00 (15 years ago)

Author:

Eric.Larour

Message:

read modis image

Location:

issm/trunk/src/m/classes/public

Files:

• mesh/meshyams.m (modified) (1 diff)
• mesh/rifts/meshprocessrifts.m (modified) (2 diffs)
• mesh/rifts/meshyamsrecreateriftsegments.m (modified) (1 diff)
• modis.m (added)
• plot/plot_rifts.m (modified) (2 diffs)

issm/trunk/src/m/classes/public/mesh/meshyams.m

@@ -55 +55 @@
 else
         md=mesh(md,domainoutline,riftoutline,resolution);
-        md=meshprocessrifts(md);
+        md=meshprocessrifts(md,domainoutline);
 end
 disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);

issm/trunk/src/m/classes/public/mesh/rifts/meshprocessrifts.m

@@ -1 @@
-function md=meshprocessrifts(md)
+function md=meshprocessrifts(md,domainoutline)
 %MESHPROCESSRIFTS - process mesh when rifts are present
 %
 %   split rifts inside mesh (rifts are defined by presence of
 %   segments inside the domain outline)
+%   if domain outline is provided, check for rifts that could touch it, and open them up.
 %
 %   Usage:
 %      md=meshprocessrifts(md)
+%
+%   Ex: 
+%      md=meshprocessrifts(md);
+%      md=meshprocessrifts(md,'DomainOutline.exp');
+%
 
-%some checks on list of arguments
-if ((nargin~=1) | (nargout~=1)),
-        help meshprocessrifts
-        error('meshprocessrifts error messagei: bad usage');
+%some checks on arguments: 
+if nargout~=1,
+        error('meshprocessrifts usage error:');
 end
+
+if nargin~=2,
+        error('meshprocessrifts usage error:');
+end
+
 
 [md.elements,md.x,md.y,md.segments,md.segmentmarkers,md.rifts]=TriMeshProcessRifts(md.elements,md.x,md.y,md.segments,md.segmentmarkers);
@@ -33 +43 @@
         md.rifts(i).tip2coordinates=[md.x(md.rifts(i).tips(2)) md.y(md.rifts(i).tips(2))];
 end
+
+%In case we have rifts that open up the domain outline, we need to open them: 
+flags=ContourToMesh(md.elements,md.x,md.y,expread(domainoutline,1),'node',0);
+found=0;
+for i=1:md.numrifts,
+        if flags(md.rifts(i).tips(1))==0,
+                found=1;
+                break;
+        end
+        if flags(md.rifts(i).tips(2))==0,
+                found=1;
+                break;
+        end
+end
+if found,
+        md=meshprocessoutsiderifts(md,domainoutline);
+end
+
+%get elements that are not correctly oriented in the correct direction:
+aires=GetAreas(md.elements,md.x,md.y);
+pos=find(aires<0);
+md.elements(pos,:)=[md.elements(pos,2) md.elements(pos,1) md.elements(pos,3)];

issm/trunk/src/m/classes/public/mesh/rifts/meshyamsrecreateriftsegments.m

@@ -7 +7 @@
                         rift=md.rifts(i);
 
-                        %find tip1 and tip2 for this rift, in the new mesh created by yams.
-                        pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
-                        tip1=md.segments(pos,1);
-                        pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
-                        tip2=md.segments(pos,1);
+                        %closed rifts first:
+                        if length(rift.tips)==2,
 
-                        %start from tip1, and build segments of this rift. 
-                        pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
-                        pos_record=[pos_record; pos];
-                        riftsegs=md.segments(pos,:);
-                        while 1,
-                                A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
-                                %find other segment that holds B.
-                                pos=find(md.segments(:,1)==B);
+                                %find tip1 and tip2 for this rift, in the new mesh created by yams.
+                                pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
+                                tip1=md.segments(pos,1);
+                                pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
+                                tip2=md.segments(pos,1);
+
+                                %start from tip1, and build segments of this rift. 
+                                pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
+                                pos_record=[pos_record; pos];
+                                riftsegs=md.segments(pos,:);
+                                while 1,
+                                        A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
+                                        %find other segment that holds B.
+                                        pos=find(md.segments(:,1)==B);
+                                        pos_record=[pos_record; pos];
+                                        riftsegs=[riftsegs; md.segments(pos,:)];
+                                        if riftsegs(end,2)==tip1, 
+                                                break;
+                                        end
+                                end
+                                md.rifts(i).segments=riftsegs;
+                                md.rifts(i).tips=[tip1 tip2];
+
+                        else
+                                %ok, this is a rift that opens up to the domain outline.  One tip is going to be 
+                                %double, the other one, single. We are going to start from the single tip, towards the two 
+                                %other doubles
+
+                                %find tip1 and tip2 for this rift, in the new mesh created by yams.
+                                pos1=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
+                                tip1=md.segments(pos1,1);
+                                pos2=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
+                                tip2=md.segments(pos2,1);
+                                if length(tip1)==2,
+                                        %swap.
+                                        temp=tip1; tip1=tip2; tip2=temp;
+                                        temp=pos1; pos1=pos2; pos2=temp;
+                                        pos=pos1;
+                                else
+                                        pos=pos1;
+                                end
+
+                                pos_record=[pos_record; pos];
+                                riftsegs=md.segments(pos,:);
+                                while 1,
+                                        A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
+                                        %find other segment that holds B.
+                                        pos=find(md.segments(:,1)==B);
+                                        pos_record=[pos_record; pos];
+                                        riftsegs=[riftsegs; md.segments(pos,:)];
+                                        if ((riftsegs(end,2)==tip2(1)) | (riftsegs(end,2)==tip2(2))), 
+                                                %figure out which tip we reached
+                                                if riftsegs(end,2)==tip2(1), index=2; else index=1; end
+                                                break;
+                                        end
+                                end
+
+                                %ok, now, we start from the other tip2, towards tip1
+                                pos=pos2(index);
                                 pos_record=[pos_record; pos];
                                 riftsegs=[riftsegs; md.segments(pos,:)];
-                                if riftsegs(end,2)==tip1, 
-                                        break;
+                                while 1,
+                                        A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
+                                        %find other segment that holds B.
+                                        pos=find(md.segments(:,1)==B);
+                                        pos_record=[pos_record; pos];
+                                        riftsegs=[riftsegs; md.segments(pos,:)];
+                                        if riftsegs(end,2)==tip1, 
+                                                break;
+                                        end
                                 end
+                                md.rifts(i).segments=riftsegs;
+                                md.rifts(i).tips=[tip1 tip2(1) tip2(2)];
+
                         end
-                        md.rifts(i).segments=riftsegs;
-                        md.rifts(i).tips=[tip1 tip2];
                 end
         end

issm/trunk/src/m/classes/public/plot/plot_rifts.m

@@ -17 +17 @@
 if isstruct(md.rifts),
         
-        for i=1:size(md.rifts,1),
+        for i=1:md.numrifts,
                 penaltypairs=md.rifts(i).penaltypairs;
 
@@ -26 +26 @@
                         x(penaltypairs(j,1))=x(penaltypairs(j,1))-normal(1)*offset;
                         y(penaltypairs(j,1))=y(penaltypairs(j,1))-normal(2)*offset;
+                end
+                if length(md.rifts(i).tips)==3,
+                        tip=md.rifts(i).tips(3);
+                        %who is tip connected to? 
+                        if isconnected(md.elements,penaltypairs(1,1),tip),
+                                normal(1)=penaltypairs(1,5);
+                                normal(2)=penaltypairs(1,6);
+                                x(tip)=x(tip)-normal(1)*offset;
+                                y(tip)=y(tip)-normal(2)*offset;
+                        end
+
+                        if isconnected(md.elements,penaltypairs(1,2),tip),
+                                normal(1)=penaltypairs(1,5);
+                                normal(2)=penaltypairs(1,6);
+                                x(tip)=x(tip)+normal(1)*offset;
+                                y(tip)=y(tip)+normal(2)*offset;
+                        end
+                        if isconnected(md.elements,penaltypairs(end,1),tip),
+                                normal(1)=penaltypairs(end,5);
+                                normal(2)=penaltypairs(end,6);
+                                x(tip)=x(tip)-normal(1)*offset;
+                                y(tip)=y(tip)-normal(2)*offset;
+                        end
+                        if isconnected(md.elements,penaltypairs(end,2),tip),
+                                normal(1)=penaltypairs(end,5);
+                                normal(2)=penaltypairs(end,6);
+                                x(tip)=x(tip)+normal(1)*offset;
+                                y(tip)=y(tip)+normal(2)*offset;
+                        end
                 end
         end