WIP other strategy using directly face integrals

include/solvers/tracer.h

@@ -298,7 +298,7 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
   assemble_system_matrix() override;
 
   /**
-   *  @brief Assembles the matrix associated with the solver
+   *  @brief Assembles the matrix associated with the solver when DG elements are used.
    */
   void
   assemble_system_matrix_dg();
@@ -309,6 +309,12 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
   void
   assemble_system_rhs() override;
 
+  /**
+   *  @brief Assembles the rhs associated with the solver when DG elements are used.
+   */
+  void
+  assemble_system_rhs_dg();
+
 
   /**
    * @brief Assemble the local matrix for a given cell.

include/solvers/tracer_assemblers.h

@@ -102,6 +102,49 @@ class TracerAssemblerCore : public TracerAssemblerBase<dim>
   std::shared_ptr<SimulationControl> simulation_control;
 };
 
+
+
+/**
+ * @brief Class that assembles the core of the Tracer equation for DG elements.
+ * This class assembles the weak form of:
+ * \f$\mathbf{u} \cdot \nabla T - D \nabla^2 =0 \f$
+ * @tparam dim An integer that denotes the number of spatial dimensions
+ *
+ * @ingroup assemblers
+ */
+
+
+template <int dim>
+class TracerAssemblerDGCore : public TracerAssemblerBase<dim>
+{
+public:
+  TracerAssemblerDGCore(std::shared_ptr<SimulationControl> simulation_control)
+    : simulation_control(simulation_control)
+  {}
+
+  /**
+   * @brief assemble_matrix Assembles the matrix
+   * @param scratch_data (see base class)
+   * @param copy_data (see base class)
+   */
+  virtual void
+  assemble_matrix(TracerScratchData<dim>    &scratch_data,
+                  StabilizedMethodsCopyData &copy_data) override;
+
+
+  /**
+   * @brief assemble_rhs Assembles the rhs
+   * @param scratch_data (see base class)
+   * @param copy_data (see base class)
+   */
+  virtual void
+  assemble_rhs(TracerScratchData<dim>    &scratch_data,
+               StabilizedMethodsCopyData &copy_data) override;
+
+  std::shared_ptr<SimulationControl> simulation_control;
+};
+
+
 /**
  * @brief Class that assembles the transient time arising from BDF time
  * integration for the Tracer equations. For example, if a BDF1 scheme is

include/solvers/tracer_scratch_data.h

@@ -89,13 +89,20 @@ class TracerScratchData
                        quadrature,
                        update_values | update_quadrature_points |
                          update_JxW_values | update_gradients | update_hessians)
+    , fe_face_values_tracer(mapping,
+                            fe_tracer,
+                            face_quadrature,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values |
+                              update_normal_vectors)
     , fe_interface_values_tracer(mapping,
                                  fe_tracer,
                                  face_quadrature,
                                  update_values | update_gradients |
                                    update_quadrature_points |
                                    update_JxW_values | update_normal_vectors)
     , fe_values_fd(mapping, fe_fd, quadrature, update_values)
+    , fe_face_values_fd(mapping, fe_fd, face_quadrature, update_values)
   {
     allocate();
   }
@@ -120,6 +127,12 @@ class TracerScratchData
                        sd.fe_values_tracer.get_quadrature(),
                        update_values | update_quadrature_points |
                          update_JxW_values | update_gradients | update_hessians)
+    , fe_face_values_tracer(sd.fe_face_values_tracer.get_mapping(),
+                            sd.fe_face_values_tracer.get_fe(),
+                            sd.fe_face_values_tracer.get_quadrature(),
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values |
+                              update_normal_vectors)
     , fe_interface_values_tracer(sd.fe_interface_values_tracer.get_mapping(),
                                  sd.fe_interface_values_tracer.get_fe(),
                                  sd.fe_interface_values_tracer.get_quadrature(),
@@ -130,6 +143,10 @@ class TracerScratchData
                    sd.fe_values_fd.get_fe(),
                    sd.fe_values_fd.get_quadrature(),
                    update_values)
+    , fe_face_values_fd(sd.fe_face_values_fd.get_mapping(),
+                        sd.fe_face_values_fd.get_fe(),
+                        sd.fe_face_values_fd.get_quadrature(),
+                        update_values)
   {
     allocate();
   }
@@ -220,6 +237,47 @@ class TracerScratchData
             this->laplacian_phi[q][k] = trace(this->hess_phi[q][k]);
           }
       }
+
+    // Reinitialize face values for all faces
+    n_faces          = cell->n_faces();
+    n_faces_q_points = fe_face_values_tracer.get_quadrature().size();
+
+    face_JxW =
+      std::vector<std::vector<double>>(n_faces,
+                                       std::vector<double>(n_faces_q_points));
+
+
+    this->phi_face = std::vector<std::vector<std::vector<double>>>(
+      n_faces,
+      std::vector<std::vector<double>>(n_faces_q_points,
+                                       std::vector<double>(n_dofs)));
+
+    this->tracer_face_value =
+      std::vector<std::vector<double>>(n_faces,
+                                       std::vector<double>(n_faces_q_points));
+
+    this->face_normal = 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_tracer.reinit(cell, face);
+        this->fe_face_values_tracer.get_function_values(
+          current_solution, this->tracer_face_value[face]);
+
+        for (unsigned int q = 0; q < n_faces_q_points; ++q)
+          {
+            this->face_normal[face][q] =
+              this->fe_face_values_tracer.normal_vector(q);
+
+            face_JxW[face][q] = fe_face_values_tracer.JxW(q);
+            for (const unsigned int k : fe_face_values_tracer.dof_indices())
+              {
+                this->phi_face[face][q][k] =
+                  this->fe_face_values_tracer.shape_value(k, q);
+              }
+          }
+      }
   }
 
   /** @brief Reinitialize the velocity, calculated by the fluid dynamics while also taking into account ALE
@@ -304,6 +362,7 @@ class TracerScratchData
 
   // FEValues for the Tracer problem
   FEValues<dim>          fe_values_tracer;
+  FEFaceValues<dim>      fe_face_values_tracer;
   FEInterfaceValues<dim> fe_interface_values_tracer;
 
   unsigned int n_dofs;
@@ -332,14 +391,30 @@ class TracerScratchData
   std::vector<std::vector<double>>         laplacian_phi;
   std::vector<std::vector<Tensor<1, dim>>> grad_phi;
 
+  // Scratch for the face boundary condition
+  std::vector<std::vector<double>> face_JxW;
+  // First vector is face number, second quadrature point, third DOF
+  std::vector<std::vector<std::vector<double>>> phi_face;
+
+  // First vector is face number, second quadrature point
+  std::vector<std::vector<double>>         tracer_face_value;
+  std::vector<std::vector<Tensor<1, dim>>> face_normal;
+
+
+  unsigned int n_faces;
+  unsigned int n_faces_q_points;
+
 
   /**
    * Scratch component for the Navier-Stokes component
    */
   FEValuesExtractors::Vector velocities;
   // This FEValues must mandatorily be instantiated for the velocity
   FEValues<dim>               fe_values_fd;
+  FEFaceValues<dim>           fe_face_values_fd;
   std::vector<Tensor<1, dim>> velocity_values;
+  // First vector is face number, second quadrature point
+  std::vector<std::vector<Tensor<1, dim>>> velocity_face_value;
 };
 
 #endif