BUG: spatial: The mesh creates elements with zero angles in the Delaunay triangulation

[luiscamilosc]

Describe your issue.

Hi, I have a code in Python that read the points of the coordinates of my nodes and then, create the mesh expected in *.vtu format. However, some of the elements with small sizes (local refinements in my previous mesh) generated triangles with zero angles inside of the mesh. Do you know how can I fix it, is you need a short meeting to show you my problem lcsc@dhigroup.com

See this link with my problem and script: https://stackoverflow.com/questions/76082920/problem-mesh-creation-vtu-files-with-zero-angles-in-the-elements-created-base

Thanks for your amazing contribution!! Regards, lcsc

Reproducing Code Example

import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import os
from scipy.spatial import Delaunay

path = r"C:\Users\\"
filename = "old_mesh.csv"
datos = np.genfromtxt(path + filename, delimiter=',')

# delete X,Y first line
if np.array_equal(datos[0], np.array([0, 0])):
    datos = np.delete(datos, 0, axis=0)

#verifiy if there are NaN values
if np.isnan(datos).any():
    print("El archivo CSV contiene valores faltantes o NaN")
    # elimina las filas con valores faltantes o NaN
    datos = datos[~np.isnan(datos).any(axis=1)]

# creation of the mesh based in the nodes

tri = Delaunay(datos)

# vtu file creation 
with open(path + os.path.basename(filename)[0:-4] + '_outputmesh.vtu', "w") as f:
    f.write('<?xml version="1.0"?>\n')
    f.write('<VTKFile type="UnstructuredGrid" version="1" byte_order="LittleEndian" header_type="UInt64">\n')
    f.write('\t<UnstructuredGrid>\n')
    f.write('\t\t<Piece NumberOfPoints="{}" NumberOfCells="{}">\n'.format(len(datos), len(tri.simplices)))
    f.write('\t\t\t<Points>\n')
    f.write('\t\t\t\t<DataArray type="Float64" Name="Points" NumberOfComponents="3" format="ascii">\n')
    for i in range(len(datos)):
        f.write('\t\t\t\t\t{} {} {}\n'.format(datos[i][0], datos[i][1], 0))
    f.write('\t\t\t\t</DataArray>\n')
    f.write('\t\t\t</Points>\n')
    f.write('\t\t\t<Cells>\n')
    f.write('\t\t\t\t<DataArray type="Int64" Name="connectivity" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(tri.simplices)):
        f.write('\t\t\t\t\t{} {} {}\n'.format(tri.simplices[i][0], tri.simplices[i][1], tri.simplices[i][2]))
    f.write('\t\t\t\t</DataArray>\n')
    f.write('\t\t\t\t<DataArray type="Int64" Name="offsets" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(tri.simplices)):
        f.write('\t\t\t\t\t{}\n'.format((i+1)*3))
    f.write('\t\t\t\t</DataArray>\n')
    f.write('\t\t\t\t<DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(tri.simplices)):
        f.write('\t\t\t\t\t5\n')
    f.write('\t\t\t\t</DataArray>\n')
    f.write('\t\t\t</Cells>\n')
    f.write('\t\t</Piece>\n')
    f.write('\t</UnstructuredGrid>\n')
    f.write('</VTKFile>')
#  plot triangle mesh
fig, ax = plt.subplots()
ax.triplot(datos[:,0], datos[:,1], tri.simplices)
ax.set_aspect('equal')
plt.show()

Error message

The code run without problems, the problem is in the output generated

SciPy/NumPy/Python version and system information

1.10.1 1.24.2 sys.version_info(major=3, minor=9, micro=6, releaselevel='final', serial=0)
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "0.29.33",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.12.1",
      "include directory": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\pip-build-env-u63ta2f1\\overlay\\Lib\\site-packages\\pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "openblas",
      "found": true,
      "version": "0.3.18",
      "detection method": "pkgconfig",
      "include directory": "c:/opt/openblas/if_32/64/include",
      "lib directory": "c:/opt/openblas/if_32/64/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= PRESCOTT MAX_THREADS=4",
      "pc file directory": "c:/opt/openblas/if_32/64/lib/pkgconfig"
    },
    "lapack": {
      "name": "openblas",
      "found": true,
      "version": "0.3.18",
      "detection method": "pkgconfig",
      "include directory": "c:/opt/openblas/if_32/64/include",
      "lib directory": "c:/opt/openblas/if_32/64/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= PRESCOTT MAX_THREADS=4",
      "pc file directory": "c:/opt/openblas/if_32/64/lib/pkgconfig"
    }
  },
  "Python Information": {
    "path": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\cibw-run-a1px0t3e\\cp39-win_amd64\\build\\venv\\Scripts\\python.exe",
    "version": "3.9"
  }
}

