function b=femval_time(x,y,expr,T,femstr)
%FEMVAL_TIME
%   B = FEMVAL_TIME(X,Y,EXPR,T,FEMSTR) returns the values B of
%   the expression EXPR for the FEM-structure with the name FEMSTR
%   at points specified by X and Y. The FEM-structure with the name
%   FEMSTR must be present in the main MATLAB work space and have
%   a time-dependent solution, that is, a field sol.u and a field
%   sol.tlist. 
%
%   The input T is the time at which the values are evaluated.
%   T does not need to be a time-step in sol.tlist, because (linear)
%   interpolation between the solutions in sol.u is used.
%
%   Note also that the mesh for the current model using FEMVAL_TIME does not
%   need to be the same as the mesh in the FEM-structure with the name
%   FEMSTR. Values of B for coordinates that are outside the mesh
%   will be NaN.
%
%   Example:
%   Given a FEM-structure fem1 in the main MATLAB work space defining
%   the variable ux on subdomain level, and containing a time dependent
%   solution for the time steps 0:0.1:1. Another time-dependent model
%   in the GUI is loaded, possibly with a different mesh and different time
%   steps. In this model, a PDE coefficient, boundary coefficient, or
%   initial value can be set to 'femval_time(x,y,'ux',t,'fem1')'. An expression
%   including femval_time can also be used directly in the expressions in the
%   post mode dialog boxes.
%
%   Cautionary:
%   FEMVAL_STAT and FEMVAL_TIME take x and y as input and uses postinterp
%   on these points. This means that for coordinates being mesh nodes
%   between elements with a discontinuity between them, postinterp will
%   just "choose" the value from one of the elements. This can be a
%   problem when trying to plot discontinuities using these functions in,
%   e.g., POSTPLOT, or entering them in the Plot Parameters dialog box.
%   However, when integrating with POSTINT, or using the Integration dialog
%   boxes, evaluation is only done in internal points, so no such coordinates
%   will be passed to FEMVAL_STAT and FEMVAL_TIME, hence this is not a problem.
%
%   See also FEMVAL_STAT.

oldfem=evalin('base',femstr);
u=interp1(oldfem.sol.tlist,oldfem.sol.u',T)';
b=postinterp(oldfem,expr,[x(:)';y(:)'],'u',u,'T',T);
b=reshape(b,size(x));