This answer extends @user21's to include different mesh densities of the inclusion along the X, Y, and Z directions.

It is important to note that the current mesher (version 12.1.1) likes to produce an isotropic mesh. One can accomplish the different mesh densities by creating a parameterized (I, J, K) structured mesh that ranges between zero and the number of elements in each direction. Then, one can rescale the coordinates from I, J, K space to user scaled coordinates.

First, let's create an isotropic structured mesh:

nx = 10; ny = 40; nz = 5;
isoMesh =
   ToElementMesh[Cuboid[{0, 0, 0}, {nx, ny, nz}], 
     "MeshOrder" -> 1, MaxCellMeasure -> 1,
     "RegionMarker" -> {{{nx, ny, nz}/2, 1}}, 
     "MeshElementType" -> TetrahedronElement];
isoMesh["Wireframe"]

Second, let's create a rescaling transformation function from I, J, K space to user scaled coordinates:

scaledToUser = 
  RescalingTransform[{{0, nx}, {0, ny}, {0, nz}}, {{-xI/2, 
     xI/2}, {-yI/2, yI/2}, {-zI/2, zI/2}}];

Now, we can create the inner mesh by simply rescaling the coordinates like so:

innerMesh = 
  ToElementMesh[
   "Coordinates" -> scaledToUser /@ isoMesh["Coordinates"], 
   "MeshElements" -> isoMesh["MeshElements"]];
innerMesh["Wireframe"]

Now, simply follow@user21's workflow with the new definition of innermesh and you can achieve different mesh densities along the X, Y, Z directions.

finalMesh[
  "Wireframe"["MeshElement" -> "MeshElements", 
    "MeshElementStyle" -> (Directive[FaceForm[#], 
              EdgeForm[]] &  /@ {Orange, Blue}), 
    PlotRange -> {All, All, {-zM, zI/2}}]]
finalMesh[
  "Wireframe"["MeshElement" -> "MeshElements", 
    "MeshElementStyle" -> (Directive[FaceForm[#], 
              EdgeForm[]] &  /@ {Orange, Blue}), 
    PlotRange -> {All, {0, yI/2}, {-zM, zI/2}}]]

Structured hex mesh approach

As I alluded to in the comments, if you want to use a structured hex mesh for the inclusion, you probably want to propagate that through the entire mesh because the current version of Mathematica does not support pyramid and wedge type elements in 3D.

Depending on the nature of the physics you're trying to solve, there often can be sharp gradients at the interface regions. In this case, your solutions can often benefit by having a boundary layer mesh (or anisotropic mesh) where you have fine element layers at the interface that grow exponentially into the domain. These types of meshes can be quite economical in terms of element count.

Workflow

Helper functions

First, we will define some helper functions to create an anisotropic mesh.

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(* Define Some Helper Functions For Structured Quad Mesh*)
pointsToMesh[data_] :=
  MeshRegion[Transpose[{data}], 
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] := 
 Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] := 
 Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] := 
 Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] := 
 Quiet@Abs@
   FindRoot[firstElmHeight[x0, xf, n, r] - fElm, {r, 1.0001, 100000}, 
     Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] := 
 N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := Module[{seg},
  seg = leftSegmentGrowth[len, n, fElm];
  flipSegment[seg]
  ]
reflectRight[pts_] := With[{rt = ReflectionTransform[{1}, {Last@pts}]},
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] := 
 With[{rt = ReflectionTransform[{-1}, {First@pts}]},
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]

Construct a tensor product grid using RegionProduct tensor product mesh

Now, we can glue a bunch of segments that have different refinement strategies together in the horizontal, vertical, and depth directions as shown in the following workflow.

(*Define parameters*)
(*Lengths*)
h = 100;(*Horizontal*)
v = 10;(*Vertical*)
d = h;(*Depth*)
(*Number of elements per segment*)
nh = 10;
nv = 10;
nd = 10;
(*Association for Clearer Region Assignment*)
reg = <|"main" -> 1, "incl" -> 2|>;
(*Create mesh segments*)
(*Horizontal segments*)
(* left segment *)
(*First element is 1/50th of seg length*)
sh = rightSegmentGrowth[h, nh, h/50];
Print["Horizontal and depth segments"]
rh = pointsToMesh@(reflectRight@reflectRight[sh] - 2 h)
(*Vertical segment*)
sv = rightSegmentGrowth[v, nv, v/50];
Print["Vertical segment"]
rv = pointsToMesh@(reflectRight@reflectRight[sv] - 2 v)
(*Create tensor product grid with RegionProduct*)
rp = RegionProduct[rh, rv, rh];
(*Show the mesh*)
Print["Highlighted RegionProduct mesh"]
HighlightMesh[rp, Style[1, Orange]]

Convert MeshRegion to ElementMesh with region markers

(*Extract Coords from RegionProduct*)
crd = MeshCoordinates[rp];
(*grab hexa element incidents RegionProduct mesh*)
inc = Delete[0] /@ MeshCells[rp, 3];
mesh = ToElementMesh["Coordinates" -> crd, 
   "MeshElements" -> {HexahedronElement[inc]}];
(*Extract bmesh*)
bmesh = ToBoundaryMesh[mesh];
(*Inclusion RegionMember Function*)
Ω3Dinclusion = Cuboid[{-h, -v, -h}, {h, v, h}];
rmf = RegionMember[Ω3Dinclusion];
regmarkerfn = If[rmf[#], reg["main"], reg["incl"]] &;
(*Get mean coordinate of each hexa for region marker assignment*)
mean = Mean /@ GetElementCoordinates[mesh["Coordinates"], #] & /@ 
    ElementIncidents[mesh["MeshElements"]] // First;
regmarkers = regmarkerfn /@ mean;
(*Create and view element mesh*)
Print["Converted Hexa Element Mesh Cutaway Drawing"]
mesh = ToElementMesh["Coordinates" -> mesh["Coordinates"], 
   "MeshElements" -> {HexahedronElement[inc, regmarkers]}];
mesh[
  "Wireframe"["MeshElement" -> "MeshElements", 
    "MeshElementStyle" -> (Directive[Opacity[0.5], FaceForm[#](*, 
              EdgeForm[]*)] &  /@ {Blue, Orange}), 
  ViewPoint -> {-1.5, 0.8, -3}, ViewVertical -> {0, 1, 0}, 
    PlotRange -> {{0, 2 h}, {0, 2 v}, {0, 2 h}}]]

Using a fully structured hex mesh, we created a fairly economical mesh (46656 hex elements) with very fine refinement at the interface.