Skip to content

Instantly share code, notes, and snippets.

@FabioLuporini

FabioLuporini/gist:14e79457d6b15823c1cd

Last active February 17, 2016 10:22
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save FabioLuporini/14e79457d6b15823c1cd to your computer and use it in GitHub Desktop.
Save FabioLuporini/14e79457d6b15823c1cd to your computer and use it in GitHub Desktop.
Download ZIP
Application of sharing elimination to a hyperelastic model
Raw
gistfile1.txt
The UFL specification of the model is available at:
https://github.com/firedrakeproject/firedrake-bench/blob/experiments/forms/firedrake_forms.py
--- Input expression in normal form ---
double (*w1)[1] = (double (*)[1])c1;
double (*w2)[1] = (double (*)[1])c2;
// Compute Jacobian
double J[4];
compute_jacobian_triangle_2d(J, coordinate_dofs);
// Compute Jacobian inverse and determinant
double K[4];
double detJ;
compute_jacobian_inverse_triangle_2d(K, detJ, J);
const double det = fabs(detJ);
static const double W1 = 0.5;
static const double FE0_D10[1][6] = {{-1.0, 1.0, 0.0, 0.0, 0.0, 0.0}};
static const double FE0_C1_D01[1][6] = {{0.0, 0.0, 0.0, -1.0, 0.0, 1.0}};
static const double FE0_C1_D10[1][6] = {{0.0, 0.0, 0.0, -1.0, 1.0, 0.0}};
static const double FE0_D01[1][6] = {{-1.0, 0.0, 1.0, 0.0, 0.0, 0.0}};
static const double FE1[1][1] = {{1.0}};
double F0 = 0.0;
double F1 = 0.0;
double F2 = 0.0;
double F3 = 0.0;
double F4 = 0.0;
double F5 = 0.0;
for (int r = 0; r < 6; r += 1)
{
F3 += (w0[r][0] * FE0_C1_D01[0][r]);
F2 += (w0[r][0] * FE0_C1_D10[0][r]);
F1 += (w0[r][0] * FE0_D01[0][r]);
F0 += (w0[r][0] * FE0_D10[0][r]);
}
for (int r = 0; r < 1; r += 1)
{
F4 += (w1[0][0] * FE1[0][0]);
F5 += (w2[0][0] * FE1[0][0]);
}
for (int j = 0; j < 6; j += 1)
{
for (int k = 0; k < 6; k += 1)
{
A[j][k] += (((((K[1] * FE0_D10[0][j]) + (K[3] * FE0_D01[0][j])) * ((((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) *
(((((((K[2] * F3) + (K[0] * F2)) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0]
* F0) + 1.0))) / 2.0)) + ((((((K[2] * F3) + (K[0] * F2)) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0))) * F4) + ((((((((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) *
((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[0]
* FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_D10[0][k]) + (K[3] *
FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) + ((((((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k]))
* ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) +
(((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_D10[0][k]) +
(K[3] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0))) * ((K[2] * F1) + (K[0] * F0) + 1.0) * F4) + (((K[3]
* F1) + (K[1] * F0)) * (((((((((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] *
FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01
[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[3] * F3) +
(K[1] * F2) + 1.0))) / 2.0)) + ((((((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) +
(((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] *
FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[3]
* F3) + (K[1] * F2) + 1.0))) / 2.0))) * F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k
])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) +
(((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_C1_D10[0][k])
+ (K[3] * FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0))) / 2.0) + (((((K[0] * FE0_C1_D10[0][k]) + (K[2] *
FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3)
+ (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (((K[0] *
FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) * 2.0))) + (((K[1] * FE0_D10[0][k])
+ (K[3] * FE0_D01[0][k])) * (((((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2)
+ 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0)) + ((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F
0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0))) * F4) + (((F5) / 2.0) *
((1.0)) * ((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F
3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0) + (((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) +
(K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0)) * 2.0))))) + (((K[0] * FE0_C1_D10[0][j]) + (K[2]
* FE0_C1_D01[0][j])) * ((((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * (((((((K[2] * F3) + (K[0] * F2)) *
((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) + ((((((K[
2] * F3) + (K[0] * F2)) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0] * F0) +
1.0))) / 2.0))) * F4) + ((((((((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) +
(((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] *
FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[2] * F1)
+ (K[0] * F0) + 1.0))) / 2.0)) + ((((((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2)))
+ (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] *
FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[2] * F1)
+ (K[0] * F0) + 1.0))) / 2.0))) * ((K[3] * F3) + (K[1] * F2) + 1.0) * F4) + (((K[2] * F3) + (K[0] * F2)) * (((((((((K[0] *
FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] *
FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) +
(K[0] * F0) + 1.0)) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) +
((((((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_C1_D10[0][k]) +
(K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] *
F1) + (K[0] * F0) + 1.0)) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.
0))) * F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0)))
+ (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] *
FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[3]
* F3) + (K[1] * F2) + 1.0))) / 2.0) + (((((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] *
F2))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k])
+ (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) *
((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) * 2.0))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((((((((
(K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0)
+ 1.0)) + (-1.0)) / 2.0)) + ((((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) +
1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0))) * F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[3] * F1) + (K[1] *
F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) /
2.0) + (((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) +
(K[0] * F0) + 1.0)) + (-1.0)) / 2.0)) * 2.0))))) + (((K[1] * FE0_C1_D10[0][j]) + (K[3] * FE0_C1_D01[0][j])) * ((((K[2] * F
3) + (K[0] * F2)) * (((((((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] *
FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01
[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] *
F0) + 1.0))) / 2.0)) + ((((((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0]
* FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D
01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0]
* F0) + 1.0))) / 2.0))) * F4) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * (((((((K[2] * F3) + (K[0] * F2))
* ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) +
((((((K[2] * F3) + (K[0] * F2)) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0]
* F0) + 1.0))) / 2.0))) * F4) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * (((((((((K[3] * F1) + (K[1] * F
0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) /
2.0)) + ((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3)
+ (K[1] * F2) + 1.0)) + (-1.0)) / 2.0))) * F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) +
(K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0) + (((((K[2] * F3)
+ (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) +
(-1.0)) / 2.0)) * 2.0))) + ((((((((((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) +
(((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] *
FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[3]
* F3) + (K[1] * F2) + 1.0))) / 2.0)) + ((((((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0)))
+ (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] *
FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[3]
* F3) + (K[1] * F2) + 1.0))) / 2.0))) * F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k
])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) +
(((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_C1_D10[0][k])
+ (K[3] * FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0))) / 2.0) + (((((K[0] * FE0_C1_D10[0][k]) + (K[2] *
FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3)
+ (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (((K[0] *
FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) * 2.0)) * ((K[3] * F3) + (K[1] * F2)
+ 1.0)))) + (((K[0] * FE0_D10[0][j]) + (K[2] * FE0_D01[0][j])) * ((((K[3] * F1) + (K[1] * F0)) * (((((((K[1] * FE0_C1_D10[
0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) *
((K[3] * F1) + (K[1] * F0))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.
0)) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) + ((((((K[1] *
FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] *
FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[3] * F3) +
(K[1] * F2) + 1.0)) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0))) *
F4) + (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * (((((((K[2] * F3) + (K[0] * F2)) * ((K[3] * F3) + (K[1] * F2) +
1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) + ((((((K[2] * F3) + (K[0] * F2)) *
((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0))) * F4) +
(((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * (((((((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) +
(((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0)) + ((((((K[2] * F3) + (K[0] * F2))
* ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0)))
* F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] *
F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0) + (((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] *
F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0)) * 2.0))) + ((((((((((K[0]
* FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] *
FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) +
(K[0] * F0) + 1.0)) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) +
((((((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_C1_D10[0][k]) +
(K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] *
F1) + (K[0] * F0) + 1.0)) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.
0))) * F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0)))
+ (((K[1] * FE0_D10[0][k]) + (K[3] * FE0_D01[0][k])) * ((K[3] * F1) + (K[1] * F0))) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] *
FE0_C1_D01[0][k])) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[1] * FE0_C1_D10[0][k]) + (K[3] * FE0_C1_D01[0][k])) * ((K[3]
* F3) + (K[1] * F2) + 1.0))) / 2.0) + (((((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] *
F2))) + (((K[0] * FE0_C1_D10[0][k]) + (K[2] * FE0_C1_D01[0][k])) * ((K[2] * F3) + (K[0] * F2))) + (((K[0] * FE0_D10[0][k])
+ (K[2] * FE0_D01[0][k])) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (((K[0] * FE0_D10[0][k]) + (K[2] * FE0_D01[0][k])) *
((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) * 2.0)) * ((K[2] * F1) + (K[0] * F0) + 1.0))))) * det * W1);
}
}
--- Expression after sharing elimination as performed by COFFEE ---
double (*w1)[1] = (double (*)[1])c1;
double (*w2)[1] = (double (*)[1])c2;
// Compute Jacobian
double J[4];
compute_jacobian_triangle_2d(J, coordinate_dofs);
// Compute Jacobian inverse and determinant
double K[4];
double detJ;
compute_jacobian_inverse_triangle_2d(K, detJ, J);
const double det = fabs(detJ);
static const double W1 = 0.5;
static const double FE0_D10[1][4] __attribute__((aligned(32))) = {{-1.0, 1.0, 0.0}};
static const double FE0_C1_D01[1][4] __attribute__((aligned(32))) = {{-1.0, 0.0, 1.0}};
static const double FE1[1][4] __attribute__((aligned(32))) = {{1.0}};
double F0 = 0.0;
double F1 = 0.0;
double F2 = 0.0;
double F3 = 0.0;
double F4 = 0.0;
double F5 = 0.0;
for (int r = 0; r < 2; r += 1)
{
F2 += (w0[r+3][0] * FE0_D10[0][r]);
F0 += (w0[r][0] * FE0_D10[0][r]);
}
for (int r = 0; r < 3; r += 1)
{
F3 += (w0[r+3][0] * FE0_C1_D01[0][r]);
F1 += (w0[r][0] * FE0_C1_D01[0][r]);
}
for (int r = 0; r < 1; r += 1)
{
F4 += (w1[0][0] * FE1[0][0]);
F5 += (w2[0][0] * FE1[0][0]);
}
double CONST_0_1_0 = ((K[3] * F3) + (K[1] * F2) + 1.0);
double CONST_0_1_1 = ((K[3] * F1) + (K[1] * F0));
double CONST_0_1_2 = ((K[2] * F3) + (K[0] * F2));
double CONST_0_1_3 = (((K[2] * F1) + (K[0] * F0) + 1.0) * F4);
double CONST_0_1_4 = ((((((((K[2] * F3) + (K[0] * F2)) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0)) + ((((((K[2] * F3) + (K[0] * F2)) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (((K[3] * F1) + (K[1] * F0)) * ((K[2] * F1) + (K[0] * F0) + 1.0))) / 2.0))) * F4);
double CONST_0_1_5 = (((((((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0)) + ((((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0))) * F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0) + (((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0)) * 2.0));
double CONST_0_1_6 = ((K[2] * F1) + (K[0] * F0) + 1.0);
double CONST_0_1_7 = (((((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0)) + ((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0))) * F4) + (((F5) / 2.0) * ((1.0)) * ((((((K[3] * F1) + (K[1] * F0)) * ((K[3] * F1) + (K[1] * F0))) + (((K[3] * F3) + (K[1] * F2) + 1.0) * ((K[3] * F3) + (K[1] * F2) + 1.0)) + (-1.0)) / 2.0) + (((((K[2] * F3) + (K[0] * F2)) * ((K[2] * F3) + (K[0] * F2))) + (((K[2] * F1) + (K[0] * F0) + 1.0) * ((K[2] * F1) + (K[0] * F0) + 1.0)) + (-1.0)) / 2.0)) * 2.0));
double CONST_0_1_8 = (((K[3] * F3) + (K[1] * F2) + 1.0) * F4);
double J_0_1_0[8] __attribute__((aligned(32))) = {0.0};
double J_0_1_1[8] __attribute__((aligned(32))) = {0.0};
double J_0_1_2[8] __attribute__((aligned(32))) = {0.0};
double J_0_1_3[8] __attribute__((aligned(32))) = {0.0};
double K_0_3_0[8] __attribute__((aligned(32))) = {0.0};
double K_0_3_1[8] __attribute__((aligned(32))) = {0.0};
double K_0_3_2[8] __attribute__((aligned(32))) = {0.0};
double K_0_4_0[8] __attribute__((aligned(32))) = {0.0};
double K_0_4_1[8] __attribute__((aligned(32))) = {0.0};
double K_0_5_0[8] __attribute__((aligned(32))) = {0.0};
double K_0_5_1[8] __attribute__((aligned(32))) = {0.0};
double K_0_6_0[8] __attribute__((aligned(32))) = {0.0};
double K_0_6_1[8] __attribute__((aligned(32))) = {0.0};
double K_0_6_2[8] __attribute__((aligned(32))) = {0.0};
double K_0_6_3[8] __attribute__((aligned(32))) = {0.0};
for (int k = 0; k < 3; k += 1)
{
J_0_1_0[k] = ((K[0] * FE0_D10[0][k]) + (K[2] * FE0_C1_D01[0][k]));
J_0_1_1[k] = ((K[1] * FE0_D10[0][k]) + (K[3] * FE0_C1_D01[0][k]));
J_0_1_2[k+3] = ((K[1] * FE0_D10[0][k]) + (K[3] * FE0_C1_D01[0][k]));
J_0_1_3[k+3] = ((K[0] * FE0_D10[0][k]) + (K[2] * FE0_C1_D01[0][k]));
}
for (int k = 0; k < 6; k += 1)
{
K_0_3_0[k] = (((J_0_1_2[k]) * CONST_0_1_2) + ((J_0_1_0[k]) * CONST_0_1_1) + ((J_0_1_3[k]) * CONST_0_1_0) + ((J_0_1_1[k]) * CONST_0_1_6));
K_0_3_1[k] = (((J_0_1_3[k]) * CONST_0_1_2) + ((J_0_1_3[k]) * CONST_0_1_2) + ((J_0_1_0[k]) * CONST_0_1_6) + ((J_0_1_0[k]) * CONST_0_1_6));
K_0_3_2[k] = (((J_0_1_1[k]) * CONST_0_1_1) + ((J_0_1_1[k]) * CONST_0_1_1) + ((J_0_1_2[k]) * CONST_0_1_0) + ((J_0_1_2[k]) * CONST_0_1_0));
K_0_4_0[k] = (((K_0_3_0[k] / 2.0)) + ((K_0_3_0[k] / 2.0)));
K_0_4_1[k] = ((K_0_3_2[k] / 2.0) + (K_0_3_1[k] / 2.0));
K_0_5_0[k] = (((((K_0_3_1[k] / 2.0)) + ((K_0_3_1[k] / 2.0))) * F4) + (((F5) / 2.0) * ((1.0)) * K_0_4_1[k] * 2.0));
K_0_5_1[k] = (((((K_0_3_2[k] / 2.0)) + ((K_0_3_2[k] / 2.0))) * F4) + (((F5) / 2.0) * ((1.0)) * K_0_4_1[k] * 2.0));
K_0_6_0[k] = (((J_0_1_2[k]) * CONST_0_1_4) + (K_0_4_0[k] * CONST_0_1_8) + (CONST_0_1_2 * K_0_5_0[k]) + ((J_0_1_3[k]) * CONST_0_1_5)) * det * W1;
K_0_6_1[k] = (((J_0_1_0[k]) * CONST_0_1_4) + (K_0_4_0[k] * CONST_0_1_3) + (CONST_0_1_1 * K_0_5_1[k]) + ((J_0_1_1[k]) * CONST_0_1_7)) * det * W1;
K_0_6_2[k] = ((CONST_0_1_2 * K_0_4_0[k] * F4) + ((J_0_1_3[k]) * CONST_0_1_4) + ((J_0_1_2[k]) * CONST_0_1_7) + (K_0_5_1[k] * CONST_0_1_0)) * det * W1;
K_0_6_3[k] = ((CONST_0_1_1 * K_0_4_0[k] * F4) + ((J_0_1_1[k]) * CONST_0_1_4) + ((J_0_1_0[k]) * CONST_0_1_5) + (K_0_5_0[k] * CONST_0_1_6)) * det * W1;
}
for (int j = 0; j < 3; j += 1)
{
for (int k = 0; k < 6; k += 1)
{
_A[0][j][k] += ((K_0_6_3[k] * J_0_1_0[j])) + ((K_0_6_1[k] * J_0_1_1[j]));
_A[0][j+3][k] += ((K_0_6_2[k] * J_0_1_2[j+3])) + ((K_0_6_0[k] * J_0_1_3[j+3]));
}
}
----------------------------------------------------------------------------------
# The following Python program contains ILP formulations for the sharing elimination algorithm
# extracted from both generic graphs and real variational formulations.
# In particular, Problem 8 derives from the "Pressure equation" used for experimentation
# in "Optimizations for quadrature representations of finite element tensors through
# automated code generation", Ølgaard and Wells (http://dl.acm.org/citation.cfm?id=1644009)
from pulp import *
from collections import defaultdict
def solver(VERTICES, EDGES, MAX):
# declare your variables
x = LpVariable.dicts('x', VERTICES, 0, 1, LpBinary)
y = LpVariable.dicts('y', [(i, j) for i, j in EDGES] + [(j, i) for i, j in EDGES], 0, 1, LpBinary)
# defines the problem
prob = LpProblem("Factorizer", LpMinimize)
# defines the constraints
for i in VERTICES:
prob += lpSum(y[(i, j)] for j in VERTICES if (i, j) in y) <= MAX[i]*x[i]
for i, j in EDGES:
prob += y[(i, j)] + y[(j, i)] == 1
# defines the objective function to maximize
prob += lpSum(x[i] for i in VERTICES)
# solve the problem
status = prob.solve(GLPK(msg=0))
LpStatus[status]
for i in prob.variables():
print (i.name, "=", i.varValue)
print "PROBLEM 1:"
VERTICES = [0, 1, 2, 3, 4, 5, 6]
EDGES = [(0, 1), (0, 3), (1, 2), (2, 3), (2, 5), (2, 6), (4, 5), (4, 6)]
MAX = {0: 2, 1: 2, 2: 4, 3: 2, 4: 2, 5: 2, 6: 2}
solver(VERTICES, EDGES, MAX)
print "******\n"
print "PROBLEM 2:"
VERTICES = range(0, 10)
EDGES = [(0, 1), (0, 2), (0, 3), (1, 4), (1, 5), (2, 6), (2, 7), (3, 8), (3, 9)]
MAX = defaultdict(int)
for i, j in EDGES:
MAX[i] += 1
MAX[j] += 1
solver(VERTICES, EDGES, MAX)
print "******\n"
print "PROBLEM 3:"
VERTICES = range(0, 11)
EDGES = [(0, 1), (0, 2), (0, 3), (1, 4), (1, 5), (2, 6), (2, 7), (3, 8), (3, 9),
(5, 10), (6, 10)]
MAX = defaultdict(int)
for i, j in EDGES:
MAX[i] += 1
MAX[j] += 1
solver(VERTICES, EDGES, MAX)
print "******\n"
print "PROBLEM 4:"
VERTICES = range(0, 13)
EDGES = [(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (6, 12), (7, 8),
(7, 9), (6, 11), (9, 10), (11, 10)]
MAX = defaultdict(int)
for i, j in EDGES:
MAX[i] += 1
MAX[j] += 1
solver(VERTICES, EDGES, MAX)
print "******\n"
print "PROBLEM 5:"
VERTICES = range(0, 12)
EDGES = [(0, 1), (0, 2), (1, 3), (2, 3), (0, 9), (9, 10), (9, 11), (0, 4), (0, 5),
(0, 6), (8, 4), (8, 5), (8, 6), (7, 4), (7, 5), (7, 6)]
MAX = defaultdict(int)
for i, j in EDGES:
MAX[i] += 1
MAX[j] += 1
solver(VERTICES, EDGES, MAX)
print "******\n"
print "PROBLEM 6:"
VERTICES = range(0, 11)
EDGES = [(0, 1), (0, 2), (1, 3), (2, 3), (0, 4), (0, 5), (0, 6), (0, 7), (8, 4),
(8, 5), (8, 6), (9, 5), (9, 6), (10, 7), (10, 4)]
MAX = defaultdict(int)
for i, j in EDGES:
MAX[i] += 1
MAX[j] += 1
solver(VERTICES, EDGES, MAX)
print "******\n"
print "PROBLEM 7:"
VERTICES = range(0, 12)
EDGES = [(0, 1), (0, 2), (1, 3), (2, 3), (0, 4), (0, 5), (0, 6), (0, 7), (8, 4),
(8, 5), (8, 6), (9, 5), (9, 6), (10, 7), (10, 4), (7, 11)]
MAX = defaultdict(int)
for i, j in EDGES:
MAX[i] += 1
MAX[j] += 1
solver(VERTICES, EDGES, MAX)
print "******\n"
print "PROBLEM 8:"
VERTICES = range(0, 7)
EDGES = [(0, 3), (0, 4), (0, 5), (0, 6), (1, 3), (1, 4), (1, 5), (1, 6),
(2, 4), (2, 5), (2, 6)]
MAX = defaultdict(int)
for i, j in EDGES:
MAX[i] += 1
MAX[j] += 1
solver(VERTICES, EDGES, MAX)
print "******\n"
----------------------------------------------------------------------------------
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment