source: issm/trunk/src/m/model/contourenvelope.m@ 9733

Last change on this file since 9733 was 9733, checked in by seroussi, 14 years ago

keep building mesh

File size: 3.6 KB
RevLine 
[2518]1function segments=contourenvelope(md,varargin)
[2517]2%CONTOURENVELOPE - build a set of segments enveloping a contour .exp
3%
4% Usage:
[2518]5% segments=contourenvelope(md,varargin)
6%
7% Example:
8% segments=contourenvelope(md,'Stream.exp');
[9641]9% segments=contourenvelope(md,md.mask.elementonfloatingice)
[2518]10% segments=contourenvelope(md);
[2517]11
12%some checks
[2518]13if nargin>2,
14 help contourenvelope
15 error('contourenvelope error message: bad usage');
[2517]16end
[2518]17if nargin==2,
[8730]18 flags=varargin{1};
19
20 if ischar(flags),
21 file=flags;
22 file=varargin{1};
23 if ~exist(file),
24 error(['contourenvelope error message: file ' file ' not found']);
25 end
26 isfile=1;
27 elseif isnumeric(flags),
28 %do nothing for now
29 isfile=0;
30 else
31 error('contourenvelope error message: second argument should a file or an elements flag');
[2518]32 end
33end
[2517]34
35%Now, build the connectivity tables for this mesh.
36%Computing connectivity
[9729]37if size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices,
[9733]38 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
[2517]39end
[9728]40if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements,
[9733]41 md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
[2517]42end
43
44%get nodes inside profile
[9728]45mesh.elementconnectivity=md.mesh.elementconnectivity;
[2518]46if nargin==2,
[8730]47 if isfile,
48 %get flag list of elements and nodes inside the contour
[9733]49 nodein=ContourToMesh(md.mesh.elements,md.x,md.y,file,'node',1);
50 elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2));
[8730]51 %modify element connectivity
52 elemout=find(~elemin);
[9728]53 mesh.elementconnectivity(elemout,:)=0;
54 mesh.elementconnectivity(find(ismember(mesh.elementconnectivity,elemout)))=0;
[8730]55 else
56 %get flag list of elements and nodes inside the contour
[9725]57 nodein=zeros(md.mesh.numberofvertices,1);
58 elemin=zeros(md.mesh.numberofelements,1);
[8730]59
60 pos=find(flags);
61 elemin(pos)=1;
[9733]62 nodein(md.mesh.elements(pos,:))=1;
[8730]63
64 %modify element connectivity
65 elemout=find(~elemin);
[9728]66 mesh.elementconnectivity(elemout,:)=0;
67 mesh.elementconnectivity(find(ismember(mesh.elementconnectivity,elemout)))=0;
[8730]68 end
[2518]69end
[2517]70
71%Find element on boundary
72%First: find elements on the boundary of the domain
[9728]73flag=mesh.elementconnectivity;
[2518]74if nargin==2,
75 flag(find(flag))=elemin(flag(find(flag)));
76end
[5602]77elementonboundary=double(prod(flag,2)==0 & sum(flag,2)>0);
[2517]78
79%Find segments on boundary
80pos=find(elementonboundary);
81num_segments=length(pos);
82segments=zeros(num_segments,3);
83count=1;
84
85for i=1:num_segments,
86 el1=pos(i);
[9728]87 els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:)));
[2517]88 if length(els2)>1,
[9733]89 flag=intersect(md.mesh.elements(els2(1),:),md.mesh.elements(els2(2),:));
90 nods1=md.mesh.elements(el1,:);
[2517]91 nods1(find(nods1==flag))=[];
92 segments(count,:)=[nods1 el1];
93
[9733]94 ord1=find(nods1(1)==md.mesh.elements(el1,:));
95 ord2=find(nods1(2)==md.mesh.elements(el1,:));
[2517]96
[8298]97 %swap segment nodes if necessary
[2517]98 if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
99 temp=segments(count,1);
100 segments(count,1)=segments(count,2);
101 segments(count,2)=temp;
102 end
103 segments(count,1:2)=fliplr(segments(count,1:2));
104 count=count+1;
105 else
[9733]106 nods1=md.mesh.elements(el1,:);
107 flag=setdiff(nods1,md.mesh.elements(els2,:));
[2517]108 for j=1:3,
109 nods=nods1; nods(j)=[];
110 if any(ismember(flag,nods)),
111 segments(count,:)=[nods el1];
[9733]112 ord1=find(nods(1)==md.mesh.elements(el1,:));
113 ord2=find(nods(2)==md.mesh.elements(el1,:));
[2517]114 if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
115 temp=segments(count,1);
116 segments(count,1)=segments(count,2);
117 segments(count,2)=temp;
118 end
119 segments(count,1:2)=fliplr(segments(count,1:2));
120 count=count+1;
121 end
122 end
123 end
124end
Note: See TracBrowser for help on using the repository browser.