[tylerjereddy]

I can't tell if there is an actual bug here just yet, but could you provide your data in a format other than a .png image? It looks like you've linked an image to the CSV file over at SO here: https://i.stack.imgur.com/hbZTX.png

But if the data is that small, it shouldn't be too much work to copy-paste it into your original comment above I think.

It may also be helpful to assert on some property of the data you believe should be true, but that is not, to facilitate with reproduction outside of some other plotting/program.

[luiscamilosc]

Dear @tylerjereddy thank you so much for your help. The data is long so I will load it here ;)

My real problem is described in the following steps:

1. I have multiple triangle meshes with global and local refinements per unstructured mesh. I need to merge them in a single file and keep the same number of nodes of all of them and the elements (cells) created as well.
2. I merge my *.csv files into one file and delete the repeat nodes to don't have problems with the new triangulation. I can't delete the input data "nodes or vertices" due to a restriction problem of my expected output.
3. Then I use your amazing package and the data worked ~99.9% well. However, the problems appear when I have zones with an existing local refinement that have triangles with less than 1 m of separation of their vertices and the Delaunay function failed, creating elements with internal angles of zero or deleting nodes(vertices) from the input, for me is a bug or restriction of this method.

I appreciate it if you can support me in identifying the bug or guide me with a good solution using your methods.

I attached: to. [ 'new_mesh.csv', 'old_mesh.csv'] #the *.cvs files with X,Y coordinates of my vertices b. Meshing methodology.pdf a presentation to show you my goal and problem c. remesh_csv_vtu_feflow.py the scripting created for the execution of the process d. the polygon mesh based (old_mesh_antamine_poly.shp) to do a comparison with the output and the new one in the area of interest.

I appreciate to don't sharing this information with third parties for a confidential thing with the company that I'm envolved Please to download the information use this link, that will be available some days: https://bluewhale.us.dhigroup.com/rs?key=DxRTefT9Y8STMSE3Dz2GkLhpyiD9jAbM&encKey=V5RigXyx2AtJHFCdAGXDNveAxD2gBXQW

[rkern]

If the file is confidential, then please delete it. This is a public website.

[luiscamilosc]

Please to download the information use this link, that will be available some days: https://bluewhale.us.dhigroup.com/rs?key=DxRTefT9Y8STMSE3Dz2GkLhpyiD9jAbM&encKey=V5RigXyx2AtJHFCdAGXDNveAxD2gBXQW

[tylerjereddy]

It looks like this might benefit from a much smaller/focused/concise reproducer, or at least one that can generate its own data programmatically rather than needing to download a bunch of stuff.

[luiscamilosc]

Dear @tylerjereddy I shared with you the files that reproduce the error. My code merges two *csv files, then, create a .vtu file to see the consistency of the mesh and also, plots the mesh too.

I shared to you the code complete, if you prefer you can use the upper version or the attached version... then, you can reproduce the complete problem...

Also, I added a short presentation that shows the problem with small triangles generated for the Delaunay method producing bad quality and deformed elements for the triangulation (for me a bug in the execution)... If you can't reproduce the error tell me to explain you better... but the problem persists when the nodes are very close and the points probably have a colinearity between them.

[tylerjereddy]

Maybe we'll take a look, thanks. My point was that it would be really nice if you could generate everyting programmatically in a small Python script without needing external file assets/downloads/me to read a presentation deck.

[tylerjereddy]

colinearity between them

There's probably not much we can do about degeneracies, sometimes changing the Qhull flags can help a little. Also, such degeneracies are almost always possible to produce creatively in a Python/NumPy-based reproducer, using i.e., cooked up example meshes with intentional colinearity, or whatever degenenarcy you may have.

[rkern]

Delaunay uses the Qhull library underneath. Here is its FAQ on the subject. There is, unfortunately, not a whole lot that can be done about collinear points and still get a triangulation out.