Ryfujk

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 has holes.

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.

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["NDSolve`FEM`"]

mesh = ToElementMesh[mreg,"MeshOrder"->1]

(* ElementMesh[{{7., 353.}, {7., 353.}}, {TriangleElement["<" 1057 ">"]}] *)



mesh["Wireframe"]

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.