function [values,grad]=FEMValueM(xy,nodes,elem,u,defaultvalue);
% ... = FEMValueM(...)
% evaluate the solution at a given set of points
%
% values=FEMValueM(xy,nodes,elem,u,defaultvalue)
% values=FEMValueM(xy,nodes,elem,u)
%
% [values,grad]=FEMValue(xy,nodes,elem,u,defaultvalue)
% [values,grad]=FEMValue(xy,nodes,elem,u)
% 
% nodes elem u  information about the mesh and the solution
% defaultvalue  value returned if a point is not in the domain
%               if not specified NaN is returned
% xy  coordinates of the points at which the solution is to be evaluated
%
% values  values of the solution at the points given by xy
% grad    values of the gradient of the solution at the points given by xy

% this function runs VERY slow
% algorithm and implementation have to be improved

if (nargin<4) usage("values=FEMValueM(xy,nodes,elem,u,defaultvalue)");endif
if (nargin<5) defaultvalue=NaN;endif
gradient=false;
if (nargout>=2) gradient=true;endif

[nx,dummy]=size(xy);
[nelem,dummy]=size(elem);

xMin=zeros(nelem,1);  xMax=zeros(nelem,1);
yMin=zeros(nelem,1);  yMax=zeros(nelem,1);
for k=1:nelem
  xMin(k)=min([nodes(elem(k,1),1),nodes(elem(k,2),1),nodes(elem(k,3),1)]);
  xMax(k)=max([nodes(elem(k,1),1),nodes(elem(k,2),1),nodes(elem(k,3),1)]);
  yMin(k)=min([nodes(elem(k,1),2),nodes(elem(k,2),2),nodes(elem(k,3),2)]);
  yMax(k)=max([nodes(elem(k,1),2),nodes(elem(k,2),2),nodes(elem(k,3),2)]);
endfor

values=ones(nx,1)*defaultvalue;
grad  =ones(nx,2)*defaultvalue;
for kxy=1:nx
  for k=1:nelem  # go through all elements
    flag=false;  % set it to true if point was in the element
    if ((xy(kxy,1)>=xMin(k)-eps)&&(xy(kxy,1)<=xMax(k)+eps)&&...
        (xy(kxy,2)>=yMin(k)-eps)&&(xy(kxy,2)<=yMax(k)+eps))  
                        # examine this element closely
      a=nodes(elem(k,2),1:2)-nodes(elem(k,1),1:2);
      b=nodes(elem(k,3),1:2)-nodes(elem(k,1),1:2);
      dx=xy(kxy,:)-nodes(elem(k,1),1:2);
      coeff=dx/[a;b];
      if (coeff(1)>=-eps&&coeff(2)>=-eps&&(coeff(1)+coeff(2)<=1+eps))
        flag=true;
        break;
      endif;
    endif; # examine this element closely
  endfor # k
  if(flag)
    u0=u(elem(k,1));
    values(kxy)=u0+coeff*[u(elem(k,2))-u0;u(elem(k,3))-u0];
    if (gradient)
      u0=u(elem(k,1));
      c=cross([a,u(elem(k,2))-u0],[b,u(elem(k,3))-u0]);
      grad(kxy,1)=-c(1)/c(3);
      grad(kxy,2)=-c(2)/c(3);
    endif
  else values(kxy)=defaultvalue;
  endif # flag
endfor # kxy
endfunction

