Basis creation fails on punctured multipatch topologies

[JochenHinz]

I have spotted a possible bug.

I am currently working with "punctured" multipatch topologies, i.e., topologies with interior interface edges across which the basis is only C^-1 continuous, instead of the usual C^0 inter-patch coupling.

What I do is, I assign differing ids to the 2 patch boundaries that are associated with the edge on either side of the two patches that share it. Usually basis creation works fine, however I am adding a minimal example that reproduces the bug.

from nutils import mesh
from nutils.topology import PatchBoundary, MultipatchTopology, Patch


def bug():

  patches = (0, 4, 2, 5), (4, 1, 5, 3), (2, 5, 6, 7), (5, 3, 7, 8), \
            (6, 7, 9, 8)

  patchverts = (0, 0), (0, 1), (1, 0), (1, 1), (0, .5), (1, .5), \
               (2, 0), (1.5, .5), (2, 1), (2.5, .5)

  domain, geom = mesh.multipatch(patches=patches, patchverts=patchverts, nelems=1)
  basis = domain.basis('spline', 2)  # works

  # assign to the first occurrence of boundary id (6, 9) the new id (13, 14)
  # to create a C^-1 continuity across the interior interface edge
  map_old_id_new = {(6, 9): (13, 14), (9, 6): (14, 13)}

  new_boundaries = []
  for patch in domain.patches:
    mybndries = []
    for bndry in patch.boundaries:
      id, dim, side, reverse, transpose = bndry.id, bndry.dim, bndry.side, bndry.reverse, bndry.transpose
      newid = map_old_id_new.get(id, id)
      if id in map_old_id_new:  # id has been mapped -> make sure the second occurrence does not get mapped as well
        del map_old_id_new[id]
        del map_old_id_new[id[::-1]]

      mybndries.append(PatchBoundary(newid, dim, side, reverse, transpose))
    new_boundaries.append(mybndries)

  new_patches = []
  for patch, bndry in zip(domain.patches, new_boundaries):
    new_patches.append(Patch(patch.topo, patch.verts, bndry))

  newdomain = MultipatchTopology(new_patches)
  basis = newdomain.basis('spline', 2)  # works

  # do the same for edge (4, 5)

  # assign to the first occurrence of boundary id (4, 5) the new id (13, 14)
  # to create a C^-1 continuity across the interior interface edge
  map_old_id_new = {(4, 5): (13, 14), (5, 4): (14, 13)}

  new_boundaries = []
  for patch in domain.patches:
    mybndries = []
    for bndry in patch.boundaries:
      id, dim, side, reverse, transpose = bndry.id, bndry.dim, bndry.side, bndry.reverse, bndry.transpose
      newid = map_old_id_new.get(id, id)
      if id in map_old_id_new:  # id has been mapped -> make sure the second occurrence does not get mapped as well
        del map_old_id_new[id]
        del map_old_id_new[id[::-1]]

      mybndries.append(PatchBoundary(newid, dim, side, reverse, transpose))
    new_boundaries.append(mybndries)

  new_patches = []
  for patch, bndry in zip(domain.patches, new_boundaries):
    new_patches.append(Patch(patch.topo, patch.verts, bndry))

  newdomain = MultipatchTopology(new_patches)
  basis = newdomain.basis('spline', 2)  # crashes


if __name__ == '__main__':
  bug()

I have also been able to identify the problem.
in the nutils.topology.MultipatchTopology.basis_spline function, stopping the interpreter in line 3134, I get commonboundarydofs[(2, 5)] = [ [6, 7, 8], [18, 19, 20] ], commonboundarydofs[(3, 5)] = [ [17, 16, 15], [29, 28, 27]] and commonboundarydofs[(5, 7)] = [ [20, 23, 26], [27, 30, 33] ].
This means that DOF 8 is coupled to {15, 20, 27}. However, before line 3140, I get merge[15] = 15.

I have been able to fix this issue by running lines 3138 and 3139 twice, i.e., performing the merging (on the same merge array) twice.

[joostvanzwieten]

#765 contains a fix. Can you test branch fix-merge-dofs?

The problem manifests itself already with the following code (added as unittest):

patches = (0, 4, 2, 5), (9, 1, 5, 3), (2, 5, 6, 7), (5, 3, 7, 8)
domain, geom = mesh.multipatch(patches=patches, nelems=1)
basis = domain.basis('spline', 1)

[JochenHinz]

I can confirm that this fixes the error.
Thanks

[joostvanzwieten]

Thanks for testing.