smooth a

fluid/fluid_problem.c

@@ -1012,6 +1012,17 @@ static void fluid_problem_dg_to_cg(FluidProblem *problem, const double *dg, doub
   }
 }
 
+static void fluid_problem_cg_to_dg(FluidProblem *problem, const double *cg, double *dg) {
+  const Mesh *mesh = problem->mesh;
+  for (int iel = 0; iel < mesh->n_elements; ++iel) {
+    const int *el = &mesh->elements[iel*N_N];
+    for (int i = 0; i< N_N; ++i) {
+      dg[iel*N_N+i] = cg[el[i]];
+    }
+  }
+}
+
+
 static void fluid_problem_dg_to_cg_grad(FluidProblem *problem, const double *dg, double *grad_cg) {
   const Mesh *mesh = problem->mesh;
   for (int i = 0; i< mesh->n_nodes*D; ++i) {
@@ -1123,10 +1134,14 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
   compute_weak_boundary_conditions(problem, solution_old, dt, all_local_vector, all_local_matrix);
   fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
   if (problem->n_fluids == 2){
-    double *grad_a_cg = malloc(sizeof(double)*mesh->n_nodes*D);
-    fluid_problem_dg_to_cg_grad(problem, problem->concentration, problem->grad_a_cg);
+    double *a_cg = malloc((problem->mesh->n_nodes)*sizeof(double));
+    double *a_dg = malloc((problem->mesh->n_elements)*3*sizeof(double));
+    fluid_problem_dg_to_cg(problem,problem->concentration,a_cg);
+    fluid_problem_cg_to_dg(problem,a_cg,a_dg);
+    fluid_problem_dg_to_cg_grad(problem, a_dg, problem->grad_a_cg);
     fluid_problem_surface_tension(problem, solution_old, problem->grad_a_cg, all_local_vector);
-    free(grad_a_cg);
+    free(a_cg);
+    free(a_dg);
   }
   for (int iel=0; iel < mesh->n_elements; ++iel) {
     double *local_vector = &all_local_vector[local_size*iel];