Surface Boundaries

General
1

Hi; I want to split my simulation into several regions.
In order to do so, I group the surface elements. That works as expected.
(Minimum working example attached)
However as soon as I try to assign the boundaries of the surface elements some name (in order to be able to utilize these connecting elements later on), I get inconsistency warnings.
In order to have a consistent quad mesh, I generate the mesh myself from a list of polygons (here called ‘mesh’) and a corresponding list of material indices (‘indsM’).
The code is a simple extension of the one given in the docs.

    from netgen.meshing import MeshPoint, Element2D, Element1D, FaceDescriptor
    from netgen.meshing import Mesh as NMMesh 

    net_mesh = NMMesh()
    net_mesh.dim = 2
    ptsI = [net_mesh.Add(MeshPoint(Pnt(*pt, 0))) for pt in mesh.pts]
    for i, m in enumerate(mat_names):
        net_mesh.SetMaterial(i+1, m)
        net_mesh.Add(FaceDescriptor(surfnr=i+1,domin=0,bc=i+1))
    for ind, pol in zip(indsM, mesh.polsI):
        net_mesh.Add(Element2D(ind+1, [ptsI[i] for i in pol[:4]]))
    net_mesh.Save('geom_m_r.vol')
    for i, j in segsI:
        net_mesh.Add(Element1D([ptsI[i], ptsI[j]]))
    net_mesh.Save('geom_m.vol')

This produces the attached *.vol files.
geom_m.vol (11.8 KB)
geom_m_r.vol (5.9 KB)
The code to produce the warning is straightforward.


from pathlib import Path
from functools import reduce
from itertools import groupby
import re
import operator

import netgen.meshing as meshing
from ngsolve import Mesh, H1, FESpace, BilinearForm


if __name__ == '__main__':
    prog = re.compile(r'.*_R(\d?)')
    def _group(m, prog=prog):
        return int(prog.match(m).group(1))//4

    p = Path(__file__).parent
    mesh2 = meshing.Mesh()
    mesh2.Load(str(p / 'geom_m_r.vol'))
    mesh = Mesh(mesh2)

    regions = [reduce(operator.add, (mesh.Materials(m) for m in g))
        for k, g in groupby(sorted(mesh.GetMaterials(), key=_group), _group)]

    fes = H1(mesh, order=1)

    print ("ndof =", fes.ndof)
    print(sum(fes.FreeDofs()))

    fes = FESpace([H1(mesh, order=1, definedon=dom) for dom in regions])

    print ("ndof =", fes.ndof)
    print(sum(fes.FreeDofs()))

    u, v = fes.TnT()
    a = BilinearForm(fes)
            
    a.Assemble()

How can I name the surface boundaries without seemingly introducing new degrees of freedom when the simulation gets split.
I also tried to create the mesh using geom2d; something of the form

    geo = CSG2d()
    for i, pol in enumerate(mesh.pols):
        ll = [el for pt, bc in zip(pol[:4], repeat('default')) for el in [tuple(pt), EI(bc=bc)]]
        rect = Solid2d(ll, mat=f'mat_{i:d}')
        geo.Add(rect)

This worked insofar as the boundaries can be assigned a name; however I do not know how to convert this geometry in a Mesh as is, meaning without introducing new elements.

I hope somebody can clarify my confusion; in addition is there an example how all this can be accomplished in 3D.
Best to all