I am a bit late for the party but it can be done in Mathematica. Here is how:

Get the file and find the bounds:

file = URLDownload[
   "http://www.aviduratas.de/wc/model-scaled-x7.stl"];
myQuartic = Import[file, "BoundaryMeshRegion"];
RegionBounds[myQuartic]

(*{{-48.5411, 48.4823}, {-49.0028, 49.0095}, {-49.0087, 49.0065}}*)

Create a set of objects - use something more intelligent here:

s = 20;
balls = RegionUnion[
   Flatten[Table[
     Ball[{i, j, k}, 3], {i, -50, 50, s}, {j, -50, 50, s}, {k, -50, 
      50, s}]]];

To show case this, we will want to find the intersections of these:

Show[myQuartic, Region[balls]]

Load and use OpenCascadeLink:

Needs["OpenCascadeLink`"]
s1 = OpenCascadeShapeImport[file];
s2 = OpenCascadeShape[balls];
s3 = OpenCascadeShapeDifference[s1, s2];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[s3];

Visualize:

bmesh["Wireframe"[
  "MeshElementStyle" -> Directive[FaceForm[LightBlue], EdgeForm[]]]]