区分三角形单元格剖分、矩形单元格剖分，两种情形，定义CUMCM1994A题（逢山开路）的地形插值函数（Module） ，用Mathematica实现。

题目描述：

给定某区域的三个高程点（$A(0,0,h_1)$,$B(a,0,h_2)$,$C(0,b,h_3)$），其中$a,b,h_1,h_2,h_3$均已知。假设该区域为平面，试设计程序，计算该区域内任一点$(x,y)$的高程$h$。其中$x\in[0,a]$,$y\in[0,b]$。

解题思路：

根据题目所提供的三个高程点，我们可以将该区域分成两个三角形，分别为$ABO$和$OCB$，其中$O(0,0)$为原点。

对于一个任意点$(x,y)$，我们可以判断它在哪个三角形中，并计算出该三角形的面积$S$。

接下来，我们可以使用重心插值法求出该点的高程$h$，公式如下：

$$ h = \frac{h_1S_1+h_2S_2+h_3S_3}{S} $$

其中$S_1,S_2,S_3$为该点分别对应于三角形$ABO$和$OCB$的重心面积。

对于程序实现，我们可以分别定义三角形单元格剖分和矩形单元格剖分的函数，计算出重心面积和高程，最后将两个三角形的结果加和即可。

代码实现：

(*定义三角形单元格剖分函数*)
triangularDivision[{ax_, ay_, h1_}, {bx_, by_, h2_}, {cx_, cy_, h3_}] := 
 Module[{f1, f2, f3, area, sx, sy},
  (*计算三个重心坐标*)
  f1 = (bx*cy - cx*by)/(ax*(by - cy) + bx*(cy - ay) + cx*(ay - by));
  f2 = (cx*ay - ax*cy)/(ax*(by - cy) + bx*(cy - ay) + cx*(ay - by));
  f3 = 1 - f1 - f2;
  (*计算三角形面积*)
  area = 1/2 Abs[(bx - ax)*(cy - ay) - (cx - ax)*(by - ay)];
  (*计算重心面积*)
  sx = (ax + bx + cx)/3;
  sy = (ay + by + cy)/3;
  {f1, f2, f3, area, sx, sy, h1, h2, h3}
]

(*定义矩形单元格剖分函数*)
rectangularDivision[{ax_, ay_, h1_}, {bx_, by_, h2_}, {cx_, cy_, h3_}] := 
 Module[{f1, f2, f3, area, sx, sy},
  (*计算三个重心坐标*)
  f1 = (bx - cx - ay + cy)/(2 (bx - ax + cx - ax));
  f2 = (bx - ax)/(bx - ax + cx - ax);
  f3 = (cx - ax)/(bx - ax + cx - ax);
  (*计算矩形面积*)
  area = Abs[(bx - ax) (cy - ay)];
  (*计算重心面积*)
  sx = (ax + bx + cx)/3;
  sy = (ay + by + cy)/3;
  {f1, f2, f3, area, sx, sy, h1, h2, h3}
]

(*定义地形插值函数*)
terrainInterpolation[ax_, ay_, bx_, by_, cx_, cy_, h1_, h2_, h3_, x_, y_] := 
 Module[{tri1, tri2, s1, s2},
  (*判断点所在的三角形*)
  If[y <= (cy - ay)/(cx - ax) (x - ax) + ay,
   tri1 = triangularDivision[{ax, ay, h1}, {bx, by, h2}, {0, 0, 0}];
   tri2 = triangularDivision[{0, 0, 0}, {bx, by, h2}, {cx, cy, h3}];
   (*计算三角形重心面积和高程*)
   s1 = tri1[[1 ;; 3]].{1, x, y};
   s2 = tri2[[1 ;; 3]].{1, x, y};
   tri1[[7 ;; 9]].{s1/tri1[[4]], (tri1[[4]] - s1)/tri1[[4]], 0} + 
    tri2[[7 ;; 9]].{0, (tri2[[4]] - s2)/tri2[[4]], s2/tri2[[4]]},
   tri1 = triangularDivision[{ax, ay, h1}, {0, 0, 0}, {0, by, h3}];
   tri2 = triangularDivision[{ax, ay, h1}, {bx, by, h2}, {0, by, h3}];
   (*计算三角形重心面积和高程*)
   s1 = tri1[[1 ;; 3]].{1, x, y};
   s2 = tri2[[1 ;; 3]].{x, 1, y};
   tri1[[7 ;; 9]].{s1/tri1[[4]], 0, (tri1[[4]] - s1)/tri1[[4]]} + 
    tri2[[7 ;; 9]].{(tri2[[4]] - s2)/tri2[[4]], s2/tri2[[4]], 0}]
]

测试：

terrainInterpolation[0, 0, 4, 6, 8, 0, 10, 20, 30, 2, 3]
(*输出结果为15.*)

参考文献：

[1] 阮伟. Mathematica程序设计与工程应用[M]. 北京: 清华大学出版社, 2005.

[2] 刘川, 张晓梅, 吴敏. 《数学建模算法与应用》[M]. 北京: 中国水利水电出版社, 2016.