OF5版本以后重心坐标与绝对坐标的相互转换的问题。

[上级]

最近在研读OpenFOAM中particle类的barycentric coordinates 和global coordinates转换矩阵的程序，遇到一点问题，现总结出来和大家讨论一下。
OF5.0以后对颗粒追踪采用了barycentric方法，目的是更好的捕捉颗粒轨迹与tetrahedra网格面的交点，四面体内任意一点$\mathbf{p}$的global coordinates 可以表示为barycentric coordinates (${\lambda }_i$)与顶点的乘积和：
$$\begin{array}{}\text{(1)}& \mathbf{p}=\sum _{i=1}^4{\lambda }_i{\mathbf{v}}_{\mathbf{i}}\end{array}$$
其中${\mathbf{v}}_{\mathbf{i}}$是四面体顶点坐标。
示意图如下：

如上图，四面体中有一点$\mathbf{p}$，则$\mathbf{p}$点的重心坐标（${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4$）的计算方法为点$\mathbf{p}$将四面体分成的四个小四面体与大四面体的体积之比:
$$\begin{array}{}\text{(2)}& {\lambda }_1=\frac{V_{\mathrm{p}\mathrm{b}\mathrm{c}\mathrm{d}}}{V_{\mathrm{a}\mathrm{b}\mathrm{c}\mathrm{d}}}=1-{\lambda }_2-{\lambda }_3-{\lambda }_4,\end{array}$$
$$\begin{array}{}\text{(3)}& {\lambda }_2=\frac{V_{\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{d}}}{V_{\mathrm{a}\mathrm{b}\mathrm{c}\mathrm{d}}}=\frac{6V_{\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{d}}}{6V_{\mathrm{a}\mathrm{b}\mathrm{c}\mathrm{d}}}=\frac{\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})}{\mathbf{a}\mathbf{b}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})};,\end{array}$$
$$\begin{array}{}\text{(4)}& {\lambda }_3=\frac{V_{\mathrm{p}\mathrm{a}\mathrm{b}\mathrm{d}}}{V_{\mathrm{a}\mathrm{b}\mathrm{c}\mathrm{d}}}=\frac{6V_{\mathrm{p}\mathrm{a}\mathrm{b}\mathrm{d}}}{6V_{\mathrm{a}\mathrm{b}\mathrm{c}\mathrm{d}}}=\frac{\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{d}\times \mathbf{a}\mathbf{b})}{\mathbf{a}\mathbf{b}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})}，\end{array}$$
$$\begin{array}{}\text{(5)}& {\lambda }_4=\frac{V_{\mathrm{p}\mathrm{a}\mathrm{b}\mathrm{c}}}{V_{\mathrm{a}\mathrm{b}\mathrm{c}\mathrm{d}}}=\frac{6V_{\mathrm{p}\mathrm{a}\mathrm{b}\mathrm{c}}}{6V_{\mathrm{a}\mathrm{b}\mathrm{c}\mathrm{d}}}=\frac{\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{b}\times \mathbf{a}\mathbf{c})}{\mathbf{a}\mathbf{b}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})}，\end{array}$$
其中 $V_{\mathrm{p}\mathrm{b}\mathrm{c}\mathrm{d}}$, $V_{\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{d}}$, $V_{\mathrm{p}\mathrm{a}\mathrm{b}\mathrm{d}}$, 和 $V_{\mathrm{p}\mathrm{a}\mathrm{b}\mathrm{c}}$ 是小四面体的体积, $\mathbf{a}\mathbf{p}=\mathbf{p}-\mathbf{a}$, $\mathbf{a}\mathbf{b}=\mathbf{b}-\mathbf{a}$, $\mathbf{a}\mathbf{c}=\mathbf{c}-\mathbf{a}$, $\mathbf{a}\mathbf{d}=\mathbf{d}-\mathbf{a}$.
定义从barycentric coordinates到global coordinates的转换矩阵$\mathbf{A}$，则
$$\begin{array}{}\text{(6)}& \begin{array}{r}\mathbf{p}=\mathbf{A}\mathit{\lambda }\left(\mathbf{p}\right)\to \left[\begin{array}{c}{x}_{\mathrm{p}}\\ y_{\mathrm{p}}\\ z_{\mathrm{p}}\end{array}\right]=\left[\begin{array}{cccc}{x}_{\mathrm{a}}& x_{\mathrm{b}}& x_{\mathrm{c}}& x_{\mathrm{d}}\\ y_{\mathrm{a}}& y_{\mathrm{b}}& y_{\mathrm{c}}& y_{\mathrm{d}}\\ z_{\mathrm{a}}& z_{\mathrm{b}}& z_{\mathrm{c}}& z_{\mathrm{d}}\end{array}\right]\left[\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ {\lambda }_3\\ {\lambda }_4\end{array}\right];.\end{array}\end{array}$$
其中${\lambda }_i$就是$\mathbf{p}$的重心坐标，矩阵$\mathbf{A}$就是由四面体各个顶点的global coordinates组成的坐标矩阵。
既然有了正向变换，同样需要逆向变换，定义一个逆向变换矩阵$\mathbf{T}=[\mathbf{b}\mathbf{d}\times \mathbf{b}\mathbf{c}\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d}\mathbf{a}\mathbf{d}\times \mathbf{a}\mathbf{b}\mathbf{a}\mathbf{b}\times \mathbf{a}\mathbf{c}]$将global coordinates转换为barycentric coordinates:
$$\begin{array}{}\text{(7)}& \mathit{\lambda }=\frac{\mathbf{a}\mathbf{p}\cdot \mathbf{T}}{\mathrm{d}\mathrm{e}\mathrm{t}\mathbf{A}}=\frac{[\mathbf{a}\mathbf{p}\cdot (\mathbf{b}\mathbf{d}\times \mathbf{b}\mathbf{c})\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{d}\times \mathbf{a}\mathbf{b})\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{b}\times \mathbf{a}\mathbf{c})]}{\mathbf{a}\mathbf{b}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})},\end{array}$$
公式$\text{(7)}$的四个组成成分其实就是上面的公式$\text{(2)}$，$\text{(3)}$，$\text{(4)}$，$\text{(5)}$。$\mathrm{d}\mathrm{e}\mathrm{t}\mathbf{A}$是矩阵$\mathbf{A}$的行列式，表示为$\mathbf{a}\mathbf{b}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})$。
至此，barycentric coordinates与global coordinates的正向变换和逆向变化都已经定义完成了，上面的分析是基于我看OF6的代码以及查阅相关的文献总结出来的。
但是有一个问题困扰了我很久了 ，那就是逆向变换矩阵的第一个元素是：
$$\begin{array}{}\text{(8)}& {\mathit{\lambda }}_1=\frac{\mathbf{a}\mathbf{p}\cdot (\mathbf{b}\mathbf{d}\times \mathbf{b}\mathbf{c})}{\mathbf{a}\mathbf{b}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})},\end{array}$$
这个和公式$\text{(2)}$的定义不同，因为分子上算的不是图中小四面体pbcd的体积，就这儿我无法理解为什么OF里要这么定义，并且我看OF5-OF8里面都是这么定义的，我觉得有点问题，所有发出来和大家讨论一下，我觉得我自己是解决不了这个问题了

