[2518] | 1 | function 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] | 13 | if nargin>2,
|
---|
| 14 | help contourenvelope
|
---|
| 15 | error('contourenvelope error message: bad usage');
|
---|
[2517] | 16 | end
|
---|
[2518] | 17 | if 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
|
---|
| 33 | end
|
---|
[2517] | 34 |
|
---|
| 35 | %Now, build the connectivity tables for this mesh.
|
---|
| 36 | %Computing connectivity
|
---|
[9725] | 37 | if size(md.nodeconnectivity,1)~=md.mesh.numberofvertices,
|
---|
| 38 | md.nodeconnectivity=NodeConnectivity(md.elements,md.mesh.numberofvertices);
|
---|
[2517] | 39 | end
|
---|
[9725] | 40 | if size(md.elementconnectivity,1)~=md.mesh.numberofelements,
|
---|
[2517] | 41 | md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
|
---|
| 42 | end
|
---|
| 43 |
|
---|
| 44 | %get nodes inside profile
|
---|
| 45 | elementconnectivity=md.elementconnectivity;
|
---|
[2518] | 46 | if nargin==2,
|
---|
[8730] | 47 | if isfile,
|
---|
| 48 | %get flag list of elements and nodes inside the contour
|
---|
| 49 | nodein=ContourToMesh(md.elements,md.x,md.y,file,'node',1);
|
---|
| 50 | elemin=(sum(nodein(md.elements),2)==size(md.elements,2));
|
---|
| 51 | %modify element connectivity
|
---|
| 52 | elemout=find(~elemin);
|
---|
| 53 | elementconnectivity(elemout,:)=0;
|
---|
| 54 | elementconnectivity(find(ismember(elementconnectivity,elemout)))=0;
|
---|
| 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;
|
---|
| 62 | nodein(md.elements(pos,:))=1;
|
---|
| 63 |
|
---|
| 64 | %modify element connectivity
|
---|
| 65 | elemout=find(~elemin);
|
---|
| 66 | elementconnectivity(elemout,:)=0;
|
---|
| 67 | elementconnectivity(find(ismember(elementconnectivity,elemout)))=0;
|
---|
| 68 | end
|
---|
[2518] | 69 | end
|
---|
[2517] | 70 |
|
---|
| 71 | %Find element on boundary
|
---|
| 72 | %First: find elements on the boundary of the domain
|
---|
| 73 | flag=elementconnectivity;
|
---|
[2518] | 74 | if nargin==2,
|
---|
| 75 | flag(find(flag))=elemin(flag(find(flag)));
|
---|
| 76 | end
|
---|
[5602] | 77 | elementonboundary=double(prod(flag,2)==0 & sum(flag,2)>0);
|
---|
[2517] | 78 |
|
---|
| 79 | %Find segments on boundary
|
---|
| 80 | pos=find(elementonboundary);
|
---|
| 81 | num_segments=length(pos);
|
---|
| 82 | segments=zeros(num_segments,3);
|
---|
| 83 | count=1;
|
---|
| 84 |
|
---|
| 85 | for i=1:num_segments,
|
---|
| 86 | el1=pos(i);
|
---|
| 87 | els2=elementconnectivity(el1,find(elementconnectivity(el1,:)));
|
---|
| 88 | if length(els2)>1,
|
---|
| 89 | flag=intersect(md.elements(els2(1),:),md.elements(els2(2),:));
|
---|
| 90 | nods1=md.elements(el1,:);
|
---|
| 91 | nods1(find(nods1==flag))=[];
|
---|
| 92 | segments(count,:)=[nods1 el1];
|
---|
| 93 |
|
---|
| 94 | ord1=find(nods1(1)==md.elements(el1,:));
|
---|
| 95 | ord2=find(nods1(2)==md.elements(el1,:));
|
---|
| 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
|
---|
| 106 | nods1=md.elements(el1,:);
|
---|
| 107 | flag=setdiff(nods1,md.elements(els2,:));
|
---|
| 108 | for j=1:3,
|
---|
| 109 | nods=nods1; nods(j)=[];
|
---|
| 110 | if any(ismember(flag,nods)),
|
---|
| 111 | segments(count,:)=[nods el1];
|
---|
| 112 | ord1=find(nods(1)==md.elements(el1,:));
|
---|
| 113 | ord2=find(nods(2)==md.elements(el1,:));
|
---|
| 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
|
---|
| 124 | end
|
---|