I was using gradient, switched to laplace.

makefile

@@ -11,7 +11,7 @@ CXXFLAGS = -std=c++23 -I$(IMGUI_DIR) -I$(IMGUI_DIR)/backends -I$(STB_DIR)
 CXXFLAGS += `pkg-config --cflags glfw3`
 CFLAGS = $(CXXFLAGS)
 LDFLAGS := -L./imgui -limgui -lGL
-LINUX_GL_LIBS = -lGL
+LINUX_GL_LIBS = -lGL -ltbb
 PKG_FLAGS := $(LINUX_GL_LIBS) `pkg-config --static --cflags --libs glfw3`
 CXXFLAGS += $(PKG_FLAGS)
 

tests/g2temp.cpp

@@ -400,9 +400,9 @@ int main() {
     bool show_another_window = false;
     ImVec4 clear_color = ImVec4(0.45f, 0.55f, 0.60f, 1.00f);
     static float f = 30.0f;
-    static int i1 = 1024;
+    static int i1 = 256;
-    static int i2 = 1024;
+    static int i2 = 256;
-    static int i3 = 480;
+    static int i3 = 4800;
     static int i4 = 8;
     static int noisemod = 42;
     static float fs = 1.0;

util/grid/grid2.hpp

@@ -394,7 +394,7 @@ public:
     }
 
     size_t getOrCreatePositionVec(const Vec2& pos, float radius = 0.0f, bool create = true) {
-        TIME_FUNCTION;
+        //TIME_FUNCTION; //called too many times and average time is less than 0.0000001 so ignore it.
         if (radius == 0.0f) {
             Vec2 gridPos = spatialGrid.worldToGrid(pos);
             auto cellIt = spatialGrid.grid.find(gridPos);
@@ -423,7 +423,7 @@ public:
     }
 
     std::vector<size_t> getPositionVecRegion(const Vec2& pos, float radius = 1.0f) const {
-        TIME_FUNCTION;
+        //TIME_FUNCTION;
         float searchRadius = (radius == 0.0f) ? std::numeric_limits<float>::epsilon() : radius;
         
         // Get candidates from spatial grid
@@ -987,7 +987,7 @@ public:
     //         if (tempIT != tempMap.end()) {
     //             Temp oldtemp = tempIT->second;
     //             tempMap.erase(id);
-    //             float newtemp = Temp::calGrad(pos, getTemps(id));
+    //             float newtemp = Temp::calLapl(pos, getTemps(id));
     //             float newtempMult = (newtemp-oldtemp.temp) * timestep;
     //             oldtemp.temp = newtempMult;
     //             tempMap.emplace(id, oldtemp);
@@ -1017,7 +1017,7 @@ public:
                     neighborTemps.emplace(Positions.at(neighborId), tempMap.at(neighborId));
                 }
             }
-            tempObj.calGrad(pos, neighborTemps);
+            tempObj.calLapl(pos, neighborTemps);
             float newtemp = tempObj.temp;
             float tempdiff = (oldtemp - newtemp) * (deltaTime / 1000);
             // if (std::abs(newtemp) < EPSILON) {

util/simblocks/temp.hpp

@@ -78,52 +78,33 @@ public:
         return num / den;
     }
     
-    void calGrad(const Vec2& testPos, const std::unordered_map<Vec2, Temp>& others) {
+    void calLapl(const Vec2& testPos, const std::unordered_map<Vec2, Temp>& others) {
-        int meh;
+        //TIME_FUNCTION;
-        //std::cout << meh++ <<  std::endl;
+        float dt = 0.032;
-        std::vector<std::pair<Vec2, float>> nearbyPoints;
+        float sumWeights = 0.0f;
-        for (const auto& [point, ctemp] : others) {
+        float sumTempWeights = 0.0f;
-            if (point.distance(testPos) <= 25) {
+        float searchRadius = 25.0f;
-                nearbyPoints.emplace_back(point, ctemp.temp);
-            }
-        }
-        //std::cout << meh++ <<  std::endl;
-        float sumX = 0,sumY = 0,sumT = 0,sumX2 = 0,sumY2 = 0,sumXY = 0,sumXT = 0, sumYT = 0;
-        int n = nearbyPoints.size();
-        //std::cout << meh++ <<  std::endl;
-        for (const auto& [point, ctemp] : nearbyPoints) {
-            float x = point.x - testPos.x;
-            float y = point.y - testPos.y;
 
-            sumX += x;
+        for (const auto& [point, tempObj] : others) {
-            sumY += y;
+            float dist = testPos.distance(point);
-            sumT += ctemp;
+            
-            sumX2 += x * x;
+            if (dist < 0.001f || dist > searchRadius) continue;
-            sumY2 += y * y;
+
-            sumXY += x * y;
+            float weight = 1.0f / (dist * dist);
-            sumXT += x * ctemp;
+
-            sumYT += y * ctemp;
+            sumTempWeights += weight * tempObj.temp;
+            sumWeights += weight;
         }
 
-        //std::cout << meh++ <<  std::endl;
+        if (sumWeights < 1e-10f) return;
-        float det = sumX2 * sumY2 - sumXY * sumXY;
-    
-        Vec2 calpoint;
-        if (std::abs(det) < 1e-10) {
-            calpoint = Vec2(0, 0); // Singular matrix, cannot solve
-        }
-        else {
-            float a = (sumXT * sumY2 - sumYT * sumXY) / det;
-            float b = (sumX2 * sumYT - sumXY * sumXT) / det;
-            calpoint =  Vec2(a, b); // ∇T = (∂T/∂x, ∂T/∂y)
-        }
-        //std::cout << meh++ <<  std::endl;
-        Vec2 closest = findClosestPoint(testPos, others);
-        Vec2 displacement = testPos - closest;
 
-        float estimatedTemp = temp + (calpoint.x * displacement.x + calpoint.y * displacement.y);
+        float equilibriumTemp = sumTempWeights / sumWeights;
-        temp = estimatedTemp;
+
+        float rate = this->diffusivity * 0.01f;
         
+        float lerpFactor = 1.0f - std::exp(-rate * dt);
+
+        this->temp += (equilibriumTemp - this->temp) * lerpFactor;
     }
     
 };