Skip to content

Instantly share code, notes, and snippets.

@wence-

wence-/gist:6ad1df6d3d68b1b300f8

Created May 5, 2014 15:45
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save wence-/6ad1df6d3d68b1b300f8 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save wence-/6ad1df6d3d68b1b300f8 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gistfile1.c
#include <math.h>
#define compute_jacobian_inverse_triangle_3d(K, det, J) \
do { const double d_0 = J[2]*J[5] - J[4]*J[3]; \
const double d_1 = J[4]*J[1] - J[0]*J[5]; \
const double d_2 = J[0]*J[3] - J[2]*J[1]; \
const double c_0 = J[0]*J[0] + J[2]*J[2] + J[4]*J[4]; \
const double c_1 = J[1]*J[1] + J[3]*J[3] + J[5]*J[5]; \
const double c_2 = J[0]*J[1] + J[2]*J[3] + J[4]*J[5]; \
const double den = c_0*c_1 - c_2*c_2; \
const double det2 = d_0*d_0 + d_1*d_1 + d_2*d_2; \
det = sqrt(det2); \
K[0] = (J[0]*c_1 - J[1]*c_2) / den; \
K[1] = (J[2]*c_1 - J[3]*c_2) / den; \
K[2] = (J[4]*c_1 - J[5]*c_2) / den; \
K[3] = (J[1]*c_0 - J[0]*c_2) / den; \
K[4] = (J[3]*c_0 - J[2]*c_2) / den; \
K[5] = (J[5]*c_0 - J[4]*c_2) / den; } while (0)
static inline void form_interior_facet_integral_0_otherwise (double A[6][6] , double **vertex_coordinates , int **cell_orientation_ , unsigned int facet_p[2] )
{
unsigned int facet_0 = facet_p[0];
unsigned int facet_1 = facet_p[1];
double **vertex_coordinates_0 = vertex_coordinates;
double **vertex_coordinates_1 = vertex_coordinates + 3;
const int cell_orientation_0 = cell_orientation_[0][0];
const int cell_orientation_1 = cell_orientation_[1][0];
// Compute Jacobian
double J_0[6];
J_0[0] = vertex_coordinates_0[1] [0] - vertex_coordinates_0[0] [0];
J_0[1] = vertex_coordinates_0[2] [0] - vertex_coordinates_0[0] [0];
J_0[2] = vertex_coordinates_0[7] [0] - vertex_coordinates_0[6] [0];
J_0[3] = vertex_coordinates_0[8] [0] - vertex_coordinates_0[6] [0];
J_0[4] = vertex_coordinates_0[13][0] - vertex_coordinates_0[12][0];
J_0[5] = vertex_coordinates_0[14][0] - vertex_coordinates_0[12][0];
// Compute Jacobian inverse and determinant
double K_0[6];
double detJ_0;
compute_jacobian_inverse_triangle_3d(K_0, detJ_0, J_0);
if (cell_orientation_0 == 1)
detJ_0 *= -1;
// Compute Jacobian
double J_1[6];
J_1[0] = vertex_coordinates_1[1] [0] - vertex_coordinates_1[0] [0];
J_1[1] = vertex_coordinates_1[2] [0] - vertex_coordinates_1[0] [0];
J_1[2] = vertex_coordinates_1[7] [0] - vertex_coordinates_1[6] [0];
J_1[3] = vertex_coordinates_1[8] [0] - vertex_coordinates_1[6] [0];
J_1[4] = vertex_coordinates_1[13][0] - vertex_coordinates_1[12][0];
J_1[5] = vertex_coordinates_1[14][0] - vertex_coordinates_1[12][0];
// Compute Jacobian inverse and determinant
double K_1[6];
double detJ_1;
compute_jacobian_inverse_triangle_3d(K_1, detJ_1, J_1);
if (cell_orientation_1 == 1)
detJ_1 *= -1;
// Facet determinant 2D in 3D (edge)
// Get vertices on edge
unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
const unsigned int v0 = edge_vertices[facet_0][0];
const unsigned int v1 = edge_vertices[facet_0][1];
// Compute scale factor (length of edge scaled by length of reference interval)
const double dx0 = vertex_coordinates_0[v1 + 0][0] - vertex_coordinates_0[v0 + 0][0];
const double dx1 = vertex_coordinates_0[v1 + 6][0] - vertex_coordinates_0[v0 + 6][0];
const double dx2 = vertex_coordinates_0[v1 + 12][0] - vertex_coordinates_0[v0 + 12][0];
const double det = sqrt(dx0*dx0 + dx1*dx1 + dx2*dx2);
// Compute facet normal for triangles in 3D
const unsigned int vertex_00 = facet_0;
// Get coordinates corresponding the vertex opposite this
// static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
const unsigned int vertex_01 = edge_vertices[facet_0][0];
const unsigned int vertex_02 = edge_vertices[facet_0][1];
// Define vectors n = (p2 - p0) and t = normalized (p2 - p1)
double n_00 = vertex_coordinates_0[vertex_02 + 0][0] - vertex_coordinates_0[vertex_00 + 0][0];
double n_01 = vertex_coordinates_0[vertex_02 + 6][0] - vertex_coordinates_0[vertex_00 + 6][0];
double n_02 = vertex_coordinates_0[vertex_02 + 12][0] - vertex_coordinates_0[vertex_00 + 12][0];
double t_00 = vertex_coordinates_0[vertex_02 + 0][0] - vertex_coordinates_0[vertex_01 + 0][0];
double t_01 = vertex_coordinates_0[vertex_02 + 6][0] - vertex_coordinates_0[vertex_01 + 6][0];
double t_02 = vertex_coordinates_0[vertex_02 + 12][0] - vertex_coordinates_0[vertex_01 + 12][0];
const double t_0_length = sqrt(t_00*t_00 + t_01*t_01 + t_02*t_02);
t_00 /= t_0_length;
t_01 /= t_0_length;
t_02 /= t_0_length;
// Subtract, the projection of (p2 - p0) onto (p2 - p1), from (p2 - p0)
const double ndott_0 = t_00*n_00 + t_01*n_01 + t_02*n_02;
n_00 -= ndott_0*t_00;
n_01 -= ndott_0*t_01;
n_02 -= ndott_0*t_02;
const double n_0_length = sqrt(n_00*n_00 + n_01*n_01 + n_02*n_02);
// Normalize
n_00 /= n_0_length;
n_01 /= n_0_length;
n_02 /= n_0_length;
// Compute facet normal for triangles in 3D
const unsigned int vertex_10 = facet_1;
// Get coordinates corresponding the vertex opposite this
// static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
const unsigned int vertex_11 = edge_vertices[facet_1][0];
const unsigned int vertex_12 = edge_vertices[facet_1][1];
// Define vectors n = (p2 - p0) and t = normalized (p2 - p1)
double n_10 = vertex_coordinates_1[vertex_12 + 0][0] - vertex_coordinates_1[vertex_10 + 0][0];
double n_11 = vertex_coordinates_1[vertex_12 + 6][0] - vertex_coordinates_1[vertex_10 + 6][0];
double n_12 = vertex_coordinates_1[vertex_12 + 12][0] - vertex_coordinates_1[vertex_10 + 12][0];
double t_10 = vertex_coordinates_1[vertex_12 + 0][0] - vertex_coordinates_1[vertex_11 + 0][0];
double t_11 = vertex_coordinates_1[vertex_12 + 6][0] - vertex_coordinates_1[vertex_11 + 6][0];
double t_12 = vertex_coordinates_1[vertex_12 + 12][0] - vertex_coordinates_1[vertex_11 + 12][0];
const double t_1_length = sqrt(t_10*t_10 + t_11*t_11 + t_12*t_12);
t_10 /= t_1_length;
t_11 /= t_1_length;
t_12 /= t_1_length;
// Subtract, the projection of (p2 - p0) onto (p2 - p1), from (p2 - p0)
const double ndott_1 = t_10*n_10 + t_11*n_11 + t_12*n_12;
n_10 -= ndott_1*t_10;
n_11 -= ndott_1*t_11;
n_12 -= ndott_1*t_12;
const double n_1_length = sqrt(n_10*n_10 + n_11*n_11 + n_12*n_12);
// Normalize
n_10 /= n_1_length;
n_11 /= n_1_length;
n_12 /= n_1_length;
// Facet area
// Cell volume of triangle in 3D
static const double W2[2] = {0.5, 0.5};
static const double FE0_f0_C1[2][3] = {{0.211324865405187, -0.211324865405187, -0.788675134594813},
{0.788675134594813, -0.788675134594813, -0.211324865405187}};
static const double FE0_f0_C0[2][3] = {{0.788675134594813, 0.211324865405187, 0.788675134594813},
{0.211324865405187, 0.788675134594813, 0.211324865405187}};
static const double FE0_f1_C0[2][3] = {{0.0, 1.0, 0.0},
{0.0, 1.0, 0.0}};
static const double FE0_f2_C0[2][3] = {{0.211324865405187, 0.788675134594813, 0.211324865405187},
{0.788675134594813, 0.211324865405187, 0.788675134594813}};
static const double FE0_f2_C1[2][3] = {{0.0, 0.0, -1.0},
{0.0, 0.0, -1.0}};
for (int ip = 0; ip < 2; ++ip)
{
for (int j = 0; j < 3; ++j)
{
for (int k = 0; k < 3; ++k)
{
A[j][k] += ((((((1.0 / detJ_0) * J_0[5] * FE0_f0_C1[ip][k]) + ((1.0 / detJ_0) * J_0[4] * FE0_f0_C0[ip][k])) * n_02) + ((((1.0 / detJ_0) * J_0[3] * FE0_f0_C1[ip][k]) + ((1.0 / detJ_0) * J_0[2] * FE0_f0_C0[ip][k])) * n_01) + ((((1.0 / detJ_0) * J_0[1] * FE0_f0_C1[ip][k]) + ((1.0 / detJ_0) * J_0[0] * FE0_f0_C0[ip][k])) * n_00)) * (((((1.0 / detJ_0) * J_0[5] * FE0_f0_C1[ip][j]) + ((1.0 / detJ_0) * J_0[4] * FE0_f0_C0[ip][j])) * n_02) + ((((1.0 / detJ_0) * J_0[3] * FE0_f0_C1[ip][j]) + ((1.0 / detJ_0) * J_0[2] * FE0_f0_C0[ip][j])) * n_01) + ((((1.0 / detJ_0) * J_0[1] * FE0_f0_C1[ip][j]) + ((1.0 / detJ_0) * J_0[0] * FE0_f0_C0[ip][j])) * n_00)) * det * W2[ip]);
/*
* }
*
* }
* for (int j = 0; j < 3; ++j)
* {
* for (int k = 0; k < 3; ++k)
* {
*/
A[(j + 3)][(k + 3)] += ((((((1.0 / detJ_1) * J_1[5] * FE0_f0_C1[ip][k]) + ((1.0 / detJ_1) * J_1[4] * FE0_f0_C0[ip][k])) * n_12) + ((((1.0 / detJ_1) * J_1[3] * FE0_f0_C1[ip][k]) + ((1.0 / detJ_1) * J_1[2] * FE0_f0_C0[ip][k])) * n_11) + ((((1.0 / detJ_1) * J_1[1] * FE0_f0_C1[ip][k]) + ((1.0 / detJ_1) * J_1[0] * FE0_f0_C0[ip][k])) * n_10)) * (((((1.0 / detJ_1) * J_1[5] * FE0_f0_C1[ip][j]) + ((1.0 / detJ_1) * J_1[4] * FE0_f0_C0[ip][j])) * n_12) + ((((1.0 / detJ_1) * J_1[3] * FE0_f0_C1[ip][j]) + ((1.0 / detJ_1) * J_1[2] * FE0_f0_C0[ip][j])) * n_11) + ((((1.0 / detJ_1) * J_1[1] * FE0_f0_C1[ip][j]) + ((1.0 / detJ_1) * J_1[0] * FE0_f0_C0[ip][j])) * n_10)) * det * W2[ip]);
}
}
}
}
void wrap_form_interior_facet_integral_0_otherwise(int start, int end,
int *arg0_0_map0_0, int *arg0_0_map1_0, double *arg1_0,
int *arg1_0_map0_0, int *arg2_0, int *arg2_0_map0_0, unsigned int *arg3_0
) {
double *arg1_0_vec[18];
int *arg2_0_vec[2];
for ( int n = start; n < end; n++ ) {
int i = n;
arg1_0_vec[0] = arg1_0 + (arg1_0_map0_0[i * 6 + 0])* 3;
arg1_0_vec[1] = arg1_0 + (arg1_0_map0_0[i * 6 + 1])* 3;
arg1_0_vec[2] = arg1_0 + (arg1_0_map0_0[i * 6 + 2])* 3;
arg1_0_vec[3] = arg1_0 + (arg1_0_map0_0[i * 6 + 3])* 3;
arg1_0_vec[4] = arg1_0 + (arg1_0_map0_0[i * 6 + 4])* 3;
arg1_0_vec[5] = arg1_0 + (arg1_0_map0_0[i * 6 + 5])* 3;
arg1_0_vec[6] = arg1_0 + (arg1_0_map0_0[i * 6 + 0])* 3 + 1;
arg1_0_vec[7] = arg1_0 + (arg1_0_map0_0[i * 6 + 1])* 3 + 1;
arg1_0_vec[8] = arg1_0 + (arg1_0_map0_0[i * 6 + 2])* 3 + 1;
arg1_0_vec[9] = arg1_0 + (arg1_0_map0_0[i * 6 + 3])* 3 + 1;
arg1_0_vec[10] = arg1_0 + (arg1_0_map0_0[i * 6 + 4])* 3 + 1;
arg1_0_vec[11] = arg1_0 + (arg1_0_map0_0[i * 6 + 5])* 3 + 1;
arg1_0_vec[12] = arg1_0 + (arg1_0_map0_0[i * 6 + 0])* 3 + 2;
arg1_0_vec[13] = arg1_0 + (arg1_0_map0_0[i * 6 + 1])* 3 + 2;
arg1_0_vec[14] = arg1_0 + (arg1_0_map0_0[i * 6 + 2])* 3 + 2;
arg1_0_vec[15] = arg1_0 + (arg1_0_map0_0[i * 6 + 3])* 3 + 2;
arg1_0_vec[16] = arg1_0 + (arg1_0_map0_0[i * 6 + 4])* 3 + 2;
arg1_0_vec[17] = arg1_0 + (arg1_0_map0_0[i * 6 + 5])* 3 + 2;
arg2_0_vec[0] = arg2_0 + (arg2_0_map0_0[i * 2 + 0])* 1;
arg2_0_vec[1] = arg2_0 + (arg2_0_map0_0[i * 2 + 1])* 1;
double buffer_arg0_0[6][6] = {{0}};
form_interior_facet_integral_0_otherwise(buffer_arg0_0, arg1_0_vec, arg2_0_vec, arg3_0 + i * 2);
}
}
int main(void)
{
int start = 0;
int end = 30;
int arg0_0_map0_0[30][6] = {{ 1, 2, 0, 9, 8, 0},
{ 1, 2, 0, 13, 12, 1},
{ 1, 2, 0, 3, 4, 2},
{ 3, 4, 2, 11, 10, 3},
{ 3, 4, 2, 5, 6, 4},
{ 5, 6, 4, 19, 18, 5},
{ 5, 6, 4, 7, 8, 6},
{ 7, 8, 6, 17, 16, 7},
{ 7, 8, 6, 9, 8, 0},
{ 9, 8, 0, 14, 15, 9},
{11, 10, 3, 21, 13, 10},
{11, 10, 3, 29, 18, 11},
{13, 12, 1, 23, 15, 12},
{13, 12, 1, 21, 13, 10},
{14, 15, 9, 25, 17, 14},
{14, 15, 9, 23, 15, 12},
{17, 16, 7, 27, 19, 16},
{17, 16, 7, 25, 17, 14},
{19, 18, 5, 29, 18, 11},
{19, 18, 5, 27, 19, 16},
{22, 20, 21, 28, 20, 29},
{22, 20, 21, 21, 13, 10},
{22, 20, 21, 24, 22, 23},
{24, 22, 23, 23, 15, 12},
{24, 22, 23, 26, 24, 25},
{26, 24, 25, 25, 17, 14},
{26, 24, 25, 28, 26, 27},
{28, 26, 27, 27, 19, 16},
{28, 26, 27, 28, 20, 29},
{28, 20, 29, 29, 18, 11}};
int arg0_0_map1_0[30][6] = {{ 1, 2, 0, 9, 8, 0},
{ 1, 2, 0, 13, 12, 1},
{ 1, 2, 0, 3, 4, 2},
{ 3, 4, 2, 11, 10, 3},
{ 3, 4, 2, 5, 6, 4},
{ 5, 6, 4, 19, 18, 5},
{ 5, 6, 4, 7, 8, 6},
{ 7, 8, 6, 17, 16, 7},
{ 7, 8, 6, 9, 8, 0},
{ 9, 8, 0, 14, 15, 9},
{11, 10, 3, 21, 13, 10},
{11, 10, 3, 29, 18, 11},
{13, 12, 1, 23, 15, 12},
{13, 12, 1, 21, 13, 10},
{14, 15, 9, 25, 17, 14},
{14, 15, 9, 23, 15, 12},
{17, 16, 7, 27, 19, 16},
{17, 16, 7, 25, 17, 14},
{19, 18, 5, 29, 18, 11},
{19, 18, 5, 27, 19, 16},
{22, 20, 21, 28, 20, 29},
{22, 20, 21, 21, 13, 10},
{22, 20, 21, 24, 22, 23},
{24, 22, 23, 23, 15, 12},
{24, 22, 23, 26, 24, 25},
{26, 24, 25, 25, 17, 14},
{26, 24, 25, 28, 26, 27},
{28, 26, 27, 27, 19, 16},
{28, 26, 27, 28, 20, 29},
{28, 20, 29, 29, 18, 11}};
double arg1_0[12][3] = {{-0.52573111, 0.85065081, 0. },
{-0.85065081, 0. , 0.52573111},
{ 0. , 0.52573111, 0.85065081},
{ 0.52573111, 0.85065081, 0. },
{ 0. , 0.52573111, -0.85065081},
{-0.85065081, 0. , -0.52573111},
{ 0.85065081, 0. , 0.52573111},
{ 0. , -0.52573111, 0.85065081},
{-0.52573111, -0.85065081, 0. },
{ 0. , -0.52573111, -0.85065081},
{ 0.85065081, 0. , -0.52573111},
{ 0.52573111, -0.85065081, 0. }};
int arg1_0_map0_0[30][6] = {{ 0, 1, 2, 0, 1, 5},
{ 0, 1, 2, 1, 2, 7},
{ 0, 1, 2, 0, 2, 3},
{ 0, 2, 3, 2, 3, 6},
{ 0, 2, 3, 0, 3, 4},
{ 0, 3, 4, 3, 4, 10},
{ 0, 3, 4, 0, 4, 5},
{ 0, 4, 5, 4, 5, 9},
{ 0, 4, 5, 0, 1, 5},
{ 0, 1, 5, 1, 5, 8},
{ 2, 3, 6, 2, 6, 7},
{ 2, 3, 6, 3, 6, 10},
{ 1, 2, 7, 1, 7, 8},
{ 1, 2, 7, 2, 6, 7},
{ 1, 5, 8, 5, 8, 9},
{ 1, 5, 8, 1, 7, 8},
{ 4, 5, 9, 4, 9, 10},
{ 4, 5, 9, 5, 8, 9},
{ 3, 4, 10, 3, 6, 10},
{ 3, 4, 10, 4, 9, 10},
{ 6, 7, 11, 6, 10, 11},
{ 6, 7, 11, 2, 6, 7},
{ 6, 7, 11, 7, 8, 11},
{ 7, 8, 11, 1, 7, 8},
{ 7, 8, 11, 8, 9, 11},
{ 8, 9, 11, 5, 8, 9},
{ 8, 9, 11, 9, 10, 11},
{ 9, 10, 11, 4, 9, 10},
{ 9, 10, 11, 6, 10, 11},
{ 6, 10, 11, 3, 6, 10}};
int arg2_0[20] = {0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0};
int arg2_0_map0_0[30][2] = {{ 0, 4},
{ 0, 6},
{ 0, 1},
{ 1, 5},
{ 1, 2},
{ 2, 9},
{ 2, 3},
{ 3, 8},
{ 3, 4},
{ 4, 7},
{ 5, 15},
{ 5, 19},
{ 6, 16},
{ 6, 15},
{ 7, 17},
{ 7, 16},
{ 8, 18},
{ 8, 17},
{ 9, 19},
{ 9, 18},
{10, 14},
{10, 15},
{10, 11},
{11, 16},
{11, 12},
{12, 17},
{12, 13},
{13, 18},
{13, 14},
{14, 19}};
unsigned int arg3_0[30][2] = {{2, 2},
{0, 2},
{1, 2},
{0, 2},
{1, 2},
{0, 2},
{1, 2},
{0, 2},
{1, 1},
{0, 2},
{1, 2},
{0, 2},
{1, 2},
{0, 1},
{0, 2},
{1, 1},
{1, 2},
{0, 1},
{1, 1},
{0, 1},
{1, 1},
{2, 0},
{0, 1},
{2, 0},
{0, 1},
{2, 0},
{0, 1},
{2, 0},
{0, 0},
{2, 0}};
wrap_form_interior_facet_integral_0_otherwise(start, end,
(int *)arg0_0_map0_0, (int *)arg0_0_map1_0,
(double *)arg1_0,
(int *)arg1_0_map0_0,
(int *)arg2_0,
(int *)arg2_0_map0_0,
(unsigned int *)arg3_0);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment