Mesh.facets_around() not working after refinement

[duarte-jfs]

I was working on the mesh refinement, trying to understand how can boundaries be kept, and came across the issue of the mesh.facets_around() not working properly after a refinement.

from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
import shapely
import shapely.affinity
from skfem import Basis, ElementTriP0, adaptive_theta
from skfem.io.meshio import from_meshio

from femwell.maxwell.waveguide import compute_modes, eval_error_estimator
from femwell.mesh import mesh_from_OrderedDict

polygons_paper = OrderedDict(
    left=shapely.LineString(((0, 0), (0, 1))),
    right=shapely.LineString(((1, 0), (1, 1))),
    top=shapely.LineString(((0, 1), (1, 1))),
    square=shapely.box(0.6, 0.6, 0.8, 0.8),
    core=shapely.box(0, 0, 0.5, 0.5),
    clad=shapely.box(0, 0, 1, 1),
)

boundaries = [
    [lambda x: x[0] == np.min(x[0]), lambda x: x[0] == np.max(x[0])],
    [lambda x: x[0] == np.min(x[0]), lambda x: x[0] == np.max(x[0])],
    [lambda x: x[0] == np.min(x[0]), lambda x: x[1] == np.max(x[1])],
    [lambda x: x[0] == np.min(x[0]), lambda x: x[1] == np.max(x[1])],
]

epsilons_paper = [
    {"core": 2.25, "clad": 1, 'square': 1},
    {"core": 8, "clad": 1, 'square': 1},
    {"core": 1, "clad": 2.25, 'square': 1},
    {"core": 1, "clad": 8, 'square': 1},
]

neff_values_paper = [1.27627404, 2.65679692, 1.387926425, 2.761465320]
num = 3

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

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]
)

check_boundary(mesh)
# print(mesh.)
mesh = refined(adaptive_theta(eval_error_estimator(modes[0].basis, modes[0].E), theta=0.5))
check_boundary(mesh)

Before refinement:

After refinement:

[kinnala]

Thanks for the report. Based on this, all named boundaries should be removed from the mesh during adaptive refine: https://github.com/kinnala/scikit-fem/blob/c6b36757a2ba3ebb048abfc9d49a9a9205155a35/skfem/mesh/mesh_tri_1.py#L389

What is check_mesh doing?

[kinnala]

Anyhow, don’t expect this to be fixed any time soon. It is a somewhat challenging indexing problem to solve in general.

[duarte-jfs]

Sorry, the check_boundary is just taking care of the plots. Getting elements coordinates and normal vectors and plotting.

Based on this, all named boundaries should be removed from the mesh during adaptive refine:

The thing is, the named boundaries are not preserved already. What is preserved is the sub domains, but the ordering may be messed up. I don't know

Anyhow, don’t expect this to be fixed any time soon. It is a somewhat challenging indexing problem to solve in general.

Yes, I can only imagine. It would be a great feature though

[kinnala]

I don't quite follow what happens in this issue. Can you provide the code of check_mesh?

[duarte-jfs]

I don't quite follow what happens in this issue. Can you provide the code of check_mesh?

def check_boundary(mesh):
    plt.figure()
    plt.scatter(*mesh.p, marker = 'o', facecolor = 'None', edgecolor = 'black', s = 15)

    # boundary = mesh.boundaries['square___clad']
    boundary = mesh.facets_around(mesh.subdomains['square'])
    facet_basis = basis0.boundary(boundary)
    facets = mesh.facets[:,boundary]
    colors = plt.cm.jet(np.linspace(0,1,facets.shape[1]))

    for i in range(facets.shape[1]):
        # plt.text(*data[:,i], f'{i}')

        x = [mesh.p[0,facets[0,i]], mesh.p[0,facets[1,i]]]
        y = [mesh.p[1,facets[0,i]], mesh.p[1,facets[1,i]]]

        norm = np.asarray([facet_basis.normals[0,i,0], 
                           facet_basis.normals[1,i,0]])

        norm = norm/np.sqrt(np.sum(np.abs(norm)**2)) * 0.1

        plt.arrow(np.mean(x), np.mean(y), dx = norm[0], dy = norm[1], color = colors[i], 
                  head_width = 0.03, 
                  head_length = 0.01,
                  length_includes_head = True)

        plt.plot(x, y, color = colors[i])
        # plt.text(np.mean(x), np.mean(y), f"{i}")

    plt.scatter(mesh.p[0,facets[0]], mesh.p[1,facets[0]], c = colors, s = 20)
    plt.scatter(mesh.p[0,facets[1]], mesh.p[1,facets[1]], c = colors, s = 20)

[kinnala]

I'm unable to reproduce this with my own minimal example, perhaps it was fixed by another issue, or then I misunderstood the issue:

from skfem import *
import numpy as np

m = MeshTri().refined(4).with_subdomains({'test': lambda x: (x[0] < 0.5) * (x[0] > 0.3) * (x[1] < 0.5) * (x[1] > 0.3)})
m = m.with_boundaries({'test': m.facets_around('test')})
m.draw(boundaries=True).show()
m = m.refined(np.arange(200))
m = m.with_boundaries({'test': m.facets_around('test')})
m.draw(boundaries=True).show()

Before refine:

After refine:

The normal vector should be pointing opposite to the red triangle.