0

I have a 3D triangulated closed surface. This surface is immersed inside a rectangular Cartesian grid. The surface is stored in a STL format. My goal is to compute the fraction of the cells which are cut by the surface. In other words, when a cell has an intersection with surface and in fact is cut by it, it is divided into an interior and exterior sub cells. My aim is to find the volume fraction of the interior sub-cell and the cell cut by the surface. For example, in the picture below which shows a 2D case for illustration purposes, I want to compute the fraction of the cyan area to the total area (which is cyan+grey).

enter image description here

Could someone help me find an efficient method/algorithm that can do this? I can do the implementation in either MATLAB or C.

afp_2008
  • 1,940
  • 1
  • 19
  • 46

3 Answers3

1

An exact method:

First consider intersecting the surface with the half-space delimited by a plane.

Some of the triangles are wholly inside the half-space and remain unchanged. Some others are wholly outside and are discarded. The remaining ones are cut in two, giving a triangle and a quadrilateral. Keep the right part and triangulate the quadrilateral if necessary.

You also need to consider the face generated as the region of the plane that lies inside the surface. For this, take all the new edges (those resulting from the section of the triangles with the plane), and chain them. I mean tie together the segments that have a common vertex (take care of numerical accuracy); in the end you will have one or more loops, forming simple polygons. Triangulate those polygons (by the ear clipping method for instance).

By this process, you obtain a new triangulated surface which describes the requested intersection. Repeat for the six planes that define a cell, and you get the intersection between the cell and the volume enclosed in the surface.

To compute its volume, you can sum the volumes of the prisms formed by every triangle and its projection on a fixed plane; make sure to compute these volumes algebraically (with a sign).

Here is a 2D analogy. The "new edges" are in green. enter image description here

As this process is time consuming, it will be costly if you need to repeat it for each cell of a whole grid. You can ease the computation by keeping intermediate results and working slice by slice.

A sort in the direction perpendicular to the first plane will accelerate the rejection of outside triangles.

  • When doing that for 3D surfaces (Sutherland-Hogdman reentrant clipping), the highlighted segments in the figure become polygons, and it is tedious to create them (need to maintain combinatorics), but since only the volume is needed, it is possible to create arbitrary triangles instead, see my article in http://alice.loria.fr/index.php/publications.html?redirect=0&Paper=CVTGraph_eg@2012 – BrunoLevy Oct 29 '15 at 17:33
0

Besides the exact solution in Yves Daoust's answer, if an approximate solution is acceptable, to have something that is reasonably efficient, I would "rasterize" the surface into a higher-definition grid, and count the fraction of micro-cell in each macro-cell. To "rasterize" the surface, I would use either the GPU [1] if performance is needed (but it is painful). If performance is not a real concern, then you can loop on each triangle, then for each micro-cell that falls in the bounding box of the triangle, test whether it has an intersection with the triangle. In the end, flood-fill the interior of the volume.

[1] https://developer.nvidia.com/content/basics-gpu-voxelization

BrunoLevy
  • 2,495
  • 17
  • 30
0

I believe you require this for the Volume of Fluid Interface reconstruction by PLIC method.The following are the steps generally followed by the VOF research community {ref: Rider, W. J., & Kothe, D. B. (1998). Reconstructing volume tracking. Journal of computational physics, 141(2), 112-152.}.

Step 1: Write/develop the equation to the interface. In PLIC, I write the interface of the form ax+by+c = 0. If you prefer quadratic/spline equations, do that.

Step 2: Get the coordinates of all the grid points. Eg: If you want to develop a grid such that x = 1:5 and y = 1:5, then use the following 

x = 1:5;
y = 1:5; 
[X,Y] = meshgrid(x,y);

** Step 3: Set up a tolerance value (I use 10^-7) to define regions inside/outside/on the interface.

For all points (X,Y) inside the fluid/liquid domain: 

a*X+b*Y+c >tol

All points (X,Y) outside the fluid/liquid domain

a*X+b*Y+c <-tol

Points on the interface: 

[-tol<=(a*X+b*Y+c)<=tol]
Vaidehi Oruganti
  • 1