function m = calculateOutputs( m )
%m = calculateOutputs( m )
%   Calculate all of the tensors for specified/actual/residual growth on
%   the A/B sides, and the rotations.  This creates the following
%   components of m:
%
%     m.outputs.specifiedstrain.A
%     m.outputs.specifiedstrain.B
%     m.outputs.actualstrain.A
%     m.outputs.actualstrain.B
%     m.outputs.residualstrain.A
%     m.outputs.residualstrain.B
%     m.outputs.rotation
%
%   Each of these is a per-FE quantity.  If there are N finite elements,
%   then the tensor fields are N*6 and the rotations field is N*3.  These
%   are all rates of growth or rotation, not absolute amounts.
%
%   Take the average of an A and a B quantity to get the value for the
%   midplane; take half of B-A to get the bending component.
%
%   To convert tensors to principal components, write for example:
%       [amounts,frames] = tensorsToComponents( m.outputs.actualstrain.A );
%   amounts will then be the principal components of actual growth on the
%   A side, listed in descending order.  If you want them listed in the
%   order parallel, perpendicular, and normal, replace the second line by:
%       [amounts,frames] = tensorsToComponents( m.outputs.actualstrain.A, m.cellFrames );
%
%   To resolve rotation vectors into rotation about the normal vector and
%   the remainder, write:
%
%       [inplane,outofplane] = splitVector( m.outputs.rotations, m.unitcellnormals );
%
%   inplane will be a column of scalars, the rotation rates around the
%   normal vectors, and outofplane will be an N*3 matrix of vectors in the
%   planes of the respective finite elements.
%
%   To convert any per-element quantity to a per-vertex quantity, call:
%
%       perVxQuantity = perFEtoperVertex( m, perFEquantity );
%
%   The per-element quantity can be an N*K matrix, where N is the number of
%   finite elements.  perVxQuantity will then be M*K where M is the number
%   of vertices.

    if isfield( m.celldata, 'Gglobal' )
        allgrowthtensors = permute( ...
                               reshape( [m.celldata.Gglobal],6,6,[]), ...
                               [2 1 3] ...
                           );
    else
        % When using externally supplied growth tensors.
        allgrowthtensors = reshape( repmat( [m.celldata.cellThermExpGlobalTensor], 6, 1 ), 6, 6, [] );
    end
    m.outputs.specifiedstrain.A = permute( sum( allgrowthtensors(:,1:3,:), 2 )/3, [3 1 2] ); 
    m.outputs.specifiedstrain.B = permute( sum( allgrowthtensors(:,4:6,:), 2 )/3, [3 1 2] ); 
    [m.outputs.actualstrain.A,m.outputs.actualstrain.B] = ...
        split6strains( 'displacementStrain' );
    [m.outputs.residualstrain.A,m.outputs.residualstrain.B] = ...
        split6strains( 'residualStrain' );
    m.outputs.rotations = getRotations( m, 'total' );

function [A,B] = split6strains( f )
    alltensors = reshape( [m.celldata.(f)], 6, 6, [] )/m.globalProps.timestep;
    A = permute( sum( alltensors(:,1:3,:), 2 )/3, [3 1 2] ); 
    B = permute( sum( alltensors(:,4:6,:), 2 )/3, [3 1 2] ); 
end
end

