Using the idea from Daniel Lichtblau here to use PowersRepresentations, we can go much faster:
Clear["Global`*"]; AbsoluteTiming[ center = {1, 3, 5}; r = 5 Sqrt[3]; nnvals = PowersRepresentations[r^2, 3, 2]; permVals = Flatten[Permutations /@ nnvals, 1]; signs = Tuples[{-1, 1}, {3}]; alltriples = Union[Flatten[Outer[Times, signs, permVals, 1], 1]]; alltriples = # + center & /@ alltriples; sideLength = r/(Sqrt[3/2]/2); distances = Outer[EuclideanDistance, alltriples, alltriples, 1]; tetras = Flatten[Table[ candidateVertices = Join[{i}, Flatten[Position[distances[[i]], n_ /; n == sideLength]]]; indices = Subsets[candidateVertices, {4}]; Part[alltriples, #] & /@ indices , {i, distances // Length}], 1]; sortedTetras = DeleteDuplicates[Sort /@ tetras]; allSameSidePositions = Flatten[Position[ Table[Tally[ DeleteCases[Flatten[Outer[EuclideanDistance, i, i, 1]], 0]][[All, 1]], {i, sortedTetras}], n_ /; n == {sideLength}]]; tetrahedrons = sortedTetras[[allSameSidePositions]]; ] Graphics3D[ConvexHullRegion /@ tetrahedrons] (*{0.028064, Null}*) And with r= = 33 Sqrt[3]
Clear["Global`*"]; AbsoluteTiming[ center = {1, 3, 5}; r = 33 Sqrt[3]; nnvals = PowersRepresentations[r^2, 3, 2]; permVals = Flatten[Permutations /@ nnvals, 1]; signs = Tuples[{-1, 1}, {3}]; alltriples = Union[Flatten[Outer[Times, signs, permVals, 1], 1]]; alltriples = # + center & /@ alltriples; sideLength = r/(Sqrt[3/2]/2); distances = Outer[EuclideanDistance, alltriples, alltriples, 1]; tetras = Flatten[Table[ candidateVertices = Join[{i}, Flatten[Position[distances[[i]], n_ /; n == sideLength]]]; indices = Subsets[candidateVertices, {4}]; Part[alltriples, #] & /@ indices , {i, distances // Length}], 1]; sortedTetras = DeleteDuplicates[Sort /@ tetras]; allSameSidePositions = Flatten[Position[ Table[Tally[ DeleteCases[Flatten[Outer[EuclideanDistance, i, i, 1]], 0]][[All, 1]], {i, sortedTetras}], n_ /; n == {sideLength}]]; tetrahedrons = sortedTetras[[allSameSidePositions]]; ] Graphics3D[ConvexHullRegion /@ tetrahedrons] (*{2.0544, Null}*) My method for deleting the candidate tetrahedrons that are not all the same side length is very ugly and can probably be greatly improved for speed and clarity:
allSameSidePositions = Flatten[Position[ Table[Tally[ DeleteCases[Flatten[Outer[EuclideanDistance, i, i, 1]], 0]][[All, 1]], {i, sortedTetras}], n_ /; n == {sideLength}]]; tetrahedrons = sortedTetras[[allSameSidePositions]]; But this produces the same output as @cvgmt and @kirma 's answers

