I'm having trouble finishing the calculations in Mathematica 13.2.1.0.
I'm trying to solve a 3D heat conduction equation (a type of partial differential equation) using the FEM package, but NDSolveValue didn't finish.
The fem mesh is created as intended.
Changing MaxStepSize, WorkingPrecision, and MeshRefinementFunction did not improve the situation.
I am not sure if the code of initial conditions and boundary conditions are appropriate.
I would appreciate any advice you could give me.

Code:

Needs["NDSolve`FEM`"];
simX = 0.05;
simY = 0.03;
gapX = 0.03;
gapY = 0.01;
gapZ = 0.01;
TotalArea = Rectangle[{-simX, -simY}, {simX, simY}];
Board1 = Rectangle[{-gapX, -gapY}, {gapX, -(gapY + gapZ)}];
Board2 = Rectangle[{-gapX, gapY}, {gapX, gapY + gapZ}];
Boards = RegionUnion[Board1, Board2];
Gap =  Rectangle[{-gapX, -gapY}, {gapX, gapY}];
Env = RegionDifference[TotalArea, Gap];

(* Mesh Generation *)
bmesh = ToBoundaryMesh[
  "Coordinates" ->
   {{-simX, -simY}, {-simX, simY}, {simX, simY}, {simX, -simY},
    {-gapX, -gapY}, {-gapX, -(gapY + gapZ)}, {gapX, -(gapY + 
        gapZ)}, {gapX, -gapY},
    {-gapX, gapY}, {-gapX, gapY + gapZ}, {gapX, gapY + gapZ}, {gapX, 
     gapY} },
  "BoundaryElements" ->
   {LineElement[{
      {1, 2}, {2, 3}, {3, 4}, {4, 1},
      {5, 6}, {6, 7}, {7, 8}, {8, 5},
      {9, 10}, {10, 11}, {11, 12}, {12, 9},
      {6, 7}, {7, 11}, {11, 10}, {10, 6}
      }]},
  "RegionHoles" ->
   {{0, -(gapY + gapZ/2)}, {0, gapY + gapZ/2}}
  ]
bmesh["Wireframe"]
mesh = ToElementMesh[bmesh, 
   MeshRefinementFunction -> 
    Function[{vertices, area}, 
     Block[{x, y}, {x, y} = Mean[vertices]; 
      If[Abs[x] < gapX*1.5 && Abs[y] < gapY*1.5, area > 0.5, 
       area > 0.9 ] 
      ]]];
mesh["Wireframe"]

c = 1006;(*(specific heat[J/kg ℃]）*)
rho = 1.166;(*(density[kg/m3]）*)
k = 0.0257; (*Thermal conductivity of air[W/(m K)]*)
eq = c*rho*D[T[x, y, t], t] == 
   k (D[T[x, y, t], {x, 2}] + D[T[x, y, t], {y, 2}]);
T0 = 30 + 273.15; 
T1 = 40 + 273.15;
ic = T[x, y, 0] == 
   Piecewise[{{T1, {x, y} \[Element] Env}, {T0, {x, y} \[Element] 
       Gap}}];
bc = NeumannValue[0, 
   x == -simX || x == simX || y == -simY || y == simY];
U = NDSolveValue[{eq, ic2}, T, {x, y} \[Element] mesh, {t, 0, 1}, 
   MaxStepSize -> 1, WorkingPrecision -> MachinePrecision];