Changeset 26849


Ignore:
Timestamp:
02/02/22 07:21:47 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: merging levelzero functions into isoline

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  
    11function 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 compatible
    4 %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_CONTOUR
    172
    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  
    11function expcontourlevelzero(md,mask,level,filename)
    2 %EXPCONTOURLEVELZERO - write an Argus file from a structure recovered from running contourlevelzero
    3 %
    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, EXPWRITE
    122
    13 contours=contourlevelzero(md,mask,level);
    14 expwrite(contours,filename);
     3error(['this function has been renamed: isoline(md,mask,''output'',''' filename ''');']);
Note: See TracChangeset for help on using the changeset viewer.