Merge pull request #12988 from fdrmrc/fix_extractors_in_localrhs_step-22

Use extractors for local_rhs in step-22

doc/news/changes/minor/20211123MarcoFeder

@@ -0,0 +1,3 @@
+Fixed: The step-22 tutorial now uses FEValuesExtractors also for the right-hand side
+<br>
+(Marco Feder, 2021/11/23)

examples/step-22/step-22.cc

@@ -25,6 +25,7 @@
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/function.h>
 #include <deal.II/base/utilities.h>
+#include <deal.II/base/tensor_function.h>
 
 #include <deal.II/lac/block_vector.h>
 #include <deal.II/lac/full_matrix.h>
@@ -42,6 +43,7 @@
 #include <deal.II/dofs/dof_renumbering.h>
 #include <deal.II/dofs/dof_tools.h>
 
+
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_system.h>
 #include <deal.II/fe/fe_values.h>
@@ -216,37 +218,43 @@ namespace Step22
 
 
   // We implement similar functions for the right hand side which for the
-  // current example is simply zero:
+  // current example is simply zero. In the current case, what we want to
+  // represent is the right hand side of the velocity equation, that is, the
+  // body forces acting on the fluid. (The right hand side of the divergence
+  // equation is always zero, since we are dealing with an incompressible
+  // medium.) As such, the right hand side is represented by a `Tensor<1,dim>`
+  // object: a dim-dimensional vector, and the right tool to implement such a
+  // function is to derive it from `TensorFunction<1,dim>`:
   template <int dim>
-  class RightHandSide : public Function<dim>
+  class RightHandSide : public TensorFunction<1, dim>
   {
   public:
     RightHandSide()
-      : Function<dim>(dim + 1)
+      : TensorFunction<1, dim>()
     {}
 
-    virtual double value(const Point<dim> & p,
-                         const unsigned int component = 0) const override;
+    virtual Tensor<1, dim> value(const Point<dim> &p) const override;
 
-    virtual void vector_value(const Point<dim> &p,
-                              Vector<double> &  value) const override;
+    virtual void value_list(const std::vector<Point<dim>> &p,
+                            std::vector<Tensor<1, dim>> &value) const override;
   };
 
 
   template <int dim>
-  double RightHandSide<dim>::value(const Point<dim> & /*p*/,
-                                   const unsigned int /*component*/) const
+  Tensor<1, dim> RightHandSide<dim>::value(const Point<dim> & /*p*/) const
   {
-    return 0;
+    return Tensor<1, dim>();
   }
 
 
   template <int dim>
-  void RightHandSide<dim>::vector_value(const Point<dim> &p,
-                                        Vector<double> &  values) const
+  void RightHandSide<dim>::value_list(const std::vector<Point<dim>> &vp,
+                                      std::vector<Tensor<1, dim>> &values) const
   {
-    for (unsigned int c = 0; c < this->n_components; ++c)
-      values(c) = RightHandSide<dim>::value(p, c);
+    for (unsigned int c = 0; c < vp.size(); ++c)
+      {
+        values[c] = RightHandSide<dim>::value(vp[c]);
+      }
   }
 
 
@@ -629,7 +637,7 @@ namespace Step22
     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 
     const RightHandSide<dim>    right_hand_side;
-    std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim + 1));
+    std::vector<Tensor<1, dim>> rhs_values(n_q_points, Tensor<1, dim>());
 
     // Next, we need two objects that work as extractors for the FEValues
     // object. Their use is explained in detail in the report on @ref
@@ -661,6 +669,7 @@ namespace Step22
     // the outer loop.
     std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
     std::vector<double>                  div_phi_u(dofs_per_cell);
+    std::vector<Tensor<1, dim>>          phi_u(dofs_per_cell);
     std::vector<double>                  phi_p(dofs_per_cell);
 
     for (const auto &cell : dof_handler.active_cell_iterators())
@@ -670,8 +679,8 @@ namespace Step22
         local_preconditioner_matrix = 0;
         local_rhs                   = 0;
 
-        right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
-                                          rhs_values);
+        right_hand_side.value_list(fe_values.get_quadrature_points(),
+                                   rhs_values);
 
         for (unsigned int q = 0; q < n_q_points; ++q)
           {
@@ -680,6 +689,7 @@ namespace Step22
                 symgrad_phi_u[k] =
                   fe_values[velocities].symmetric_gradient(k, q);
                 div_phi_u[k] = fe_values[velocities].divergence(k, q);
+                phi_u[k]     = fe_values[velocities].value(k, q);
                 phi_p[k]     = fe_values[pressure].value(k, q);
               }
 
@@ -729,22 +739,17 @@ namespace Step22
                 // is overloaded for symmetric tensors, yielding the scalar
                 // product between the two tensors.
                 //
-                // For the right-hand side we use the fact that the shape
-                // functions are only non-zero in one component (because our
-                // elements are primitive).  Instead of multiplying the tensor
-                // representing the dim+1 values of shape function i with the
-                // whole right-hand side vector, we only look at the only
-                // non-zero component. The function
-                // FiniteElement::system_to_component_index will return
-                // which component this shape function lives in (0=x velocity,
-                // 1=y velocity, 2=pressure in 2d), which we use to pick out
-                // the correct component of the right-hand side vector to
-                // multiply with.
-                const unsigned int component_i =
-                  fe.system_to_component_index(i).first;
-                local_rhs(i) += (fe_values.shape_value(i, q)   // (phi_u_i(x_q)
-                                 * rhs_values[q](component_i)) // * f(x_q))
-                                * fe_values.JxW(q);            // * dx
+                // For the right-hand side, we need to multiply the (vector of)
+                // velocity shape functions with the vector of body force
+                // right-hand side components, both evaluated at the current
+                // quadrature point. We have implemented the body forces as a
+                // `TensorFunction<1,dim>`, so its values at quadrature points
+                // are already tensors for which the application of `operator*`
+                // against the velocity components of the shape function results
+                // in the dot product, as intended.
+                local_rhs(i) += phi_u[i]            // phi_u_i(x_q)
+                                * rhs_values[q]     // * f(x_q)
+                                * fe_values.JxW(q); // * dx
               }
           }