function [tri,tet,pen]=bods(x,xa,y,ya,z,za,F,Fa)


The computational array is body-centered-cubic.  The physical array
need not be strictly bcc; it can be bent bcc, but not bent so far as
to change the topology.  We have
to be able to find x,y,z (& F) on the bases of the six pyramids that 
surround each of the central nodes at xa,ya,za (& Fa).


% function [tri,tet,pen]=bods(x,xa,y,ya,z,za,F,Fa)
%
% x,y,z,F is data on a simple cubic array
% xa,ya,za,Fa is data at the centers
% cols of tri are [xxxyyyzzz]' near F=0
% cols of tet are [xxxxyyyyzzzz]' near F=0
% cols of pen are [xxxxxyyyyyzzzzz]' near F=0


Just the help message for bods().  The surface patches for a five-sided
pyramid are triangles, tetragons and pentagons.  The returned data can 
easily be reshaped as rows or cols of x,y,z.


if nargin<8 help bods; return; end

% perturb accidental zeros

N=find(F==0); F(N)=eps*(2*round(rand(size(N)))-1);
N=find(Fa==0); Fa(N)=eps*(2*round(rand(size(N)))-1);


Finding the "accidental zeros" at vertices is an important trick.  Once
they have been replaced by arbitrarily small nonzeros, there are no special
cases in the task of finding zeros on edges.  The accidental zeros almost
never occur in measured data - it's the test cases that often do us in.


% wesnba indices for x,y,z,F in cubes


Leftover from happy days at Thinking Machines, tending the old connection machine.


N=size(F); NN=reshape(1:prod(N),N);


(x,y,z) were generated by matlab's meshgrid().  I don't remember how 
that worked any more, but it doesn't matter!  The arrays were
x(i,j,k), y(i,j,k) and z(i,j,k) with i=1:N(1), j=1:N(2) and k=1:N(3), 
The array F(i,j,k) was the function F(x(i,j,k),y(i,j,k),z(i,j,k)).  I
think the i-axis is the x-axis in function style (left to right), but 
the algorithms don't depend on that.  Call the i-axis the (west to east)-
axis, the j-axis the (south to north)-axis and the k-axis the (below 
to above)-axis.  (It's news-axes gone crazy.)  NN is just the natural 
numbers in the same kind of array



NN(i,j,k)=i+N(1)(j-1)+N(1)N(2)(k-1)


(Anybody remember ravel in APL?)  Incidentally, Matlab enumerates
in Fortran-style, so the first col of the first of N(3) arrays it would
display on a page is [1:N(1)]'  That would put west at the top of the arrays,
south at the left, and we would stack later arrays on top of
earlier ones. I could be wrong about that, but it wouldn't matter.


NN=reshape(NN(1:end-1,1:end-1,1:end-1),1,prod(N-1)); % [b-sw corners;


NN(1,1,1)=1 is as far west, as far south, as low and as small as it gets,
so these are indices of F(:) at lower southwest corners of the cubes that
contain the central nodes Xa,Ya,Za (& Fa).  They are arranged in a row, 
ready to have the other seven rows placed beneath them.


NN=[NN;N(1)+NN]; NN=[NN;1+NN]; NN=[NN;N(1)*N(2)+NN]; %  to ;a-ne corners]

row:   1      2      3      4      5      6      7      8
    [ b-sw ; b-nw ; b-se ; b-ne ; a-sw ; a-nw ; a-se ; a-ne ] 


xa=xa(:)'; ya=ya(:)'; za=za(:)'; Fa=Fa(:)';   % central node data in rows


In variants of this where xa, ya, za, and Fa were not input arguments, data
at the central node would be: xc=0.125*sum(x(NN)), yc=0.125*sum(y(NN)), and
zc=0.125sum(z(NN)), with Fc=Fof(xc,yc,zc) (function style) or (mean value)
Fc=0.125*sum(F(NN)).  Incidentally, I would use 'c' instead of 'a' if I were
to rewrite this.  And don't forget to check for Fc==0.


N=find(abs(sum(sign([Fa;F(NN)])))<8);      % this hack is worth doing 


If the values of F at the nine nodes of a bcc-cell all have the same sign
(never zero) F=0 can't be there.


NN=NN(:,N); xa=xa(N); ya=ya(N); za=za(N); Fa=Fa(N);


Keep only the cols of NN and elements of xa,ya,za,Fa that matter.

      
% define WESNBA pyramids (rows of NN for bases)

P={[5 7 3 1],[2 4 8 6],[1 2 6 5],[3 7 8 4],[1 3 4 2],[5 6 8 7]};  % a cell array

  {  south  ,  north  ,   west  ,  east   ,  below  ,  above  }

The orientation of nodes on faces is right-handed relative to the inward normal.
pyras() uses that when it encodes info about which edges have sign changes. Now
this could be important: any vertex of a tetrahedron can be called the apex
(in the first row), and tetras() expects the other three to be in right-handed 
order, relative to the normal that points toward the apex.  Feeding pyras() and 
tetras() mixed right- and left-handed bases will probably result in trash, but 
maybe not. (Sometime when I am back in Providence, I'll have a look to see how 
hard this can be kicked before it breaks.)


tri=[]; tet=[]; pen=[];                   % let pyras() do all the work


Start with nothing.


for k=1:6
   N=NN(P{k},:);
   [p3,p4,p5]=pyras([xa;x(N)],[ya;y(N)],[za;z(N)],[Fa;F(N)]);
   tri=[tri,p3]; tet=[tet,p4]; pen=[pen,p5];
end


Append the cols of 9, 12 and 15 coordinates at the vertices of patches.