faster smooth porosity

fluid/fluid_problem.c

@@ -1152,7 +1152,7 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
   double L2 = pow2(0.002);
   int *nodes = NULL;
   int *stack = NULL;
-  double *nodes_w = NULL;
+  double *nodes_w = malloc(sizeof(double)*mesh->n_nodes);
   int *neighbours = malloc(sizeof(int)*mesh->n_elements*(D+1));
   for (int i = 0; i< mesh->n_elements; ++i){
     touched_el_flag[i] = 0;
@@ -1169,12 +1169,12 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
   double *volume = problem->particle_volume;
   for (int i = 0; i < mesh->n_nodes; ++i){
     problem->porosity[i] = 0;
+    nodes_w[i] = -1;
   }
   for (int i = 0; i < problem->n_particles; ++i) {
     if(problem->particle_element_id[i] == -1) continue;
     double *xp = problem->particle_position+i*D;
     vectorClear(nodes);
-    vectorClear(nodes_w);
     vectorClear(stack);
     int p = 0;
     double wtot = 0;
@@ -1185,25 +1185,29 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
       int close = 0;
       for (int j = 0; j < N_N; ++j) {
         int jnode = mesh->elements[elid*N_N+j];
+        int found = nodes_w[jnode] != -1;
+        /*for (int l = 0; l < vectorSize(nodes); ++l){
+          if (nodes[l] == jnode) {
+            found = 1;
+            break;
+          }
+        }
+        */
+        if (found) continue;
         double *x = mesh->x + jnode*3;
         double d2 = 0;
         for (int k = 0; k < D; ++k)
           d2 += pow2(x[k]-xp[k]);
+        *vectorPush(&nodes) = jnode;
         if (d2 < L2) {
           close = 1;
-          int found = 0;
-          for (int l = 0; l < vectorSize(nodes); ++l){
-            if (nodes[l] == jnode) {
-              found = 1;
-              break;
-            }
-          }
-          if (found) continue;
-          *vectorPush(&nodes) = jnode;
           double w = pow2(sqrt(L2)-sqrt(d2))*problem->node_volume[jnode];
-          *vectorPush(&nodes_w) = w;
+          nodes_w[jnode] = w;
           wtot += w;
         }
+        else {
+          nodes_w[jnode] = 0;
+        }
       }
       if (close) {
         for (int j = 0; j < N_N; ++j) {
@@ -1217,7 +1221,8 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
     }
     if (vectorSize(nodes) >= 3) {
       for (size_t j = 0; j < vectorSize(nodes); ++j) {
-        problem->porosity[nodes[j]] += volume[i]*nodes_w[j]/wtot;
+        problem->porosity[nodes[j]] += volume[i]*nodes_w[nodes[j]]/wtot;
+        nodes_w[nodes[j]] = -1;
       }
     }
     else {
@@ -1242,8 +1247,8 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
   }
   free(neighbours);
   vectorFree(nodes);
-  vectorFree(nodes_w);
   vectorFree(stack);
+  free(nodes_w);
 }
 
 FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type, int drag_in_stab)