Implement local mesh refinement in GMSH, trying to remove degenerate mesh elements.

Add Solution parameter `min_mesh_quality` defaulting to 5e-7 which determined a mesh element
quality threshold under which elements are refined by redoing the whole mesh with locally modified mesh size fields.

Set `mesh_optimizer` default to Netgen in ElmerSolution

Author: Tuomas Myllari

klayout_package/python/kqcircuits/simulations/export/elmer/elmer_solution.py

@@ -48,13 +48,18 @@ class ElmerSolution(Solution):
             For example, if the dictionary is {'substrate*': 10, 'substrate*&vacuum': [2, 5], 'global_max': 100}, then
             the maximal mesh element length is 10 inside the substrates and 2 on region which is less than 5 units away
             from any substrate-vacuum interface. Outside these regions, the mesh element size can increase up to 100.
-        mesh_optimizer: Dictionary to determine mesh optimization, or None (default) to ignore optimization. The
-            dictionary can contain keywords 'method', 'force', 'niter' and 'dimTags'. See Gmsh manual
-            (gmsh.model.mesh.optimize) for details. The default value for 'method' is 'Netgen'.
+        mesh_optimizer: Dictionary to determine mesh optimization, or None to ignore optimization. The dictionary can
+            contain keywords 'method', 'force', 'niter' and 'dimTags'. See Gmsh manual (gmsh.model.mesh.optimize) for
+            details. The default value is {'method': 'Netgen'}.
         vtu_output: Output vtu files to view fields in Paraview.
                     Turning this off will make the simulations slightly faster
         save_elmer_data: Save the full Elmer model after simulation. This can be used to restart the simulation
                          or extract result field values as a post-processing step.
+        min_mesh_quality: If the initial Gmsh mesh contains elements with quality below this limit, a local mesh
+            refinement around those elements is applied. This causes the whole mesh to be recomputed and can therefore
+            cause performance issues if set too high. The aim of the remeshing is to prevent fatal errors due
+            to degenerate elements in Elmer. The default value of 5e-7 is determined experimentally and might need
+            adjustment. Setting this to 0 disables the feature.
 
         linear_system_method: Method for solving the FEM linear system of equations in Elmer. For iterative methods use
                 "GCR", "bicgstab" or any other iterative solver mentioned in ElmerSolver manual section 4.3.1.
@@ -90,9 +95,10 @@ class ElmerSolution(Solution):
     is_axisymmetric: bool = False
     mesh_levels: int = 1
     mesh_size: dict = field(default_factory=dict)
-    mesh_optimizer: dict | None = None
+    mesh_optimizer: dict | None = field(default_factory=lambda: {"method": "Netgen"})
     vtu_output: bool = True
     save_elmer_data: bool = False
+    min_mesh_quality: float = 5e-7
 
     linear_system_method: str = "GCR"
     convergence_tolerance: float = 1.0e-9
@@ -118,6 +124,8 @@ def __post_init__(self):
         if self.mg_smoothing_iterations is None:
             defaults = {"sgs": 1, "wjacobi": 4, "cjacobi": 4, "cg": 8}
             object.__setattr__(self, "mg_smoothing_iterations", defaults.get(self.mg_smoother.lower(), 1))
+        if self.mesh_optimizer == {}:
+            object.__setattr__(self, "mesh_optimizer", {"method": "Netgen"})
 
 
 @dataclass(kw_only=True, frozen=True)

klayout_package/python/scripts/simulations/circular_capacitor_c_deembed_sim.py

@@ -74,10 +74,7 @@
     "n": 64,
 }
 
-solution = ElmerCapacitanceSolution(
-    mesh_size=refine_metal_edges(2.0, 0.5),
-    mesh_optimizer={},
-)
+solution = ElmerCapacitanceSolution(mesh_size=refine_metal_edges(2.0, 0.5))
 
 # Prepare output directory
 dir_path = create_or_empty_tmp_directory(Path(__file__).stem + var_str + "_output")

klayout_package/python/scripts/simulations/circular_capacitor_epr_sim.py

@@ -97,7 +97,7 @@
 
 solution = ElmerEPR3DSolution(
     mesh_size=refine_metal_edges(2.0, 0.5),
-    mesh_optimizer={},
+    mesh_optimizer={"method": "Netgen", "force": True},
     voltage_excitations=voltage_excitations,
 )
 

klayout_package/python/scripts/simulations/circular_transmon_single_island_epr_sim.py

