vecgrad:=proc()
 local i,v,q,co,h,vq,hq,hinv,dia,vdia,ddia;
 description ` Vektorgradient:  v dyadisch nabla`;
 v:=args[1];
 q:=args[2];
 h:=[1,1,1];
 if nargs>2 then 
   co:=args[3];
   if type(co,`=`)then
     if rhs(co)=cylindrical then
       h:=[1,q[1],1] fi;
     if rhs(co)=spherical then
       h:=[1,q[1],q[1]*sin(q[2])] fi
   else h:=co
   fi
 fi;
 vdia:=array(sparse,1..3,1..3):
 for i to 3 do vdia[i,i]:=v[i] od;
 vq:=jacobian(v,q);
 hq:=jacobian(h,q);
 hinv:=array(sparse,1..3,1..3):
 for i to 3 do hinv[i,i]:=1/h[i] od;
 dia:=evalm(hq&*hinv&*v);
 ddia:=array(sparse,1..3,1..3):
 for i to 3 do ddia[i,i]:=dia[i] od;
 evalm((vq-hinv&*transpose(hq)&*vdia
 +ddia)&*hinv);
 end: 
 tensdiv:=proc()
 local i,l,t,q,co,h,hinv,hq,vol,tz,tzeile,di;
 description ` Tensordivergenz:  t punkt nabla`;
 t:=args[1];
 q:=args[2];
 h:=[1,1,1];
 if nargs>2 then 
   co:=args[3];
   if type(co,`=`)then
     if rhs(co)=cylindrical then
       h:=[1,q[1],1] fi;
     if rhs(co)=spherical then
       h:=[1,q[1],q[1]*sin(q[2])] fi
   else h:=co
   fi
 fi;
 hinv:=array(sparse,1..3,1..3):
 for i to 3 do hinv[i,i]:=1/h[i] od;
 vol:=h[1]*h[2]*h[3];
 tz:=evalm(t&*hinv*vol);
 hq:=jacobian(h,q);
 for i to 3 do
     di[i]:=diverge(row(tz,i),q)/vol +1/h[i]*
 add((t[l,i]*hq[i,l]-t[l,l]*hq[l,i])/h[l],
 l=1..3)
 od;
 vector([di[1],di[2],di[3]])
 end: