source: issm/trunk/src/m/mesh/findsegments.m@ 13395

Last change on this file since 13395 was 13395, checked in by Mathieu Morlighem, 12 years ago

merged trunk-jpl and trunk for revision 13393

File size: 2.6 KB
RevLine 
[3115]1function segments=findsegments(md,varargin)
[1346]2%FINDSEGMENTS - build segments model field
3%
[3115]4% Optional inputs:
[9728]5% 'mesh.elementconnectivity'
[3115]6%
[1346]7% Usage:
[3115]8% segments=findsegments(md,varargin);
[1346]9
[3115]10%get options
11options=pairoptions(varargin{:});
12
13%Get connectivity
[9728]14mesh.elementconnectivity=getfieldvalue(options,'mesh.elementconnectivity',md.mesh.elementconnectivity);
[3115]15
16%Now, build the connectivity tables for this mesh if not correclty done
[9728]17if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements,
18 if exist(options,'mesh.elementconnectivity'),
19 error(' ''mesh.elementconnectivity'' option does not have thge right size.');
[3115]20 else
[9733]21 mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
[3115]22 end
[2082]23end
[1346]24
25%Recreate the segments
[9728]26elementonboundary=double(mesh.elementconnectivity(:,3)==0);
[1346]27pos=find(elementonboundary);
28num_segments=length(pos);
29segments=zeros(num_segments,3);
[1347]30count=1;
31
[2612]32%loop over the segments
[1346]33for i=1:num_segments,
[2612]34
35 %get current element on boundary
[1346]36 el1=pos(i);
[2612]37
38 %get elements connected to el1
[9728]39 els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:)));
[2612]40
41 %el1 is connected to 2 other elements
[1347]42 if length(els2)>1,
[2612]43
44 %get nodes of el1
[9733]45 nods1=md.mesh.elements(el1,:);
[2612]46
47 %find the common vertices to the two elements connected to el1 (1 or 2)
[9733]48 flag=intersect(md.mesh.elements(els2(1),:),md.mesh.elements(els2(2),:));
[2612]49
50 %get the vertices on the boundary and build segment
51 nods1(find(ismember(nods1,flag)))=[];
[1347]52 segments(count,:)=[nods1 el1];
[1346]53
[8298]54 %swap segment nodes if necessary
[9733]55 ord1=find(nods1(1)==md.mesh.elements(el1,:));
56 ord2=find(nods1(2)==md.mesh.elements(el1,:));
[1347]57 if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
58 temp=segments(count,1);
59 segments(count,1)=segments(count,2);
60 segments(count,2)=temp;
61 end
62 segments(count,1:2)=fliplr(segments(count,1:2));
63 count=count+1;
[2612]64
65 %el1 is connected to only one element
[1347]66 else
[2612]67 %get nodes of el1
[9733]68 nods1=md.mesh.elements(el1,:);
[2612]69
70 %find the vertex the el1 to not share with els2
[9733]71 flag=setdiff(nods1,md.mesh.elements(els2,:));
[2612]72
[1347]73 for j=1:3,
74 nods=nods1; nods(j)=[];
75 if any(ismember(flag,nods)),
[2612]76
[1347]77 segments(count,:)=[nods el1];
[2612]78
[8298]79 %swap segment nodes if necessary
[9733]80 ord1=find(nods(1)==md.mesh.elements(el1,:));
81 ord2=find(nods(2)==md.mesh.elements(el1,:));
[1347]82 if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
83 temp=segments(count,1);
84 segments(count,1)=segments(count,2);
85 segments(count,2)=temp;
86 end
87 segments(count,1:2)=fliplr(segments(count,1:2));
88 count=count+1;
89 end
90 end
[1346]91 end
92end
Note: See TracBrowser for help on using the repository browser.