Duplicate points

[griff10000]

I have found that dmsh sometimes generates duplicate points, which makes it difficult to use with meshless methods as this can result in singular matrices. An example is given below which produces three duplicate points.

import dmsh
import optimesh

geo = dmsh.Union([dmsh.Rectangle(-1.0, +1.0,-0.2, +0.2),
                  dmsh.Rectangle(+0.5, +1.0,-0.6, +1.0),
                  dmsh.Circle([-0.5, 0.0], 0.4)])
geo = dmsh.Difference(geo, dmsh.Circle([+0.5, 0.2], 0.2))
r = 0.1   
pts, cells = dmsh.generate(geo, edge_size=r, tol=1.0e-10)
# optionally optimize the mesh
pts, cells = optimesh.optimize_points_cells(pts, cells, "CVT (full)",
                                            1.0e-10, 100)

dmsh.helpers.show(pts, cells, geo)

I note there is a helper function unique_rows() that may help, but is there a function call that can delete duplicate points and also the associated cells?

[nschloe]

Can you try and simplify this more? E.g., is optimesh required for this? How about the tolerance? A simpler mesh?

In any case, duplicate points are bad and should indeed now happen.

Edit: Also, which points are duplicated?

[griff10000]

I have tried various values for r, setting tol=1.0e-14 and eliminating optimesh, but still get duplicate points.

The following function removes duplicates from the pts array, but that still leaves the associated cells to be removed.

def testForDups(pts) :
  # ref: https://stackoverflow.com/questions/31097247/remove-duplicate-rows-of-a-numpy-array
  new_array = [tuple(row) for row in pts]
  pts_unique = np.unique(new_array,axis=0)

  return pts_unique

And the following function returns the indices of duplicate rows,

def findDups(pts) :
  # Ref: https://stackoverflow.com/questions/48099479/how-do-you-find-and-save-duplicated-rows-in-a-numpy-array
  unq, count = np.unique(pts, axis=0, return_counts=True)
  repeated_groups = unq[count > 1]

  pts_dupIndices = []
  for repeated_group in repeated_groups:
      repeated_idx = np.argwhere(np.all(pts == repeated_group, axis=1))
      pts_dupIndices.append(repeated_idx.ravel())
  return(np.array(pts_dupIndices))

for the above example, it returns,

[[ 1 13]
 [ 2 16]]

Which indicates that rows 1 and 13 of pts are duplicated and rows 2 and 16 are also duplicated.

[griff10000]

Further experimentation shows that in the above example, the duplicate points are generated due to the union of two objects which have part of their boundaries exactly overlapping. By moving the overlapping boundary part of one object so that it lies within the other object boundary, duplicate points are avoided. Therefore, changing the above example to,

import dmsh
import optimesh

geo = dmsh.Union([dmsh.Rectangle(-1.0, +0.9,-0.2, +0.2),
                  dmsh.Rectangle(+0.5, +1.0,-0.6, +1.0),
                  dmsh.Circle([-0.5, 0.0], 0.4)])
geo = dmsh.Difference(geo, dmsh.Circle([+0.5, 0.2], 0.2))
r = 0.1   
pts, cells = dmsh.generate(geo, edge_size=r, tol=1.0e-10)
# optionally optimize the mesh
pts, cells = optimesh.optimize_points_cells(pts, cells, "CVT (full)",
                                            1.0e-10, 100)

dmsh.helpers.show(pts, cells, geo)

avoids this particular problem. However, this approach is unlikely to solve all problems.

[gdmcbain]

Yeah, I have another example which doesn't have overlapping boundaries but still has duplicate points. The geometry is that of DFG benchmark 2D-2, as reused in the FEniCS tutorial. Looking to redo that in pure Python, without mshr (gdmcbain/fenics-tuto-in-skfem#5), I tried dmsh but with current main this generates 1670 points of which only 1662 are unique.

import dmsh
import numpy as np

geo = dmsh.Difference(
    dmsh.Rectangle(0.0, 2.2, 0.0, 0.41),
    dmsh.Circle([.2, .2], 0.05)
)

points, triangles = dmsh.generate(geo, 0.025, tol=1e-9)

tmp = np.ascontiguousarray(points)  # test from skfem.Mesh._validate
if points.shape[0] != np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1])).shape[0]:
    print(points.shape[0], np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1])).shape[0])

[gdmcbain]

I should add that this example was working (no duplicate points) last year; I only noticed it broken a couple of weeks ago. It still works with

pip install -U dmsh==0.2.2 meshplex==0.13.6

(Not pinning meshplex raises AttributeError: 'MeshTri' object has no attribute 'node_coords'.)

[nschloe]

Intermediate note: This happens because there are points which are not part of any cell:

is_part_of_cell = np.zeros(len(points), dtype=bool)
is_part_of_cell[cells.flat] = True
assert np.all(is_part_of_cell)

Kicking out those points should work for you while I'm working on a solution. (Note that the cells array has to be updated, too.)

[nschloe]

Fixed now.