I have the following three functions $\lambda_i$ and $q_i$:
λ1 = (a x + b y + c)/.First@Solve[{a x1 + b y1 + c == 1, a x2 + b y2 + c == 0,
a x3 + b y3 + c == 0}, {a, b, c}];
λ2 = (a x + b y + c)/.First@Solve[{a x1 + b y1 + c == 0, a x2 + b y2 + c == 1,
a x3 + b y3 + c == 0}, {a, b, c}];
λ3 = (a x + b y + c)/.First@Solve[{a x1 + b y1 + c == 0, a x2 + b y2 + c == 0,
a x3 + b y3 + c == 1}, {a, b, c}];
q1 = 1;
q2 = x - (x1+x2+x3)/3;
q3 = y - (y1+y2+y3)/3;
Then we have 7 functions $\phi_i$ defined as follows
φ1 = λ1 (2 λ1 - 1) + 3 λ1 λ2 λ3;
φ2 = λ2 (2 λ2 - 1) + 3 λ1 λ2 λ3;
φ3 = λ3 (2 λ3 - 1) + 3 λ1 λ2 λ3;
φ4 = 4 λ2 λ3 - 12 λ1 λ2 λ3;
φ5 = 4 λ3 λ1 - 12 λ1 λ2 λ3;
φ6 = 4 λ1 λ2 - 12 λ1 λ2 λ3;
φ7 = 27 λ1 λ2 λ3;
I am interested in integrals in this form $$ \int \limits_{\text{triangle}} \phi_i \: \mathrm{d} S\\ \int \limits_{\text{triangle}} q_i \phi_j \: \mathrm{d} S\\ \int \limits_{\text{triangle}} q_i \partial_x \phi_j \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} q_i \partial_y \phi_j \: \mathrm{d} S\\ \int \limits_{\text{triangle}} \phi_i \phi_j \: \mathrm{d} S\\ \int \limits_{\text{triangle}} \phi_i \partial_x \phi_j \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} \phi_i \partial_y \phi_j \: \mathrm{d} S\\ \int \limits_{\text{triangle}} \phi_i \phi_j \partial_x \phi_k \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} \phi_i \phi_j \partial_y \phi_k \: \mathrm{d} S\\ \int \limits_{\text{triangle}} \phi_i \partial_x \phi_j \partial_x \phi_k \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} \phi_i \partial_x \phi_j \partial_y \phi_k \: \mathrm{d} S, \quad \int \limits_{\text{triangle}} \phi_i \partial_y \phi_j \partial_y \phi_k \: \mathrm{d} S\\ $$ where "triangle" is defined by its vertices $(x_1, y_1)$, $(x_2, y_2)$, $(x_3, y_3)$
So I tried it with $\phi_i, \phi_j$. It is always proportional to the are of the triangle, so I am interested only in the numeric constant in front of it (integral = something*expression for the area)
Δ =
1/2 Abs[x2 y1 - x3 y1 - x1 y2 + x3 y2 + x1 y3 - x2 y3];
resij = {};
For[i = 1, i <= 7, i++,
For[j = 1, j <= 7, j++,
expr = ToExpression[ToString[φ] <> ToString[i]]*
ToExpression[ToString[φ] <> ToString[j]] // Simplify;
res = 1/Δ Integrate[
expr, {x, y} ∈
Triangle[{{x1, y1}, {x2, y2}, {x3, y3}}]] // Simplify //
Expand;
Print["(" <> ToString[i] <> "," <> ToString[j] <> ") " <>
ToString[N@res]]
AppendTo[resij, res];
]]
Although it takes a very long time to compute, it gives the desired result:
(1,1) 0.0321429
(1,2) -0.00674603
(1,3) -0.00674603
(1,4) -0.00793651
...
2520 Partition[resij, 7] // MatrixForm
$$ \begin{pmatrix} 81 & -17 &-17 & -20 & 36 & 36 & 27 \\ -17 & 81 & -17 & 36 &-20 & 36 & 27 \\ -17 & -17 & 81 & 36 & 36 & -20 & 27 \\ -20 & 36 & 36 &208 &-16 &-16 & 108 \\ 36 & -20 & 36 & -16 & 208 &-16 & 108 \\ 36 & 36 & -20 & -16 & -16 & 208 & 108 \\ 27 & 27 & 37 & 108 & 108 & 108 & 729 \end{pmatrix} $$
Then I tried it with $\phi_i \phi_j \partial_x \phi_k$, which is always in form $$ \frac{\Delta}{|\Delta|} \left( a y_1 + b y_2 + c y_3 \right) $$ where $\Delta$ is signed triangle area. So I tried:
Δ =
1/2 Cross[{x2 - x1, y2 - y1, 0}, {x3 - x1, y3 - y1, 0}][[3]];
resijdxk = {};
For[i = 1, i <= 7, i++,
For[j = 1, j <= 7, j++,
For[k = 1, k <= 7, k++,
expr = (Abs[Δ]/Δ)^-1 ToExpression[
ToString[φ] <> ToString[i]]*
ToExpression[ToString[φ] <> ToString[j]]*
D[ToExpression[ToString[φ] <> ToString[k]], x] //
Simplify;
res = Integrate[
expr, {x, y} ∈
Triangle[{{x1, y1}, {x2, y2}, {x3, y3}}]] // Simplify //
Expand;
AppendTo[
resijdxk, {Coefficient[res, y1], Coefficient[res, y2],
Coefficient[res, y3]}];
Print["(" <> ToString[i] <> "," <> ToString[j] <> "," <>
ToString[k] <> ") " <> "(" <>
ToString[N@Coefficient[res, y1]] <> "," <>
ToString[N@Coefficient[res, y2]] <> "," <>
ToString[N@Coefficient[res, y3]] <> ")"]
]]]
(1,1,1) (0.,0.0309524,-0.0309524)
(1,1,2) (0.00892857,-0.00297619,-0.00595238)
...
which gets the result, however, one in 40 seconds. It would take 7*7*7*40 seconds = 3 hours and 48 minutes to compute.
What is a faster way to compute this? I have a suspicion that this is just a bunch of spagetti code I wrote and it can be somehow optimalized, but I don't know how.