function out=element_tetra4(task,nodes,coef)
% common stuff
% reference tetra: x>0 sum(x)<1
% nodes coordinates:
%tet_ref_nodes_xyz=[0 1 0 0 
%                   0 0 1 0
%                   0 0 0 1];
% basis functions
% 1-sum(x), x(1),x(2),x(3)
f=zeros(2,2,2);
f1=f; f1(2,1,1)=-1;f1(1,2,1)=-1;f1(1,1,2)=-1;f1(2,2,2)=1;
f2=f; f2(1,2,2)=1;
f3=f; f3(2,1,2)=1;
f4=f; f4(2,2,1)=1;
fun_tetra={f1,f2,f3,f4};
for i=1:3, % isoparametric mapping
    map_tetra{i}=node(i,1)*f1+node(i,2)*f2+node(i,3)*f3+node(i,4)*f4;
end
%lambda=1000; mu=1;
%c_e3d = elast_coeff(lambda,mu,3);
switch task
    case 'stiff'
        % compute local stiffness matrix
        % 1 point gaussian enough, grads of shape functions constant
        K=stiff_pol(c_e3d,fun_tetra,map_tetra,[1;1;1]/4,1/6);
        out=reshape(K,[12,12]);
    otherwise
        error('bad function')
end