计算 GMSH 地理脚本内的二次拟合系数答案

【问题标题】：Calculating coefficients of a quadratic fit inside GMSH geo script计算 GMSH 地理脚本内的二次拟合系数
【发布时间】：2018-04-29 17:22:07
【问题描述】：

我遇到了一个问题，需要计算要在 Gmsh .geo 文件中使用的矩阵的系数。我使用 Matlab 创建了一个小脚本来求解 2x2 矩阵。代码如下；

clear all;clc;close all 
R = 0.509515; % m
R1 = 0.559467; % m
R2 = 0.579594; % m
L = 0;
L1 = -0.08568;
L2 = 0.13974;
y = [R1-R;R2-R];
x = [L1^2 L1;L2^2 L2];
c = x\y;
x1 = linspace(-0.08568,0.13974,20);

从上面的小代码中，我需要使用矩阵c中的系数，我将其作为c1和c2手动转移到以下geo文件中

R = 0.509515; //% m
R11 = 0.10; //% m
R1 = 0.559467; //% m
R2 = 0.579594; //% m
L = 0;
L1 = -0.08568; 
L2 = 0.13974;
c1 = 4.811;
c2 = -0.1708; 

nPoints = 20; // Number of discretization points (top-right point of the inlet region)
For i In {1 : nPoints}
  x = L1 + (L2-L1)*i/(nPoints + 1);
  pList[i] = newp;
  Point(pList[i]) = {x,
                ( c1*x^2+c2*x+R ),
                0};
EndFor
Spline(1) = {20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1};

请注意，在地理文件中，我无法使用 linspace，因为我不知道 GMSH 中是否存在 linspace 语法。我的目标是让geo文件求解矩阵以拟合二次多项式。

【问题讨论】：

标签： matlab linear-algebra mesh

【解决方案1】：

最简单的选择就是在 GMSH 内求解一个 2 x 2 线性方程组。

通过添加以下计算而不是直接插入 c1 和 c2 的值，会将工作从 Matlab 转移到 GMSH。注意，因为它是一个 2 x 2 系统，所以我没有考虑数值稳定性；但是，对于这样一个实际应用的小系统，无论如何它应该不重要。

// Calculating solution of Ax=y
// A = [L1^2 L1; L2^2 L2]; y = [R1-R; R2-R]; x = [c1;c2]
// Ainv = 1/det(A)*[L2 -L1; -L2^2 L1^2]
// x = Ainv*y
det = L1*L1*L2-L1*L2*L2;
If (Abs(det) < 1E-10) //Yeah, I am lazy.
    Printf("Ill-conditioned fit required. Will crash.");
EndIf
factor = 1./det;
c1 = factor*(L2*(R1-R) - L1*(R2-R));
c2 = factor*(-L2*L2*(R1-R) +L1*L1*(R2-R));

【讨论】：