Add Capability of Newton Raphson Functions to deal with History Variables

For several extensions to be added to the library in the future, history variables are necessary (e.g. accumulated plastic strain, accumulated fatigue effect, ...). Therefore, the Newton-Raphson functions phase_field.fem.solver.pf.newton_raphson and phase_field.fem.solver.equilibrium.newton_raphson need the capability to deal with this. Since the various extensions of the library might have varying number of history variables, this must be done in a general way.

To be precise, this means that (maybe a struct of) history variables at the Gauss points must be added to the input as well as to the output, similar to the currently implemented history variable for the phase-field (optimally, this should be integrated into the new history variable struct).

This might look like:

function [p_field, res_0, res, history_vars_new] = newton_raphson(...
    sys, irrev,...
    strain_en_undgr,...
    p_field, p_field_old, label_el,...
    penalty_irrev, penalty_recov,...
    history_vars_old,...
    args...
)

Accordingly, the history variables need to be passed into the Fortan routines. For the field problems which do not need history variables, zeros need to be passed in.