0
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

This his how we you can compute the requestest integrals on the standard triangle {{0,0},{1,0},{0,1}} as follows:

We start by defining the basis functions

λ1 = x;
λ2 = y;
λ3 = 1 - x - y;
q1 = 1;
q2 = x - 1/3;
q3 = y - 1/3;
φ = {λ1 (2 λ1 - 1) + 3 λ1 λ2 λ3, λ2 (2 λ2 - 1) + 3 λ1 λ2 λ3, λ3 (2 λ3 - 1) + 3 λ1 λ2 λ3, 4 λ2 λ3 - 12 λ1 λ2 λ3, 4 λ3 λ1 - 12 λ1 λ2 λ3, 4 λ1 λ2 - 12 λ1 λ2 λ3, 27 λ1 λ2 λ3};
Dφ = D[φ, {{x, y}, 1}];
Φ = {Transpose[{φ}], Dφ, {{q1}, {q2}, {q3}}};

Memoized way to compute integrals of monomials:

ClearAll[int];
mem : int[{a_, b_}] := mem=Integrate[Integrate[x^a y^b, {y, 0, 1 - x}], {x, 0, 1}]

The working horse function; it takes a integer list (expecting to contain1, 2, and 3, only), reads the respective tensors from Φ and computes the integrals of their tensor products.

getIntegrals[idx_?VectorQ] := Module[{a},
  a = 
    Outer[TensorProduct, 
     Sequence @@ Table[Φ[[i]], {i, idx}], 1];
  Developer`ToPackedArray[
   Map[
    Total[
      KeyValueMap[
       {key, value} \[Function] int[key] value,
       Association@CoefficientRules[#, {x, y}]
       ]
      ] &,
    a,
    {ArrayDepth[a]}
    ]
   ]
  ]

Computing all requested integrals on the standard triangle:

results = {
    getIntegrals[{1}],
    getIntegrals[{3, 1}],
    getIntegrals[{3, 2}],
    getIntegrals[{1, 1}],
    getIntegrals[{1, 2}],
    getIntegrals[{1, 1, 2}],
    getIntegrals[{1, 2, 2}]
    }; // AbsoluteTiming // First

1.11174

Here result[[i]] contains the integrals in the i-th row of the OP. For example, the mass matrix (fourth row) can be found as follows:

Flatten[results[[4]], {{1, 3}, {2, 4}}] // TeXForm

$$ \left( \begin{array}{ccccccc} \frac{9}{560} & -\frac{17}{5040} & -\frac{17}{5040} & -\frac{1}{252} & \frac{1}{140} & \frac{1}{140} & \frac{3}{560} \\ -\frac{17}{5040} & \frac{9}{560} & -\frac{17}{5040} & \frac{1}{140} & -\frac{1}{252} & \frac{1}{140} & \frac{3}{560} \\ -\frac{17}{5040} & -\frac{17}{5040} & \frac{9}{560} & \frac{1}{140} & \frac{1}{140} & -\frac{1}{252} & \frac{3}{560} \\ -\frac{1}{252} & \frac{1}{140} & \frac{1}{140} & \frac{13}{315} & -\frac{1}{315} & -\frac{1}{315} & \frac{3}{140} \\ \frac{1}{140} & -\frac{1}{252} & \frac{1}{140} & -\frac{1}{315} & \frac{13}{315} & -\frac{1}{315} & \frac{3}{140} \\ \frac{1}{140} & \frac{1}{140} & -\frac{1}{252} & -\frac{1}{315} & -\frac{1}{315} & \frac{13}{315} & \frac{3}{140} \\ \frac{3}{560} & \frac{3}{560} & \frac{3}{560} & \frac{3}{140} & \frac{3}{140} & \frac{3}{140} & \frac{81}{560} \\ \end{array} \right) $$

$\endgroup$
  • $\begingroup$ The first solution looks fine (giving the desired result, up to a constant). However, I don't get the second solution. I'm looking for coefficients that appear in front of the x1, x2, x3 or y1, y2, y3 when I divide that integral by the sign of signed integral area. Your results just give one number per integral. Other than that, I think that Outer[Times, phi, phi, phi] should have the last one D[phi, x] or D[phi, y], but that's just detail. I'm more interested in those constants, rather than the pure value of the integral over the basic simplex. $\endgroup$ – user16320 Jan 9 '18 at 17:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.