WIP

chaos-polymtl · Sep 24, 2024 · b41fac6 · b41fac6
1 parent 4d9751c
commit b41fac6
Show file tree

Hide file tree

Showing 3 changed files with 74 additions and 20 deletions.
diff --git a/include/solvers/tracer_scratch_data.h b/include/solvers/tracer_scratch_data.h
@@ -262,6 +262,7 @@ class TracerScratchData
     for (const auto face : cell->face_indices())
       {
         fe_face_values_tracer.reinit(cell, face);
+        const auto cell_neighbor = cell->neighbor(face);
         this->fe_face_values_tracer.get_function_values(
           current_solution, this->tracer_face_value[face]);
 
@@ -323,6 +324,17 @@ class TracerScratchData
         velocity_values[q] += drift_velocity_tensor;
       }
 
+
+    // To fix for ALE
+    this->velocity_face_value = std::vector<std::vector<Tensor<1, dim>>>(
+      n_faces, std::vector<Tensor<1, dim>>(n_faces_q_points));
+    for (const auto face : cell->face_indices())
+      {
+        fe_face_values_fd.reinit(cell, face);
+        this->fe_face_values_fd[velocities].get_function_values(
+          current_solution, this->velocity_face_value[face]);
+      }
+
     if (!ale.enabled())
       return;
 

diff --git a/source/solvers/tracer.cc b/source/solvers/tracer.cc
@@ -489,19 +489,19 @@ Tracer<dim>::assemble_system_rhs_dg()
 
       std::vector<double> average_value(q_points.size());
       fe_iv.get_average_of_function_values(evaluation_point, average_value);
-      Tensor<1, dim> beta;
-      //  beta[0] = 1;
-      //  for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
-      //    {
-      //      const double beta_dot_n = beta * normals[qpoint];
-      //      for (unsigned int i = 0; i < n_dofs; ++i)
-      //        for (unsigned int j = 0; j < n_dofs; ++j)
-      //          copy_data_face.cell_matrix(i, j) -=
-      //            fe_iv.jump_in_shape_values(i, qpoint) // [\phi_i]
-      //            * average_value[qpoint]               //
-      //            * beta_dot_n                          // (\beta .n)
-      //            * JxW[qpoint];                        // dx
-      //    }
+      Tensor<1, dim> beta =
+        scratch_data.for (unsigned int qpoint = 0; qpoint < q_points.size();
+                          ++qpoint)
+      {
+        const double beta_dot_n = beta * normals[qpoint];
+        for (unsigned int i = 0; i < n_dofs; ++i)
+          for (unsigned int j = 0; j < n_dofs; ++j)
+            copy_data_face.cell_matrix(i, j) -=
+              fe_iv.jump_in_shape_values(i, qpoint) // [\phi_i]
+              * average_value[qpoint]               //
+              * beta_dot_n                          // (\beta .n)
+              * JxW[qpoint];                        // dx
+      }
     };
 
 
@@ -902,10 +902,10 @@ Tracer<dim>::postprocess_tracer_flow_rate(const VectorType &current_solution_fd)
                            normal_vector_tracer) *
                         fe_face_values_tracer.JxW(q);
                     } // end loop on quadrature points
-                }     // end face is a boundary face
-            }         // end loop on faces
-        }             // end condition cell at boundary
-    }                 // end loop on cells
+                } // end face is a boundary face
+            } // end loop on faces
+        } // end condition cell at boundary
+    } // end loop on cells
 
 
   // Sum across all cores

diff --git a/source/solvers/tracer_assemblers.cc b/source/solvers/tracer_assemblers.cc
@@ -264,9 +264,29 @@ TracerAssemblerDGCore<dim>::assemble_matrix(
     } // end loop on quadrature points
 
   // Loop over the faces to assemble the flux term
-  // for (const auto face : scratch_data.->face_indices())
-  //  {
-  //  }
+  for (unsigned int f = 0; f < scratch_data.n_faces; ++f)
+    {
+      for (unsigned int q = 0; q < scratch_data.n_faces_q_points; ++q)
+        {
+          const double flux = scratch_data.velocity_face_value[f][q] *
+                              scratch_data.face_normal[f][q];
+          // If the flux is negative the the velocity is outbound since
+          // deal.II normales point outwards In that case assemble the flux
+          // if (flux > 0)
+          {
+            for (unsigned int i = 0; i < n_dofs; ++i)
+              {
+                const auto phi_T_i = scratch_data.phi_face[f][q][i];
+                for (unsigned int j = 0; j < n_dofs; ++j)
+                  {
+                    const auto phi_T_j = scratch_data.phi_face[f][q][j];
+                    local_matrix(i, j) +=
+                      flux * phi_T_i * phi_T_j * scratch_data.face_JxW[f][q];
+                  }
+              }
+          }
+        }
+    }
 }
 
 
@@ -311,6 +331,28 @@ TracerAssemblerDGCore<dim>::assemble_rhs(TracerScratchData<dim> &scratch_data,
                           JxW;
         }
     } // end loop on quadrature points
+
+  // Loop over the faces to assemble the flux term
+  for (unsigned int f = 0; f < scratch_data.n_faces; ++f)
+    {
+      for (unsigned int q = 0; q < scratch_data.n_faces_q_points; ++q)
+        {
+          const double flux = scratch_data.velocity_face_value[f][q] *
+                              scratch_data.face_normal[f][q];
+          // If the flux is negative the the velocity is outbound since
+          // deal.II normales point outwards In that case assemble the flux
+          // if (flux > 0)
+          {
+            for (unsigned int i = 0; i < n_dofs; ++i)
+              {
+                const auto phi_T_i      = scratch_data.phi_face[f][q][i];
+                const auto tracer_value = scratch_data.tracer_face_value[f][q];
+                local_rhs(i) -=
+                  flux * phi_T_i * tracer_value * scratch_data.face_JxW[f][q];
+              }
+          }
+        }
+    }
 }
 template class TracerAssemblerDGCore<2>;
 template class TracerAssemblerDGCore<3>;