In a simple example I try to solve the heat equation using NDSolve and Method->"FiniteElement". I know that NDSolve gives the solution as an interpolationfunction inside an ElementMesh.

First the mesh is created

Needs["NDSolve`FEM`"]
disc = ImplicitRegion[  x^2 + y^2 <= 1 && 0 <= z <= 1/10, {x, y , z}];
mesh = DiscretizeRegion[disc, MaxCellMeasure -> {"Volume" -> Pi/10}];
MeshCells[mesh, 3] // Length  (*36279*)

The mesh has around 36000 elements.

Simulation inside this mesh

U = NDSolveValue[{ 0 ==    Laplacian[u[ x, y, z ] , {x, y, z }]+ 50NeumannValue[1,  x >= 0 && y >= 0 && z == 1/10     ]-  1     NeumannValue[u[ x, y, z] - 20, z == 1/10] },u, 
Element[{x, y, z }, mesh ] , Method -> "FiniteElement"] 
DensityPlot[U[x, y, 1/10],Element[{x, y}, ImplicitRegion[x^2 + y^2 <= 1, {x, y}]],ColorFunction -> (ColorData["TemperatureMap"][#1] &),MaxRecursion -> 4]

gives a feasible solution

with considerably refined mesh

U["ElementMesh"]
(*TetrahedronElement["<" 160748 ">"]*) 

consisting of now 160748 elements, that means 440% increase of meshsize!

My question: Is it possible to keep the initial mesh mesh unchanged(fixed) by NDSolve?

Thanks!