Preprocessor purge round 11: GRAVFFT

By splitting gravity implementations into separate translation units.

meson.build

@@ -350,6 +350,7 @@ conf_data.set('APPID', app_id)
 conf_data.set('APPDATA', get_option('app_data'))
 conf_data.set('APPVENDOR', get_option('app_vendor'))
 conf_data.set('LUACONSOLE', lua_variant != 'none' ? 'true' : 'false')
+conf_data.set('GRAVFFT', enable_gravfft ? 'true' : 'false')
 
 data_files = []
 

src/Config.template.h

@@ -3,7 +3,6 @@
 
 // Boolean macros (defined / not defined), would be great to get rid of them all.
 #mesondefine NOHTTP
-#mesondefine GRAVFFT
 #mesondefine RENDERER
 #mesondefine FONTEDITOR
 #mesondefine BETA
@@ -16,6 +15,7 @@
 #mesondefine MACOSX
 #mesondefine X86
 
+constexpr bool GRAVFFT                  = @GRAVFFT@;
 constexpr bool LUACONSOLE               = @LUACONSOLE@;
 constexpr bool ALLOW_FAKE_NEWER_VERSION = @ALLOW_FAKE_NEWER_VERSION@;
 constexpr bool USE_UPDATESERVER         = @USE_UPDATESERVER@;

src/config/font/meson.build

@@ -2,7 +2,6 @@ font_conf_data = conf_data
 font_conf_data.set('FONTEDITOR', true)
 font_conf_data.set('RENDERER', false)
 font_conf_data.set('NOHTTP', true)
