From be723b7edeebe4430c5df701ec467c058e573334 Mon Sep 17 00:00:00 2001
From: Zhuoran Wang <zhrwang@rams.colostate.edu>
Date: Mon, 16 Sep 2019 10:41:18 -0600
Subject: [PATCH] Add Darcy velocity.

---
 examples/step-61/step-61.cc | 60 +++++++++++++++++++++++++++++++++++--
 1 file changed, 57 insertions(+), 3 deletions(-)

diff --git a/examples/step-61/step-61.cc b/examples/step-61/step-61.cc
index 2ee5436a25..0ea6c6dc19 100644
--- a/examples/step-61/step-61.cc
+++ b/examples/step-61/step-61.cc
@@ -43,6 +43,7 @@
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_dg_vector.h>
 #include <deal.II/fe/fe_system.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/fe_face.h>
@@ -64,7 +65,7 @@
 // program into its own namespace:
 namespace Step61
 {
-  using namespace dealii;
+    using namespace dealii;
 
   // @sect3{The WGDarcyEquation class template}
 
@@ -96,6 +97,9 @@ namespace Step61
 
     FESystem<dim>   fe;
     DoFHandler<dim> dof_handler;
+      
+    FE_DGRaviartThomas<dim> fe_dgrt;
+    DoFHandler<dim>         dof_handler_dgrt;
 
     AffineConstraints<double> constraints;
 
@@ -104,6 +108,8 @@ namespace Step61
 
     Vector<double> solution;
     Vector<double> system_rhs;
+      
+    Vector<double>   darcy_velocity;
   };
 
 
@@ -254,7 +260,10 @@ namespace Step61
   template <int dim>
   WGDarcyEquation<dim>::WGDarcyEquation(const unsigned int degree)
     : fe(FE_DGQ<dim>(degree), 1, FE_FaceQ<dim>(degree), 1)
-    , dof_handler(triangulation)
+    , dof_handler(triangulation),
+    
+    fe_dgrt(0),
+    dof_handler_dgrt (triangulation)
 
   {}
 
@@ -293,12 +302,15 @@ namespace Step61
   void WGDarcyEquation<dim>::setup_system()
   {
     dof_handler.distribute_dofs(fe);
+    dof_handler_dgrt.distribute_dofs (fe_dgrt);
 
     std::cout << "   Number of pressure degrees of freedom: "
               << dof_handler.n_dofs() << std::endl;
 
     solution.reinit(dof_handler.n_dofs());
     system_rhs.reinit(dof_handler.n_dofs());
+      
+    darcy_velocity.reinit (dof_handler_dgrt.n_dofs());
 
     {
       constraints.clear();
@@ -683,9 +695,19 @@ namespace Step61
                                         update_values | update_normal_vectors |
                                           update_quadrature_points |
                                           update_JxW_values);
+      
+    FEValues<dim> fe_values_dgrt (fe_dgrt, quadrature_formula,
+                                  update_values   | update_gradients |
+                                  update_quadrature_points | update_JxW_values);
+      
+    FEFaceValues<dim> fe_face_values_dgrt (fe_dgrt, face_quadrature_formula,
+                                           update_values   | update_normal_vectors |
+                                           update_quadrature_points |
+                                           update_JxW_values);
 
     const unsigned int dofs_per_cell    = fe.dofs_per_cell;
     const unsigned int dofs_per_cell_rt = fe_rt.dofs_per_cell;
+    const unsigned int dofs_per_cell_dgrt = fe_dgrt.dofs_per_cell;
 
     const unsigned int n_q_points    = fe_values.get_quadrature().size();
     const unsigned int n_q_points_rt = fe_values_rt.get_quadrature().size();
@@ -696,6 +718,7 @@ namespace Step61
 
 
     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices_dgrt(dofs_per_cell_dgrt);
 
     FullMatrix<double> cell_matrix_M(dofs_per_cell_rt, dofs_per_cell_rt);
     FullMatrix<double> cell_matrix_G(dofs_per_cell_rt, dofs_per_cell);
@@ -717,6 +740,7 @@ namespace Step61
     const FEValuesExtractors::Scalar pressure(dim);
     const FEValuesExtractors::Scalar interior(0);
     const FEValuesExtractors::Scalar face(1);
+    const FEValuesExtractors::Vector velocities_dgrt(0);
 
     const ExactVelocity<dim> exact_velocity;
 
@@ -736,9 +760,14 @@ namespace Step61
     // components. Then, we multiply all these coefficients and call them beta.
     // The numerical velocity is the product of beta and the basis functions of
     // the Raviart-Thomas space.
-    for (const auto &cell : dof_handler.active_cell_iterators())
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end(),
+    cell_dgrt = dof_handler_dgrt.begin_active();
+    for (; cell!=endc; ++cell,++cell_dgrt)
       {
         fe_values.reinit(cell);
+        fe_values_dgrt.reinit(cell_dgrt);
 
         const typename Triangulation<dim>::active_cell_iterator cell_rt = cell;
         fe_values_rt.reinit(cell_rt);
@@ -832,6 +861,15 @@ namespace Step61
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               cell_velocity(k) +=
                 -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
+          
+        // We compute Darcy velocity.
+        // This is same as cell_velocity but used to graph Darcy velocity.
+        cell_dgrt->get_dof_indices(local_dof_indices_dgrt);
+        for (unsigned int k = 0; k<dofs_per_cell_dgrt;++k)
+          for (unsigned int j = 0; j<dofs_per_cell_dgrt; ++j)
+            for (unsigned int i = 0; i<dofs_per_cell;++i)
+              darcy_velocity(local_dof_indices_dgrt[k]) +=
+                -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j));
 
         // Now, we can calculate the numerical velocity at each quadrature point
         // and compute the $L_2$ error on each cell.
@@ -948,6 +986,22 @@ namespace Step61
       std::ofstream face_output("Pressure_Face.vtu");
       data_out_faces.write_vtu(face_output);
     }
+    // Output Darcy velocity vectors.
+    {
+      std::vector<std::string> solution_names (dim, "velocity");
+      std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      data_component_interpretation
+       (dim, DataComponentInterpretation::component_is_part_of_vector);
+          
+      DataOut<dim> data_out_dgrt;
+      data_out_dgrt.attach_dof_handler (dof_handler_dgrt);
+      data_out_dgrt.add_data_vector (darcy_velocity, solution_names,
+                                     DataOut<dim>::type_dof_data,
+                                     data_component_interpretation);
+      data_out_dgrt.build_patches (fe_dgrt.degree);
+      std::ofstream dgrt_output ("Darcy_velocity.vtk");
+      data_out_dgrt.write_vtk (dgrt_output);
+    }
   }
 
 
-- 
2.39.5