function [N,D,E,R,Q] = driver_getfem_mod (h_getfem, nrefmax, h_utrue);
% inputs:
% h_getfem = string with function name, for instance, h_getfem = 'getfem_cos':
if (nargin < 1)
  h_getfem = 'getfem';
end;
% nrefmax = maximum number of refinements, for instantce, nrefmax = 6:
if (nargin < 2)
  nrefmax = 6;
end;
% h_utrue = string with expression for true solution as function of the
% independent variables, for instance,
%   h_utrue = '(cos(pi*x/2))^2*(cos(pi*y/2))^2'
if (nargin < 3)
  h_utrue = [];
end;

% reference solution on finest mesh:
fem_ref = feval(h_getfem,nrefmax);

if (isempty(h_utrue)) % if no true solution available, then
  % obtain mass matrix for finest mesh:
  fema = fem_ref;
  fema.equ.c = 0;
  fema.equ.a = 1;
  fema.xmesh = meshextend (fema);
  [Mass,L,M,N] = assemble (fema);
  p_ref = fema.mesh.p;
  clear fema L M N; % clear to save memory

  % interpolate ref. sol. to ref. mesh:
  u_ref = postinterp(fem_ref, 'u', p_ref);
end;

% preallocate vectors:
N = repmat (NaN, [nrefmax 1]);
D = repmat (NaN, [nrefmax 1]);
E = repmat (NaN, [nrefmax 1]);
R = repmat (NaN, [nrefmax 1]);
Q = repmat (NaN, [nrefmax 1]);

clear fem_ref; % clear to save memory

for nref = nrefmax-1 : -1 : 0
  fem = feval(h_getfem,nref);
  N(nref+1) = size(fem.mesh.e,2);
  D(nref+1) = length(fem.sol.u);

  if (isempty(h_utrue)) % if no true solution available, then
    % interpolate sol. to ref. mesh:
    u_int = postinterp(fem, 'u', p_ref);

    % compute error as a column vector:
    e = u_ref(:) - u_int(:);

    % compute L2-norm of nodal error:
    E(nref+1) = sqrt(e'*Mass*e);
  else
    str_expression = ['(u-(',h_utrue,'))^2'];
    I=postint(fem,str_expression, ...
           'unit','', ...
           'dl',[1]);
    E(nref+1) = sqrt(I);
  end;
end;

R(2:end) = E(1:nrefmax-1) ./ E(2:nrefmax);
Q(2:end) = log2(R(2:end));

fprintf(' r       N       D            E            R            Q\n');
for nref = 0 : nrefmax-1
  fprintf('%2d %7d %7d %12.4e %12.4e %12.4e\n', ...
          nref, N(nref+1), D(nref+1), E(nref+1), R(nref+1), Q(nref+1));
end;