-font_conf_data.set('GRAVFFT', false)
 configure_file(
 	input: config_template,
 	output: 'Config.h',

src/config/powder/meson.build

@@ -2,7 +2,6 @@ powder_conf_data = conf_data
 powder_conf_data.set('FONTEDITOR', false)
 powder_conf_data.set('RENDERER', false)
 powder_conf_data.set('NOHTTP', not enable_http)
-powder_conf_data.set('GRAVFFT', enable_gravfft)
 configure_file(
 	input: config_template,
 	output: 'Config.h',

src/config/render/meson.build

@@ -2,7 +2,6 @@ render_conf_data = conf_data
 render_conf_data.set('FONTEDITOR', false)
 render_conf_data.set('RENDERER', true)
 render_conf_data.set('NOHTTP', true)
-render_conf_data.set('GRAVFFT', false)
 configure_file(
 	input: config_template,
 	output: 'Config.h',

src/gui/game/IntroText.h

@@ -45,9 +45,10 @@ inline ByteString IntroText()
 	{
 		sb << " LUACONSOLE";
 	}
-#ifdef GRAVFFT
+	if constexpr (GRAVFFT)
-	sb << " GRAVFFT";
+	{
-#endif
+		sb << " GRAVFFT";
+	}
 #ifdef REALISTIC
 	sb << " REALISTIC";
 #endif

src/simulation/FftGravity.cpp

@@ -0,0 +1,202 @@
+#include "Gravity.h"
+#include "Misc.h"
+#include <cstring>
+#include <cmath>
+#include <fftw3.h>
+
+struct FftGravity : public Gravity
+{
+	bool grav_fft_status = false;
+	float *th_ptgravx = nullptr;
+	float *th_ptgravy = nullptr;
+	float *th_gravmapbig = nullptr;
+	float *th_gravxbig = nullptr;
+	float *th_gravybig = nullptr;
+
+	fftwf_complex *th_ptgravxt, *th_ptgravyt, *th_gravmapbigt, *th_gravxbigt, *th_gravybigt;
+	fftwf_plan plan_gravmap, plan_gravx_inverse, plan_gravy_inverse;
+
+	void grav_fft_init();
+	void grav_fft_cleanup();
+
+	FftGravity() : Gravity(CtorTag{})
+	{
+	}
+
+	~FftGravity();
+};
+
+FftGravity::~FftGravity()
+{
+	stop_grav_async();
+	grav_fft_cleanup();
+}
+
+void FftGravity::grav_fft_init()
+{
+	int xblock2 = XCELLS*2;
+	int yblock2 = YCELLS*2;
+	int fft_tsize = (xblock2/2+1)*yblock2;
+	float distance, scaleFactor;
+	fftwf_plan plan_ptgravx, plan_ptgravy;
+	if (grav_fft_status) return;
+
+	//use fftw malloc function to ensure arrays are aligned, to get better performance
+	th_ptgravx = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
+	th_ptgravy = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
+	th_ptgravxt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
+	th_ptgravyt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
+	th_gravmapbig = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
+	th_gravmapbigt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
+	th_gravxbig = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
+	th_gravybig = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
+	th_gravxbigt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
+	th_gravybigt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
+
+	//select best algorithm, could use FFTW_PATIENT or FFTW_EXHAUSTIVE but that increases the time taken to plan, and I don't see much increase in execution speed
+	plan_ptgravx = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_ptgravx, th_ptgravxt, FFTW_MEASURE);
+	plan_ptgravy = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_ptgravy, th_ptgravyt, FFTW_MEASURE);
+	plan_gravmap = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_gravmapbig, th_gravmapbigt, FFTW_MEASURE);
+	plan_gravx_inverse = fftwf_plan_dft_c2r_2d(yblock2, xblock2, th_gravxbigt, th_gravxbig, FFTW_MEASURE);
+	plan_gravy_inverse = fftwf_plan_dft_c2r_2d(yblock2, xblock2, th_gravybigt, th_gravybig, FFTW_MEASURE);
+
+	//NCELL*4 is size of data array, scaling needed because FFTW calculates an unnormalized DFT
+	scaleFactor = -float(M_GRAV)/(NCELL*4);
+	//calculate velocity map caused by a point mass
+	for (int y = 0; y < yblock2; y++)
+	{
+		for (int x = 0; x < xblock2; x++)
+		{
+			if (x == XCELLS && y == YCELLS)
+				continue;
+			distance = hypotf(float(x-XCELLS), float(y-YCELLS));
+			th_ptgravx[y * xblock2 + x] = scaleFactor * (x - XCELLS) / powf(distance, 3);
+			th_ptgravy[y * xblock2 + x] = scaleFactor * (y - YCELLS) / powf(distance, 3);
+		}
+	}
+	th_ptgravx[yblock2 * xblock2 / 2 + xblock2 / 2] = 0.0f;
+	th_ptgravy[yblock2 * xblock2 / 2 + xblock2 / 2] = 0.0f;
+
+	//transform point mass velocity maps
+	fftwf_execute(plan_ptgravx);
+	fftwf_execute(plan_ptgravy);
+	fftwf_destroy_plan(plan_ptgravx);
+	fftwf_destroy_plan(plan_ptgravy);
+	fftwf_free(th_ptgravx);
+	fftwf_free(th_ptgravy);
+
+	//clear padded gravmap
+	memset(th_gravmapbig, 0, xblock2 * yblock2 * sizeof(float));
+
+	grav_fft_status = true;
+}
+
+void FftGravity::grav_fft_cleanup()
+{
+	if (!grav_fft_status) return;
+	fftwf_free(th_ptgravxt);
+	fftwf_free(th_ptgravyt);
+	fftwf_free(th_gravmapbig);
+	fftwf_free(th_gravmapbigt);
+	fftwf_free(th_gravxbig);
+	fftwf_free(th_gravybig);
+	fftwf_free(th_gravxbigt);
+	fftwf_free(th_gravybigt);
+	fftwf_destroy_plan(plan_gravmap);
+	fftwf_destroy_plan(plan_gravx_inverse);
+	fftwf_destroy_plan(plan_gravy_inverse);
+	grav_fft_status = false;
+}
+
+void Gravity::get_result()
+{
+	std::swap(gravy, th_gravy);
+	std::swap(gravx, th_gravx);
+	std::swap(gravp, th_gravp);
+}
+
+void Gravity::update_grav()
+{
+	auto *fftGravity = static_cast<FftGravity *>(this);
+	if (!fftGravity->grav_fft_status)
+		fftGravity->grav_fft_init();
+
+	auto *th_gravmapbig = fftGravity->th_gravmapbig;
+	auto *th_gravxbig = fftGravity->th_gravxbig;
+	auto *th_gravybig = fftGravity->th_gravybig;
+	auto *th_ptgravxt = fftGravity->th_ptgravxt;
+	auto *th_ptgravyt = fftGravity->th_ptgravyt;
+	auto *th_gravmapbigt = fftGravity->th_gravmapbigt;
+	auto *th_gravxbigt = fftGravity->th_gravxbigt;
+	auto *th_gravybigt = fftGravity->th_gravybigt;
+	auto &plan_gravmap = fftGravity->plan_gravmap;
+	auto &plan_gravx_inverse = fftGravity->plan_gravx_inverse;
+	auto &plan_gravy_inverse = fftGravity->plan_gravy_inverse;
+
+	int xblock2 = XCELLS*2, yblock2 = YCELLS*2;
+	int fft_tsize = (xblock2/2+1)*yblock2;
+	float mr, mc, pr, pc, gr, gc;
+	if (memcmp(th_ogravmap, th_gravmap, sizeof(float)*NCELL) != 0)
+	{
+		th_gravchanged = 1;
+
+		membwand(th_gravmap, gravmask, NCELL*sizeof(float), NCELL*sizeof(unsigned));
+		//copy gravmap into padded gravmap array
+		for (int y = 0; y < YCELLS; y++)
+		{
+			for (int x = 0; x < XCELLS; x++)
+			{
+				th_gravmapbig[(y+YCELLS)*xblock2+XCELLS+x] = th_gravmap[y*XCELLS+x];
+			}
+		}
+		//transform gravmap
+		fftwf_execute(plan_gravmap);
+		//do convolution (multiply the complex numbers)
+		for (int i = 0; i < fft_tsize; i++)
+		{
+			mr = th_gravmapbigt[i][0];
+			mc = th_gravmapbigt[i][1];
+			pr = th_ptgravxt[i][0];
+			pc = th_ptgravxt[i][1];
+			gr = mr*pr-mc*pc;
+			gc = mr*pc+mc*pr;
+			th_gravxbigt[i][0] = gr;
+			th_gravxbigt[i][1] = gc;
+			pr = th_ptgravyt[i][0];
+			pc = th_ptgravyt[i][1];
+			gr = mr*pr-mc*pc;
+			gc = mr*pc+mc*pr;
+			th_gravybigt[i][0] = gr;
+			th_gravybigt[i][1] = gc;
+		}
+		//inverse transform, and copy from padded arrays into normal velocity maps
+		fftwf_execute(plan_gravx_inverse);
+		fftwf_execute(plan_gravy_inverse);
+		for (int y = 0; y < YCELLS; y++)
+		{
+			for (int x = 0; x < XCELLS; x++)
+			{
+				th_gravx[y*XCELLS+x] = th_gravxbig[y*xblock2+x];
+				th_gravy[y*XCELLS+x] = th_gravybig[y*xblock2+x];
+				th_gravp[y*XCELLS+x] = hypotf(th_gravxbig[y*xblock2+x], th_gravybig[y*xblock2+x]);
+			}
+		}
+	}
+	else
+	{
+		th_gravchanged = 0;
+	}
+
+	// Copy th_ogravmap into th_gravmap (doesn't matter what th_ogravmap is afterwards)
+	std::swap(th_gravmap, th_ogravmap);
+}
+
+GravityPtr Gravity::Create()
+{
+	return GravityPtr(new FftGravity());
+}
+
+void GravityDeleter::operator ()(Gravity *ptr) const
+{
+	delete static_cast<FftGravity *>(ptr);
+}

