function dtscale = dtscale3D()

% compute inscribed circle diameter as characteristic for grid to choose timestep

Globals3D;

vmask1 = find(abs(r + s + t + 3) < NODETOL);
vmask2 = find(abs(r - 1) < NODETOL);
vmask3 = find(abs(s - 1) < NODETOL);
vmask4 = find(abs(t - 1) < NODETOL);
vmask = [vmask1; vmask2; vmask3; vmask4];
vx = x(vmask(:),:); vy = y(vmask(:),:); vz = z(vmask(:),:);

vectors1 = [vx(2,:) - vx(1,:); vy(2,:) - vy(1,:); vz(2,:) - vz(1,:)];
vectors2 = [vx(3,:) - vx(1,:); vy(3,:) - vy(1,:); vz(3,:) - vz(1,:)];
vectors3 = [vx(4,:) - vx(1,:); vy(4,:) - vy(1,:); vz(4,:) - vz(1,:)];
volumes = abs(dot(vectors1,cross(vectors2,vectors3,1),1)) / 6;

face_def = [1, 2, 3; 1, 2, 4; 2, 3, 4; 1, 3, 4];
areas = zeros(N_faces,K);
for i = 1:N_faces
    face = face_def(i,:);
    vectors1 = [vx(face(2),:) - vx(face(1),:); vy(face(2),:) - vy(face(1),:); vz(face(2),:) - vz(face(1),:)];
    vectors2 = [vx(face(3),:) - vx(face(1),:); vy(face(3),:) - vy(face(1),:); vz(face(3),:) - vz(face(1),:)];
    cross_products = cross(vectors1,vectors2,1);
    norms = sqrt(cross_products(1,:) .^ 2 + cross_products(2,:) .^ 2 + cross_products(3,:) .^ 2);
    areas(i,:) = 0.5 * norms;
end

dtscale = 6 * volumes ./ sum(areas,1);

end