function [deldotdel, delphi] = rbdeldotdel(cfg)
%
% [deldotdel, delphi]=rbdeldotdel(cfg)
%
% Compute deldotdel=<grad(phi_i).grad(phi_j)>, where <> means spatial integration
% inside elements, "." means dot-product, grad() means gradience, phi means
% linear basis function in a tetrahedron. For a linear function phi,
% grad(phi) is a constant across the element.
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     cfg: the redbird data structure
%
% output:
%     deldotdel: precomputed deldotdel=<grad(phi_i).grad(phi_j)>. For each
%         element, deldotdel is a 4x4 symmetric matrix - to store this data
%         efficiently; we only store the upper triangule of the matrix per
%         element - this gives a Ne x 10 matrix (Ne is the number of tets);
%         for each row: [d11, d12, d13, d14, d22, d23, d24, d33, d34, d44]
%     delphi: gradient of the basis functions (grad(phi_i)) in each
%         tetrahedral element; a 3x4 matrix for each element with a total
%         dimension is 3 x 4 x Ne, where
%            3 - gradient direction, for x, y and z
%            4 - which basis functions - for node 1, 2, 3 and 4
%            Ne - number of element
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details
%
% -- this function is part of Redbird-m toolbox
%

no = cfg.node;
el = cfg.elem(:, 1:4);

no = reshape(no(el', :)', 3, 4, size(el, 1));

delphi = zeros(3, 4, size(el, 1));

col = [4 2 3 2
       3 1 4 3
       2 4 1 4
       1 3 2 1];

for coord = 1:3
    idx = 1:3;
    idx(coord) = [];
    for i = 1:4
        delphi(coord, i, :) = squeeze(((no(idx(1), col(i, 1), :) - no(idx(1), col(i, 2), :)) .* (no(idx(2), col(i, 3), :) - no(idx(2), col(i, 4), :)) - ...
                                       (no(idx(1), col(i, 3), :) - no(idx(1), col(i, 4), :)) .* (no(idx(2), col(i, 1), :) - no(idx(2), col(i, 2), :)))) ./ (cfg.evol(:) * 6);
    end
end

deldotdel = zeros(size(el, 1), 10);
count = 1;

for i = 1:4
    for j = i:4
        deldotdel(:, count) = sum(squeeze(delphi(:, i, :) .* delphi(:, j, :)), 1);
        count = count + 1;
    end
end
deldotdel = deldotdel .* repmat(cfg.evol(:), 1, 10);
