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

Last change on this file was 25836, checked in by Mathieu Morlighem, 4 years ago

merged trunk-jpl and trunk for revision 25834

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