[李东岳]

简直是大作！！我一直在想研究下这方面工作没倒出空来。我先留言学习一下，在讨论

感谢分享！！！！！！！

[李东岳]

首先跟随你的思想：按照理论，应该是：
$$\begin{array}{}\text{(9)}& \mathit{\lambda }=\frac{\mathbf{a}\mathbf{p}\cdot \mathbf{T}}{\mathrm{d}\mathrm{e}\mathrm{t}\mathbf{A}}=\frac{[\mathbf{b}\mathbf{p}\cdot (\mathbf{b}\mathbf{d}\times \mathbf{b}\mathbf{c})\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{d}\times \mathbf{a}\mathbf{b})\mathbf{a}\mathbf{p}\cdot (\mathbf{a}\mathbf{b}\times \mathbf{a}\mathbf{c})]}{\mathbf{a}\mathbf{b}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})},\end{array}$$
同时
$$\begin{array}{}\text{(10)}& \mathbf{b}\mathbf{p}\ne \mathbf{a}\mathbf{p}\end{array}$$
所以代码跟理论对不上？

[上级]

@东岳 是的李老师，我的疑惑就在这儿，$\mathbf{b}\mathbf{p}\ne \mathbf{a}\mathbf{p}$，和代码对不上，细读代码发现不一样。

[李东岳]

问题解决了么？

[上级]

仔细读了代码，现在还剩最后一点小问题，其实说到底还是原来的问题，barycentric追踪是通过bary displacement元素的正负来判断颗粒运动轨迹与四面体网格面的交点的，例如这个图形

