From 9068920de3dd61b4c58fd407597fd53e1408dd88 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tam=C3=A1s=20B=C3=A1lint=20Misius?= Date: Sun, 8 Jan 2023 13:10:16 +0100 Subject: [PATCH] Preprocessor purge round 11: GRAVFFT By splitting gravity implementations into separate translation units. --- meson.build | 1 + src/Config.template.h | 2 +- src/config/font/meson.build | 1 - src/config/powder/meson.build | 1 - src/config/render/meson.build | 1 - src/gui/game/IntroText.h | 7 +- src/simulation/FftGravity.cpp | 202 +++++++++++++++++++++++++++++ src/simulation/Gravity.cpp | 223 +------------------------------- src/simulation/Gravity.h | 38 ++---- src/simulation/GravityPtr.h | 9 ++ src/simulation/PlainGravity.cpp | 53 ++++++++ src/simulation/Simulation.cpp | 3 +- src/simulation/Simulation.h | 3 +- src/simulation/meson.build | 7 + 14 files changed, 294 insertions(+), 257 deletions(-) create mode 100644 src/simulation/FftGravity.cpp create mode 100644 src/simulation/GravityPtr.h create mode 100644 src/simulation/PlainGravity.cpp diff --git a/meson.build b/meson.build index 825bb41b3..bbc3a7016 100644 --- a/meson.build +++ b/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 = [] diff --git a/src/Config.template.h b/src/Config.template.h index 69a7d9727..9a4df1616 100644 --- a/src/Config.template.h +++ b/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@; diff --git a/src/config/font/meson.build b/src/config/font/meson.build index 2a84c60f4..05098e993 100644 --- a/src/config/font/meson.build +++ b/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', diff --git a/src/config/powder/meson.build b/src/config/powder/meson.build index 86d580840..3a4f94559 100644 --- a/src/config/powder/meson.build +++ b/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', diff --git a/src/config/render/meson.build b/src/config/render/meson.build index 94adcd954..f0b6e64f6 100644 --- a/src/config/render/meson.build +++ b/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', diff --git a/src/gui/game/IntroText.h b/src/gui/game/IntroText.h index 22454d1c1..99b3758d6 100644 --- a/src/gui/game/IntroText.h +++ b/src/gui/game/IntroText.h @@ -45,9 +45,10 @@ inline ByteString IntroText() { sb << " LUACONSOLE"; } -#ifdef GRAVFFT - sb << " GRAVFFT"; -#endif + if constexpr (GRAVFFT) + { + sb << " GRAVFFT"; + } #ifdef REALISTIC sb << " REALISTIC"; #endif diff --git a/src/simulation/FftGravity.cpp b/src/simulation/FftGravity.cpp new file mode 100644 index 000000000..b2f8d1ec6 --- /dev/null +++ b/src/simulation/FftGravity.cpp @@ -0,0 +1,202 @@ +#include "Gravity.h" +#include "Misc.h" +#include +#include +#include + +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(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); + th_ptgravy = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); + th_ptgravxt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); + th_ptgravyt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); + th_gravmapbig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); + th_gravmapbigt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); + th_gravxbig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); + th_gravybig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); + th_gravxbigt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); + th_gravybigt = reinterpret_cast(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(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(ptr); +} diff --git a/src/simulation/Gravity.cpp b/src/simulation/Gravity.cpp index aa7b2ee6f..b83b9f6d1 100644 --- a/src/simulation/Gravity.cpp +++ b/src/simulation/Gravity.cpp @@ -9,9 +9,7 @@ #include "Simulation.h" #include "SimulationData.h" -#define GRAV_DIFF - -Gravity::Gravity() +Gravity::Gravity(CtorTag) { // 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(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_ptgravy = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_ptgravxt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); - th_ptgravyt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); - th_gravmapbig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_gravmapbigt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); - th_gravxbig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_gravybig = reinterpret_cast(fftwf_malloc(xblock2 * yblock2 * sizeof(float))); - th_gravxbigt = reinterpret_cast(fftwf_malloc(fft_tsize * sizeof(fftwf_complex))); - th_gravybigt = reinterpret_cast(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); - std::swap(gravx, th_gravx); - std::swap(gravp, th_gravp); -#endif + get_result(); } 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 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 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; diff --git a/src/simulation/Gravity.h b/src/simulation/Gravity.h index b5c3f330d..de97a5a8b 100644 --- a/src/simulation/Gravity.h +++ b/src/simulation/Gravity.h @@ -1,20 +1,17 @@ #pragma once #include "Config.h" +#include "GravityPtr.h" #include #include #include - -#ifdef GRAVFFT -#include -#endif +#include 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(); - - -#ifdef GRAVFFT - void grav_fft_init(); - void grav_fft_cleanup(); -#endif + + struct CtorTag // Please use Gravity::Create(). + { + }; 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(); - ~Gravity(); + static GravityPtr Create(); }; diff --git a/src/simulation/GravityPtr.h b/src/simulation/GravityPtr.h new file mode 100644 index 000000000..b58e096e3 --- /dev/null +++ b/src/simulation/GravityPtr.h @@ -0,0 +1,9 @@ +#pragma once +#include + +class Gravity; +struct GravityDeleter +{ + void operator ()(Gravity *ptr) const; +}; +using GravityPtr = std::unique_ptr; diff --git a/src/simulation/PlainGravity.cpp b/src/simulation/PlainGravity.cpp new file mode 100644 index 000000000..c728913d4 --- /dev/null +++ b/src/simulation/PlainGravity.cpp @@ -0,0 +1,53 @@ +#include "Gravity.h" +#include "Misc.h" +#include +#include + +// 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; +} diff --git a/src/simulation/Simulation.cpp b/src/simulation/Simulation.cpp index 027fcd6ae..b11addb4b 100644 --- a/src/simulation/Simulation.cpp +++ b/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 diff --git a/src/simulation/Simulation.h b/src/simulation/Simulation.h index b52626b56..23f878d0a 100644 --- a/src/simulation/Simulation.h +++ b/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 signs; diff --git a/src/simulation/meson.build b/src/simulation/meson.build index adb49630d..7f82344d8 100644 --- a/src/simulation/meson.build +++ b/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')