@@ -55,10 +55,7 @@
 }
 
 
-solution = ElmerEPR3DSolution(
-    mesh_size=refine_metal_edges(3.0, 0.5),
-    mesh_optimizer={},
-)
+solution = ElmerEPR3DSolution(mesh_size=refine_metal_edges(3.0, 0.5))
 
 # Prepare output directory
 dir_path = create_or_empty_tmp_directory(Path(__file__).stem + "_output")

klayout_package/python/scripts/simulations/double_pads_epr_sim.py

@@ -82,7 +82,6 @@
         # Set opposite voltages [1, -1] on the islands and smaller voltage 0.1 on coupler
         voltage_excitations=[1, 0.1, -1],
         mesh_size=refine_metal_edges(5.0, 0.5),
-        mesh_optimizer={},
         p_element_order=2,
     )
 # Get layout

klayout_package/python/scripts/simulations/elmer/gmsh_helpers.py

@@ -159,7 +159,7 @@ def produce_mesh(json_data: dict[str, Any], msh_file: Path) -> None:
     # Set meshing
     mesh_size = json_data.get("mesh_size", {})
     workflow = json_data.get("workflow", {})
-    set_meshing(mesh_size, new_tags, workflow)
+    mesh_field_ids, global_max_mesh_size = set_meshing(mesh_size, new_tags, workflow)
 
     # Remove layers without material
     for name, data in layers.items():
@@ -222,6 +222,34 @@ def produce_mesh(json_data: dict[str, Any], msh_file: Path) -> None:
     gmsh.model.mesh.generate(3)
 
     optimize_mesh(json_data.get("mesh_optimizer"))
+
+    min_quality = json_data["min_mesh_quality"]
+    if min_quality > 0:
+        deg_elements = get_elements_by_quality(min_quality)
+        if len(deg_elements) > 0:
+            logging.warning(f"Found {len(deg_elements)} poor elements. Starting iterative refinement")
+
+            refined_mesh_field_ids = mesh_field_ids.copy()
+            max_iter = 5
+            for i in range(max_iter):
+                print_mesh_quality_metrics()
+                logging.warning(f"Iterative refinement {i + 1}: Remeshing...")
+
+                refined_mesh_field_ids = remesh_with_local_box_refinement(
+                    deg_elements, refined_mesh_field_ids, global_max_mesh_size
+                )
+                optimize_mesh(json_data.get("mesh_optimizer"))
+
+                deg_elements = get_elements_by_quality(min_quality)
+
+                if len(deg_elements) == 0:
+                    logging.warning(f"Iterative refinement {i + 1}: Completed succesfully.")
+                    break
+                logging.warning(f"Iterative refinement {i + 1}: {len(deg_elements)} poor elements still left")
+            else:
+                logging.warning(f"Iterative refinement: Failed to improve the mesh in {max_iter} iterations")
+            print_mesh_quality_metrics()
+
     gmsh.write(str(msh_file))
 
     # Open mesh viewer
@@ -231,13 +259,86 @@ def produce_mesh(json_data: dict[str, Any], msh_file: Path) -> None:
     gmsh.finalize()
 
 
