How to mesh a region with inclusions touching the boundary

I cannot take credit for the following code, but it allows one to join several boundary meshes together.

Needs["NDSolve`FEM`"]
(* Code to join multiple boundary meshes *)
ClearAll[validInputQ]
validInputQ[bm1_, bm2_] := 
 BoundaryElementMeshQ[bm1] && 
  BoundaryElementMeshQ[
   bm2] && (bm1["EmbeddingDimension"] === 
    bm2["EmbeddingDimension"]) && (bm1["MeshOrder"] === 
    bm2["MeshOrder"] === 1)
BoundaryElementMeshJoin[bm1_, bm2_, 
   opts : OptionsPattern[ToBoundaryMesh]] /; validInputQ[bm1, bm2] := 
 Module[{c1, c2, nc1, newBCEle, newPEle, eleTypes, markers}, 
  c1 = bm1["Coordinates"];
  c2 = bm2["Coordinates"];
  nc1 = Length[c1];
  newBCEle = bm2["BoundaryElements"];
  eleTypes = Head /@ newBCEle;
  If[ElementMarkersQ[newBCEle], markers = ElementMarkers[newBCEle], 
   markers = Sequence[]];
  newBCEle = 
   MapThread[#1[##2] &, {eleTypes, ElementIncidents[newBCEle] + nc1, 
     markers}];
  newPEle = bm2["PointElements"];
  eleTypes = Head /@ newPEle;
  If[ElementMarkersQ[newPEle], markers = ElementMarkers[newPEle], 
   markers = Sequence[]];
  newPEle = 
   MapThread[#1[##2] &, {eleTypes, ElementIncidents[newPEle] + nc1, 
     markers}];
  ToBoundaryMesh["Coordinates" -> Join[c1, c2], 
   "BoundaryElements" -> Flatten[{bm1["BoundaryElements"], newBCEle}],
    "PointElements" -> Flatten[{bm1["PointElements"], newPEle}], opts]]
BoundaryElementMeshJoin[r1_, r2_, r3__] := 
  BoundaryElementMeshJoin[BoundaryElementMeshJoin[r1, r2], r3];

Now, we can create a domain region and inclusion region boundary meshes and join them together to create a multiregion element mesh.

(* Code to create and join several distinct boundary meshes *)
(* Create Regions *)
srdom = Rectangle[{0, 0}, {10, 5}];
nrdom = ToNumericalRegion[srdom];
srdisk = Disk[{3, 2}, 3/2];
srrect = Rectangle[{9, 4}, {10, 5}];
(* Estimate Region Bounds of Full Domain *)
symbolicBounds = RegionBounds[srdom];
(* Create Boundary Meshes for Each Region *)
(bm1 = ToBoundaryMesh[srdom, symbolicBounds])["Wireframe"];
(bm2 = ToBoundaryMesh[srdisk, symbolicBounds])["Wireframe"];
(bm3 = ToBoundaryMesh[srrect, symbolicBounds])["Wireframe"];
(* Join Boundary Meshes *)
(bm = BoundaryElementMeshJoin[bm1, bm2])["Wireframe"];
(bm = BoundaryElementMeshJoin[bm, bm3])["Wireframe"];
SetNumericalRegionElementMesh[nrdom, bm];
meshTriangle = 
  ToElementMesh[nrdom, 
   "RegionMarker" -> {{{1, 1}, 1, 0.5}, {{3, 2}, 2, 0.5}, {{9.5, 4.5},
       3, 0.5}}];
meshTriangle[
 "Wireframe"[
  "MeshElementStyle" -> {FaceForm[Red], FaceForm[Green], 
    FaceForm[Yellow]}, ImageSize -> Medium]]

Now, we can solve your PDE using a piecewise function and the region markers for the inclusions and domain defined above.

c = Evaluate[Piecewise[{{DiagonalMatrix@{100, 20}, ElementMarker == 2},
     {DiagonalMatrix@{100, 20}, ElementMarker == 3},
     {DiagonalMatrix@{3, 2}, True}}]];
pde = Inactive[Div][c.Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0;
bc = {DirichletCondition[u[x, y] == 0, x == 0], 
   DirichletCondition[u[x, y] == 100, x == 10]};
usol = NDSolveValue[{pde, bc}, u, {x, y} \[Element] meshTriangle];
ContourPlot[usol[x, y], Element[{x, y}, meshTriangle], 
 AspectRatio -> Automatic]

Let's make the disk a resistor.

c = Evaluate[
   Piecewise[{{DiagonalMatrix@{0.0003, 0.0002}, ElementMarker == 2},
     {DiagonalMatrix@{100, 20}, ElementMarker == 3},
     {DiagonalMatrix@{3, 2}, True}}]];
pde = Inactive[Div][c.Inactive[Grad][u[x, y], {x, y}], {x, y}] == 0;
bc = {DirichletCondition[u[x, y] == 0, x == 0], 
   DirichletCondition[u[x, y] == 100, x == 10]};
usol = NDSolveValue[{pde, bc}, u, {x, y} \[Element] meshTriangle];
ContourPlot[usol[x, y], Element[{x, y}, meshTriangle], 
 PlotPoints -> All, AspectRatio -> Automatic]

It seems to be behaving as expected.

This is an extentsion to Tim Laska's answer. I have added the BoundaryElementMeshJoin (and a few other Boolean operations) for boundary element meshes into the FEMAddOns paclet. The installation of the paclet is now very easy since the installation can be done via the FEMAddOnsInstall resource function.

Install and load the paclet:

ResourceFunction["FEMAddOnsInstall"][]
Needs["FEMAddOns`"]
(* Paclet[FEMAddOns, 1.3.2] *)

The run the above code as usual:

(*Code to create and join several distinct boundary meshes*)(*Create \
Regions*)srdom = Rectangle[{0, 0}, {10, 5}];
nrdom = ToNumericalRegion[srdom];
srdisk = Disk[{3, 2}, 3/2];
srrect = Rectangle[{9, 4}, {10, 5}];
(*Estimate Region Bounds of Full Domain*)

symbolicBounds = RegionBounds[srdom];
(*Create Boundary Meshes for Each Region*)
(bm1 = 
    ToBoundaryMesh[srdom, symbolicBounds])["Wireframe"];
(bm2 = ToBoundaryMesh[srdisk, symbolicBounds])["Wireframe"];
(bm3 = ToBoundaryMesh[srrect, symbolicBounds])["Wireframe"];
(*Join Boundary Meshes*)
(bm = 
    BoundaryElementMeshJoin[bm1, bm2, bm3])["Wireframe"];
SetNumericalRegionElementMesh[nrdom, bm];
meshTriangle = 
  ToElementMesh[nrdom, 
   "RegionMarker" -> {{{1, 1}, 1, 0.5}, {{3, 2}, 2, 0.5}, {{9.5, 4.5},
       3, 0.5}}];
meshTriangle[
 "Wireframe"[
  "MeshElementStyle" -> {FaceForm[Red], FaceForm[Green], 
    FaceForm[Yellow]}, ImageSize -> Medium]]