How to use mesh refinement and ensure named boundaries remain?

[duarte-jfs]

I was playing around with mesh refinement tutorial (https://helgegehring.github.io/femwell/photonics/examples/refinement.html) and I noticed that along the iteration a warning pops up as

Named boundaries invalidated by a call to Mesh.refined()

This seems to be issue from skfem, and it can cause problems when we want to have later reference to zthose named boundaries for some integrals. Does anyone know a solution for it?

[duarte-jfs]

Looks like an old issue (kinnala/scikit-fem#104 (comment)) but not general solution seems to have been given (?)

[HelgeGehring]

Hmm, I thought it was fixed in kinnala/scikit-fem#1024 🤔

The strange thing is that it doesn't happen directly in the beginning? Also, it should change the mode if it suddenly stops being defined?
We should probably start by understanding why it happens exactly then? And if the boundaries are still defined correctly?

[duarte-jfs]

Sorry for the long delay in replying, I got wrapped up in other things.
A lot of the code has gone over my head, but here's my understanding of it:

That issue (kinnala/scikit-fem#1024) was resolved the preservation of subdomains (are subdomains like the inside of square while the boundaries are the line of the rectangle?). And there is a very pythonic change in the code in skfem/mesh/mesh_tri_1.py:

This suggests me that it only the subdomains that were preserved rather than preserving boundaries. Indeed, the subdomains are preserved. For instance, the example below (from the femwell docs) shows that:

mesh = from_meshio(
    mesh_from_OrderedDict(polygons_paper, {}, default_resolution_max=0.1, filename="mesh.msh")
)

prev = 0
for i in range(10):
    basis0 = Basis(mesh, ElementTriP0())
    epsilon = basis0.zeros()
    for subdomain, e in epsilons_paper[num].items():
        epsilon[basis0.get_dofs(elements=subdomain)] = e

    modes = compute_modes(
        basis0, epsilon, wavelength=1.5, num_modes=1, order=1, metallic_boundaries=boundaries[num]
    )
    # modes[0].show(modes[0].E.real, direction="x")
    # print(f"Error: {modes[0].n_eff - neff_values_paper[num]:.5g} | {np.abs(modes[0].n_eff-prev):.5g}")
    
    # prev = modes[0].n_eff

    print('=============================================')
    print(i)
    print(mesh._subdomains.keys())
    print(mesh.boundaries.keys() if mesh.boundaries is not None else None)
    

    print(f'old mesh has boundaries:{True if mesh.boundaries is not None else False}')
    mesh = mesh.refined(adaptive_theta(eval_error_estimator(modes[0].basis, 
                                                            modes[0].E), 
                                       theta=0.5))
    print(f'New mesh has boundaries:{True if mesh.boundaries is not None else False}')
    print('============================================')
    # break
    # mesh.draw().show()

This outputs:

=============================================
0
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
dict_keys(['left', 'right', 'top', 'core___clad'])
old mesh has boundaries:True
New mesh has boundaries:True
============================================
=============================================
1
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
dict_keys(['left', 'bottom', 'right', 'top'])
old mesh has boundaries:True
New mesh has boundaries:True
============================================
=============================================
2
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
dict_keys(['left', 'bottom', 'right', 'top'])
old mesh has boundaries:True
New mesh has boundaries:True
============================================
=============================================
3
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
dict_keys(['left', 'bottom', 'right', 'top'])
old mesh has boundaries:True
New mesh has boundaries:True
============================================
=============================================
4
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
dict_keys(['left', 'bottom', 'right', 'top'])
old mesh has boundaries:True
New mesh has boundaries:True
============================================
=============================================
5
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
dict_keys(['left', 'bottom', 'right', 'top'])
old mesh has boundaries:True
New mesh has boundaries:True
============================================
=============================================
6
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
dict_keys(['left', 'bottom', 'right', 'top'])
old mesh has boundaries:True
New mesh has boundaries:True
============================================
Named boundaries invalidated by a call to Mesh.refined()
=============================================
7
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
dict_keys(['left', 'bottom', 'right', 'top'])
old mesh has boundaries:True
New mesh has boundaries:False
============================================
=============================================
8
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
None
old mesh has boundaries:False
New mesh has boundaries:False
============================================
=============================================
9
dict_keys(['core', 'clad', 'gmsh:bounding_entities'])
None
old mesh has boundaries:False
New mesh has boundaries:False
============================================

My understanding for now is flawed because I do not understand why after the first iteration it has boundaries at all since the ._adaptive() runs and the replace(..., _boundaries = None) is evaluated.

The strange thing is that it doesn't happen directly in the beginning?

I'm confused as well

Also, it should change the mode if it suddenly stops being defined?

I don't it should since is only changing our references to it and it is only the boundaries that are defined. In fact, the basis0.get_dofs() continues to work just fine... I imagine this is a consequence of the preservation of _subdomains?

I'll update here if I find something more