stuck in mesh.flip_until_delaunay()

[kayarre]

If you're having problems optimizing a mesh, remember to

• [x] specify your optimesh version (optimesh -v): 0.6.3
• [x] specify what command you were running: quasi_newton_uniform_full
• [x] attach the mesh (or even better: a small mesh that reproduces the problem): portion of code and mesh attached.

I finagled to the point that I can get it cvt to start running but it appears stuck in the flip_until_delaunay() region.

I know There could be an issue with the gradient and normals, especially since I am not super confident I have the derivatives correct or the correct catches near the poles.

Maybe could you answer what flip_until_delaunay does? and maybe how to short circuit it?

import numpy as np  # from numpy import mgrid, empty, sin, pi
import pandas as pd
from scipy import special
import meshio
import pooch
import os
import optimesh

class star_object(object):
    def __init__(self, file_path=None):  # , res=10):
        # res = 4 if res < 4 else res  # ternary
        # self.radius = 1.0
        # self.center = [0.0, 0.0, 0.0]
        # self.thetaResolution = int(res)
        # self.phiResolution = int(res)
        # self.startTheta = 0.0
        # self.endTheta = 360.0
        # self.startPhi = 0.0
        # self.endPhi = 180.0
        # self.LatLongTessellation = False
        # self.output = vtk.vtkPolyData()
        # self.tol = 1.0E-8
        if file_path != None:
            self.file_path = file_path
            self.read_file()
        else:
            self.file_path = "/home/krs/code/python/Tools/mesh/c109-20001.anm"
        # self.read_file()

    def read_file(self, file_path=None):
        if file_path != None:
            self.file_path = file_path

        with open(self.file_path, mode="r") as f:
            data = pd.read_csv(f, sep="\s+", names=["n", "m", "a", "aj"])
            # print(data.head())
        self.n = data["n"].to_numpy()
        self.m = data["m"].to_numpy()
        # print(n.shape[0])
        self.coeff = np.empty((self.n.shape[0]), dtype=complex)
        self.coeff.real = data["a"].to_numpy()
        self.coeff.imag = data["aj"].to_numpy()
        # print(coeff[0])
        self.set_norm_coeff()
        # normalizing the coefficients
        # self.set_coeff_to_norm()

    def set_norm_coeff(self):
        self.norm_coeff = self.coeff / self.coeff[0]

    def set_coeff_to_norm(self):
        self.coeff = self.norm_coeff

    def cart_to_spherical_np(self, xyz):
        # physics notation
        ptsnew = np.zeros(xyz.shape)
        xy = xyz[0] ** 2 + xyz[1] ** 2  # r

        ptsnew[0] = np.sqrt(xy + xyz[2] ** 2)
        ptsnew[1] = np.arctan2(  # theta
            np.sqrt(xy), xyz[2]
        )  # for elevation angle defined from Z-axis down
        # ptsnew[:, 1] = np.arctan2(  # theta
        #     xyz[:, 2], np.sqrt(xy)
        # )  # for elevation angle defined from XY-plane up
        # print(ptsnew[1])
        # quit()
        ptsnew[2] = np.arctan2(xyz[1], xyz[0])  # phi

        return ptsnew

    # scalar
    def get_radius(self, theta, phi):
        radius = np.zeros(theta.shape, dtype=np.complex)
        for idx in range(self.n.shape[0]):
            radius += self.coeff[idx] * special.sph_harm(
                self.m[idx], self.n[idx], theta, phi
            )
        return radius.real

    def get_derivatives(self, r, theta, phi):
        radius = np.zeros(theta.shape, dtype=np.complex)
        r_theta = np.zeros(theta.shape, dtype=np.complex)
        r_phi = np.zeros(theta.shape, dtype=np.complex)
        Pnm_phi = np.zeros(theta.shape, dtype=np.complex)
        for idx in range(self.n.shape[0]):

            r_theta += (
                1.0j
                * self.m[idx]
                * self.coeff[idx]
                * special.sph_harm(self.m[idx], self.n[idx], theta, phi)
            )

            f_nm = np.sqrt(
                ((2.0 * self.n[idx] + 1) * special.factorial(self.n[idx] - self.m[idx]))
                / (4.0 * np.pi * (special.factorial(self.n[idx] + self.m[idx])))
            )

            pre_coef = self.coeff[idx] * f_nm * np.exp(1.0j * self.m[idx] * theta)

            abs_m = abs(self.m[idx])

            if self.n[idx] == 0 and self.m[idx] == 0:
                # r_phi += pre_coef * 0.0
                pass
            elif self.n[idx] == -self.m[idx]:
                Pnm_phi = (
                    0.5
                    * np.power(-1, abs_m)
                    / (2.0 * special.factorial(2.0 * self.n[idx] - 1))
                    * special.lpmv(self.n[idx] - 1, self.n[idx], np.cos(phi))
                )
            elif self.n[idx] == self.m[idx]:
                Pnm_phi = self.n[idx] * special.lpmv(
                    self.n[idx] - 1, self.n[idx], np.cos(phi)
                )
            elif self.m[idx] == 0:
                Pnm_phi = -special.lpmv(1, self.n[idx], np.cos(phi))
            else:
                Pnm_phi = 0.5 * (
                    (self.n[idx] + self.m[idx])
                    * (self.n[idx] - self.m[idx] + 1)
                    * special.lpmv(self.m[idx] - 1, self.n[idx], np.cos(phi))
                    - special.lpmv(self.m[idx] + 1, self.n[idx], np.cos(phi))
                )

            phi_has_nan = np.isnan(pre_coef).any()
            if phi_has_nan:
                raise ValueError(" output has nans")
            array_has_nan = np.isnan(Pnm_phi).any()
            if array_has_nan:
                raise ValueError(" output has nans")
            r_phi += pre_coef * Pnm_phi

        S = np.absolute(
            r
            * (
                (r_theta ** 2.0)
                + (r_phi ** 2.0) * (np.sin(phi) ** 2.0)
                + (r ** 2.0) * np.sin(phi)
            )
            ** (0.5)
        )
        t_idx = np.where(S < 0.0000001)
        array_has_nan = np.isnan(S).any()
        if array_has_nan:
            raise ValueError("output has nans")

        n_x = r * (
            r_theta * np.sin(theta)
            - r_phi * np.sin(phi) * np.cos(phi) * np.cos(theta)
            + r * np.square(np.sin(phi)) * np.cos(theta)
        )  # / S

        n_y = r * (
            -r_theta * np.cos(theta)
            - r_phi * np.sin(phi) * np.cos(phi) * np.sin(theta)
            + r * np.square(np.sin(phi)) * np.sin(theta)
        )  # / S

        # print(n_y)

        n_z = r * (
            r_phi * np.square(np.sin(phi)) + r * np.cos(phi) * np.sin(phi)
        )  # / S

        coords = np.stack((n_x.real, n_y.real, n_z.real), axis=0)

        array_has_nan = np.isnan(coords).any()
        if array_has_nan:
            raise ValueError("output has nans")
        test_norm = np.sqrt(np.sum(np.square(coords), axis=0))
        south_pole = np.where(phi == np.pi)
        north_pole = np.where(phi == 0.0)
        coords[:, south_pole[0]] = np.array([[0.0], [0.0], [-1.0]])
        coords[:, north_pole[0]] = np.array([[0.0], [0.0], [1.0]])
        test_norm[south_pole] = 1.0
        test_norm[north_pole] = 1.0
        return coords   / test_norm

