density normalization work for pressurecell

docs/src/architecture/fluidsim/fluidsimnotes.md

@@ -33,4 +33,11 @@ edge boundaries
 extension of the "warm up" idea
  Run minified solver that just performs velocity step (no density)
  when advecting, check if density is at the advection position
- if there is density, pull it along as well and fully activate the empty chunk
+ if there is density, pull it along as well and fully activate the empty chunk
+
+==current frame==
+newest u,v,w depend on density (for grav)
+neighboring divergences depend on newest u,v,w
+pressure field depends on neighboring divergences
+final u,v,w depends on pressure field
+final density depends on final u,v,w

src/main/c/includes/fluid/queue/chunk.h

@@ -57,6 +57,16 @@ static char INTEREST_MODIFIER_DIMS[CHUNK_MAX_INTEREST_LEVEL+1] = {
  */
 #define INTEREST(tree,level,x,y,z) tree[level][x/INTEREST_MODIFIER_DIMS[level]*INTEREST_MODIFIER_DIMS[level]*INTEREST_MODIFIER_DIMS[level]+(y/INTEREST_MODIFIER_DIMS[level])*INTEREST_MODIFIER_DIMS[level]+(z/INTEREST_MODIFIER_DIMS[level])]
 
+/**
+ * The data for this chunk that is specific to the pressure cell method
+ */
+typedef struct {
+    /**
+     * The sum of density
+     */
+    double densitySum;
+} PressureCellData;
+
 /**
  * A chunk
 */
@@ -141,11 +151,6 @@ typedef struct {
      */
     int simLOD;
 
-    /**
-     * Octree that stores whether a location is of interest or not
-     */
-    char ** interestTree;
-
     /**
      * The convergence of this chunk
      */
@@ -155,6 +160,11 @@ typedef struct {
      * The number of iterations this chunk took to project
      */
     int projectionIterations;
+
+    /**
+     * The data for pressure cell work in particular
+     */
+    PressureCellData pressureCellData;
     
 } Chunk;
 

src/main/c/src/fluid/queue/chunk.c

@@ -22,13 +22,13 @@ LIBRARY_API Chunk * chunk_create(){
      * 8x8x8
      * 16x16x16
      */
-    rVal->interestTree = (char **)calloc(1,sizeof(char *) * 6);
+    // rVal->interestTree = (char **)calloc(1,sizeof(char *) * 6);
-    //allocating extra data for ghost cells even though we will not evaluate the ghost cells
+    // //allocating extra data for ghost cells even though we will not evaluate the ghost cells
-    rVal->interestTree[0] = (char *)calloc(3*3*3,sizeof(char));
+    // rVal->interestTree[0] = (char *)calloc(3*3*3,sizeof(char));
-    rVal->interestTree[1] = (char *)calloc(4*4*4,sizeof(char));
+    // rVal->interestTree[1] = (char *)calloc(4*4*4,sizeof(char));
-    rVal->interestTree[2] = (char *)calloc(6*6*6,sizeof(char));
+    // rVal->interestTree[2] = (char *)calloc(6*6*6,sizeof(char));
-    rVal->interestTree[3] = (char *)calloc(10*10*10,sizeof(char));
+    // rVal->interestTree[3] = (char *)calloc(10*10*10,sizeof(char));
-    rVal->interestTree[4] = (char *)calloc(18*18*18,sizeof(char));
+    // rVal->interestTree[4] = (char *)calloc(18*18*18,sizeof(char));
 
     return rVal;
 }

src/main/c/src/fluid/sim/pressurecell/bounds.c

@@ -96,52 +96,51 @@ LIBRARY_API void pressurecell_update_bounds(Environment * environment, Chunk * c
 /**
  * Updates the interest tree for this chunk
  */
-LIBRARY_API void pressurecell_update_interest(Environment * environment, Chunk * chunk){
+// LIBRARY_API void pressurecell_update_interest(Environment * environment, Chunk * chunk){
-    int level, x, y, z;
+//     int level, x, y, z;
-    int dim;
+//     int dim;
-    char ** interestTree = chunk->interestTree;
+//     float * densityArr = chunk->d[CENTER_LOC];
-    float * densityArr = chunk->d[CENTER_LOC];
+//     float * densitSrcArr = chunk->d0[CENTER_LOC];
-    float * densitSrcArr = chunk->d0[CENTER_LOC];
+//     // for(level = 0; level < CHUNK_MAX_INTEREST_LEVEL; level++){
-    // for(level = 0; level < CHUNK_MAX_INTEREST_LEVEL; level++){
+//     //     dim = pow(2,level);
-    //     dim = pow(2,level);
+//     //     for(x = 0; x < dim; x++){
-    //     for(x = 0; x < dim; x++){
+//     //         for(y = 0; y < dim; y++){
-    //         for(y = 0; y < dim; y++){
+//     //             for(z = 0; z < dim; z++){
-    //             for(z = 0; z < dim; z++){
+//     //                 if(
-    //                 if(
+//     //                     densityArr[IX(x,y,z)] > 0 ||
-    //                     densityArr[IX(x,y,z)] > 0 ||
+//     //                     densityArr[IX(x+1,y,z)] > 0 ||
-    //                     densityArr[IX(x+1,y,z)] > 0 ||
+//     //                     densityArr[IX(x-1,y,z)] > 0 ||
-    //                     densityArr[IX(x-1,y,z)] > 0 ||
+//     //                     densityArr[IX(x,y+1,z)] > 0 ||
-    //                     densityArr[IX(x,y+1,z)] > 0 ||
+//     //                     densityArr[IX(x,y-1,z)] > 0 ||
-    //                     densityArr[IX(x,y-1,z)] > 0 ||
+//     //                     densityArr[IX(x,y,z+1)] > 0 ||
-    //                     densityArr[IX(x,y,z+1)] > 0 ||
+//     //                     densityArr[IX(x,y,z-1)] > 0
-    //                     densityArr[IX(x,y,z-1)] > 0
+//     //                     ){
-    //                     ){
 
-    //                 }
+//     //                 }
-    //             }
+//     //             }
-    //         }
+//     //         }
-    //     }
+//     //     }
-    // }
+//     // }
-    for(x = 1; x < DIM-1; x++){
+//     for(x = 1; x < DIM-1; x++){
-        for(y = 1; y < DIM-1; y++){
+//         for(y = 1; y < DIM-1; y++){
-            for(z = 1; z < DIM-1; z++){
+//             for(z = 1; z < DIM-1; z++){
-                if(
+//                 if(
-                    densityArr[IX(x,y,z)] > MIN_FLUID_VALUE ||
+//                     densityArr[IX(x,y,z)] > MIN_FLUID_VALUE ||
-                    densityArr[IX(x+1,y,z)] > MIN_FLUID_VALUE ||
+//                     densityArr[IX(x+1,y,z)] > MIN_FLUID_VALUE ||
-                    densityArr[IX(x-1,y,z)] > MIN_FLUID_VALUE ||
+//                     densityArr[IX(x-1,y,z)] > MIN_FLUID_VALUE ||
-                    densityArr[IX(x,y+1,z)] > MIN_FLUID_VALUE ||
+//                     densityArr[IX(x,y+1,z)] > MIN_FLUID_VALUE ||
-                    densityArr[IX(x,y-1,z)] > MIN_FLUID_VALUE ||
+//                     densityArr[IX(x,y-1,z)] > MIN_FLUID_VALUE ||
-                    densityArr[IX(x,y,z+1)] > MIN_FLUID_VALUE ||
+//                     densityArr[IX(x,y,z+1)] > MIN_FLUID_VALUE ||
-                    densityArr[IX(x,y,z-1)] > MIN_FLUID_VALUE ||
+//                     densityArr[IX(x,y,z-1)] > MIN_FLUID_VALUE ||
-                    densitSrcArr[IX(x,y,z)] > MIN_FLUID_VALUE
+//                     densitSrcArr[IX(x,y,z)] > MIN_FLUID_VALUE
-                ){
+//                 ){
-                    INTEREST(interestTree,0,x,y,z) = 1;
+//                     INTEREST(interestTree,0,x,y,z) = 1;
-                    INTEREST(interestTree,1,x,y,z) = 1;
+//                     INTEREST(interestTree,1,x,y,z) = 1;
-                    INTEREST(interestTree,2,x,y,z) = 1;
+//                     INTEREST(interestTree,2,x,y,z) = 1;
-                    INTEREST(interestTree,3,x,y,z) = 1;
+//                     INTEREST(interestTree,3,x,y,z) = 1;
-                    INTEREST(interestTree,4,x,y,z) = 1;
+//                     INTEREST(interestTree,4,x,y,z) = 1;
-                }
+//                 }
-            }
+//             }
-        }
+//         }
-    }
+//     }
-}
+// }

src/main/c/src/fluid/sim/pressurecell/normalization.c

@@ -1,6 +1,8 @@
 
 
 #include "fluid/sim/pressurecell/normalization.h"
+#include "fluid/queue/chunkmask.h"
+#include "fluid/queue/chunk.h"
 
 
 
@@ -9,21 +11,49 @@
  * Calculates the expected density and pressure
  */
 LIBRARY_API void fluid_pressurecell_calculate_expected_intake(Environment * env, Chunk * chunk){
-
+    int x, y, z;
+    double sum;
+    for(x = 1; x < DIM-1; x++){
+        for(y = 1; y < DIM-1; y++){
+            for(z = 1; z < DIM-1; z++){
+                sum = sum + chunk->d[CENTER_LOC][IX(x,y,z)] + chunk->d0[CENTER_LOC][IX(x,y,z)];
+            }
+        }
+    }
+    chunk->pressureCellData.densitySum = sum;
 }
 
 /**
  * Calculates the ratio to normalize the chunk by
  */
 LIBRARY_API void fluid_pressurecell_calculate_normalization_ratio(Environment * env, Chunk * chunk){
-
+    int x, y, z;
+    double sum;
+    for(x = 1; x < DIM-1; x++){
+        for(y = 1; y < DIM-1; y++){
+            for(z = 1; z < DIM-1; z++){
+                sum = sum + chunk->d[CENTER_LOC][IX(x,y,z)];
+            }
+        }
+    }
+    double expected = chunk->pressureCellData.densitySum;
+    double normalizationRatio = expected / sum;
+    chunk->pressureCellData.densitySum = normalizationRatio;
 }
 
 /**
  * Normalizes the chunk
  */
 LIBRARY_API void fluid_pressurecell_normalize_chunk(Environment * env, Chunk * chunk){
-
+    int x, y, z;
+    double ratio = chunk->pressureCellData.densitySum;
+    for(x = 1; x < DIM-1; x++){
+        for(y = 1; y < DIM-1; y++){
+            for(z = 1; z < DIM-1; z++){
+                chunk->d[CENTER_LOC][IX(x,y,z)] = chunk->d[CENTER_LOC][IX(x,y,z)] * ratio;
+            }
+        }
+    }
 }
 
 

src/main/c/src/fluid/sim/pressurecell/pressurecell.c

@@ -40,6 +40,7 @@ LIBRARY_API void fluid_pressurecell_simulate(
 
     for(int i = 0; i < numChunks; i++){
         Chunk * currentChunk = chunks[i];
+        fluid_pressurecell_calculate_expected_intake(environment,currentChunk);
         pressurecell_update_bounds(environment,currentChunk);
         // pressurecell_update_interest(environment,currentChunk);
     }
@@ -120,6 +121,11 @@ LIBRARY_API void fluid_pressurecell_simulate(
     //
     for(int i = 0; i < numChunks; i++){
         Chunk * currentChunk = chunks[i];
+        fluid_pressurecell_calculate_normalization_ratio(environment,currentChunk);
+    }
+    for(int i = 0; i < numChunks; i++){
+        Chunk * currentChunk = chunks[i];
+        fluid_pressurecell_normalize_chunk(environment,currentChunk);
         fluid_pressurecell_clearArr(currentChunk->d0[CENTER_LOC]);
         fluid_pressurecell_clearArr(currentChunk->u0[CENTER_LOC]);
         fluid_pressurecell_clearArr(currentChunk->v0[CENTER_LOC]);