function A=dudvElem(vol,G,r,s)
  % vol=simplex.volume(q); % volume of d-simplex 
  % G=simplex.GradBaCo(q); % gradient of barycentric coordinates of d-simplex (d+1)-by-n array
  
  A=vol*G(:,s)*G(:,r)';
end
