Changeset 3250

Timestamp:

03/10/10 16:27:21 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added rift on border capability

File:

• issm/trunk/src/m/classes/public/bamg.m (modified) (1 diff)

issm/trunk/src/m/classes/public/bamg.m

@@ -62 +62 @@
 bamg_geometry.SubDomains=zeros(0,3);
 if exist(options,'domain'),
+
+        %Check that file exists
         domainfile=getfieldvalueerr(options,'domain');
         if ~exist(domainfile,'file')
                 error(['bamg error message: file ' domainfile ' not found ']);
-        else
-                domain=expread(domainfile);
-
-                %build geometry
-                count=0;
-
-                for i=1:length(domain),
-
-                        %if current profile is closed
-                        if (domain(i).x(1)==domain(i).x(end) & domain(i).y(1)==domain(i).y(end)),
-                                nods=domain(i).nods-1; %the domain are closed 1=end;
-                                bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
-                                bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  1*ones(nods,1)]];
-                                if i>1,
-                                        %if closed : hole
-                                        clockwise=-1;
-                                        bamg_geometry.SubDomains=[2 count+1 clockwise 1];
-
+        end
+
+        %Build geometry 
+        domain=expread(domainfile);
+        count=0;
+        for i=1:length(domain),
+
+                %Check that the domain is closed
+                if (domain(i).x(1)~=domain(i).x(end) | domain(i).y(1)~=domain(i).y(end)),
+                        error('bamg error message: all contours provided in ''domain'' should be closed');
+                end
+                %Checks that all holes are INSIDE the principla domain outline
+                if i>1,
+                        flags=ContourToNodes(domain(i).x,domain(i).y,domain(1),0);
+                        if any(~flags),
+                                error('bamg error message: All holes should be stricly inside the principal domain');
+                        end
+                end
+
+                %Add all points to bamg_geometry
+                nods=domain(i).nods-1; %the domain are closed 1=end;
+                bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(1:nods) domain(i).y(1:nods) ones(nods,1)]];
+                bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods count+1])  1*ones(nods,1)]];
+
+                %update counter
+                count=count+nods;
+        end
+
+        %take care of rifts
+        if exist(options,'rifts'),
+
+                %Check that file exists
+                riftfile=getfieldvalueerr(options,'rifts');
+                if ~exist(riftfile,'file')
+                        error(['bamg error message: file ' riftfile ' not found ']);
+                end
+                rift=expread(riftfile);
+
+                for i=1:length(rift),
+
+                        %detect wether all points of the rift are inside the domain
+                        flags=ContourToNodes(rift(i).x,rift(i).y,domain(1),0);
+                        if ~flags,
+                                error('one Rift has all his points outside of the domain outline'),
+
+                        elseif any(~flags),
+                                %We LOTS of work to do
+                                disp('Rift tip outside of or on the domain has been detected and is being processed...');
+
+                                %check that only one point is outsie (for now)
+                                if sum(~flags)~=1,
+                                        error('bamg error message: only one point outside of the domain is supported yet');
                                 end
 
-                        elseif i>1
-                                %rift
-                                nods=domain(i).nods-1;
-                                bamg_geometry.Vertices=[bamg_geometry.Vertices; [domain(i).x(:) domain(i).y(:) ones(nods+1,1)]];
+                                %Move tip outside to the first position
+                                if flags(1)==0,
+                                        %OK, first point is outside (do nothing),
+                                elseif (flags(end)==0),
+                                        rift(i).x=flipud(rift(i).x);
+                                        rift(i).y=flipud(rift(i).y);
+                                else
+                                        error('bamg error message: only a rift tip can be outside of the domain');
+                                end
+
+                                %Get cordinate of intersection point
+                                x1=rift(i).x(1); y1=rift(i).y(1);
+                                x2=rift(i).x(2); y2=rift(i).y(2);
+                                for j=1:length(domain(1).x)-1;
+                                        if SegIntersect([x1 y1; x2 y2],[domain(1).x(j) domain(1).y(j); domain(1).x(j+1) domain(1).y(j+1)]),
+
+                                                %Get position of the two nodes of the edge in domain
+                                                i1=j;
+                                                i2=mod(j+1,domain(1).nods-1);
+
+                                                %rift is crossing edge [i1 i2] of the domain
+                                                %Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
+                                                x3=domain(1).x(i1); y3=domain(1).y(i1);
+                                                x4=domain(1).x(i2); y4=domain(1).y(i2);
+                                                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 intersection point to Vertices
+                                                bamg_geometry.Vertices=[bamg_geometry.Vertices; x y 1];
+                                                count=count+1;
+
+                                                %Decompose the crossing edge in 2 subedges
+                                                pos=find(bamg_geometry.Edges(:,1)==i1 & bamg_geometry.Edges(:,2)==i2);
+                                                if isempty(pos) error('bamg error message: a problem occured...'); end
+                                                bamg_geometry.Edges=[bamg_geometry.Edges(1:pos-1,:);...
+                                                        bamg_geometry.Edges(pos,1) count                           bamg_geometry.Edges(pos,3);...
+                                                        count                      bamg_geometry.Edges(pos,2)   bamg_geometry.Edges(pos,3);...
+                                                        bamg_geometry.Edges(pos+1:end,:)];
+
+                                                %OK, no we can add our own rift
+                                                nods=rift(i).nods-1;
+                                                bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]];
+                                                bamg_geometry.Edges=[bamg_geometry.Edges;...
+                                                        count  count+1  2 ;...
+                                                        [transp(count+1:count+nods-1) transp([count+2:count+nods])  2*ones(nods-1,1)]];
+                                                count=count+nods;
+
+                                                disp(' done'); break;
+                                        end
+                                end
+                        else
+                                nods=rift(i).nods-1;
+                                bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(:) rift(i).y(:) ones(nods+1,1)]];
                                 bamg_geometry.Edges=[bamg_geometry.Edges; [transp(count+1:count+nods) transp([count+2:count+nods+1])  2*ones(nods,1)]];
-
-                        else
-                                error('bamg error message: the first domain should be closed');
+                                count=count+nods+1;
                         end
-                        count=count+nods;
-                end
-                if exist(options,'hVertices'),
-                        bamg_geometry.hVertices=getfieldvalueerr(options,'hVertices');
-                        if length(bamg_geometry.hVertices)~=nods,
-                                error(['hVertices option should be of size [' num2str(nods) ',1]']);
-                        end
-                end
-                bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
-                bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
-                bamg_geometry.NumSubDomains=size(bamg_geometry.SubDomains,1);
-        end 
+                end
+        end
+
+        if exist(options,'hVertices'),
+                bamg_geometry.hVertices=getfieldvalueerr(options,'hVertices');
+                if length(bamg_geometry.hVertices)~=nods,
+                        error(['hVertices option should be of size [' num2str(nods) ',1]']);
+                end
+        end
+
+        %update other fields of bamg_geometry
+        bamg_geometry.NumVertices=size(bamg_geometry.Vertices,1);
+        bamg_geometry.NumEdges=size(bamg_geometry.Edges,1);
+        bamg_geometry.NumSubDomains=size(bamg_geometry.SubDomains,1);
 end
 %}}}