Yet another variation of documenting the bilinear form in step-22.

diff --git a/examples/step-22/step-22.cc b/examples/step-22/step-22.cc
index 712d9eda862b86886ff485843c8a7ef31a5cf616..6a49ce7a17c828088576d96be1dfd5735ddcb43e 100644 (file)
--- a/examples/step-22/step-22.cc
+++ b/examples/step-22/step-22.cc
@@ -687,24 +687,52 @@ namespace Step22
                 phi_p[k]     = fe_values[pressure].value(k, q);
               }
 
+            // Now finally for the bilinear forms of both the system matrix and
+            // the matrix we use for the preconditioner. Recall that the
+            // formulas for these two are
+            // @f{align*}{
+            //   A_{ij} &= a(\varphi_i,\varphi_j)
+            //   \\     &= \underbrace{2(\varepsilon(\varphi_{i,\textbf{u}}),
+            //                           \varepsilon(\varphi_{j,\textbf{u}}))_{\Omega}}
+            //                        _{(1)}
+            //           \;
+            //             \underbrace{- (\textrm{div}\; \varphi_{i,\textbf{u}},
+            //                            \varphi_{j,p})_{\Omega}}
+            //                        _{(2)}
+            //           \;
+            //             \underbrace{- (\varphi_{i,p},
+            //                            \textrm{div}\;
+            //                            \varphi_{j,\textbf{u}})_{\Omega}}
+            //                        _{(3)}
+            // @f}
+            // and
+            // @f{align*}{
+            //   M_{ij} &= \underbrace{(\varphi_{i,p},
+            //                          \varphi_{j,p})_{\Omega}}
+            //                        _{(4)},
+            // @f}
+            // respectively, where $\varphi_{i,\textbf{u}}$ and $\varphi_{i,p}$
+            // are the velocity and pressure components of the $i$th shape
+            // function. The various terms above are then easily recognized in
+            // the following implementation:
             for (unsigned int i = 0; i < dofs_per_cell; ++i)
               {
                 for (unsigned int j = 0; j <= i; ++j)
                   {
                     local_matrix(i, j) +=
-                      // (2 * (grad^s phi_u_i(x_q) * grad^s phi_u_j(x_q))
-                      (2 * (symgrad_phi_u[i] * symgrad_phi_u[j])
-                       // - div phi_u_i(x_q) * phi_p_j(x_q)
-                       - div_phi_u[i] * phi_p[j]
-                       // - phi_p_i(x_q) * div phi_u_j(x_q))
-                       - phi_p[i] * div_phi_u[j]) //
-                      * fe_values.JxW(q);         // * dx
+                      (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) // (1)
+                       - div_phi_u[i] * phi_p[j]                 // (2)
+                       - phi_p[i] * div_phi_u[j])                // (3)
+                      * fe_values.JxW(q);                        // * dx
 
                     local_preconditioner_matrix(i, j) +=
-                      (phi_p[i] * phi_p[j]) // (phi_p_i(x_q) * phi_p_j(x_q))
+                      (phi_p[i] * phi_p[j]) // (4)
                       * fe_values.JxW(q);   // * dx
                   }
-
+                // Note that in the implementation of (1) above, `operator*`
+                // 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
@@ -716,7 +744,6 @@ namespace Step22
                 // 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)
@@ -725,10 +752,6 @@ namespace Step22
               }
           }
 
-        // Note that operator* is overloaded for symmetric tensors,
-        // yielding the scalar product between the two tensors in the first
-        // line of the local matrix contribution.
-
         // Before we can write the local data into the global matrix (and
         // simultaneously use the AffineConstraints object to apply
         // Dirichlet boundary conditions and eliminate hanging node constraints,