# Holes in ElementMesh with ToElementMesh of ImplicitRegion

I am trying to plot a function in a region below a level curve of the function and within a cell. I have been doing this by calculating an ElementMesh using ImplicitRegion and ToElementMesh, but the result is inaccurate.

Here is the cell (it's just a square),

cell = Parallelogram[{-0.5, -0.5}, {{1., 0.}, {0., 1.}}];
Graphics[{Transparent, EdgeForm[Thick], cell}]


and the function,

f[kx_, ky_, n_] :=
Sort[Eigenvalues[{{(-1. + kx)^2 + (-1. + ky)^2, -0.23, 0., -0.23,
0.12, 0., 0., 0.,
0.}, {-0.23, (-1. + kx)^2 + (0. + ky)^2, -0.23, 0.12, -0.23,
0.12, 0., 0., 0.}, {0., -0.23, (-1. + kx)^2 + (1. + ky)^2, 0.,
0.12, -0.23, 0., 0., 0.}, {-0.23, 0.12,
0., (0. + kx)^2 + (-1. + ky)^2, -0.23, 0., -0.23, 0.12,
0.}, {0.12, -0.23,
0.12, -0.23, (0. + kx)^2 + (0. + ky)^2, -0.23, 0.12, -0.23,
0.12}, {0., 0.12, -0.23, 0., -0.23, (0. + kx)^2 + (1. + ky)^2,
0., 0.12, -0.23}, {0., 0., 0., -0.23, 0.12,
0., (1. + kx)^2 + (-1. + ky)^2, -0.23, 0.}, {0., 0., 0.,
0.12, -0.23,
0.12, -0.23, (1. + kx)^2 + (0. + ky)^2, -0.23}, {0., 0., 0.,
0., 0.12, -0.23, 0., -0.23, (1. + kx)^2 + (1. + ky)^2}}]][[
n]];
Plot3D[f[x, y, 4], {x, y} \[Element] cell, PlotPoints -> 50]


and what the region should look like,

isovalue = 1.29897233417072;
ContourPlot[f[x, y, 4], {x, y} \[Element] cell,
Contours -> {isovalue}, ColorFunction -> GrayLevel,
PlotPoints -> 100]


This is what I have tried

reg = ToElementMesh[
ImplicitRegion[
f[x, y, 4] < isovalue && {x, y} \[Element] cell, {x, y}],
"MaxBoundaryCellMeasure" -> 0.01, MeshQualityGoal -> 1,
PerformanceGoal -> "Quality", MaxCellMeasure -> 0.01,
"BoundaryMeshGenerator" -> "Continuation"];
RegionPlot[reg]


The region is no more accurate when I decrease MaxCellMeasure or MaxBoundaryCellMeasure. I also tried the solution suggested here.

Another approach is:

reg = ToElementMesh[
ImplicitRegion[
f[x, y, 4] < isovalue && {x, y} \[Element] cell, {x, y}],
"MaxBoundaryCellMeasure" -> 0.01, MeshQualityGoal -> 1,
PerformanceGoal -> "Quality", MaxCellMeasure -> 0.01,
"BoundaryMeshGenerator" -> {"RegionPlot", "SamplePoints" -> 41}];

reg["Wireframe"]


One thing to be a bit careful about is the question if the holes intersect the boundary. From the mesh it does not look like it but the math might say it.

I hope I interpreted your question correctly that you want a more accurate ElementMesh representation of the region.

First we create a high quality Graphics of the region of interest.

isovalue = 1.29897233417072;
(* Add some margins to plot range to get connected region. *)
tolerance = 0.05;
plot = ContourPlot[
f[x, y, 4],
{x, y} ∈ Cuboid[{-0.5, -0.5} - tolerance, {0.5, 0.5} + tolerance],
Contours -> {isovalue},
ColorFunction -> GrayLevel,
(* We need high quality plot for ImageMesh later. *)
PlotPoints -> 200,
Frame -> None
]


Create MeshRegion from Graphics object.

mreg = ImageMesh[ColorNegate[plot]]


And convert it to ElementMesh.

Needs["NDSolveFEM"]
mesh = ToElementMesh[mreg,"MeshOrder"->1]
(* ElementMesh[{{7., 353.}, {7., 353.}}, {TriangleElement["<" 1057 ">"]}] *)

mesh["Wireframe"]