#include <cmath>
void ComputeMembrane(
   const double *xa, // [input 3x1] World position of vertex A
   const double *xb, // [input 3x1] World position of vertex B
   const double *xc, // [input 3x1] World position of vertex C
   const double *Xa, // [input 2x1] Material position of vertex A
   const double *Xb, // [input 2x1] Material position of vertex B
   const double *Xc, // [input 2x1] Material position of vertex C
   double e,         // [input 1x1] Young's modulus -pascals
   double nu,        // [input 1x1] Poisson's ratio 0.3 or 0.2
   const double *P,  // [input 2x3] Projection matrix (column major)
   const double *Q,  // [input 2x2] Polar decomposition of \bar{F} (column major)
   double *W,        // [output 1x1] Membrane potential energy
   double *f,        // [output 9x1] Membrane force vector
   double *K)        // [output 9x9] Membrane stiffness matrix
{
double xax = xa[0];
double xay = xa[1];
double xaz = xa[2];
double xbx = xb[0];
double xby = xb[1];
double xbz = xb[2];
double xcx = xc[0];
double xcy = xc[1];
double xcz = xc[2];
double Xax = Xa[0];
double Xay = Xa[1];
double Xbx = Xb[0];
double Xby = Xb[1];
double Xcx = Xc[0];
double Xcy = Xc[1];
// P and Q are stored as column-major matrices
double P00 = P[0];
double P10 = P[1];
double P01 = P[2];
double P11 = P[3];
double P02 = P[4];
double P12 = P[5];
double Q00 = Q[0];
double Q10 = Q[1];
double Q01 = Q[2];
double Q11 = Q[3];
double t7 = Xax * Xby - Xax * Xcy - Xbx * Xay + Xcx * Xay + Xbx * Xcy - Xcx * Xby;
double t8 = t7 / 2;
double t9 = 1 + nu;
double t12 = e / t9 / 2;
double t15 = P00 * Q00 + P10 * Q10;
double t17 = 0.1e1 / t7;
double t18 = (xbx - xax) * t17;
double t19 = Xcy - Xay;
double t22 = (xcx - xax) * t17;
double t23 = -Xby + Xay;
double t25 = t18 * t19 + t22 * t23;
double t26 = t15 * t25;
double t29 = P01 * Q00 + P11 * Q10;
double t31 = (xby - xay) * t17;
double t34 = (xcy - xay) * t17;
double t36 = t31 * t19 + t34 * t23;
double t37 = t29 * t36;
double t40 = P02 * Q00 + P12 * Q10;
double t42 = (xbz - xaz) * t17;
double t45 = (xcz - xaz) * t17;
double t47 = t42 * t19 + t45 * t23;
double t48 = t40 * t47;
double t49 = t26 + t37 + t48 - 1;
double t50 = t49 * t49;
double t53 = P00 * Q01 + P10 * Q11;
double t57 = P01 * Q01 + P11 * Q11;
double t61 = P02 * Q01 + P12 * Q11;
double t63 = t53 * t25 + t57 * t36 + t61 * t47;
double t64 = t63 * t63;
double t65 = -Xcx + Xax;
double t67 = Xbx - Xax;
double t69 = t18 * t65 + t22 * t67;
double t73 = t31 * t65 + t34 * t67;
double t77 = t42 * t65 + t45 * t67;
double t79 = t15 * t69 + t29 * t73 + t40 * t77;
double t80 = t79 * t79;
double t81 = t53 * t69;
double t82 = t57 * t73;
double t83 = t61 * t77;
double t84 = t81 + t82 + t83 - 1;
double t85 = t84 * t84;
double t88 = e * nu;
double t89 = 1 / t9;
double t92 = 1 / (1 - 2 * nu);
double t93 = t89 * t92;
double t94 = t26 + t37 + t48 - 2 + t81 + t82 + t83;
double t95 = t94 * t94;
double W00 = t8 * (t12 * (t50 + t64 + t80 + t85) + t88 * t93 * t95 / 2);
double t100 = t49 * t15;
double t101 = t17 * t19;
double t102 = t17 * t23;
double t103 = -t101 - t102;
double t105 = t63 * t53;
double t107 = t79 * t15;
double t108 = t17 * t65;
double t109 = t17 * t67;
double t110 = -t108 - t109;
double t112 = t84 * t53;
double t117 = t88 * t89;
double t118 = t92 * t94;
double t121 = t15 * t103 + t53 * t110;
double f01 = -t8 * (t12 * (2 * t100 * t103 + 2 * t105 * t103 + 2 * t107 * t110 + 2 * t112 * t110) + t117 * t118 * t121);
double t126 = t49 * t29;
double t128 = t63 * t57;
double t130 = t79 * t29;
double t132 = t84 * t57;
double t139 = t29 * t103 + t57 * t110;
double f02 = -t8 * (t12 * (2 * t126 * t103 + 2 * t128 * t103 + 2 * t130 * t110 + 2 * t132 * t110) + t117 * t118 * t139);
double t144 = t49 * t40;
double t146 = t63 * t61;
double t148 = t79 * t40;
double t150 = t84 * t61;
double t157 = t40 * t103 + t61 * t110;
double f03 = -t8 * (t12 * (2 * t144 * t103 + 2 * t146 * t103 + 2 * t148 * t110 + 2 * t150 * t110) + t117 * t118 * t157);
double t169 = t15 * t17;
double t171 = t53 * t17;
double t173 = t169 * t19 + t171 * t65;
double f04 = -t8 * (t12 * (2 * t100 * t101 + 2 * t105 * t101 + 2 * t107 * t108 + 2 * t112 * t108) + t117 * t118 * t173);
double t185 = t29 * t17;
double t187 = t57 * t17;
double t189 = t185 * t19 + t187 * t65;
double f05 = -t8 * (t12 * (2 * t126 * t101 + 2 * t128 * t101 + 2 * t130 * t108 + 2 * t132 * t108) + t117 * t118 * t189);
double t201 = t40 * t17;
double t203 = t61 * t17;
double t205 = t201 * t19 + t203 * t65;
double f06 = -t8 * (t12 * (2 * t144 * t101 + 2 * t146 * t101 + 2 * t148 * t108 + 2 * t150 * t108) + t117 * t118 * t205);
double t219 = t169 * t23 + t171 * t67;
double f07 = -t8 * (t12 * (2 * t100 * t102 + 2 * t105 * t102 + 2 * t107 * t109 + 2 * t112 * t109) + t117 * t118 * t219);
double t233 = t185 * t23 + t187 * t67;
double f08 = -t8 * (t12 * (2 * t126 * t102 + 2 * t128 * t102 + 2 * t130 * t109 + 2 * t132 * t109) + t117 * t118 * t233);
double t247 = t201 * t23 + t203 * t67;
double f09 = -t8 * (t12 * (2 * t144 * t102 + 2 * t146 * t102 + 2 * t148 * t109 + 2 * t150 * t109) + t117 * t118 * t247);
double t252 = t15 * t15;
double t253 = t103 * t103;
double t255 = t53 * t53;
double t257 = t110 * t110;
double t263 = t121 * t121;
double K0101 = t8 * (t12 * (2 * t252 * t253 + 2 * t252 * t257 + 2 * t255 * t253 + 2 * t255 * t257) + t88 * t93 * t263);
double K0102 = t8 * (t12 * (2 * t29 * t253 * t15 + 2 * t29 * t257 * t15 + 2 * t57 * t253 * t53 + 2 * t57 * t257 * t53) + t117 * t92 * t139 * t121);
double t282 = t40 * t253;
double t284 = t61 * t253;
double t286 = t40 * t257;
double t288 = t61 * t257;
double t293 = t92 * t157;
double K0103 = t8 * (t12 * (2 * t282 * t15 + 2 * t286 * t15 + 2 * t284 * t53 + 2 * t288 * t53) + t117 * t293 * t121);
double t297 = t252 * t17;
double t298 = t19 * t103;
double t300 = t255 * t17;
double t302 = t65 * t110;
double t308 = t92 * t173;
double K0104 = t8 * (t12 * (2 * t297 * t298 + 2 * t297 * t302 + 2 * t300 * t298 + 2 * t300 * t302) + t117 * t308 * t121);
double t313 = t19 * t15 * t103;
double t316 = t19 * t53 * t103;
double t319 = t65 * t15 * t110;
double t322 = t65 * t53 * t110;
double t326 = t12 * (2 * t185 * t313 + 2 * t185 * t319 + 2 * t187 * t316 + 2 * t187 * t322);
double t327 = t92 * t189;
double K0105 = t8 * (t117 * t327 * t121 + t326);
double t337 = t12 * (2 * t201 * t313 + 2 * t201 * t319 + 2 * t203 * t316 + 2 * t203 * t322);
double t338 = t92 * t205;
double K0106 = t8 * (t117 * t338 * t121 + t337);
double t342 = t23 * t103;
double t345 = t67 * t110;
double t351 = t92 * t219;
double K0107 = t8 * (t12 * (2 * t297 * t342 + 2 * t297 * t345 + 2 * t300 * t342 + 2 * t300 * t345) + t117 * t351 * t121);
double t355 = t23 * t15;
double t356 = t355 * t103;
double t358 = t23 * t53;
double t359 = t358 * t103;
double t361 = t67 * t15;
double t362 = t361 * t110;
double t364 = t67 * t53;
double t365 = t364 * t110;
double t369 = t12 * (2 * t185 * t356 + 2 * t185 * t362 + 2 * t187 * t359 + 2 * t187 * t365);
double t370 = t92 * t233;
double K0108 = t8 * (t117 * t370 * t121 + t369);
double t380 = t12 * (2 * t201 * t356 + 2 * t201 * t362 + 2 * t203 * t359 + 2 * t203 * t365);
double t381 = t92 * t247;
double K0109 = t8 * (t117 * t381 * t121 + t380);
double K0201 = K0102;
double t385 = t29 * t29;
double t387 = t57 * t57;
double t394 = t139 * t139;
double K0202 = t8 * (t12 * (2 * t385 * t253 + 2 * t387 * t253 + 2 * t385 * t257 + 2 * t387 * t257) + t88 * t93 * t394);
double K0203 = t8 * (t12 * (2 * t282 * t29 + 2 * t284 * t57 + 2 * t286 * t29 + 2 * t288 * t57) + t117 * t293 * t139);
double K0204 = t8 * (t117 * t308 * t139 + t326);
double t411 = t385 * t17;
double t413 = t387 * t17;
double K0205 = t8 * (t12 * (2 * t411 * t298 + 2 * t413 * t298 + 2 * t411 * t302 + 2 * t413 * t302) + t117 * t327 * t139);
double t437 = t12 * (2 * t201 * t19 * t29 * t103 + 2 * t203 * t19 * t57 * t103 + 2 * t201 * t65 * t29 * t110 + 2 * t203 * t65 * t57 * t110);
double K0206 = t8 * (t117 * t338 * t139 + t437);
double K0207 = t8 * (t117 * t351 * t139 + t369);
double K0208 = t8 * (t12 * (2 * t411 * t342 + 2 * t413 * t342 + 2 * t411 * t345 + 2 * t413 * t345) + t117 * t370 * t139);
double t454 = t23 * t29;
double t457 = t23 * t57;
double t460 = t67 * t29;
double t463 = t67 * t57;
double t468 = t12 * (2 * t201 * t454 * t103 + 2 * t203 * t457 * t103 + 2 * t201 * t460 * t110 + 2 * t203 * t463 * t110);
double K0209 = t8 * (t117 * t381 * t139 + t468);
double K0301 = K0103;
double K0302 = K0203;
double t472 = t40 * t40;
double t474 = t61 * t61;
double t481 = t157 * t157;
double K0303 = t8 * (t12 * (2 * t472 * t253 + 2 * t474 * t253 + 2 * t472 * t257 + 2 * t474 * t257) + t88 * t93 * t481);
double K0304 = t8 * (t117 * t308 * t157 + t337);
double K0305 = t8 * (t117 * t327 * t157 + t437);
double t491 = t472 * t17;
double t493 = t474 * t17;
double K0306 = t8 * (t12 * (2 * t491 * t298 + 2 * t493 * t298 + 2 * t491 * t302 + 2 * t493 * t302) + t117 * t338 * t157);
double K0307 = t8 * (t117 * t351 * t157 + t380);
double K0308 = t8 * (t117 * t370 * t157 + t468);
double K0309 = t8 * (t12 * (2 * t491 * t342 + 2 * t493 * t342 + 2 * t491 * t345 + 2 * t493 * t345) + t117 * t381 * t157);
double K0401 = K0104;
double K0402 = K0204;
double K0403 = K0304;
double t519 = t7 * t7;
double t520 = 0.1e1 / t519;
double t521 = t252 * t520;
double t522 = t19 * t19;
double t524 = t255 * t520;
double t526 = t65 * t65;
double t532 = t173 * t173;
double K0404 = t8 * (t12 * (2 * t521 * t522 + 2 * t521 * t526 + 2 * t524 * t522 + 2 * t524 * t526) + t88 * t93 * t532);
double t536 = t29 * t520;
double t537 = t522 * t15;
double t539 = t57 * t520;
double t540 = t522 * t53;
double t542 = t526 * t15;
double t544 = t526 * t53;
double K0405 = t8 * (t12 * (2 * t536 * t537 + 2 * t536 * t542 + 2 * t539 * t540 + 2 * t539 * t544) + t117 * t327 * t173);
double t552 = t40 * t520;
double t554 = t61 * t520;
double K0406 = t8 * (t12 * (2 * t552 * t537 + 2 * t554 * t540 + 2 * t552 * t542 + 2 * t554 * t544) + t117 * t338 * t173);
double t564 = t23 * t19;
double t567 = t67 * t65;
double K0407 = t8 * (t12 * (2 * t521 * t564 + 2 * t521 * t567 + 2 * t524 * t564 + 2 * t524 * t567) + t117 * t351 * t173);
double t576 = t355 * t19;
double t578 = t358 * t19;
double t580 = t361 * t65;
double t582 = t364 * t65;
double t586 = t12 * (2 * t536 * t576 + 2 * t536 * t580 + 2 * t539 * t578 + 2 * t539 * t582);
double K0408 = t8 * (t117 * t370 * t173 + t586);
double t596 = t12 * (2 * t552 * t576 + 2 * t552 * t580 + 2 * t554 * t578 + 2 * t554 * t582);
double K0409 = t8 * (t117 * t381 * t173 + t596);
double K0501 = K0105;
double K0502 = K0205;
double K0503 = K0305;
double K0504 = K0405;
double t600 = t385 * t520;
double t602 = t387 * t520;
double t609 = t189 * t189;
double K0505 = t8 * (t12 * (2 * t600 * t522 + 2 * t602 * t522 + 2 * t600 * t526 + 2 * t602 * t526) + t88 * t93 * t609);
double K0506 = t8 * (t12 * (2 * t552 * t522 * t29 + 2 * t552 * t526 * t29 + 2 * t554 * t522 * t57 + 2 * t554 * t526 * t57) + t117 * t338 * t189);
double K0507 = t8 * (t117 * t351 * t189 + t586);
double K0508 = t8 * (t12 * (2 * t600 * t564 + 2 * t602 * t564 + 2 * t600 * t567 + 2 * t602 * t567) + t117 * t370 * t189);
double t650 = t12 * (2 * t552 * t454 * t19 + 2 * t554 * t457 * t19 + 2 * t552 * t460 * t65 + 2 * t554 * t463 * t65);
double K0509 = t8 * (t117 * t381 * t189 + t650);
double K0601 = K0106;
double K0602 = K0206;
double K0603 = K0306;
double K0604 = K0406;
double K0605 = K0506;
double t654 = t472 * t520;
double t656 = t474 * t520;
double t663 = t205 * t205;
double K0606 = t8 * (t12 * (2 * t654 * t522 + 2 * t656 * t522 + 2 * t654 * t526 + 2 * t656 * t526) + t88 * t93 * t663);
double K0607 = t8 * (t117 * t351 * t205 + t596);
double K0608 = t8 * (t117 * t370 * t205 + t650);
double K0609 = t8 * (t12 * (2 * t654 * t564 + 2 * t656 * t564 + 2 * t654 * t567 + 2 * t656 * t567) + t117 * t381 * t205);
double K0701 = K0107;
double K0702 = K0207;
double K0703 = K0307;
double K0704 = K0407;
double K0705 = K0507;
double K0706 = K0607;
double t683 = t23 * t23;
double t686 = t67 * t67;
double t692 = t219 * t219;
double K0707 = t8 * (t12 * (2 * t521 * t683 + 2 * t521 * t686 + 2 * t524 * t683 + 2 * t524 * t686) + t88 * t93 * t692);
double t696 = t683 * t15;
double t698 = t683 * t53;
double t700 = t686 * t15;
double t702 = t686 * t53;
double K0708 = t8 * (t12 * (2 * t536 * t696 + 2 * t536 * t700 + 2 * t539 * t698 + 2 * t539 * t702) + t117 * t370 * t219);
double K0709 = t8 * (t12 * (2 * t552 * t696 + 2 * t552 * t700 + 2 * t554 * t698 + 2 * t554 * t702) + t117 * t381 * t219);
double K0801 = K0108;
double K0802 = K0208;
double K0803 = K0308;
double K0804 = K0408;
double K0805 = K0508;
double K0806 = K0608;
double K0807 = K0708;
double t727 = t233 * t233;
double K0808 = t8 * (t12 * (2 * t600 * t683 + 2 * t600 * t686 + 2 * t602 * t683 + 2 * t602 * t686) + t88 * t93 * t727);
double K0809 = t8 * (t12 * (2 * t552 * t683 * t29 + 2 * t552 * t686 * t29 + 2 * t554 * t683 * t57 + 2 * t554 * t686 * t57) + t117 * t381 * t233);
double K0901 = K0109;
double K0902 = K0209;
double K0903 = K0309;
double K0904 = K0409;
double K0905 = K0509;
double K0906 = K0609;
double K0907 = K0709;
double K0908 = K0809;
double t752 = t247 * t247;
double K0909 = t8 * (t12 * (2 * t654 * t683 + 2 * t654 * t686 + 2 * t656 * t683 + 2 * t656 * t686) + t88 * t93 * t752);
W[0]=W00;
f[0]=f01; f[1]=f02; f[2]=f03; f[3]=f04; f[4]=f05; f[5]=f06; f[6]=f07; f[7]=f08; f[8]=f09; 
K[ 0*9+ 0]=K0101; K[ 0*9+ 1]=K0102; K[ 0*9+ 2]=K0103; K[ 0*9+ 3]=K0104; K[ 0*9+ 4]=K0105; K[ 0*9+ 5]=K0106; K[ 0*9+ 6]=K0107; K[ 0*9+ 7]=K0108; K[ 0*9+ 8]=K0109; 
K[ 1*9+ 0]=K0201; K[ 1*9+ 1]=K0202; K[ 1*9+ 2]=K0203; K[ 1*9+ 3]=K0204; K[ 1*9+ 4]=K0205; K[ 1*9+ 5]=K0206; K[ 1*9+ 6]=K0207; K[ 1*9+ 7]=K0208; K[ 1*9+ 8]=K0209; 
K[ 2*9+ 0]=K0301; K[ 2*9+ 1]=K0302; K[ 2*9+ 2]=K0303; K[ 2*9+ 3]=K0304; K[ 2*9+ 4]=K0305; K[ 2*9+ 5]=K0306; K[ 2*9+ 6]=K0307; K[ 2*9+ 7]=K0308; K[ 2*9+ 8]=K0309; 
K[ 3*9+ 0]=K0401; K[ 3*9+ 1]=K0402; K[ 3*9+ 2]=K0403; K[ 3*9+ 3]=K0404; K[ 3*9+ 4]=K0405; K[ 3*9+ 5]=K0406; K[ 3*9+ 6]=K0407; K[ 3*9+ 7]=K0408; K[ 3*9+ 8]=K0409; 
K[ 4*9+ 0]=K0501; K[ 4*9+ 1]=K0502; K[ 4*9+ 2]=K0503; K[ 4*9+ 3]=K0504; K[ 4*9+ 4]=K0505; K[ 4*9+ 5]=K0506; K[ 4*9+ 6]=K0507; K[ 4*9+ 7]=K0508; K[ 4*9+ 8]=K0509; 
K[ 5*9+ 0]=K0601; K[ 5*9+ 1]=K0602; K[ 5*9+ 2]=K0603; K[ 5*9+ 3]=K0604; K[ 5*9+ 4]=K0605; K[ 5*9+ 5]=K0606; K[ 5*9+ 6]=K0607; K[ 5*9+ 7]=K0608; K[ 5*9+ 8]=K0609; 
K[ 6*9+ 0]=K0701; K[ 6*9+ 1]=K0702; K[ 6*9+ 2]=K0703; K[ 6*9+ 3]=K0704; K[ 6*9+ 4]=K0705; K[ 6*9+ 5]=K0706; K[ 6*9+ 6]=K0707; K[ 6*9+ 7]=K0708; K[ 6*9+ 8]=K0709; 
K[ 7*9+ 0]=K0801; K[ 7*9+ 1]=K0802; K[ 7*9+ 2]=K0803; K[ 7*9+ 3]=K0804; K[ 7*9+ 4]=K0805; K[ 7*9+ 5]=K0806; K[ 7*9+ 6]=K0807; K[ 7*9+ 7]=K0808; K[ 7*9+ 8]=K0809; 
K[ 8*9+ 0]=K0901; K[ 8*9+ 1]=K0902; K[ 8*9+ 2]=K0903; K[ 8*9+ 3]=K0904; K[ 8*9+ 4]=K0905; K[ 8*9+ 5]=K0906; K[ 8*9+ 6]=K0907; K[ 8*9+ 7]=K0908; K[ 8*9+ 8]=K0909; 
}