generate_periodic_mesh usage with points and faces

[dthillaithevan]

Hi,

In cases where it's not possible to define the surfaces mathematically, is it possible to supply the points and triangles (faces) of a surface mesh to the generate_periodic_mesh function? I've tried to supply a .vtu as the first argument to the function but that didn't seem to work.

import numpy as np
import pygalmesh

#Generate points and faces defining 3D surface mesh
points, verts = ....
cells = [("triangle", faces)]
meshio.write_points_cells('test.vtu',points,cells)

mesh = pygalmesh.generate_periodic_mesh(
    'test.vtu',
    [0, 0, 0, 1, 1, 1],
)

Do you know if something like this would be possible?

Thanks!

[nschloe]

There is generate_volume_mesh_from_surface_mesh(), but it doesn't support periodic meshes.

[dthillaithevan]

Ok, thanks.

Do you know if there is any way to modify the eval function used in the periodic Schwarz example to accept a grid of points rather than a single coordinate?

At the moment the code appears to be looping over the 3D coordinates, x, one at a time which is slowing down the code significantly.

It would be great instead of eval accepting a single [3,] coordinate, it could accept a grid of [n,3] coordinates (n is the total number of points in the domain), then it would be possible to use vectorized functions within eval to compute the resulting surface.

[nschloe]

It would be great instead of eval accepting a single [3,] coordinate, it could accept a grid of [n,3] coordinates (n is the total number of points in the domain), then it would be possible to use vectorized functions within eval to compute the resulting surface.

I couldn't agree more.

At the moment the code appears to be looping over the 3D coordinates, x, one at a time which is slowing down the code significantly.

Indeed it does, and there's currently no other way. This is the relevant piece of code. The problem is that CGAL needs to have the points supplied individually.

I've opened an issue on CGAL on this (https://github.com/CGAL/cgal/issues/3874) where you could chime in, but I suppose we shouldn't hold our breaths for the change to happen.

[dthillaithevan]

Thanks for your quick response and for opening the issue on CGAL.

[lrineau]

(I am the main author/develop of the CGAL 3D mesh algorithm used by pygalmesh.)

The theoretical algorithm behind the CGAL implementation of the meshing process is intrinsically sequential. Its proof needs that: at each step, the point that is inserted in the mesh is the one with the biggest empty sphere centered on it.

But, the CGAL implementation can also be run in a parallel mode, using multi-threading, and gain that way almost as much as the number of physical CPU cores available in your computer. That is the best we can do.

... Well, almost! There are two modifications we could do:

• For an implicit function, like a distance function, the mesh domain class needs to compute the intersection between a segment and the iso-level of the function, and it proceeds by bisection along the segment. It could operate differently, using an higher level of dividing the segment into smaller parts. Or, if the distance function allows for first derivation, it could use a Newton algorithm to speed up the computaton.
• Each time a new vertex is inserted in the mesh, there is a lot of new Delaunay triangles that are created around the new vertex, and they are all processed sequentially, in a loop. They could be processed in parallel. That would mean that the mesh domain class would receive a bucket of segments queries to intersect with the surface, instead of just one at a time.

Would that help?

[nschloe]

Thanks @lrineau for the insights!

It appears that all boundary steps can be performed at the same time, so that'd be an easy improvement.

But, the CGAL implementation can also be run in a parallel mode, using multi-threading,

How does the logic work here? Don't you have to find one point after another, defeating all possibility for parallelization?

Would that help?

It depends on how much time is spent on each task. How many domain function evaluations have to be done for the vertex insertion vs. the boundary steps, on average?