Inconsistency with plotchannels

mwillette

Hello,
I've been doing GlaDS hydrology runs and trying to plot the edges directly connected to the grounding line. I've selected all edges in md.mesh.edges connected to a grounding line node, and the nodes that make those connections (found with a loop first, then verified manually). When I make a plot highlighting these nodes, it appears correctly.

When I use plotchannels to plot these edges only, some appear missing and there are some additional edges that are two nodes away from the grounding line, but I'm certain the missing edges are in the set and the extraneous edges are not.
channels=md5M.results.TransientSolution(10).ChannelDischarge; % selecting a time when all discharge is 0
channels(GLedgenums)=25; % making it so only the grounding line edges are visible
plotchannels(md5M,channels,'min',1,'colormap',themap,'max',25,'log',10)

Then, when I try to plot a subset of the grounding line edges,
plotchannels(md5M,abs(md5M.results.TransientSolution(i).ChannelDischarge(topedges)),...
different edges that haven't appeared in either of the previous plots show up.

Does anyone know what might be happening here? Do md.mesh.edges and md.results.TransientSolution.ChannelDischarge have different numbering systems?

Thanks,
Michael

mathieumorlighem

mwillette Hi Michael

the short answer is YES, they have a different numbering system... sorry! Here is how they are created by GlaDS and plotchannel (which is not the same ordering as the one from bamg):

%Maximum number of edges
maxnbf = 3*md.mesh.numberofelements;
%Initialize intermediaries
edges = zeros(maxnbf,3);
exchange = zeros(maxnbf,1);
%Chaining algorithm
head_minv = -1*ones(md.mesh.numberofvertices,1);
next_face = zeros(maxnbf,1);
%Initialize number of faces
nbf = 0;
for i=1:md.mesh.numberofelements
   for j=1:3
      %Get the two indices of the edge number j of the ith triangle
      v1 = md.mesh.elements(i,j);
      if(j==3)
         v2 = md.mesh.elements(i,1);
      else
         v2 = md.mesh.elements(i,j+1);
      end
      %sort
      if(v2<v1)
         v3=v2; v2=v1; v1=v3;
      end
      %This edge a priori has not been processed yet
      exists = false;
      %Go through all processed faces connected to v1 and check whether we have seen this edge yet
      e=head_minv(v1);
      while(e~=-1)
         if(edges(e,2)==v2)
            exists = true;
            break;
         end
         e=next_face(e);
      end
      %If this edge is new, add it to the lists
      if(~exists)
         %Update edges
         edges(nbf+1,1) = v1; %vertex 1
         edges(nbf+1,2) = v2; %vertex 2
         edges(nbf+1,3) = i;  %element 1 (ignore second one)
         if(v1~=md.mesh.elements(i,j)) exchange(nbf+1)=1; end
         %Update chain
         next_face(nbf+1) = head_minv(v1);
         head_minv(v1)    = nbf+1;
         %Increase number of faces
         nbf=nbf+1;
      end
   end
end

edges = edges(1:nbf,:);
pos = find(exchange);
v3 = edges(pos,1);
edges(pos,1) = edges(pos,2);
edges(pos,2) = v3;

Let me know if you have questions!
Cheers
Mathieu

mwillette

Hi Mathieu,

Thanks for the quick response, this should be all I need to sort out the channels.

All the best,
Michael