Index: /issm/trunk/src/m/classes/public/bamg.m =================================================================== --- /issm/trunk/src/m/classes/public/bamg.m (revision 3249) +++ /issm/trunk/src/m/classes/public/bamg.m (revision 3250) @@ -62,48 +62,133 @@ 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 %}}}