Filling, Printout3D with RegionPlot3D

I modified your code slightly, setting BoxRatios to Automatic so the on-screen display matches the true dimensions, and naming the plot ms.

 ms = Block[{f = 3, r = 6, base = 5, l = 15, c = 9, holes, ms}, holes = Table[((x - c Cos[a])^2 + (y - c Sin[a])^2) > 9, {a, 0, 5 Pi/3, Pi/3}] /. List -> And; RegionPlot3D[ z < base + (x^2 + y^2)/(4 f) && holes, {x, -l, l}, {y, -l, l}, {z, 0, base + r^2/(4 f)}, BoxRatios -> Automatic] ] 

First discretize your region, producing a MeshRegion object that comprises the boundary surface. It is 2-dimensional, as you stated, since it is just the boundary surface of the solid region.

dms = DiscretizeGraphics[ms] 

Using the same data (the vertices and face indices in the MeshRegion), change the head to BoundaryMeshRegion to get the solid version.

BoundaryMeshRegion[MeshCoordinates[dms], MeshCells[dms, 2]] 

STL only supports faces (as opposed to solids), so I'd be surprised if converting to BoundaryMeshRegion and triangulating fixed things.

Printout3D does some repair with the option Method -> "PerformModelRepair" (I believe this is the default setting). In your case, your model comes back repaired:

 FindMeshDefects[RepairMesh[DiscretizeGraphics[ms]], All, "Cell"] 

<|"FlippedFaces" -> {}, "HoleEdges" -> {}, "TinyFaces" -> {}, "SingularVertices" -> {}, "DanglingEdges" -> {}, "SingularEdges" -> {}, "TinyComponents" -> {}, "TJunctionEdges" -> {}, "IsolatedVertices" -> {}, "OverlappingFaces" -> {}|> 

As for the external software, it's hard to know what's going wrong, but it sounds like something is preventing it from inferring a solid. Which software are you using?

Perhaps a different technique would give a better result.

The highest fidelity method is probably subtracting away the cylinders and paraboliod from the cuboid. We can do this with mesh based boolean operations in 11.2+:

δ = .5; fill = BoundaryDiscretizeGraphics[Cuboid[{-15, -15, 0}, {15, 15, 8}]]; cyls = Table[With[{cx = c Cos[a], cy = c Sin[a]}, Cylinder[{{cx, cy, -1}, {cx, cy, 9}}, 3]], {a, 0, 5 Pi/3, Pi/3}] dc = RegionUnion @@ (BoundaryDiscretizeRegion[#, MaxCellMeasure -> .1δ] & /@ cyls); parab = ImplicitRegion[z > base + (x^2 + y^2)/(4 f), {{x, -l, l}, {y, -l, l}, {z, 0, 8}}]; db = BoundaryDiscretizeRegion[parab, MaxCellMeasure -> δ]; model = RegionDifference[RegionDifference[fill, dc], db] 

Note the weird artifacts in the paraboloid. I think this is just an issue with rendering multi-cell polygons, as Normal[Show[mesh]] renders just fine.

This model is also has no defects:

Values @ FindMeshDefects[RepairMesh[model], All, "Cell"] 

{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}} 

Give this one a try and see if it works.