You can't load in multiple geometries seperately in 3D to my knowledge. From my current experience, the only way to properly manage multiple domains (or cell-zones) in 3D within the PDE toolbox is via the
geometryFromMesh(model,nodes,elements,groupsID)
tool. Here groupsID is a list of values the same length as the number of elements with values (in your case) of 1 and 2 depending on which cell a specific element is in.
The downside of this is you need to have a prexisting mesh made elsewhere rather than using MATLAB inbuilt meshing software. My personal recommendation to do this would be the opensource GMSH software. This can import your STL file and create a mesh (linear, quadratic etc) as well as optimise the mesh via a whole host of algorithms.
Depending on how your PDE system is setup, the "Coherence;" command in GMSH could be really useful - this forces a coherent mesh with no overlapping elements allowing MATLAB to interpret your geometry as a single coherent system with distinct geometrical functions. However, if you're interested in analytically distinct regions this step is unimportant (I note this here because it took a long time for me to figure this bit out).
This principle extends to more domains - my current has 15 geometrical regions enclosed in a solid box. The resullting mesh is simply a solid box but I can still identify and manipulate the various geometrical features inside the box.
