function D = IsotropicStiffnessMatrix( K, nu, nuNormal )
%D = IsotropicStiffnessMatrix( K, nu, nuNormal )
%   Compute the 6*6 stiffness matrix for an isotropic material with bulk
%   modulus K and Poisson's ratio nu.  If K, nu, and nuNormal are vectors
%   of the same length N > 1, D is returned as an 6*6*N matrix.  nuNormal
%   defaults to nu (and in fact is not used anyway).

    if nargin < 3
        nuNormal = nu;
    end
    
    lambda = 3*K.*nu./(1+nu);  % Lame's first parameter
  % mu = lambda.*(1/(2*nu) - 1);
    mu = 3*K.*(1-2*nu)./(2*(1+nu));  % Shear modulus, Lame's second parameter.
    a = lambda + mu + mu;
    b = lambda;
    c = mu;
    if length(K)==1
        D = [ a b b 0 0 0;
              b a b 0 0 0;
              b b a 0 0 0;
              0 0 0 c 0 0;
              0 0 0 0 c 0;
              0 0 0 0 0 c ];
    else
        a = reshape(a,1,1,[]);
        b = reshape(b,1,1,[]);
        c = reshape(c,1,1,[]);
        z = zeros(size(c));
        D = [ 
                [ [ [a b b];[b a b];[b b a] ], zeros(3,3,length(K)) ];
                [ zeros(3,3,length(K)) [ [c z z];[z c z];[z z c] ] ]
            ];
    end
end
