Finding the arithmetic mean of the nodal solution is PDE toolbox

Hi,

I am simulating a nonlinear system of PDE using PDE toolbox. I am using 2D rectangle geometry and I have to find the arithmetic mean of the nodal solution of the previous nodes and substitute it in next node for calculating the film thickness. Is there any possible way to calculate it?

model = createpde;
%% Bearing Parameters
epsi    = 1.6962;

cl      = 31.8*10^-6;                 % Lubrication gap clearance (m)
N       = 55500;                    % Rotations per minute (rpm)

l       = 38.1*10^-3;               % Axial bearing length (m)
r       = 19.05*10^-3;              % Bearing radius (m)
p_a     = 1*10^5;                   % Ambient Pressure (Pa)
mu      = 1.836 * 10^-5;            % Dynamic Viscosity (Ns/m^2)
Omg     = (2*pi*N)/60;              % Angular speed of Journal (Rad/s)
rc      = (r/cl)^2;
A       = (6*mu*Omg*rc / p_a);      % Bearing number
%% Foil Parameters
t_T     = 0.2032 * 10^-3;           % Thickness of top foil (m)
t_B     = 0.1016 * 10^-3;           % Thickness of bump foil (m)
E       = 2.07*10^11;                 % Young's modulus of top and bump foil (Pa)
nu      = 0.3;                     % Poisson's ratio
s_b     = 4.572*10^-3;              % Bump foil pitch (m)
l_o     = 1.778*10^-3;              % Bump half length (m)

%% Model Geometry
rec     = [3,4,0, 2*pi, 2*pi, 0, l/r, l/r, 0, 0]';
g = decsg(rec);
geometryFromEdges(model,g);

%% Mesh
generateMesh(model,'GeometricOrder','linear','Hmax',0.1);

%% Assigning the Boundary conditions with all edges equal to ambient pressure (p_a)
applyBoundaryCondition(model,'dirichlet','Edge',1:model.Geometry.NumEdges,'u',1);

%% Assigning the Initial condition as pressure equal to ambient pressure initially

setInitialConditions(model,1.785);


%% Assign the coefficients of the Non-Linear 2-D Reynolds equation for compressible fluids to generic PDE eq

al = (2*p_a*s_b/(cl*E))*((l_o/t_B)^3)*(1-nu^2);                                         % foil flexibility per unit area

h =@(location,state) 1 + epsi.*cos(location.x) + al.*(state.u - 1);                                 % film thickness

h_grad =@(location,state) -epsi.*sin(location.x) + al.*state.ux;                   % Gradient of foil thickness 


f =@(location,state) A.*( h(location,state).* state.ux + state.u .* h_grad(location,state) ); % 1st term in the Reynold eq's rhs                            % 2st term in the Reynold eq's rhs

c = @(location,state) -( state.u.* h(location,state).^3 ); 
                                      
specifyCoefficients(model,'m',0,...
                          'd',0,...
                          'c',c,...
                          'a',0,...
                          'f',f);

%% Solve PDE
R = solvepde(model);