+def get_element_qualities(element_type: int | None = None) -> dict[int, float]:
+    """Get element quality metrics in dictionary format ElementTag -> Quality.
+    `element_type` can be used to filter by element type, for example tetras=4"""
+    qualities = {}
+    elemTypes, elemTags, _ = gmsh.model.mesh.getElements(3)
+    for etype, tags in zip(elemTypes, elemTags):
+        if element_type and element_type != etype:
+            continue
+
+        tag_qualities = gmsh.model.mesh.getElementQualities(tags)
+        qualities.update(dict(zip(tags, tag_qualities)))
+
+    return qualities
+
+
+def get_elements_by_quality(threshold: float, element_type: int | None = None) -> list[int]:
+    """Get Tags of Elements with quality below given threshold."""
+    return sorted([e for e, q in get_element_qualities(element_type).items() if q < threshold])
+
+
+def remesh_with_local_box_refinement(
+    refined_elements: list[int],
+    existing_mesh_fields: list,
+    global_max: float,
+    mesh_size_factor: float = 1.0,
+) -> list:
+    """Applies mesh refinement trying to improve quality of target elements. Mesh is refined in a box around each
+    target element with a mesh size equal to average edge length of the element.
+
+    Args:
+        refined_elements: List of elements to refine.
+        existing_mesh_fields: GMSH mesh field list of the current mesh
+        global_max: Global maximum mesh size
+        mesh_size_factor: Mesh size in the refinement box relative to the average edge length of the element.
+
+    Returns:
+        Gmsh mesh size fields for the refined mesh (including the original fields)
+    """
+    box_fields = []
+    for e in refined_elements:
+        nodeTags = gmsh.model.mesh.getElement(e)[1]
+        node_coords = [gmsh.model.mesh.getNode(n)[0] for n in nodeTags]
+
+        average_dist = sum(coord_dist(p1, p2) for p1, p2 in itertools.combinations(node_coords, 2)) / 6
+        # Set mesh size inside refinement box as average edge length times a factor
+        f_box = gmsh.model.mesh.field.add("Box")
+        mesh_size_inside = mesh_size_factor * average_dist
+        gmsh.model.mesh.field.setNumber(f_box, "VIn", mesh_size_inside)
+        gmsh.model.mesh.field.setNumber(f_box, "VOut", global_max)
+
+        gmsh.model.mesh.field.setNumber(f_box, "Thickness", (global_max - mesh_size_inside) / 2)
+        for coord, values in zip(["X", "Y", "Z"], zip(*node_coords)):
+            gmsh.model.mesh.field.setNumber(f_box, coord + "Min", min(values))
+            gmsh.model.mesh.field.setNumber(f_box, coord + "Max", max(values))
+
+        box_fields.append(f_box)
+
+    all_mesh_fields = existing_mesh_fields + box_fields
+    background_field_id = gmsh.model.mesh.field.add("Min")
+    gmsh.model.mesh.field.setNumbers(background_field_id, "FieldsList", all_mesh_fields)
+    gmsh.model.mesh.field.setAsBackgroundMesh(background_field_id)
+
+    gmsh.model.mesh.clear()
+    gmsh.model.mesh.generate(3)
+
+    return all_mesh_fields
+
+
+def print_mesh_quality_metrics(element_type: int | None = None):
+    """Prints minimum and average element quality metrics"""
+    qualities = np.array(list(get_element_qualities(element_type).values()))
+    logging.warning(f"Min quality: {float(np.min(qualities)):.4g}. Mean quality: {float(np.average(qualities)):.4g}")
+
+
 def optimize_mesh(mesh_optimizer: dict | None) -> None:
     """Optimize the mesh if the mesh_optimizer is a dictionary. Ignore mesh optimization if mesh_optimizer is None."""
     if mesh_optimizer is None:
         return
     try:
-        # Try optimizing with Netgen as the default method
-        gmsh.model.mesh.optimize(**{"method": "Netgen", **mesh_optimizer})
+        gmsh.model.mesh.optimize(**mesh_optimizer)
     except Exception as error:  # pylint: disable=broad-except
         logging.warning(f"Mesh optimization failed: {error}")
 
@@ -414,13 +515,18 @@ def get_recursive_children(dim_tags: Iterable[DimTag], include_parent: bool = Fa
 
 def set_meshing(
     mesh_size: dict[str, float | list[float]], layer_dts: dict[str, list[DimTag]], workflow: dict[str, Any]
-):
+) -> tuple[list, float]:
     """Applies mesh refinement and sets meshing options
 
     Args:
         mesh_size: Dictionary to determine mesh refinement
         layer_dts: dictionary with layer names as keys and lists of dim-tags as values
         workflow: Parameters for simulation workflow
+
+    Returns:
+        Tuple of
+            1. List of mesh field ids used in the refinement
+            2. global max mesh size
     """
     # Find global maximum mesh element length
     mesh_global_max_size = mesh_size.pop("global_max", 0)
@@ -459,6 +565,7 @@ def set_meshing(
     n_threads_dict = workflow["sbatch_parameters"] if "sbatch_parameters" in workflow else workflow
     gmsh_n_threads = int(n_threads_dict.get("gmsh_n_threads", 1))
     set_meshing_options(mesh_field_ids, mesh_global_max_size, gmsh_n_threads)
+    return mesh_field_ids, mesh_global_max_size
 
 
 def set_meshing_options(mesh_field_ids: list[int], max_size: float, n_threads: int) -> None:

klayout_package/python/scripts/simulations/swissmon_epr_sim.py

@@ -92,7 +92,6 @@
             voltage_excitations=[1.0, 0.01],
             # poor mesh for faster computation
             mesh_size=refine_metal_edges(5.0, 1.0),
-            mesh_optimizer={},
         ),
     )