function cmat = elasticityC3D(E,nu)
% ELASTICITYC3D Calculate c-matrix for isotropic elastic solid
%   This convenience function is included in the toolbox for the use in
%   several examples. This undocumented function may be removed in a future
%   release.
%
%   This function uses the linearized small-displacement assumption for an
%   isotropic material. Hence, there are only two independent elasticity
%   constants, Young's modulus (E) and Poisson's ratio (nu), which are used
%   to construct c-matrix in PDE Toolbox format.

%   Copyright 2014-2016 The MathWorks, Inc.
if ~isa(E,'function_handle') && ~isa(nu,'function_handle') % Retain original workflow for numeric E and nu
    G = E/(2*(1+nu));
    c1 = E*(1-nu)/((1+nu)*(1-2*nu));
    c12 = c1*nu/(1-nu);
    C11 = [c1 0 G 0 0 G];
    C12 = [0 G 0 c12 0 0 0 0 0];
    C22 = [G 0 c1 0 0 G];
    C13 = [0 0 G 0 0 0 c12 0 0];
    C23 = [0 0 0 0 0 G 0 c12 0];
    C33 = [G 0 G 0 0 c1];
    cmat  = [C11 C12 C22 C13 C23 C33]';
else % E or nu or both are function handles
    if ~isa(E,'function_handle') % Convert numeric E to function handle
        E = @(l,s) E.*ones(1,length(l.x));
    end

    if ~isa(nu,'function_handle') % Convert numeric nu to function handle
        nu = @(l,s) nu.*ones(1,length(l.x));
    end

    cmat = @(l,s) pde.internal.FESolverUtilities.elasticityFcnHandleC3D(l,s,E,nu); %wrapper function to construct c-matrix as function handle
end

end
