You have this same idiom in many places:
    std::vector<double> volume_fractions(this->n_compositional_fields() + 1);
    compute_volume_fractions( composition, volume_fractions);
How about you let compute_volume_fractions() return the volume_fraction as a returned value, rather than through the second argument? Then you could just change this to 
    std::vector<double> volume_fractions = compute_volume_fractions( composition, volume_fractions);
which presumably would also be more efficient.