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["NDSolveFEM"]
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!

If you use ToElementMesh versus DiscretizeRegion and look at the "Wireframe", then you will see that your mesh looks pretty ugly.

mesh["Wireframe"] You can make a prettier mesh by starting with a geometry that has more cube-like aspect ratio and then scale the z-coordinate and remesh as shown below:

Needs["NDSolveFEM"]
(* Stretch Disc Region for better aspect ratio *)
disc = ImplicitRegion[x^2 + y^2 <= 1 && 0 <= z <= 1, {x, y, z}];
mesh0 = ToElementMesh[disc]
mesh0["Wireframe"]
(* Scale z-coords and create new ElementMesh *)
crd = mesh0["Coordinates"];
crd[[All, 3]] = crd[[All, 3]]/10;
mesh = ToElementMesh["Coordinates" -> crd,
"MeshElements" -> mesh0["MeshElements"],
"BoundaryElements" -> mesh0["BoundaryElements"],
"PointElements" -> mesh0["PointElements"]];
mesh["Wireframe"] Now, you can perform your NDSolveValue in the following way and it should default to the FEM method.

U = NDSolveValue[{-Laplacian[
u[x, y, z], {x, y, z}] == +50 NeumannValue[1,
x >= 0 && y >= 0 && z == 1/10] -
1 NeumannValue[u[x, y, z] - 20, z == 1/10]}, u,
Element[{x, y, z}, mesh]]
DensityPlot[U[x, y, 1/10],
Element[{x, y}, ImplicitRegion[x^2 + y^2 <= 1, {x, y}]],
ColorFunction -> (ColorData["TemperatureMap"][#1] &)] Now, you can verify that the solution required about 10x fewer elements than the original approach.

U["ElementMesh"]
(* ElementMesh[{{-1., 1.}, {-1., 1.}, {0.,
0.1}}, {TetrahedronElement["<" 16103 ">"]}] *)


# Speed Improvement

In 3D problems, it is often useful to extract boundary meshes for visualizations. Below, I show a method to extract the top surface for the DensityPlot. The overhead in extracting the top surface is dwarfed by the efficiency gains in the DensityPlot function. On my machine, the following code was 350x faster than the code in the OP an 16.7x faster than the code above.

AbsoluteTiming[
(* Stretch Disc Region for better aspect ratio *)
disc = ImplicitRegion[x^2 + y^2 <= 1 && 0 <= z <= 1, {x, y, z}];
mesh0 = ToElementMesh[disc];
(* Scale z-coords and create new ElementMesh *)
crd = mesh0["Coordinates"];
crd[[All, 3]] = crd[[All, 3]]/10;
mesh = ToElementMesh["Coordinates" -> crd,
"MeshElements" -> mesh0["MeshElements"],
"BoundaryElements" -> mesh0["BoundaryElements"],
"PointElements" -> mesh0["PointElements"]];
U = NDSolveValue[{-Laplacian[
u[x, y, z], {x, y, z}] == +50 NeumannValue[1,
x >= 0 && y >= 0 && z == 1/10] -
1 NeumannValue[u[x, y, z] - 20, z == 1/10]}, u,
Element[{x, y, z}, mesh]];
(* Extract Top Surface as 2D Mesh *)
(* Select Top Surface by Normal *)
mask = (0.9999 < -{0, 0, 1}.#) & /@ mesh0["BoundaryNormals"][];
topelmi =
(* Extract unique incidents from top surface *)
unique = DeleteDuplicates@Flatten[topelmi];
(* Create an associate to renumber connectivity *)
crd2d = mesh0["Coordinates"][[All, 1 ;; 2]];
crd2dtop = crd2d[[unique]];
newincidents = MapAt[crddict, topelmi, {All, All}];
mesh2d =
ToElementMesh["Coordinates" -> crd2dtop,
"MeshElements" -> {TriangleElement[newincidents]}];
(* Plot using new 2D mesh *)
DensityPlot[U[x, y, 1/10], {x, y} \[Element] mesh2d,
ColorFunction -> (ColorData["TemperatureMap"][#1] &)]
] • Thank you for your answer. I observed the same effect, using MeshTools to create easy structured mesh (as you did) , that the mesh keeps unchanged, as desired. – Ulrich Neumann Sep 27 '19 at 14:33
• Technically, my mesh is unstructured. On my machine, the number of mesh elements changes each time I mesh. Unstructured meshers have an easier time with domains with aspect ratios near unity. Structured or extruded meshes are the way to go for this geometry. If I saw support for a wedge element, I probably would have extruded a Disk. Also, I added a section to speed up the DensityPlot that you may be interested in since plotting seems to be the throughput limiter. – Tim Laska Sep 28 '19 at 12:22
• Thanks again, very interesting! Your basic trick meshing well proportioned geometries and scale the mesh afterwards is what I was looking for. I've to elaborate your new remarks... – Ulrich Neumann Sep 28 '19 at 12:36