function segments=findsegments(md,varargin) %FINDSEGMENTS - build segments model field % % Optional inputs: % 'mesh.elementconnectivity' % % Usage: % segments=findsegments(md,varargin); %get options options=pairoptions(varargin{:}); %Get connectivity mesh.elementconnectivity=getfieldvalue(options,'mesh.elementconnectivity',md.mesh.elementconnectivity); %Now, build the connectivity tables for this mesh if not correclty done if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements, if exist(options,'mesh.elementconnectivity'), error(' ''mesh.elementconnectivity'' option does not have thge right size.'); else mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); end end %Recreate the segments elementonboundary=double(mesh.elementconnectivity(:,3)==0); pos=find(elementonboundary); num_segments=length(pos); segments=zeros(num_segments,3); count=1; %loop over the segments for i=1:num_segments, %get current element on boundary el1=pos(i); %get elements connected to el1 els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:))); %el1 is connected to 2 other elements if length(els2)>1, %get nodes of el1 nods1=md.mesh.elements(el1,:); %find the common vertices to the two elements connected to el1 (1 or 2) flag=intersect(md.mesh.elements(els2(1),:),md.mesh.elements(els2(2),:)); %get the vertices on the boundary and build segment nods1(find(ismember(nods1,flag)))=[]; segments(count,:)=[nods1 el1]; %swap segment nodes if necessary ord1=find(nods1(1)==md.mesh.elements(el1,:)); ord2=find(nods1(2)==md.mesh.elements(el1,:)); if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), temp=segments(count,1); segments(count,1)=segments(count,2); segments(count,2)=temp; end segments(count,1:2)=fliplr(segments(count,1:2)); count=count+1; %el1 is connected to only one element else %get nodes of el1 nods1=md.mesh.elements(el1,:); %find the vertex the el1 to not share with els2 flag=setdiff(nods1,md.mesh.elements(els2,:)); for j=1:3, nods=nods1; nods(j)=[]; if any(ismember(flag,nods)), segments(count,:)=[nods el1]; %swap segment nodes if necessary ord1=find(nods(1)==md.mesh.elements(el1,:)); ord2=find(nods(2)==md.mesh.elements(el1,:)); if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ), temp=segments(count,1); segments(count,1)=segments(count,2); segments(count,2)=temp; end segments(count,1:2)=fliplr(segments(count,1:2)); count=count+1; end end end end