src/simulation/Gravity.cpp

@@ -9,9 +9,7 @@
 #include "Simulation.h"
 #include "SimulationData.h"
 
-#define GRAV_DIFF
+Gravity::Gravity(CtorTag)
-
-Gravity::Gravity()
 {
 	// Allocate full size Gravmaps
 	unsigned int size = NCELL;
@@ -30,9 +28,6 @@ Gravity::Gravity()
 Gravity::~Gravity()
 {
 	stop_grav_async();
-#ifdef GRAVFFT
-	grav_fft_cleanup();
-#endif
 
 	delete[] th_ogravmap;
 	delete[] th_gravmap;
@@ -58,84 +53,6 @@ void Gravity::Clear()
 	ignoreNextResult = true;
 }
 
-#ifdef GRAVFFT
-void Gravity::grav_fft_init()
-{
-	int xblock2 = XCELLS*2;
-	int yblock2 = YCELLS*2;
-	int fft_tsize = (xblock2/2+1)*yblock2;
-	float distance, scaleFactor;
-	fftwf_plan plan_ptgravx, plan_ptgravy;
-	if (grav_fft_status) return;
-
-	//use fftw malloc function to ensure arrays are aligned, to get better performance
-	th_ptgravx = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
-	th_ptgravy = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
-	th_ptgravxt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
-	th_ptgravyt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
-	th_gravmapbig = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
-	th_gravmapbigt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
-	th_gravxbig = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
-	th_gravybig = reinterpret_cast<float*>(fftwf_malloc(xblock2 * yblock2 * sizeof(float)));
-	th_gravxbigt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
-	th_gravybigt = reinterpret_cast<fftwf_complex*>(fftwf_malloc(fft_tsize * sizeof(fftwf_complex)));
-
-	//select best algorithm, could use FFTW_PATIENT or FFTW_EXHAUSTIVE but that increases the time taken to plan, and I don't see much increase in execution speed
-	plan_ptgravx = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_ptgravx, th_ptgravxt, FFTW_MEASURE);
-	plan_ptgravy = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_ptgravy, th_ptgravyt, FFTW_MEASURE);
-	plan_gravmap = fftwf_plan_dft_r2c_2d(yblock2, xblock2, th_gravmapbig, th_gravmapbigt, FFTW_MEASURE);
-	plan_gravx_inverse = fftwf_plan_dft_c2r_2d(yblock2, xblock2, th_gravxbigt, th_gravxbig, FFTW_MEASURE);
-	plan_gravy_inverse = fftwf_plan_dft_c2r_2d(yblock2, xblock2, th_gravybigt, th_gravybig, FFTW_MEASURE);
-
-	//NCELL*4 is size of data array, scaling needed because FFTW calculates an unnormalized DFT
-	scaleFactor = -float(M_GRAV)/(NCELL*4);
-	//calculate velocity map caused by a point mass
-	for (int y = 0; y < yblock2; y++)
-	{
-		for (int x = 0; x < xblock2; x++)
-		{
-			if (x == XCELLS && y == YCELLS)
-				continue;
-			distance = sqrtf(pow(x-XCELLS, 2.0f) + pow(y-YCELLS, 2.0f));
-			th_ptgravx[y * xblock2 + x] = scaleFactor * (x - XCELLS) / pow(distance, 3);
-			th_ptgravy[y * xblock2 + x] = scaleFactor * (y - YCELLS) / pow(distance, 3);
-		}
-	}
-	th_ptgravx[yblock2 * xblock2 / 2 + xblock2 / 2] = 0.0f;
-	th_ptgravy[yblock2 * xblock2 / 2 + xblock2 / 2] = 0.0f;
-
-	//transform point mass velocity maps
-	fftwf_execute(plan_ptgravx);
-	fftwf_execute(plan_ptgravy);
-	fftwf_destroy_plan(plan_ptgravx);
-	fftwf_destroy_plan(plan_ptgravy);
-	fftwf_free(th_ptgravx);
-	fftwf_free(th_ptgravy);
-
-	//clear padded gravmap
-	memset(th_gravmapbig, 0, xblock2 * yblock2 * sizeof(float));
-
-	grav_fft_status = true;
-}
-
-void Gravity::grav_fft_cleanup()
-{
-	if (!grav_fft_status) return;
-	fftwf_free(th_ptgravxt);
-	fftwf_free(th_ptgravyt);
-	fftwf_free(th_gravmapbig);
-	fftwf_free(th_gravmapbigt);
-	fftwf_free(th_gravxbig);
-	fftwf_free(th_gravybig);
-	fftwf_free(th_gravxbigt);
-	fftwf_free(th_gravybigt);
-	fftwf_destroy_plan(plan_gravmap);
-	fftwf_destroy_plan(plan_gravx_inverse);
-	fftwf_destroy_plan(plan_gravy_inverse);
-	grav_fft_status = false;
-}
-#endif
-
 void Gravity::gravity_update_async()
 {
 	int result;
@@ -153,16 +70,8 @@ void Gravity::gravity_update_async()
 			{
 				if (th_gravchanged && !ignoreNextResult)
 				{
-#if !defined(GRAVFFT) && defined(GRAV_DIFF)
-					memcpy(gravy, th_gravy, NCELL*sizeof(float));
-					memcpy(gravx, th_gravx, NCELL*sizeof(float));
-					memcpy(gravp, th_gravp, NCELL*sizeof(float));
-#else
 					// Copy thread gravity maps into this one
-					std::swap(gravy, th_gravy);
+					get_result();
-					std::swap(gravx, th_gravx);
-					std::swap(gravp, th_gravp);
-#endif
 				}
 				ignoreNextResult = false;
 
@@ -194,11 +103,6 @@ void Gravity::update_grav_async()
 	std::fill(&th_gravx[0], &th_gravx[0] + NCELL, 0.0f);
 	std::fill(&th_gravp[0], &th_gravp[0] + NCELL, 0.0f);
 
-#ifdef GRAVFFT
-	if (!grav_fft_status)
-		grav_fft_init();
-#endif
-
 	std::unique_lock<std::mutex> l(gravmutex);
 	while (!thread_done)
 	{
@@ -255,129 +159,6 @@ void Gravity::stop_grav_async()
 	std::fill(&gravmap[0], &gravmap[0] + NCELL, 0.0f);
 }
 
-#ifdef GRAVFFT
-void Gravity::update_grav()
-{
-	int xblock2 = XCELLS*2, yblock2 = YCELLS*2;
-	int fft_tsize = (xblock2/2+1)*yblock2;
-	float mr, mc, pr, pc, gr, gc;
-	if (memcmp(th_ogravmap, th_gravmap, sizeof(float)*NCELL) != 0)
-	{
-		th_gravchanged = 1;
-
-		membwand(th_gravmap, gravmask, NCELL*sizeof(float), NCELL*sizeof(unsigned));
-		//copy gravmap into padded gravmap array
-		for (int y = 0; y < YCELLS; y++)
-		{
-			for (int x = 0; x < XCELLS; x++)
-			{
-				th_gravmapbig[(y+YCELLS)*xblock2+XCELLS+x] = th_gravmap[y*XCELLS+x];
-			}
-		}
-		//transform gravmap
-		fftwf_execute(plan_gravmap);
-		//do convolution (multiply the complex numbers)
-		for (int i = 0; i < fft_tsize; i++)
-		{
-			mr = th_gravmapbigt[i][0];
-			mc = th_gravmapbigt[i][1];
-			pr = th_ptgravxt[i][0];
-			pc = th_ptgravxt[i][1];
-			gr = mr*pr-mc*pc;
-			gc = mr*pc+mc*pr;
-			th_gravxbigt[i][0] = gr;
-			th_gravxbigt[i][1] = gc;
-			pr = th_ptgravyt[i][0];
-			pc = th_ptgravyt[i][1];
-			gr = mr*pr-mc*pc;
-			gc = mr*pc+mc*pr;
-			th_gravybigt[i][0] = gr;
-			th_gravybigt[i][1] = gc;
-		}
-		//inverse transform, and copy from padded arrays into normal velocity maps
-		fftwf_execute(plan_gravx_inverse);
-		fftwf_execute(plan_gravy_inverse);
-		for (int y = 0; y < YCELLS; y++)
-		{
-			for (int x = 0; x < XCELLS; x++)
-			{
-				th_gravx[y*XCELLS+x] = th_gravxbig[y*xblock2+x];
-				th_gravy[y*XCELLS+x] = th_gravybig[y*xblock2+x];
-				th_gravp[y*XCELLS+x] = sqrtf(pow(th_gravxbig[y*xblock2+x],2)+pow(th_gravybig[y*xblock2+x],2));
-			}
-		}
-	}
-	else
-	{
-		th_gravchanged = 0;
-	}
-
-	// Copy th_ogravmap into th_gravmap (doesn't matter what th_ogravmap is afterwards)
-	std::swap(th_gravmap, th_ogravmap);
-}
-
-#else
-// gravity without fast Fourier transforms
-
-void Gravity::update_grav(void)
-{
-	int x, y, i, j;
-	float val, distance;
-	th_gravchanged = 0;
-#ifndef GRAV_DIFF
-	//Find any changed cells
-	for (i=0; i<YCELLS; i++)
-	{
-		if(changed)
-			break;
-		for (j=0; j<XCELLS; j++)
-		{
-			if(th_ogravmap[i*XCELLS+j]!=th_gravmap[i*XCELLS+j]){
-				changed = 1;
-				break;
-			}
-		}
-	}
-	if(!changed)
-		goto fin;
-	memset(th_gravy, 0, NCELL*sizeof(float));
-	memset(th_gravx, 0, NCELL*sizeof(float));
-#endif
-	th_gravchanged = 1;
-	membwand(th_gravmap, gravmask, NCELL*sizeof(float), NCELL*sizeof(unsigned));
-	for (i = 0; i < YCELLS; i++) {
-		for (j = 0; j < XCELLS; j++) {
-#ifdef GRAV_DIFF
-			if (th_ogravmap[i*XCELLS+j] != th_gravmap[i*XCELLS+j])
-			{
-#else
-			if (th_gravmap[i*XCELLS+j] > 0.0001f || th_gravmap[i*XCELLS+j]<-0.0001f) //Only calculate with populated or changed cells.
-			{
-#endif
-				for (y = 0; y < YCELLS; y++) {
-					for (x = 0; x < XCELLS; x++) {
-						if (x == j && y == i)//Ensure it doesn't calculate with itself
-							continue;
-						distance = sqrt(pow(j - x, 2.0f) + pow(i - y, 2.0f));
-#ifdef GRAV_DIFF
-						val = th_gravmap[i*XCELLS+j] - th_ogravmap[i*XCELLS+j];
-#else
-						val = th_gravmap[i*XCELLS+j];
-#endif
-						th_gravx[y*XCELLS+x] += M_GRAV * val * (j - x) / pow(distance, 3.0f);
-						th_gravy[y*XCELLS+x] += M_GRAV * val * (i - y) / pow(distance, 3.0f);
-						th_gravp[y*XCELLS+x] += M_GRAV * val / pow(distance, 2.0f);
-					}
-				}
-			}
-		}
-	}
-	memcpy(th_ogravmap, th_gravmap, NCELL*sizeof(float));
-}
-#endif
-
-
-
 bool Gravity::grav_mask_r(int x, int y, char checkmap[YCELLS][XCELLS], char shape[YCELLS][XCELLS])
 {
 	int x1, x2;

src/simulation/Gravity.h

@@ -1,20 +1,17 @@
 #pragma once
 #include "Config.h"
+#include "GravityPtr.h"
 
 #include <thread>
 #include <mutex>
 #include <condition_variable>
-
+#include <memory>
-#ifdef GRAVFFT
-#include <fftw3.h>
-#endif
 
 class Simulation;
 
 class Gravity
 {
-private:
+protected:
-
 	bool enabled = false;
 
 	// Maps to be processed by the gravity thread
@@ -33,18 +30,6 @@ private:
 	int gravthread_done = 0;
 	bool ignoreNextResult = false;
 
-#ifdef GRAVFFT
-	bool grav_fft_status = false;
-	float *th_ptgravx = nullptr;
-	float *th_ptgravy = nullptr;
-	float *th_gravmapbig = nullptr;
-	float *th_gravxbig = nullptr;
-	float *th_gravybig = nullptr;
-
-	fftwf_complex *th_ptgravxt, *th_ptgravyt, *th_gravmapbigt, *th_gravxbigt, *th_gravybigt;
-	fftwf_plan plan_gravmap, plan_gravx_inverse, plan_gravy_inverse;
-#endif
-
 	struct mask_el {
 		char *shape;
 		char shapeout;
@@ -56,15 +41,17 @@ private:
 	void mask_free(mask_el *c_mask_el);
 
 	void update_grav();
+	void get_result();
 	void update_grav_async();
-
+	
-
+	struct CtorTag // Please use Gravity::Create().
-#ifdef GRAVFFT
+	{
-	void grav_fft_init();
+	};
-	void grav_fft_cleanup();
-#endif
 
 public:
+	Gravity(CtorTag);
+	~Gravity();
+
 	//Maps to be used by the main thread
 	float *gravmap = nullptr;
 	float *gravp = nullptr;
@@ -84,6 +71,5 @@ public:
 	void stop_grav_async();
 	void gravity_mask();
 
-	Gravity();
+	static GravityPtr Create();
-	~Gravity();
 };

src/simulation/GravityPtr.h

@@ -0,0 +1,9 @@
+#pragma once
+#include <memory>
+
+class Gravity;
+struct GravityDeleter
+{
+	void operator ()(Gravity *ptr) const;
+};
+using GravityPtr = std::unique_ptr<Gravity, GravityDeleter>;

src/simulation/PlainGravity.cpp

@@ -0,0 +1,53 @@
+#include "Gravity.h"
+#include "Misc.h"
+#include <cmath>
+#include <cstring>
+
+// gravity without fast Fourier transforms
+
+void Gravity::get_result()
+{
+	memcpy(gravy, th_gravy, NCELL*sizeof(float));
+	memcpy(gravx, th_gravx, NCELL*sizeof(float));
+	memcpy(gravp, th_gravp, NCELL*sizeof(float));
+}
+
+void Gravity::update_grav(void)
+{
+	th_gravchanged = 1;
+	membwand(th_gravmap, gravmask, NCELL*sizeof(float), NCELL*sizeof(unsigned));
+	for (int i = 0; i < YCELLS; i++)
+	{
+		for (int j = 0; j < XCELLS; j++)
+		{
+			if (th_ogravmap[i*XCELLS+j] != th_gravmap[i*XCELLS+j])
+			{
+				for (int y = 0; y < YCELLS; y++)
+				{
+					for (int x = 0; x < XCELLS; x++)
+					{
+						if (x == j && y == i)//Ensure it doesn't calculate with itself
+							continue;
+						auto distance = hypotf(j - x, i - y);
+						float val;
+						val = th_gravmap[i*XCELLS+j] - th_ogravmap[i*XCELLS+j];
+						th_gravx[y*XCELLS+x] += M_GRAV * val * (j - x) / powf(distance, 3.0f);
+						th_gravy[y*XCELLS+x] += M_GRAV * val * (i - y) / powf(distance, 3.0f);
+						th_gravp[y*XCELLS+x] += M_GRAV * val / powf(distance, 2.0f);
+					}
+				}
+			}
+		}
+	}
+	memcpy(th_ogravmap, th_gravmap, NCELL*sizeof(float));
+}
+
+GravityPtr Gravity::Create()
+{
+	return GravityPtr(new Gravity(CtorTag{}));
+}
+
+void GravityDeleter::operator ()(Gravity *ptr) const
+{
+	delete ptr;
+}

src/simulation/Simulation.cpp

@@ -5198,7 +5198,6 @@ void Simulation::AfterSim()
 
 Simulation::~Simulation()
 {
-	delete grav;
 	delete air;
 }
 
@@ -5240,7 +5239,7 @@ Simulation::Simulation():
 	elementRecount = true;
 
 	//Create and attach gravity simulation
-	grav = new Gravity();
+	grav = Gravity::Create();
 	//Give air sim references to our data
 	grav->bmap = bmap;
 	//Gravity sim gives us maps to use

src/simulation/Simulation.h

@@ -15,6 +15,7 @@
 #include "BuiltinGOL.h"
 #include "MenuSection.h"
 #include "CoordStack.h"
+#include "GravityPtr.h"
 
 #include "Element.h"
 
@@ -37,7 +38,7 @@ class Simulation
 {
 public:
 
-	Gravity * grav;
+	GravityPtr grav;
 	Air * air;
 
 	std::vector<sign> signs;

src/simulation/meson.build

@@ -19,3 +19,10 @@ subdir('simtools')
 
 powder_files += simulation_files
 render_files += simulation_files
+
+if enable_gravfft
+	powder_files += files('FftGravity.cpp')
+else
+	powder_files += files('PlainGravity.cpp')
+endif
+render_files += files('PlainGravity.cpp')