Interface-fitted polyhedra mesh based on Cartesian mesh

科研笔记

Algorithm

Given a cube domain $D$ which include the interface $\Gamma$ with level set function $\phi$, and a mesh size $h$.

Generate the Cartesian mesh $[node, elem]$ with size $h$. Here $node$ is a $N\times 3$ matrix, $node(i,:)$ are the three coordinates of $i$-th node; $elem$ is a $NC\times 8$ matrix, and $elem(j,:)$ are the row indices of the eight vertex nodes of $j$-th element in $node$ matrix, and represents $j$-th small cube.

Find the interface small cubes. Here the interface elements inlcude:
1. exist at least two vertex nodes $p$ and $q$ with opposite sign, namely $\phi(p)\phi(q)<0$.
2. exist at least three vertex nodes on the interface, namely the value of $\phi$ on these three vertex nodes are 0.

Classify the small cubes by exterior cube,interiror cube and interface cubes，and label them by 1, -1 and 0, respectively.

Take all of the interface cubes which composite a new hex mesh $H_{cutted}$, and find all the boundary faces of this new hex mesh. Notice that, these boundary faces are square faces, and they will be the faces of the polyhedra mesh we will construct to fit the interface.

Next let us try to find the triangle faces of the polyhedra mesh near the interface by Delaunay triangulation algorithm.

Firstly, find the cut edge and the cut points by interface $\Gamma$. Find all the edges $edge$ which is a $NE\times 2$ matrix, and $edge(k,:)$ are the row indices of the two end nodes of $k$-th edge. For a edge $(p, q)$, if $\phi(p)\phi(q) < 0$, we take it as cutted edge by interface $\Gamma$. Then one can find the cutted points of all cutted edge by bisection method.

Seccondly, find the auxiliary points. For a interface cube, if it has at least one square face whos two opposite vertices on the interface by lucky coincidence, we call the center point of this square face as auxiliary points, and put all the cutted points, all the auxiliary points and all vertex nodes of all interface cubes together as a new point set $P$. Then construct Delaunay mesh $DT$ on $P$. Here we add the auxiliary points to make sure the interface can be recovered certainly.

Next, we need to get rid of some unnecessary tets in $DT$, and just keep the tets in every interface cubes. Since $(node,elem)$ is a cartesian mesh, so we can number every small cubes in these hex mesh. More inmportantly, given a points in domain $D$, we can easily known wich small cubes it blong to. For a tet in $DT$, we compute its centroid point, and find the small cube element it blong to. If the small cube element is not a interface cubes, then get rid of this tet from $DT$.

The above process can't guarantee to get rid of all unnecessary tets in $DT$. So we need an another criterion, if four vertex nodes of a tet are in different cubes, we remove it from $DT$. One just need to compare the difference of the maximum and the minimum of the coordinates of the four vertex nodes on the same dimension, if the difference is bigger than size $h$, this tet is unnecessary.

Now we have a tet mesh whos tets are all in interface cubes, and we still call this tet mesh as $DT$ for convenience. Remember that, given a tet in $DT$, we know which small interface cube it blong to.

Next we split the tets in $DT$ to two part: exterior tet set $DT_E$ and interior tet set $DT_{I}$. For a tet in $DT$, if the minimum of the $\phi$ value of the four vertex nodes is -1, we put it into $DE_I$. And

For $DT_I$, find the triangle faces whos two neighbor tets on different small interface cubes, and the triangle faces whos three vertex nodes are on interface which are the interface triangles and their direction point outside of the interface. Pay attention that, here we don't find the boundary triangle faces of $DT_I$ which are not interface triangles. Recall that, we have find the boundary square faces when we find all the interface cubes.

For $DT_E$, we do the same thing like on $DT_I$.

In most cases, every interface cube is cutted into two part and every one is a polyhedron. And we need to set a unique ID for these two part. For the interior part, we set its ID by the small interface cube index, e.g., $j$; for the polyhedron in the exterior part, we set its ID by $j+NC$. So when you find the square or triangle faces, remember record the polyhedron ID it blong to.

Of course the IDs of the polyhedra mesh we get are not continuous. One can make it continuous easily.

Currently, there are some complex cases which can't deal with by above algorithm. For example, the small interface cube is split into more parts than two. We need do more work to make the algorithm more robust. One key problem, how to reconstruct the three parts if a small interface cube is cutted into three or more than three parts by interface.