从隐式函数生成网格

Question

我想使用Python（或Python 3）从隐式函数生成体积（3D）网格：

def func(x,y,z):
    q = 0.25
    mu = q/(1.+q)
    return -(1-mu)*pow(x*x + y*y + z*z,-1./2.)  - mu*pow(pow(x-1,2) + y*y + z*z,-1./2.) - 0.5*(pow(x-mu,2) + y*y) + 1.9023266381531847

此函数具有复杂的等值面，但是我想将曲面限制在x = -0.615和x = 1.4之间，并且y = -0.6，y = 0.6，但z方向没有限制，但是有趣的部分在z = +/- 1之间。

我已经尝试过pygalmesh，但无法获得他们的示例来适应我的功能。它使我的Python内核崩溃而无输出。是否可以获取pygalmesh来执行此操作？如果没有，有什么更好的方法？

Answer 1

仅作记录，没有输出对我来说不会崩溃：

import numpy
import pygalmesh


class GrimReaper(pygalmesh.DomainBase):
    def __init__(self):
        super().__init__()

    def eval(self, x):
        q = 0.25
        mu = q / (1.0 + q)
        x, y, z = x
        return (
            -(1 - mu) / numpy.sqrt(x ** 2 + y ** 2 + z ** 2)
            - mu / numpy.sqrt((x - 1) ** 2 + y ** 2 + z ** 2)
            - 0.5 * ((x - mu) ** 2 + y ** 2)
            + 1.9023266381531847
        )

    def get_bounding_sphere_squared_radius(self):
        return 2.0


d = GrimReaper()
mesh = pygalmesh.generate_mesh(d, cell_size=0.1)

mesh.write("out.xmf")

Inserting protection balls...
  refine_balls = true
  min_balls_radius = 0
  min_balls_weight = 0
insert_corners() done. Nb of points in triangulation: 0
insert_balls_on_edges() done. Nb of points in triangulation: 0
refine_balls() done. Nb of points in triangulation: 0

construct initial points (nb_points: 12)
s.py:17: RuntimeWarning: divide by zero encountered in double_scalars
  + 1.9023266381531847
12/12 initial point(s) found...
Start surface scan...Scanning triangulation for bad facets (sequential) - number of finite facets = 50...
Number of bad facets: 0
scanning edges (lazy)
scanning vertices (lazy)
end scan. [Bad facets:0]

Refining Surface...
Legend of the following line: (#vertices,#steps,#facets to refine,#tets to refine)
(12,0,0,0)
Total refining surface time: 1.90735e-05s

Start volume scan...Scanning triangulation for bad cells (sequential)... 20 cells scanned, done.
Number of bad cells: 1
end scan. [Bad tets:1]

Refining...
Legend of the following line: (#vertices,#steps,#facets to refine,#tets to refine)
(23,11,0,70) (18839.3 vertices/s)Segmentation fault (core dumped)

没有术语- 0.5 * ((x - mu) ** 2 + y ** 2)，它可以正常工作。

segfault指向CGAL中的一个问题。可能在此处提交错误很有用。