3d implicit surface optimization is hard to follow due numpy construction

[kayarre]

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

• [x] specify your optimesh version (optimesh -v), optimesh 0.6.3
• [x] specify what command you were running: optimesh.cvt.quasi_newton_uniform_blocks, and subsequently functions in helpers.py
• [x] attach the mesh (or even better: a small mesh that reproduces the problem).

I am trying to debug the following error:

/home/krs/anaconda3/envs/meshio/lib/python3.7/site-packages/optimesh/helpers.py:141: RuntimeWarning: invalid value encountered in multiply
  mesh.node_coords -= (grad * (fval / grad_dot_grad)[:, None]).T
Traceback (most recent call last):
  File "gen_sphere.py", line 579, in <module>
    main()
  File "gen_sphere.py", line 568, in main
    implicit_surface=star_object(fname),
  File "/home/krs/anaconda3/envs/meshio/lib/python3.7/site-packages/optimesh/cvt/_block_diagonal.py", line 99, in quasi_newton_uniform_blocks
    method_name="Centroidal Voronoi Tesselation (CVT), uniform density, block-diagonal variant"
  File "/home/krs/anaconda3/envs/meshio/lib/python3.7/site-packages/optimesh/helpers.py", line 141, in runner
    mesh.node_coords -= (grad * (fval / grad_dot_grad)[:, None]).T
ValueError: operands could not be broadcast together with shapes (1002,3) (3,1002) (1002,3) 

Basically there is some inconsistency with knowing the size of the numpy arrays, having a hard time figuring out where the issue is.

here is the producing code, I need to simplify, notice you can download the files input files from: ftp://ftp.nist.gov/pub/bfrl/garbocz/Particle-shape-database/

import numpy as np
import pandas as pd
from scipy import special
import meshio
import pooch
import os
import optimesh
import meshzoo

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]
    # )  # 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
    ptsnew[:, 2] = np.arctan2(xyz[:, 1], xyz[:, 0])  # phi
    return ptsnew

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])

    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
        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)
        for idx in range(self.n.shape[0]):
            r_theta += (
                self.m[idx]
                * 1.0j
                * self.coeff[idx]
                * special.sph_harm(self.m[idx], self.n[idx], theta, phi)
            )
            f_nm = (
                ((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])))
            ) ** (0.5)
            sin_phi = np.sin(phi)
            pre_coef = np.divide(
                -self.coeff[idx] * f_nm,
                sin_phi,
                out=np.zeros(phi.shape, dtype=np.complex),
                where=sin_phi != 0,
            )
            r_phi += (
                pre_coef
                * (
                    (
                        (self.n[idx] + 1)
                        * np.cos(phi)
                        * special.lpmv(self.m[idx], self.n[idx], np.cos(phi))
                    )
                    - (
                        (self.n[idx] - self.m[idx] + 1)
                        * special.lpmv(self.m[idx], self.n[idx] + 1, np.cos(phi))
                    )
                )
            ) * np.exp(1.0j * self.m[idx] * theta)

        # S = r * (
        #     r_theta ** 2.0 + r_phi ** 2.0 * np.sin(phi) ** 2.0 + r ** 2.0 * np.sin(phi)
        # ) ** (0.5)

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

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

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

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

        return coords

    # scalar
    def f(self, x):
        sph_coord = self.Cart_to_Spherical_np(x)
        # note the switch to math notation
        radius = self.get_radius(sph_coord[2], sph_coord[1])
        f_res = radius ** 2 - (x[0] ** 2 + x[1] ** 2 + x[2] ** 2)
        return f_res

    def grad(self, x):
        if x.shape[1] != 3:
            if x.shape[0] == 3:
                print("input is ambiguous")
            print("there is a weird issue with the grad operator")
            x = x.T
        sph_coord = self.Cart_to_Spherical_np(x)

        radius = self.get_radius(sph_coord[:, 2], sph_coord[:, 1])

        der_ = self.get_derivatives(radius, sph_coord[:, 2], sph_coord[:, 1])
        # print(der_.shape)
        # quit()
        return der_

    def get_file_path(self):
        return self.file_path

    def set_file_path(self, file_path_str):
        self.file_path = file_path_str

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

def convert_nodes(nodes, star_data):
    sph = Cart_to_Spherical_np(nodes)

    r = np.zeros(sph[:, 0].shape, dtype=complex)
    for idx in range(star_data.n.shape[0]):

        # shift phi to be 0 <= phi <= pi
        r += star_data.coeff[idx] * special.sph_harm(
            star_data.m[idx], star_data.n[idx], sph[:, 2], sph[:, 1] + np.pi / 2.0
        )
    r_real = r.real
    # print(np.abs(r), r.real)
    # scalar * [nx3] * [n,] = [nx3]
    nodes *= r_real[:, None]

    return nodes

class Sphere:
    def f(self, x):
        return 1.0 - (x[0] ** 2 + x[1] ** 2 + x[2] ** 2)

    def grad(self, x):
        return -2 * x

def main():

    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()

    # points, cells = uv_sphere(20)
    points, cells = meshzoo.icosa_sphere(10)
    conv_pts = convert_nodes(points, sand) / 100.0
    # print(points, cells)
    tri_cells = [("triangle", cells)]

    mesh = meshio.write_points_cells("test.vtk", conv_pts, tri_cells)
    print("made sand surface")

    test = star_object(fname)
    print(points.shape)
    s_pts = test.f(points)
    s_grad = test.grad(points)
    print(s_pts, s_grad.shape)
    # quit()
    # You can use all methods in optimesh:
    # points, cells = optimesh.cpt.fixed_point_uniform(
    # points, cells = optimesh.odt.fixed_point_uniform(
    points_opt, cells_opt = optimesh.cvt.quasi_newton_uniform_blocks(
        conv_pts,
        cells,
        1.0e-2,
        10,
        verbose=True,
        implicit_surface=star_object(fname),
        # 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")

if __name__ == "__main__":
    main()

[kayarre]

Alrighty, So I got it past this issue by determining that the implicit function expects [3xn] shaped input and output for the f and grad operator. I was able to remedy the issue in my code.

I think now my implicit function is degenerate at the poles.