从隐式函数生成网格

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)即可正常工作。

段错误指向 CGAL 中的一个问题。也许在那里提交错误很有用。