Changeset 26849
- Timestamp:
- 02/02/22 07:21:47 (3 years ago)
- Location:
- issm/trunk-jpl/src/m/exp
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/exp/contourlevelzero.m
r26792 r26849 1 1 function contours=contourlevelzero(md,mask,level,varargin) 2 %CONTOURLEVELZERO - figure out the zero level (or offset thereof, specified by the level value)3 % of a vectorial mask, and vectorialize it into an exp or shp compatible4 %structure.5 %6 % Usage:7 % contours=contourlevelzero(md,mask,level)8 %9 % Example:10 % contours=contourlevelzero(md, md.results.TransientSolution(end).MaskOceanLevelset,0);11 % contours=contourlevelzero(md, md.results.TransientSolution(end).MaskOceanLevelset,0,'output','vector');12 %13 % Supported options:14 % 'output': either 'matrix' (with NaN delimiters) or 'struct' (default)15 %16 % See also: PLOT_CONTOUR17 2 18 %Process options 19 options = pairoptions(varargin{:}); 20 21 %process data 22 if dimension(md.mesh)==3, 23 % error('contourlevelzero error message: routine not supported for 3d meshes, project on a layer'); 24 x = md.mesh.x2d; 25 y = md.mesh.y2d; 26 index=md.mesh.elements2d; 27 else 28 x=md.mesh.x; 29 y=md.mesh.y; 30 index=md.mesh.elements; 31 end 32 33 if isprop(md.mesh,'z'), 34 z=md.mesh.z; 35 else 36 z=zeros(md.mesh.numberofvertices,1); 37 end 38 39 if isempty(mask), error('mask provided is empty'); end 40 if dimension(md.mesh)==3, 41 if length(mask)~=md.mesh.numberofvertices2d, error('mask provided should be specified at the vertices of the mesh'); end 42 else 43 if length(mask)~=md.mesh.numberofvertices, error('mask provided should be specified at the vertices of the mesh'); end 44 end 45 46 %initialization of some variables 47 numberofelements=size(index,1); 48 elementslist=1:numberofelements; 49 c=[]; 50 h=[]; 51 52 %get unique edges in mesh 53 %1: list of edges 54 edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])]; 55 %2: find unique edges 56 [edges,I,J]=unique(sort(edges,2),'rows'); 57 %3: unique edge numbers 58 vec=J; 59 %4: unique edges numbers in each triangle (2 triangles sharing the same edge will have 60 % the same edge number) 61 edges_tria=[vec(elementslist), vec(elementslist+numberofelements), vec(elementslist+2*numberofelements)]; 62 63 %segments [nodes1 nodes2] 64 Seg1=index(:,[1 2]); 65 Seg2=index(:,[2 3]); 66 Seg3=index(:,[3 1]); 67 68 %segment numbers [1;4;6;...] 69 Seg1_num=edges_tria(:,1); 70 Seg2_num=edges_tria(:,2); 71 Seg3_num=edges_tria(:,3); 72 73 %value of data on each tips of the segments 74 Data1=mask(Seg1); 75 Data2=mask(Seg2); 76 Data3=mask(Seg3); 77 78 %get the ranges for each segment 79 Range1=sort(Data1,2); 80 Range2=sort(Data2,2); 81 Range3=sort(Data3,2); 82 83 %find the segments that contain this value 84 pos1=(Range1(:,1)<level & Range1(:,2)>=level); 85 pos2=(Range2(:,1)<level & Range2(:,2)>=level); 86 pos3=(Range3(:,1)<level & Range3(:,2)>=level); 87 88 %get elements 89 poselem12=(pos1 & pos2); 90 poselem13=(pos1 & pos3); 91 poselem23=(pos2 & pos3); 92 poselem=find(poselem12 | poselem13 | poselem23); 93 numelems=length(poselem); 94 95 %if no element has been flagged, skip to the next level 96 if numelems==0, 97 warning('contourlevelzero warning message: no elements found with corresponding level value in mask'); 98 contours=struct([]); 99 return; 100 end 101 102 %go through the elements and build the coordinates for each segment (1 by element) 103 x1=zeros(numelems,1); 104 x2=zeros(numelems,1); 105 y1=zeros(numelems,1); 106 y2=zeros(numelems,1); 107 z1=zeros(numelems,1); 108 z2=zeros(numelems,1); 109 110 edge_l=zeros(numelems,2); 111 112 for j=1:numelems, 113 114 weight1=(level-Data1(poselem(j),1))/(Data1(poselem(j),2)-Data1(poselem(j),1)); 115 weight2=(level-Data2(poselem(j),1))/(Data2(poselem(j),2)-Data2(poselem(j),1)); 116 weight3=(level-Data3(poselem(j),1))/(Data3(poselem(j),2)-Data3(poselem(j),1)); 117 118 if poselem12(poselem(j)); 119 120 x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1))); 121 x2(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1))); 122 y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1))); 123 y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1))); 124 z1(j)=z(Seg1(poselem(j),1))+weight1*(z(Seg1(poselem(j),2))-z(Seg1(poselem(j),1))); 125 z2(j)=z(Seg2(poselem(j),1))+weight2*(z(Seg2(poselem(j),2))-z(Seg2(poselem(j),1))); 126 127 edge_l(j,1)=Seg1_num(poselem(j)); 128 edge_l(j,2)=Seg2_num(poselem(j)); 129 130 elseif poselem13(poselem(j)), 131 132 x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1))); 133 x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1))); 134 y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1))); 135 y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1))); 136 z1(j)=z(Seg1(poselem(j),1))+weight1*(z(Seg1(poselem(j),2))-z(Seg1(poselem(j),1))); 137 z2(j)=z(Seg3(poselem(j),1))+weight3*(z(Seg3(poselem(j),2))-z(Seg3(poselem(j),1))); 138 139 edge_l(j,1)=Seg1_num(poselem(j)); 140 edge_l(j,2)=Seg3_num(poselem(j)); 141 142 elseif poselem23(poselem(j)), 143 144 x1(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1))); 145 x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1))); 146 y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1))); 147 y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1))); 148 z1(j)=z(Seg2(poselem(j),1))+weight2*(z(Seg2(poselem(j),2))-z(Seg2(poselem(j),1))); 149 z2(j)=z(Seg3(poselem(j),1))+weight3*(z(Seg3(poselem(j),2))-z(Seg3(poselem(j),1))); 150 151 edge_l(j,1)=Seg2_num(poselem(j)); 152 edge_l(j,2)=Seg3_num(poselem(j)); 153 else 154 %it shoud not go here 155 end 156 end 157 158 %now that we have the segments, we must try to connect them... 159 160 %loop over the subcontours 161 contours=struct([]); 162 163 while ~isempty(edge_l), 164 165 %take the right edge of the second segment and connect it to the next segments if any 166 e1=edge_l(1,1); e2=edge_l(1,2); 167 xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)]; zc=[z1(1);z2(1)]; 168 169 170 %erase the lines corresponding to this edge 171 edge_l(1,:)=[]; 172 x1(1)=[]; x2(1)=[]; 173 y1(1)=[]; y2(1)=[]; 174 z1(1)=[]; z2(1)=[]; 175 176 [ro1,co1]=find(edge_l==e1); 177 178 while ~isempty(ro1) 179 180 if co1==1, 181 xc=[x2(ro1);xc]; yc=[y2(ro1);yc];zc=[z2(ro1);zc]; 182 183 %next edge: 184 e1=edge_l(ro1,2); 185 186 else 187 xc=[x1(ro1);xc]; yc=[y1(ro1);yc];zc=[z1(ro1);zc]; 188 189 %next edge: 190 e1=edge_l(ro1,1); 191 end 192 193 %erase the lines of this 194 edge_l(ro1,:)=[]; 195 x1(ro1)=[]; x2(ro1)=[]; 196 y1(ro1)=[]; y2(ro1)=[]; 197 z1(ro1)=[]; z2(ro1)=[]; 198 199 %next connection 200 [ro1,co1]=find(edge_l==e1); 201 end 202 203 %same thing the other way (to the right) 204 [ro2,co2]=find(edge_l==e2); 205 206 while ~isempty(ro2) 207 208 if co2==1, 209 xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];zc=[zc;z2(ro2)]; 210 211 %next edge: 212 e2=edge_l(ro2,2); 213 else 214 xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)]; zc=[zc;z1(ro2)]; 215 216 %next edge: 217 e2=edge_l(ro2,1); 218 end 219 220 %erase the lines of this 221 edge_l(ro2,:)=[]; 222 x1(ro2)=[]; x2(ro2)=[]; 223 y1(ro2)=[]; y2(ro2)=[]; 224 z1(ro2)=[]; z2(ro2)=[]; 225 226 %next connection 227 [ro2,co2]=find(edge_l==e2); 228 end 229 230 %save xc,yc contour: 231 contours(end+1).x=xc; 232 contours(end).y=yc; 233 contours(end).z=zc; 234 contours(end).name=''; 235 contours(end).nods=length(xc); 236 contours(end).density=1; 237 contours(end).closed=0; 238 239 end 240 241 %process output 242 if strcmp(getfieldvalue(options,'output',1),'matrix') 243 x = cell2mat(cellfun(@(x) [x;NaN],{contours(:).x},'UniformOutput',0)'); 244 y = cell2mat(cellfun(@(x) [x;NaN],{contours(:).y},'UniformOutput',0)'); 245 246 contours = [x y]; 247 end 3 error('this function has been renamed: A = isoline(md,mask);'); -
issm/trunk-jpl/src/m/exp/expcontourlevelzero.m
r26202 r26849 1 1 function expcontourlevelzero(md,mask,level,filename) 2 %EXPCONTOURLEVELZERO - write an Argus file from a structure recovered from running contourlevelzero3 %4 % Usage:5 % expcontourlevelzero(md,mask,level,filename)6 %7 % Example:8 % expcontourlevelzero(md,md.geometry.thickness,10, 'Level0.exp');9 % expcontourlevelzero(md,md.mask.ocean_levelset,0, 'Level0.exp');10 %11 % See also CONTOURLEVELZERO, EXPWRITE12 2 13 contours=contourlevelzero(md,mask,level); 14 expwrite(contours,filename); 3 error(['this function has been renamed: isoline(md,mask,''output'',''' filename ''');']);
Note:
See TracChangeset
for help on using the changeset viewer.