![\[ u_t = L(u),\qquad\mathrm{in\ }\Omega \]](form_0.png)
Co-Workers:
移动网格方法的框架是选取一个参考区域，构造参考区域和实际问题的物理区域 之间的某种映射，将参考区域上的均匀网格映射成为物理区域上的非均匀网格。 一个被称为控制函数的量，作为我们构造这个区域间映射的参量，从而，通过调 整控制函数，我们就可以操纵网格。
移动网格方法中的一个主要问题就是网格会发生缠绕。我们研究的方法主要是使 用调和映射来构造网格映射，调和映射的优秀性质保证了网格在理论上的存在唯 一性。进一步我们引进一个迭代过程来实现调和映射，消除数值原因可能造成的 网格缠绕。网格移动的过程和方程求解完全分开。
对于解函数从旧网格到新网格上的更新，我们构造了一个来源于解法的全局插值 方法，使得网格的移动和方程原有的求解格式完美配合，但是在实现上又是完全 分开的，有利于代码重用，使得我们能够开发出来一个通用的网格移动程序包。 这样的整体插值方案可以看出事实上是求解了一个人工对流方程，在不可压流体 和间断Gerlerkin方法中的运用得到了新颖的结果。
由于很多问题的奇异性往往出现在边界上，所以边界上的网格的移动是一个比较 重要的问题。以前的方法常常是将边界的网格和内部的网格分开来处理。我们提 出了一个新的能量泛函，在极小化这个能量泛函的情况下，构造了一个调和映射 的推广，使用这个映射来构造参考区域和物理区域之间的网格变换，使得边界上 的网格和内部网格耦合移动，网格的质量得到进一步的提高，并能够顺利使用到 三维的情况。
加上合适的初值和边值。假设对于这个问题已经有了一个现成的求解格式为
![\[ u_h^{(n+1)} = u_h^{(n)} + \Delta t L_h(u_h^{(n)}) \]](form_1.png)
我们将不考虑此问题的求解方法的细节，仅仅是对其求解方法做一些必要的假定，并给出一些记号 以便与网格移动部分的算法相结合。问题的解
所处的空间记为
，求解方法将在有限维 空间
中来对进行逼近，空间依赖于区域
（我们假设这是一个多面体） 上的一个剖分，也就是说， 给定区域上满足某些特定要求的一个剖分，空间就可以根据事先确定的方式唯一 的给出。对于我们的数值例子具体的来说，需要先给定一个所谓的“逻辑区域”
（对应 于“逻辑区域”的名称，我们将叫做“物理区域”），在 上有一个确定的单纯形剖分（我们用到的剖分需要能够使节点自由移动而不破坏网格 的结构，象单纯形剖分，二维的任意四边形剖分，三维的任意六面体剖分，对于象矩形的剖分这样的 节点不具有自由度的剖分，就目前我们看来，是不能进行网格移动自适应的）记为
， 其节点记为
，单元 记为
，对于每个从区域到的同胚
，我们可以得到物理区域 上的剖分
，以
作为其节点，并且具有剖分 的节点间的拓扑关系。由于是一个同胚，从而有
， 于是我们可以保证确实是给出了一个剖分。相对应的，我们把这两个区域上的两套 剖分分别叫做“逻辑剖分”和”物理剖分“。移动网格所要研究的就是给出一个比较恰当的同胚 ，来作为区域变换函数，从而得到计算区域上的逻辑剖分和物理区域上的物理剖分之间的一个 对应。我们将这个同胚离散为一个分片线性的映射（如果是二维上的任意四边形剖分，这将会是一个 分片双线性映射，对于三维的任意六面体剖分，则是一个分片的三线性映射，这个离散方式随着要求 和单元的变化可以做不同的选取，但至少要保证在单元的边界上能使单元保持形状，然后我们要求 这个映射必须是只依赖于节点的位置的，也就是说，如果给定所有的节点的位置，那么这个映射是 唯一确定的，因为至少在目前，我们的网格移动最终是要依靠节点的移动来实现的，是否有其他的 实现网格移动的途径，我们还没有考虑）。 由于物理网格和区域之间的变换将会依赖于时间
，我们在以后的讨论中会 将它们都加上一个脚标，成为
和
，表示是时刻的物理网格和区域 变换，并用
表示的节点。相应的，我们将时刻的解空间也记为 
。于是我们在这个格式的基础上加入移动网格的部分后整个求解的算法的轮廓如下：
计算控制函数 
；
![\[ -\nabla_x \left( M \nabla_x \xi \right) = 0; \]](form_21.png)
；
；
调和映射(Harmonic Mapping)的理论在数学中是一个相对比较新的概念。1944年， Fuller 给出了调和映射的定义。直到1964年，Eell和Tompson's建立了调和映射的 奠基性理论工作以后，调和映射才在几何中得到了比较充分的重视。Hamilton， Schoen，Yau（丘成桐）等建立了比较完善的存在性和唯一性的结果。关于调和映 射的研究一方面是理论上的关于其存在性(Existence)，唯一性(Uniqueness)， 正则性(Regularity)的探讨，另一方面是关于调和映射在不同的数学领域中的应 用。我们只需要用到关于其存在性和物理意义的基本内容。由于这个概念首先来 自于 Rieman 几何学，我们需要用到一些几何的概念来描述调和映射的理论。给 定两个 
维Rieman流形 
和 
，他 们上面关于局部坐标系(Local Coordinates) 
和 
的 Rieman 度量分别为 
和 
。对给定的映射 
，我们定义一个能量密度(Energy Density)：
![\[ e(\xi) = \displaystyle\frac{1}{2} G^{ij} g_{\alpha \beta}{\displaystyle\frac {\partial \xi^\alpha} {\partial x^i}}{\displaystyle\frac {\partial \xi^\beta} {\partial x^j}} \]](form_32.png)
其中，我们引进标准求和约定（同一个指标如果同时在上标和下标中出现，则表 示对该指标从 
到 进行求和）。从而，映射 
的能量是能量密度的积分
![\[ E(\xi) = \int_\Omega e(\xi) dx. \]](form_35.png)
如果 
，能量就小于无穷。如果 是能量泛函 的极小点，那么 就是调和映射。能量泛函的Euler-Lagrange方程是
![\[ \Delta_{L-B} \xi^k + G^{ij} \Gamma^k_{\alpha \beta} {\displaystyle\frac {\partial \xi^\alpha} {\partial x^i}}{\displaystyle\frac {\partial \xi^\beta} {\partial x^j}}= 0 \]](form_37.png)
其中， 
是 Laplace-Beltrimi 算子， 
是 上的第二类 Christoffel 符号。定理保证，如果
的曲率是非正的；
关于 上的度量是凸的；
那么，对于任意的从 到 的同胚 
， 都存在唯一的一个调和映射 从 到 满足 
。这是 在 Dirichlet 边界条件下的存在唯一性的结果。Hamilton 还给出了在 Neumann 边界条件 和混合边界条件下的存在唯一性结果。在我们的应用中，逻辑区域 上的度量采用的就是 Euclid 度量，所以第一个条件自然是满足的，我们将逻辑区域选 为凸区域就是为了满足第二个条件。这导致我们只能够给出单连通区域上的数值例子， 对于多连通区域，只能够期望通过区域分解或者是其他特殊技巧来处理。
在 Eells 和 Thompson's 的研究中，给出了对任意给定的同胚，构造出收敛到调和 映射的同伦的方法，那就是通过解热传导方程
![\[ \xi_t = \Delta_{L-B} \xi^k + G^{ij} \Gamma^k_{\alpha \beta} {\displaystyle\frac {\partial \xi^\alpha} {\partial x^i}}{\displaystyle\frac {\partial \xi^\beta} {\partial x^j}} \]](form_43.png)
用 
作为初值。这个方法得到的同伦类将会一致的收敛到 调和映射，并且连同 
的梯度也会一致的收敛。由于对于我们 的特殊情形，也就是 上的度量是欧氏的情形， 
，上面的 Euler 方程成为一个线性椭圆型方程，我们不 需要用到这个方法。
我们来讨论具体的实现方法。首先我们将物理区域上事先给定一个一致剖分，它的 节点是
。在选定了逻辑区域以后，我们先确定在边界上的映射的情况。一般的， 我们选择的逻辑区域和物理区域有相同的边界结构，也就是说，对物理区域的每一个顶点，棱，和面 来说，逻辑区域的边界上都有相对应的元素，并且，这些元素之间的归属关系也都保持（但并非一定 要这样），从而我们可以方便的事先人工的给出边界上的映射。我们将这个边界上的映射记为给定 的一个同胚
在边界上的限制，从而通过解 Possion 方程组
![\[ \begin{array}{rcl} \Delta \xi^k &=& 0 \\ \xi^k|_{\partial \Omega} &=& \phi^k|_{\partial \Omega} \end{array} \]](form_49.png)
其中
，我们就可以得到一个逻辑网格。这个方程的解将上的Euclid度量变成上 的诱导度量
。我们可以看到， 这个度量在逻辑区域上导致了一个非均匀的网格，但是这个网格在单元上具有将诱导度量平均化的 作用，所以，我们可以将逻辑区域上的逻辑网格看作为物理区域上的物理网格的一个代替品，而 不会有太大的偏差。通过对相同的物理区域，不同的逻辑区域上生成的逻辑网格的比较，我们可以看出， 逻辑区域 的选取对于逻辑网格平均化诱导度量的影响是非常小的。这个事实说明，我们对逻辑区域的选取 并不会使逻辑网格对物理区域上的物理网格的逼近程度受到很大的影响，从而可以比较自由的选择 逻辑区域。我们选取不同的逻辑区域进行的计算结果没有明显差别的事实也说明，逻辑区域的选取 并不会对计算的效果造成明显的影响。 Huang 对这个问题也有一些讨论，并有相似的 观点。
逻辑区域是凸区域的要求使得我们的方法的应用范围受到很多限制，因为在欧氏度量下，凸区域 只能是单连通区域。我们现在仅仅能够考虑使用区域分解方法来对付多连通区域，我们可以在那些 事先知道不会有大幅度的网格调整的位置将区域分割成几个单连通区域，然后在每个单连通区域上 使用我们现在的移动网格方法，但是需要事先知道问题的一些性质，并且在区域之间 的边界上会引进一些人工的因素，使得实现上比较麻烦一些。另一个可能的方法是在逻辑区域 上引进一个非欧度量，使得多连通区域的边界也是凸的，但是这样会导致 Euler 方程不能够化为 线性形式。
，而将物理区域上的度量记为
。 由于在二维的情形，调和映射是比较特殊的，我们还要详细讨论。由于控制函数的概念由一个 函数扩展为一个矩阵，Cao，Huang和Russell在{Huang.6}中研究了网格加密的情况和控制函数的特征 向量、特征值之间的关系（二维情形，但是可以推广到高维的情况）。对于控制函数是一个函数乘以 一个单位矩阵的情形，事实上等价于是一个函数的情形，我们可以看到网格是在控制函数比较大的 区域被加密，在控制函数比较小的区域变得稀疏，在各个方向上的变化情况都是基本相同的。如果 控制函数是矩阵，并且特征值不相等，那么在不同的特征方向，网格的变化将是不同的。在特征值 大的方向，网格会比较的稠密，特征值小的方向上，网格比较稀疏，从而，从发展的观点来看，对于 特征值变大的方向，网格会加密，特征值变小的方向，网格会变稀疏。网格变化的速度和特征值 变化的速度是正相关的。
我们在选取控制函数的时候，一般的希望它仅仅依赖于解和物理区域上的坐标，也就是说，希望能够 具有形式
。在这样的情况下，可以保证网格的求解过程是收敛的。在一些研究 中，有作者甚至将取得和区域变换本身有关系，在这样的情况下，网格的收敛是没有理论的 保证的。
在以前的大部分研究中，控制函数被取为
![\[ M = \sqrt{1 + \alpha |{\vec u}|^2 + \beta |\nabla {\vec u}|^2 + \gamma |\nabla^2 {\vec u}|^2} \]](form_55.png)
这样的形式，其中
，
，
是正的参数。在一维情形下，如果
， 
，
，这个控制函数将会导致弧长坐标系，这可能是对于一维问题最令人钟情 的自适应网格了。一般的来说，在函数值比较大，函数的梯度比较大，甚至是二次导数比较大的位置 将会是函数本身比较奇异，或者函数的变化比较大的地方，所以，这个形式的控制函数可以满足一大批 问题的需要。
我们尝试了使用后验误差估计来构造控制函数的方法，对于变分不等式，最优控制问题都获得了不错 的计算结果。
通过解构造出来的控制函数在解本身比较奇异的情况下，一般的也比较奇异，数值结果表明，对 控制函数进行磨光可以极大的提高计算的精度。事实上，磨光是为了保证网格在尺度的过度上比较 连续。1980年代的一些文献中，开始意识到磨光的重要性，并且有很多关于磨光方法的研究。这些 磨光方法都是一维的方法，对于高维的问题，尤其是对于有限元网格，不好直接推广使用。而且现在 看来，这些方法都没有普遍的意义，已经很少有人再使用，而且网格移动对磨光本身并没有太严格 的要求，仅仅是要求磨光能够保持控制函数的主要轮廓，单调性这样一些主要的性质，我们就不具体讨论 这些磨光的方法了。尽管调和映射本身也具有一定的磨光的效果，但我们还是需要对控制函数进行 一些磨光的操作。我们要求的磨光方法希望能够保持控制函数的一些主要特点，最好是易于实现， 就认为是比较合适的了。控制函数在离散的计算出来以后，常常是一个分片常数的函数，我们使用的 磨光方法是将它从分片常数的空间投影到分片线性的空间中，然后再插值到分片常数的空间中，反复 做这个操作就能够将控制函数有效的磨光。这个投影算子
在分片常数函数
上实现方式为
![\[ \pi_h (f_h)|_{X^i} := {\displaystyle\frac {\sum_{T^j: X^i {\mathrm{\ is\ the\ vertex\ of\ }} T^j} f_j |T^j|} {\sum_{T^j: X^i {\mathrm{\ is\ the\ vertex\ of\ }} T^j} |T^j|}} \]](form_64.png)
其中，
表示单元
的体积。有结果表明，这个投影算子是依测度收敛的。插值算子 
在分片线性函数
上的实现方式为
![\[ \pi^{-1}_h (f_h)|_{T^i} := {\displaystyle\frac {\sum_{X^j: X^j {\mathrm{\ is\ the\ vertex\ of\ }} T^i} f_j} {\sum_{X^j: X^j {\mathrm{\ is\ the\ vertex\ of\ }} T^i} 1}} \]](form_69.png)
其中，分母事实上就是空间的维数
加上1，整个表达式是单元的各个节点上的值的平均值。
，我们考虑如下的问题：
![\[ \begin{array} {rcl} {\vec u}_t &=& {\vec u}_0 \\ {\vec u}|_{\partial \Omega} &=& t \times {\vec u}_0|_{\partial \Omega} \\ {\vec u}|_{t=0} &=& {\vec 0} \end{array} \]](form_72.png)
这是一个我们想要求解的问题的特例。它的初值是
，当将这个问题求解到
时，就是原问题 的初始状态。我们可以将这个问题的初始物理网格取为均匀网格，然后将它求解到， 并且在求解的过程中使用移动网格，就可以获得原问题的初值的物理网格。之所以要通过 这样一个同伦来实现这个过程，主要是因为对于比较奇异的初值，它的控制函数将也会比较奇异，我们 直接从均匀网格开始去获得很奇异的控制函数的网格，将会有比较大的网格移动步长，造成计算上 的一些困难，因为尽管理论上来说，解将会是个同胚，但是数值上，如果步长太大，还是可能造成网格的 缠绕，这正是我们要引进一个迭代过程来实现网格移动的原因。事实上我们可以通过如下的观点来看待这个同伦： 即对于从初值得到的控制函数， 我们通过从0变到1，逐步计算控制函数
的网格，最后获得控制函数的网格。 在解这个问题时的时间步长可以取地比较大，因为事实上我们不需要求解时间向前的方程（在时刻 ，直接代入
作为解即可）。一般的来说，取时间步长为
左右就能够对付 非常奇异的初值了。
解在物理网格和网格上的解函数
，我们 通过事先给定的方式计算出控制函数，在网格，求解椭圆型方程组
![\[ \label{2.5.1} \begin{array} {rcl} {\displaystyle\frac {\partial}{\partial x^i}} \left( G^{ij} {\displaystyle\frac {\partial \xi^k} {\partial x^j}} \right) &=& 0 \\ \xi|_{\partial \Omega} &=& \xi_b \end{array} \]](form_79.png)
在弱解空间
中的解函数
，事实上，我们首先得到的是物理网格的 节点在调和映射下的像
，我们得到它和逻辑网格之间的差别
，这是在逻辑网格上的一个分片线性的向量场
，利用 逻辑网格到物理网格上的变换的一阶导数，我们将它插值为物理网格上的一个分片线性的向量场，
![\[ \label{2.5.2} \delta x = \nabla_\xi x \delta \xi \]](form_85.png)
其中，
在单元
上的值通过恒等式
![\[ \label{2.5.3} \begin{array}{l} \left( X^{n_{T^i,2}}_t - X^{n_{T^i,1}}_t \cdots X^{n_{T^i,n+1}}_t - X^{n_{T^i,1}}_t \right) \nabla_\xi x |_{T_c^i} \\= \left( \Xi^{n_{T^i,2}}_t - \Xi^{n_{T^i,1}}_t \cdots \Xi^{n_{T^i,n+1}}_t - \Xi^{n_{T^i,1}}_t \right) \end{array} \]](form_88.png)
得到，其中
表示
的第j个顶点的序号。从而，我们得到物理区域上的每个节点的 移动方向
，我们将物理网格的节点更新为
![\[ \label{2.5.4} X^i_t := X^i_t + \eta \delta X^{i,*} \]](form_92.png)
其中，
是网格移动的步长，理论上来说，取为
是最佳的，但是为了避免网格的缠绕问题， 我们需要将它取得小一些。在我们的实际计算中，我们的网格移动步长的取法如下：对于单元 来说，存在一个正数![$\eta^i \in (0,1]$](form_95.png)
，使得网格移动的步长为
时，单元的定向不 发生改变，事实上，用
表示关于
的代数方程
![\[ \label{2.5.5} \det \left( X^{n_{T^i,1}}_t + \lambda \delta X^{n_{T^i,1},*} \ X^{n_{T^i,2}}_t + \lambda \delta X^{n_{T^i,2},*} \ \cdots \ X^{n_{T^i,n+1}}_t + \lambda \delta X^{n_{T^i,n+1},*} \right) = 0 \]](form_99.png)
正解的集合，那么
。然后我们取
作为网格移动的步长， 在前面乘以一个是为了避免数值原因引起的网格缠绕，我们在计算中取
。 我们是通过看
是否已经足够小来判断网格是否已经足够好。
如果我们使用有限差分法，可以通过等式
![\[\label{2.5.6} \begin{array}{rcl} - {\displaystyle\frac {\partial x^l} {\partial \mu}} &=& {\displaystyle\frac {\partial \xi^k} {\partial \mu}} {\displaystyle\frac {\partial x^l} {\partial \xi^k}} \\ &=& {\displaystyle\frac {\partial}{\partial x^i}} G^{ij} {\displaystyle\frac {\partial \xi^k} {\partial x^j}} {\displaystyle\frac {\partial x^l} {\partial \xi^k}} + G^{ij} {\displaystyle\frac {\partial^2 \xi^k} {\partial x^i \partial x^j}} {\displaystyle\frac {\partial x^l} {\partial \xi^k}} \\ &=& {\displaystyle\frac {\partial}{\partial x^i}} G^{ij} \delta^l_j - G^{ij} {\displaystyle\frac {\partial \xi^k} {\partial x^\beta}} {\displaystyle\frac {\partial^2 x^\beta} {\partial \xi^\gamma \partial \xi^\delta}} {\displaystyle\frac {\partial \xi^\gamma} {\partial x^i}} {\displaystyle\frac {\partial \xi^\delta} {\partial x^j}} {\displaystyle\frac {\partial x^l} {\partial \xi^k}} \\ &=& {\displaystyle\frac {\partial}{\partial x^i}} G^{il} - G^{ij} {\displaystyle\frac {\partial^2 x^\beta} {\partial \xi^\gamma \partial \xi^\delta}} {\displaystyle\frac {\partial \xi^\gamma} {\partial x^i}} {\displaystyle\frac {\partial \xi^\delta} {\partial x^j}} \delta^l_\beta \\ &=& {\displaystyle\frac {\partial} {\partial x^i}} G^{il} - G^{ij} {\displaystyle\frac {\partial^2 x^l} {\partial \xi^\gamma \partial \xi^\delta}} {\displaystyle\frac {\partial \xi^\gamma} {\partial x^i}} {\displaystyle\frac {\partial \xi^\delta} {\partial x^j}} \end{array} \]](form_104.png)
来计算网格移动的向量场，这其中需要用到恒等式
![\[ {\displaystyle\frac {\partial^2 \xi^k} {\partial x^i \partial x^j}} = - {\displaystyle\frac {\partial \xi^k} {\partial x^\beta}} {\displaystyle\frac {\partial^2 x^\beta} {\partial \xi^\gamma \partial \xi^\delta}} {\displaystyle\frac {\partial \xi^\gamma} {\partial x^i}} {\displaystyle\frac {\partial \xi^\delta} {\partial x^j}}. \]](form_105.png)
这个形式，对于差分格式是比较方便的。当然，这个表达式也可以用于有限元方法计算网格 移动的方向。令
，有
![\[\label{2.5.7} {\displaystyle\frac {\partial x^l} {\partial \mu}}= {\tilde G}^{\gamma \delta} {\displaystyle\frac {\partial^2 x^l} {\partial \xi^\gamma \partial \xi^\delta}} - {\displaystyle\frac {\partial}{\partial x^i}} G^{il}, \]](form_107.png)
写成弱形式
![\[\label{2.5.8} \begin{array}{rl} & \left. \displaystyle\int_{\Omega_c} {\displaystyle\frac {\partial x^l} {\partial \mu}} v d\xi \right. \\ =& \left. \displaystyle\int_{\Omega_c} \left\{{\tilde G}^{\gamma \delta} {\displaystyle\frac {\partial^2 x^l} {\partial \xi^\gamma \partial \xi^\delta}} - {\displaystyle\frac {\partial}{\partial x^i}} G^{il} \right\} v d\xi \right. \\ =& - \left. \displaystyle\int_{\Omega_c} \left\{{\tilde G}^{\gamma \delta} {\displaystyle\frac {\partial x^l} {\partial \xi^\gamma}} {\displaystyle\frac {\partial v} {\partial \xi^\delta}} - G^{il} {\displaystyle\frac {\partial \xi^k} {\partial x^i}} {\displaystyle\frac {\partial v} {\partial \xi^k}} \right\} d\xi, \right. \\ \end{array} \]](form_108.png)
其中检验函数
。离散求解就得到网格节点移动的方向。这个格式的 好处在于求解的系数矩阵是逻辑网格上的质量矩阵，如果使用预优迭代法求解的话，只需要构造一次 预优矩阵，可以节约很多的计算量。如果进行一次质量集中，更是可以完全避免求解线性方程组。
及其上的解函数
，通过上面的操作又得到了 一个新的网格
，我们试图获得解函数在新网格上的表示。目前，我们还仅仅考虑 一个一般的插值方法，并没有将解函数限制在某个具有特殊要求的空间中。但是这样的问题是需要 考虑的，比如对于不可压流体的速度，需要满足连续性方程，或者是守恒律方程，解需要满足一定 的守恒条件， 关于这个问题，我们正在进行研究。这个插值的问题，是一直困扰移动网格方法的 问题。一般的多项式插值已经被考虑过了，线性的插值方法总是带来一个系统误差，事实上是解函数 被逐渐的抹平了，但是高阶的插值会带来振动。我们现在使用的方法主要有以下特点：（1）保持 解函数的
范数；（2）不会带来振动。我们主要的想法是将网格的移动看作为一个人工的流动， 而在这个流动之下，保持解函数不变。这样，事实上是引进了从旧网格到新网格上的一个同伦，我们 跟随这个同伦，将解函数逐渐过渡过来。用基函数将表示出来为
，
是
的基函数。两个网格之间的不同可以表示为分片线性的向量场
，
是线性基函数。我们使用节点之间的线性变化
![\[\label{2.6.1} X^{i}(\tau)=X^i_t + \tau (X^{i,*}_t - X^i_t),\tau \in [0,1] \]](form_119.png)
获得这个同伦。伴随着节点位置的变化，基函数也将发生相应的变化，成为
，由于 基函数完全被网格所确定，因此可以将它记为
，其中
![\[\label{2.6.2} x(\tau)= x_0 + \tau (x_1 - x_0),\tau \in [0,1] \]](form_122.png)
其中
是物理区域上的恒等映射，
将物理区域中旧网格单元
上的点映射到新网格 的相应单元上，并且保持点的面积坐标，于是
恰好是网格之间的差异向量场。 从而，在同伦路径上，我们有一系列的的近似 
，为了使变换以后的的表示形式能够尽量保持不变， 我们要求
在离散空间里的弱意义下满足
![\[\label{2.6.3} {\displaystyle\frac {\partial u_h} {\partial \tau}} = 0 \]](form_129.png)
也就是
![\[\label{2.6.4} \begin{array}{rcl} 0 &=& \displaystyle\int_\Omega {\displaystyle\frac {\partial u_h} {\partial \tau}} \omega dx \\ &=& \displaystyle\int_\Omega \left\{ {\displaystyle\frac {\partial u_i} {\partial \tau}} \phi^i + u_i {\displaystyle\frac {\partial \phi^i} {\partial \tau}} \right\} \omega dx \end{array} \]](form_130.png)
其中
是检验函数，
![\[\label{2.6.5} \begin{array}{rcl} {\displaystyle\frac {\partial \phi^i} {\partial \tau}} &=& \nabla \phi^i \cdot {\frac {d x} {d\tau}} \\ &=& \nabla \phi^i \cdot (x_1(x) - x_0(x)) \\ &=& - \nabla \phi^i \cdot \delta x \end{array} \]](form_132.png)
于是我们最后需要解方程
![\[\label{2.6.6} \int_\Omega {\displaystyle\frac {\partial u_i} {\partial \tau}} \phi^i \omega dx = \int_\Omega u_i \nabla \phi^i \cdot \delta x \omega dx \]](form_133.png)
我们来分析一下这个方程的求解本质上的意义。在网格的发生移动时，我们如果保持 节点的函数值，那么我们就引进了一个变化量
![\[\label{2.6.7} \begin{array}{rcl} \dot{u}_h &:=& u_i {\displaystyle\frac {\partial \phi^i} {\partial \tau}}\\ &=& - u_i \nabla \phi^i \cdot \delta x \end{array} \]](form_134.png)
将这个变化量在中的部分丢弃，保留下和垂直的部分。(2.6.6})就是将 的变化量取为
中和垂直的部分。
可以看出，这两个条件都是必要条件。第一个条件是为了保持几何特性，而第二个条件是为了 强制的避免产生边界网格上的网格缠绕。同时，我们将会看到，在合理的将这两个条件具体化 了以后，这样的条件也能够唯一的确定物理区域到逻辑区域的变换。
和
分别表示 物理区域和逻辑区域相应的边，我们来考虑如下的从
到 的映射
![\[ K = \{\xi_b \in C^0|\xi_b:\partial \Omega \rightarrow \partial \Omega_c,\xi_b(\Gamma_i) =\Gamma_{c,i},\mathrm{on\ every\ }\Gamma_i\mathrm{,\ it\ is\ segmentwise\ linear}\} \]](form_139.png)
明显的
是一个开集，它可以表示为可列个闭集的并集
![\[ K = \bigcup_{M=1}^\infty K_M \]](form_141.png)
其中，
![\[ K_M = \{\xi_b \in K | {\displaystyle\frac 1 M} \leq \mathrm{ess\ inf} \xi^{'}_b,\mathrm{ess\ sup}\xi^{'}_b \leq M\}. \]](form_142.png)
根据Eells和Thompson's的结果，Dirichlet边界条件下的从物理区域到逻辑区域的调和映射是存在 唯一的，我们将边界条件为
的唯一的调和映射记为
，于是考虑 优化问题
![\[ \label{3.2.1} \begin{array}{rl} \min & E(P(\xi_b))\\ s.t. & \xi_b \in K_M \end{array} \]](form_145.png)
由于能量泛函
是凸的，对
，我们有
![\[ \begin{array}{rcl} \displaystyle\frac{1}{2} \left\{E(P(\xi^{(1)}_b)) + E(P(\xi^{(2)}_b)) \right\}& \ge & E\left(\displaystyle\frac{1}{2} \left\{P(\xi^{(1)}_b) + P(\xi^{(2)}_b) \right\}\right) \\ & \ge & E\left(P\left(\displaystyle\frac{1}{2} (\xi^{(1)}_b + \xi^{(2)}_b)\right)\right) \end{array} \]](form_148.png)
易见，
，从而，泛函在
在
中也是凸的，再由于是闭集，可知优化问题(3.2.1})存在唯一解。我们将算子
表示为 一个约束，问题(3.2.1})具有如下的形式
![\[ \label{3.2.2} \begin{array} {rl} \min & \displaystyle\sum_k \displaystyle\int_\Omega G^{ij} {\displaystyle\frac {\partial {\xi^k}} {\partial {x^i}}} {\displaystyle\frac {\partial {\xi^k}} {\partial {x^j}}} dx \\ [4mm] \mathrm{s.t.} & \xi_b \in K_M \end{array} \]](form_153.png)
假设单元的尺度为
，机器的精度为
，那么，我们可以将
取为
，其中
是一个给定的小正整数。这是我们在给定的机器 精度下能够获得的最好结果了。
为了能够有效的求解这个优化问题，我们需要详细分析约束的性质。第一个约束和第二个约束明显 是线性等式约束，但是第三个约束比较复杂一些，首先，它要求将物理区域的某条边映射到逻辑区域 的某条边所在的直线上，这是一个线性等式约束，其次，它要求这个直线到直线的映射的导数在 某个给定的范围内，这是一个线性不等式约束。这个优化问题的非线性来源于这个线性不等式约束。 需要说明的是，如果不选定一个，将问题的解集限制在中的话，就不能够保证解在 中。即使物理区域上的度量是单位度量，如果区域不是凸的，解就可能不是一个同胚。 但根据我们的经验，如果物理区域是凸集的话，控制函数又不是非常奇异的话，非线性约束 （事实上可以化为线性不等式约束）将不是积极约束，从而这个优化问题将转化为一个线性方程组的求解 问题。尽管这个优化问题是一个标准的二次规划问题，有很多针对这个问题的现成的优化算法， 但由于我们对求解这个优化问题的效率有相当的要求，根据我们的数值实验，这些现成的算法 并不能满足我们在 效率上的要求，我们将只给出不等式约束不是积极约束的算例。由于这个优化问题是一个边界控制 问题，正是我们目前所感兴趣的问题之一，我们现在正在研究该问题的理论和算法，希望能够发展 一个针对这类特殊问题的，比一般的二次规划算法更好的算法。在我们给出的算例中，为了使 不等式约束是不积极的约束，我们将物理区域取为凸区域。事实上，即使物理区域是凸区域，也不能 确保不等式约束不是积极约束，但是对控制函数的选择范围应该是相当的广泛的，在我们的数值 例子中，控制函数已经比较奇异，却没有出现不等式约束是积极约束的情况。
到
，将边界节点 排在后面，为
到
。首先，目标被离散为
![\[\label{3.3.1} \sum_k \int_\Omega G^{ij} {\displaystyle\frac {\partial {\lambda^\alpha}} {\partial {x^i}}} {\displaystyle\frac {\partial {\lambda^\beta}} {\partial {x^j}}} dx \xi^k_\alpha \xi^k_\beta \]](form_162.png)
对于边界上的约束，我们将它分为两个部分。一个是对一个边界上的节点
来说，它被映射到了 一条给定的直线上，这是一个线性约束，为
![\[\label{3.3.4} \Xi^i \cdot {\mathbf n}^i = b^i \]](form_164.png)
其中，
是逻辑区域的边界在某条边上的法向，
是一个事先确定的数，事实上 和共同唯一的确定了将要映到的这条边的几何状态，由于我们事先 就确定了将要映到哪条边上，所以和是已知的。另一个部分是关于边 界上的相邻的结点之间的距离的约束，对于边界上一对相邻的节点和
，我们要求它们 之间的距离既不能太大，也不能太小，从而这个约束是
![\[\label{3.3.5} {\displaystyle\frac 1 M} \le {\displaystyle\frac {| \Xi^i - \Xi^j |} {| X^i - X^j|}} \le M \]](form_168.png)
通过一些处理，这个约束可以化为线性的不等式约束。用矩阵形式将上述的离散表示出来，令
![\[ H = \left(\int_\Omega G^{ij} {\displaystyle\frac {\partial {\lambda^\alpha}} {\partial {x^i}}} {\displaystyle\frac {\partial {\lambda^\beta}} {\partial {x^j}}} dx \right)_{\begin{array}{c} 1 \leq \alpha \leq N\\ 1 \leq \beta \leq N \end{array}} \]](form_169.png)
用
，
，
和
分别表示物理区域和 逻辑区域，内部和边界上的节点的坐标排成的向量，则目标的矩阵表示为
![\[ \left(\begin{array}{cc} \Xi^{1,T} & \Xi^{2,T} \end{array}\right) \left(\begin{array}{cc} H & 0 \\ 0 & H \end{array} \right)\left(\begin{array}{c} \Xi^1 \\ \Xi^2 \end{array}\right), \]](form_174.png)
线性约束中的等式部分具有形式
![\[ \left(\begin{array}{cc} A_1 & A_2 \end{array}\right) \left(\begin{array}{c} \Xi^1_{bound} \\ \Xi^2_{bound} \end{array}\right) = b, \]](form_175.png)
不等式部分具有形式
![\[ \left(\begin{array}{cc} B_1 & B_2 \end{array}\right) \left(\begin{array}{c} \Xi^1_{bound} \\ \Xi^2_{bound} \end{array}\right) \ge c, \]](form_176.png)
这个不等式部分具有线性形式是二维问题的独有结果，对于三维的问题这个约束不可能化为 线性形式，这是我们只宣称在二维获得结果的原因。下面，我们假设这个不等式约束不是 积极约束，于是通过引进Lagrange乘子，我们可以将这个优化问题化为一个解线性方程组 的问题
![\[\label{3.3.6} \left(\begin{array}{ccccc} H_{11} & H_{12} & 0 & 0 & 0 \\ H_{21} & H_{22} & 0 & 0 & A'_1 \\ 0 & 0 & H_{11} & H_{12} & 0 \\ 0 & 0 & H_{21} & H_{22} & A'_2 \\ 0 & A_1 & 0 & A_2 & 0 \end{array}\right) \left(\begin{array}{c} \Xi^1_{inner} \\ \Xi^1_{bound} \\ \Xi^2_{inner} \\ \Xi^2_{bound} \\ \lambda_3 \end{array}\right) = \left(\begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ b \end{array}\right) \]](form_177.png)
我们看到，由于最主要的非零元素都集中在矩阵
中，从而，我们 可以采用迭代法来比较有效的求解这个方程组。如果仅仅是用抛弃不等式约束来求解的情形， 这个方法是可以被推广到三维的。
![\[\label{3.4.1} \begin{array}{rl} \min & \displaystyle\sum_k \displaystyle\int_\Omega \sum_i \left( {\displaystyle\frac {\partial \xi^k} {\partial x^i}} \right)^2 dx \\ s.t. & \xi_b \in K_M \end{array} \]](form_179.png)
来得到。在我们的算例中，我们由于只能将物理区域取为凸区域，并将逻辑区域取为和物理区域 相同，因此没有将这个方案实现，但是我们给出了一个例子来看这样的构造方式和边界上是给定 的映射的方式给出的网格有什么不一样。这个优化问题的离散和(3.2.2})相似。
边界上网格点的移动：和仅仅移动内部的节点不同，我们现在要移动边界上的节点。尽管 在逻辑区域上，网格的移动方向是沿着边界的（事实上，由于数值误差的原因，也不是严格的沿着 边界的方向的），但是将这个移动方向插值到了物理区域上时，物理 区域上的网格移动方向可能并不是沿着边界的。因此为了保证边界上的节点被移动以后，还在边界 上，我们需要将物理区域上的网格移动方向在边界上的这些值向边界上作一个投影。
边界上不等式约束的处理：不等式约束本来具有(3.3.5})的形式，但是由于我们事先确定 了边的映射方式，于是我们根据边界映射的方式，首先可以确定相邻的两个点
, 
的像在边
上的顺序 是
比较靠近顶点
，
比较靠近顶点
，于是我们可以确定有
，再通过考虑边的倾斜度，比如说可 以保证
，那 么我们就选取作为方向，将这个边界段上的约束写为
![\[ {\displaystyle\frac 1 M} \le {\displaystyle\frac {\sqrt {1 + \alpha^{-2}}} {|X_1 - X_2|}} (\xi_2 - \xi_1) \le M. \]](form_189.png)
![\[\label{3.6.1} {\displaystyle\frac {\partial} {\partial x}}\left({\displaystyle\frac 1 M} {\displaystyle\frac {\partial \xi} {\partial x}}\right) = 0 \]](form_190.png)
边值为
, 
，就是直接寻找等分布点
![\[ x_0 = X_0 < X_1 < \cdots < X_{N-1} < X_N = x_1 \]](form_193.png)
使得
。因此，对于将边界作为 一维的方法，可以直接得到边界点的分布状况，然后就可以使用我们前面的方法调整内部节点的分布。 从而，这个方法可以改进的方面就是边界上的控制函数如何通过内部的控制函数投影得到。我们看到 至少有两个方案可以选择：一个是将控制函数的解析表达式投影到边界上，得到边界上的控制函数的 解析表达式，然后在计算时使用这个表达式计算控制函数的值；另一个方案是将控制函数在边界单元 上计算出来的值投影到单元位于边界的边上，作为边界上的控制函数的值。我们比较倾向于第二种 方案，因为它的计算用到了解在内部节点上的信息，蕴涵着一定的边界和内部耦合的因素，而第一 种方案仅仅根据解在边界节点上的信息计算控制函数，边界和内部是完全独立的。我们举一个例子 来看一看这两个方案的不同之处。假设控制函数具有形式
![\[ \label{3.6.2} M = {\sqrt {1 + |\nabla u_h |^2}} \]](form_195.png)
其中，离散解是分片线性函数，在一个三角形单元上，在节点
上的取值是
，那么，在第一种方案下，在边界上的控制函数的值是
![\[ {\sqrt {1 + \left( {\displaystyle\frac {u_1 - u_2} {| X_1 - X_2|}} \right)^2}}, \]](form_198.png)
而在第二种方案下，在边界上的控制函数的值是
![\[ {\sqrt {1 + \left({\displaystyle\frac {\partial u_h} {\partial x}}\right)^2 + \left({\displaystyle\frac {\partial {u_h}} {\partial y}}\right)^2}}, \]](form_199.png)
其中
![\[ \left({\displaystyle\frac {\partial {u_h}} {\partial x}}\right)^2 + \left({\displaystyle\frac {\partial {u_h}} {\partial y}}\right)^2 \]](form_200.png)
等于以
为顶点的三维空间中的三角形的面积和单元面积的比值。 在这样的情况下，我们使用的算法的步骤是：
我们提到，这样的方法在{Huang.2}中实现过，但是由于文献中细节不是很清楚，我们不 知道其实现是否用的是上面的方式。
1.4.7