Incompressibility constraint

The incompressibility constraint is defined at each cell center:

\[\frac{1}{J} \dif{}{\gx} \left( \jhx \ux \right) + \frac{1}{J} \dif{}{\gy} \left( \jhy \uy \right) + \frac{1}{J} \dif{}{\gz} \left( \jhz \uz \right) = 0.\]
src/logging/divergence.c
57const double hx_xm = HXXF(i  );
58const double hx_xp = HXXF(i+1);
59const double jd_xm = JDXF(i  );
60const double jd_x0 = JDXC(i  );
61const double jd_xp = JDXF(i+1);
62const double ux_xm = UX(i  , j  );
63const double ux_xp = UX(i+1, j  );
64const double uy_ym = UY(i  , j  );
65const double uy_yp = UY(i  , j+1);
66const double div = 1. / jd_x0 * (
67    - jd_xm / hx_xm * ux_xm + jd_xp / hx_xp * ux_xp
68    - jd_x0 / hy    * uy_ym + jd_x0 / hy    * uy_yp
69);
src/logging/divergence.c
79const double hx_xm = HXXF(i  );
80const double hx_xp = HXXF(i+1);
81const double jd_xm = JDXF(i  );
82const double jd_x0 = JDXC(i  );
83const double jd_xp = JDXF(i+1);
84const double ux_xm = UX(i  , j  , k  );
85const double ux_xp = UX(i+1, j  , k  );
86const double uy_ym = UY(i  , j  , k  );
87const double uy_yp = UY(i  , j+1, k  );
88const double uz_zm = UZ(i  , j  , k  );
89const double uz_zp = UZ(i  , j  , k+1);
90const double div = 1. / jd_x0 * (
91    - jd_xm / hx_xm * ux_xm + jd_xp / hx_xp * ux_xp
92    - jd_x0 / hy    * uy_ym + jd_x0 / hy    * uy_yp
93    - jd_x0 / hz    * uz_zm + jd_x0 / hz    * uz_zp
94);
src/fluid/compute_potential.c
483const double hx_xm = HXXF(i  );
484const double hx_xp = HXXF(i+1);
485const double jd_xm = JDXF(i  );
486const double jd_x0 = JDXC(i  );
487const double jd_xp = JDXF(i+1);
488const double ux_xm = UX(i  , j  );
489const double ux_xp = UX(i+1, j  );
490const double uy_ym = UY(i  , j  );
491const double uy_yp = UY(i  , j+1);
492const double div = 1. / jd_x0 * (
493    - jd_xm / hx_xm * ux_xm + jd_xp / hx_xp * ux_xp
494    - jd_x0 / hy    * uy_ym + jd_x0 / hy    * uy_yp
495);
src/fluid/compute_potential.c
532const double hx_xm = HXXF(i  );
533const double hx_xp = HXXF(i+1);
534const double jd_xm = JDXF(i  );
535const double jd_x0 = JDXC(i  );
536const double jd_xp = JDXF(i+1);
537const double ux_xm = UX(i  , j  , k  );
538const double ux_xp = UX(i+1, j  , k  );
539const double uy_ym = UY(i  , j  , k  );
540const double uy_yp = UY(i  , j+1, k  );
541const double uz_zm = UZ(i  , j  , k  );
542const double uz_zp = UZ(i  , j  , k+1);
543const double div = 1. / jd_x0 * (
544    - jd_xm / hx_xm * ux_xm + jd_xp / hx_xp * ux_xp
545    - jd_x0 / hy    * uy_ym + jd_x0 / hy    * uy_yp
546    - jd_x0 / hz    * uz_zm + jd_x0 / hz    * uz_zp
547);