def Cart_to_Spherical_np(xyz):
    # physics notation
    ptsnew = np.zeros(xyz.shape)  # np.hstack((xyz, np.zeros(xyz.shape)))
    xy = xyz[0] ** 2 + xyz[1] ** 2  # r
    ptsnew[0] = np.sqrt(xy + xyz[2] ** 2)
    ptsnew[1] = np.arctan2(  # theta
        np.sqrt(xy), xyz[2]
    ) 
    ptsnew[2] = np.arctan2(xyz[1], xyz[0])  # phi
    return ptsnew

def convert_nodes(nodes, star_data):
    sph = Cart_to_Spherical_np(nodes.T)
    r = np.zeros(sph.shape[-1], dtype=np.complex)
    for idx in range(star_data.n.shape[0]):
        r += star_data.coeff[idx] * special.sph_harm(
            star_data.m[idx], star_data.n[idx], sph[2], sph[1]
        )
    r_real = r.real
    nodes *= r_real[:, None]

    return nodes

class call_back_test(object):
    def __init__(self, max_num_steps):
        self.avg_quality = np.empty((max_num_steps + 1))

    def test_callback(self, k, mesh):
        self.avg_quality[k] = np.average(mesh.cell_quality)
        print(k)
        return

def setup_pooch():
    # Define the Pooch exactly the same (urls is None by default)
    GOODBOY = pooch.create(
        path=pooch.os_cache("particles"),
        base_url="ftp://ftp.nist.gov/pub/bfrl/garbocz/Particle-shape-database/",
        version="0.0.1",
        version_dev="master",
        registry=None,
    )
    # If custom URLs are present in the registry file, they will be set automatically
    GOODBOY.load_registry(os.path.join(os.path.dirname(__file__), "registry.txt"))

    return GOODBOY

points, cells = meshzoo.icosa_sphere(10)

pooch_particles = setup_pooch()
fname = pooch_particles.fetch("C109-sand/C109-20002.anm")
print(fname)
sand = star_object()
sand.set_file_path(fname)
sand.read_file()

conv_pts = convert_nodes(points, sand)

max_steps = 10
q_check = call_back_test(max_num_steps=max_steps)
points_opt, cells_opt = optimesh.cvt.quasi_newton_uniform_full(
    conv_pts,
    cells,
    1.0e-2,
    max_steps,
    omega=0.9,
    verbose=True,
    callback=q_check.test_callback,
    uniform_density=True,
    implicit_surface=star_object(fname),
    implicit_surface_tol=1.0e-8
    # step_filename_format="out{:03d}.vtk"
)

tri_cells_opt = [("triangle", cells_opt)]

mesh = meshio.write_points_cells("test_opt.vtk", points_opt, tri_cells_opt)
print("optimized sand surface")

test.zip

[nschloe]

You can expect anyone to go through so many lines of code. I'll close this, please reopen with a MINIMAL example.