颗粒$\mathbf{p}$运动到$\mathbf{n}$，轨迹$\mathbf{p}\mathbf{n}$与面$\underset{\_}{\mathrm{b}\mathrm{c}\mathrm{d}}$相交于点$\mathbf{s}$，那么矢量$\mathbf{p}\mathbf{n}$的barycentric displacement的计算方法为：
$$\begin{array}{}\text{(11)}& \begin{array}{rlr}{\mathit{\lambda }}_{\mathbf{p}\to \mathbf{n}}& =\mathit{\lambda }\left(\mathbf{n}\right)-\mathit{\lambda }\left(\mathbf{p}\right)=\frac{\mathbf{a}\mathbf{n}\cdot \mathbf{T}}{\mathrm{d}\mathrm{e}\mathrm{t}\mathbf{A}}-\frac{\mathbf{a}\mathbf{p}\cdot \mathbf{T}}{\mathrm{d}\mathrm{e}\mathrm{t}\mathbf{A}}=\frac{\mathbf{p}\mathbf{n}\cdot \mathbf{T}}{\mathrm{d}\mathrm{e}\mathrm{t}\mathbf{A}}& =\frac{[\mathbf{p}\mathbf{n}\cdot (\mathbf{b}\mathbf{d}\times \mathbf{b}\mathbf{c})\mathbf{p}\mathbf{n}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})\mathbf{p}\mathbf{n}\cdot (\mathbf{a}\mathbf{d}\times \mathbf{a}\mathbf{b})\mathbf{p}\mathbf{n}\cdot (\mathbf{a}\mathbf{b}\times \mathbf{a}\mathbf{c})]}{\mathbf{a}\mathbf{b}\cdot (\mathbf{a}\mathbf{c}\times \mathbf{a}\mathbf{d})};.\end{array}\end{array}$$
此时barycentric dispalcement的符号为(-,+++)，可以判断轨迹与面$\underset{\_}{\mathrm{b}\mathrm{c}\mathrm{d}}$相交，同根根据$\mathbf{p}$和$\mathbf{n}$的barycentric coordinates可以求出点$\mathbf{s}$的barycentric coordinates，再转换为点$\mathbf{s}的$gobal coordinates，然后进如下一个四面体网格追踪。
上面我只是简述了颗粒在一个四面体网格内的追踪过程，从计算barycentric dispalcement的角度来看，逆向变换矩阵$\mathbf{T}$是没有问题的，但是问题还是在于从barycentric coordinates to global coordinates的过程，矩阵$\mathbf{T}$的第一个元素有问题。如果把$\text{???}$变为${\mathit{\lambda }}_1=1-{\mathit{\lambda }}_2-{\mathit{\lambda }}_3-{\mathit{\lambda }}_4$就好了，但是程序里没有这么写，说到底还是原来的问题 如果李老师有时间的话可不可以读一下particle类，我觉得您可以帮我释疑

[上级]

上面那个公式引用是公式（8）

[BlookCFD]

（7）这个公式确实是错的，但是particle里面都用的是(11)这个公式，用来计算“位移”的barycentric坐标，而不是点的barycentric坐标，OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H有计算点的barycentric坐标，你可以比较下。

[上级]

@BlookCFD 谢谢您的指点，我有看了一眼程序，我理解错了，OF里面是通过计算barycentric displacement来获得barycentric coordinates的，我一直被src/lagrangian/basic/particle/particle.C里面的一句代码给误导了，现贴出来：

particle.C 1062-1096 OpenFOAM6
void Foam::particle::correctAfterInteractionListReferral(const label celli)
 {
     // Get the position from the barycentric data
     const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
 
     // Create some arbitrary topology for the supplied cell
     celli_ = celli;
     tetFacei_ = mesh_.cells()[celli_][0];
     tetPti_ = 1;
     facei_ = -1;
 
     // Get the reverse transform and directly set the coordinates from the
     // position. This isn't likely to be correct; the particle is probably not
     // in this tet. It will, however, generate the correct vector when the
     // position method is called. A referred particle should never be tracked,
     // so this approximate topology is good enough. By using the nearby cell we
     // minimize the error associated with the incorrect topology.
     coordinates_ = barycentric(1, 0, 0, 0);
     if (mesh_.moving())
     {
         Pair<vector> centre;
         FixedList<scalar, 4> detA;
         FixedList<barycentricTensor, 3> T;
         movingTetReverseTransform(0, centre, detA, T);
         coordinates_ += (pos - centre[0]) & T[0]/detA[0];
     }
     else
     {
         vector centre;
         scalar detA;
         barycentricTensor T;
         stationaryTetReverseTransform(centre, detA, T);
         coordinates_ += (pos - centre) & T/detA;
     }
 }

其中的这句代码

coordinates_ += (pos - centre) & T/detA;

其实前面已经先定义了

coordinates_ = barycentric(1, 0, 0, 0);

后面计算单点的barycentric coordinates其实还是用的barycentric displacement的概念来转换的。

困扰了我两周的问题终于解决了，现在异常开心，在此对@东岳 和@BlookCFD 表示万分感谢。