]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Finish writing documentation.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 25 Apr 2011 00:00:59 +0000 (00:00 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 25 Apr 2011 00:00:59 +0000 (00:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@23644 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-46/step-46.cc

index bf2e6c8202823adfa39f215b6b7e2078b3ed65f6..dd9e9e2f6c4c6352644e5be490ae5feb049eb89e 100644 (file)
@@ -249,11 +249,25 @@ RightHandSide<dim>::vector_value (const Point<dim> &p,
 
                                  // @sect3{The <code>FluidStructureProblem</code> implementation}
 
+                                 // @sect4{Constructors and helper functions}
+
                                 // Let's now get to the implementation of the
                                 // primary class of this program. The first
-                                // few functions are the constructor and
-
-
+                                // few functions are the constructor and the
+                                // helper functions that can be used to
+                                // determine which part of the domain a cell
+                                // is in. Given the discussion of these
+                                // topics in the introduction, their
+                                // implementation is rather obvious. In the
+                                // constructor, note that we have to
+                                // construct the hp::FECollection object from
+                                // the base elements for Stokes and
+                                // elasticity; using the
+                                // hp::FECollection::push_back function
+                                // assigns them spots zero and one in this
+                                // collection, an order that we have to
+                                // remember and use consistently in the rest
+                                // of the program.
 template <int dim>
 FluidStructureProblem<dim>::
 FluidStructureProblem (const unsigned int stokes_degree,
@@ -298,13 +312,32 @@ cell_is_in_solid_domain (const typename hp::DoFHandler<dim>::cell_iterator &cell
 }
 
 
-
+                                 // @sect4{Meshes and assigning subdomains}
+
+                                // The next pair of functions deals with
+                                // generating a mesh and making sure all
+                                // flags that denote subdomains are
+                                // correct. <code>make_grid</code>, as
+                                // discussed in the introduction, generates
+                                // an $8\times 8$ mesh (or an $8\times
+                                // 8\times 8$ mesh in 3d) to make sure that
+                                // each coarse mesh cell is completely within
+                                // one of the subdomains. After generating
+                                // this mesh, we loop over its boundary and
+                                // set the boundary indicator to one at the
+                                // top boundary, the only place where we set
+                                // nonzero Dirichlet boundary
+                                // conditions. After this, we loop again over
+                                // all cells to set the material indicator
+                                // &mdash; used to denote which part of the
+                                // domain we are in, to either the fluid or
+                                // solid indicator.
 template <int dim>
 void
 FluidStructureProblem<dim>::make_grid ()
 {
-// not quite what we want...  
   GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1);
+
   for (typename Triangulation<dim>::active_cell_iterator
         cell = triangulation.begin_active();
        cell != triangulation.end(); ++cell)
@@ -331,7 +364,24 @@ FluidStructureProblem<dim>::make_grid ()
 }
 
 
-
+                                // The second part of this pair of functions
+                                // determines which finite element to use on
+                                // each cell. Above we have set the material
+                                // indicator for each coarse mesh cell, and
+                                // as mentioned in the introduction, this
+                                // information is inherited from mother to
+                                // child cell upon mesh refinement.
+                                //
+                                // In other words, whenever we have refined
+                                // (or created) the mesh, we can rely on the
+                                // material indicators to be a correct
+                                // description of which part of the domain a
+                                // cell is in. We then use this to set the
+                                // active FE index of the cell to the
+                                // corresponding element of the
+                                // hp::FECollection member variable of this
+                                // class: zero for fluid cells, one for solid
+                                // cells.
 template <int dim>
 void
 FluidStructureProblem<dim>::set_active_fe_indices ()
@@ -339,18 +389,34 @@ FluidStructureProblem<dim>::set_active_fe_indices ()
   for (typename hp::DoFHandler<dim>::active_cell_iterator
          cell = dof_handler.begin_active();
        cell != dof_handler.end(); ++cell)
-    if (cell_is_in_fluid_domain(cell))
-      cell->set_active_fe_index (0);
-    else if (cell_is_in_solid_domain(cell))
-      cell->set_active_fe_index (1);
-    else
-      Assert (false, ExcNotImplemented());
+    {
+      if (cell_is_in_fluid_domain(cell))
+       cell->set_active_fe_index (0);
+      else if (cell_is_in_solid_domain(cell))
+       cell->set_active_fe_index (1);
+      else
+       Assert (false, ExcNotImplemented());
+    }
 }
 
 
-
+                                 // @sect4{<code>FluidStructureProblem::setup_dofs</code>}
+
+                                // The next step is to setup the data
+                                // structures for the linear system. To this
+                                // end, we first have to set the active FE
+                                // indices with the function immediately
+                                // above, then distribute degrees of freedom,
+                                // and then determine constraints on the
+                                // linear system. The latter includes hanging
+                                // node constraints as usual, but also the
+                                // inhomogenous boundary values at the top
+                                // fluid boundary, and zero boundary values
+                                // along the perimeter of the solid
+                                // subdomain.
 template <int dim>
-void FluidStructureProblem<dim>::setup_dofs ()
+void
+FluidStructureProblem<dim>::setup_dofs ()
 {
   set_active_fe_indices ();
   dof_handler.distribute_dofs (fe_collection);
@@ -368,6 +434,7 @@ void FluidStructureProblem<dim>::setup_dofs ()
                                              StokesBoundaryValues<dim>(),
                                              constraints,
                                              velocity_mask);
+
     std::vector<bool> elasticity_mask (dim+1+dim, false);
     for (unsigned int d=dim+1; d<dim+1+dim; ++d)
       elasticity_mask[d] = true;
@@ -378,8 +445,12 @@ void FluidStructureProblem<dim>::setup_dofs ()
                                              elasticity_mask);
   }
 
-                                  // make sure velocity is zero at
-                                  // the interface
+                                  // There are more constraints we have to
+                                  // handle, though: we have to make sure
+                                  // that the velocity is zero at the
+                                  // interface between fluid and solid. The
+                                  // following piece of code was already
+                                  // presented in the introduction:
   {
     std::vector<unsigned int> local_face_dof_indices (stokes_fe.dofs_per_face);
     for (typename hp::DoFHandler<dim>::active_cell_iterator
@@ -397,14 +468,9 @@ void FluidStructureProblem<dim>::setup_dofs ()
                face_is_on_interface = true;
              else if (cell->neighbor(f)->has_children() == true)
                {
-                                                  // neighbor does
-                                                  // have
-                                                  // children. see if
-                                                  // any of the cells
-                                                  // on the other
-                                                  // side are elastic
                  for (unsigned int sf=0; sf<cell->face(f)->n_children(); ++sf)
-                   if (cell_is_in_solid_domain (cell->neighbor_child_on_subface(f, sf)))
+                   if (cell_is_in_solid_domain (cell->neighbor_child_on_subface
+                                                (f, sf)))
                      {
                        face_is_on_interface = true;
                        break;
@@ -420,7 +486,11 @@ void FluidStructureProblem<dim>::setup_dofs ()
            }
   }
 
-
+                                  // At the end of all this, we can declare
+                                  // to the constraints object that we now
+                                  // have all constraints ready to go and
+                                  // that the object can rebuild its internal
+                                  // data structures for better efficiency:
   constraints.close ();
 
   std::cout << "   Number of active cells: "
@@ -430,6 +500,10 @@ void FluidStructureProblem<dim>::setup_dofs ()
             << dof_handler.n_dofs()
             << std::endl;
 
+                                  // The rest of this function is standard:
+                                  // Create a sparsity pattern and use it to
+                                  // initialize the matrix; then also set
+                                  // vectors to their correct sizes.
   {
     CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
                                         dof_handler.n_dofs());
@@ -446,6 +520,20 @@ void FluidStructureProblem<dim>::setup_dofs ()
 
 
 
+                                 // @sect4{<code>FluidStructureProblem::assemble_system</code>}
+
+                                // Following is the central function of this
+                                // program: the one that assembles the linear
+                                // system. It has a long section of setting
+                                // up auxiliary functions at the beginning:
+                                // from creating the quadrature formulas and
+                                // setting up the FEValues, FEFaceValues and
+                                // FESubfaceValues objects necessary to
+                                // integrate the cell terms as well as the
+                                // interface terms for the case where cells
+                                // along the interface come together at same
+                                // size or with differing levels of
+                                // refinement...
 template <int dim>
 void FluidStructureProblem<dim>::assemble_system ()
 {
@@ -465,42 +553,50 @@ void FluidStructureProblem<dim>::assemble_system ()
                                  update_JxW_values |
                                  update_gradients);
 
-  const QGauss<dim-1> face_quadrature(std::max (stokes_degree+2,
-                                               elasticity_degree+2));
+  const QGauss<dim-1> common_face_quadrature(std::max (stokes_degree+2,
+                                                      elasticity_degree+2));
 
   FEFaceValues<dim>    stokes_fe_face_values (stokes_fe,
-                                             face_quadrature,
+                                             common_face_quadrature,
                                              update_JxW_values |
                                              update_normal_vectors |
                                              update_gradients);
   FEFaceValues<dim>    elasticity_fe_face_values (elasticity_fe,
-                                                 face_quadrature,
+                                                 common_face_quadrature,
                                                  update_values);
   FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
-                                                face_quadrature,
+                                                common_face_quadrature,
                                                 update_JxW_values |
                                                 update_normal_vectors |
                                                 update_gradients);
   FESubfaceValues<dim> elasticity_fe_subface_values (elasticity_fe,
-                                                    face_quadrature,
+                                                    common_face_quadrature,
                                                     update_values);
 
-  const unsigned int   stokes_dofs_per_cell     = stokes_fe.dofs_per_cell;
-  const unsigned int   elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell;
+                                  // ...to objects that are needed to
+                                  // describe the local contributions to the
+                                  // global linear system...
+  const unsigned int        stokes_dofs_per_cell     = stokes_fe.dofs_per_cell;
+  const unsigned int        elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell;
 
-  FullMatrix<double>   local_matrix;
-  FullMatrix<double>   local_interface_matrix (elasticity_dofs_per_cell,
-                                              stokes_dofs_per_cell);
-  Vector<double>       local_rhs;
+  FullMatrix<double>        local_matrix;
+  FullMatrix<double>        local_interface_matrix (elasticity_dofs_per_cell,
+                                                   stokes_dofs_per_cell);
+  Vector<double>            local_rhs;
 
   std::vector<unsigned int> local_dof_indices;
   std::vector<unsigned int> neighbor_dof_indices (stokes_dofs_per_cell);
 
-  const RightHandSide<dim>         right_hand_side;
+  const RightHandSide<dim>  right_hand_side;
 
-  const FEValuesExtractors::Vector velocities (0);
-  const FEValuesExtractors::Scalar pressure (dim);
-  const FEValuesExtractors::Vector displacements (dim+1);
+                                  // ...to variables that allow us to extract
+                                  // certain components of the shape
+                                  // functions and cache their values rather
+                                  // than having to recompute them at every
+                                  // quadrature point:
+  const FEValuesExtractors::Vector     velocities (0);
+  const FEValuesExtractors::Scalar     pressure (dim);
+  const FEValuesExtractors::Vector     displacements (dim+1);
 
   std::vector<SymmetricTensor<2,dim> > stokes_phi_grads_u (stokes_dofs_per_cell);
   std::vector<double>                  stokes_div_phi_u   (stokes_dofs_per_cell);
@@ -508,8 +604,14 @@ void FluidStructureProblem<dim>::assemble_system ()
 
   std::vector<Tensor<2,dim> >          elasticity_phi_grad (elasticity_dofs_per_cell);
   std::vector<double>                  elasticity_phi_div  (elasticity_dofs_per_cell);
-  std::vector<Tensor<1,dim> >          elasticity_phi  (elasticity_dofs_per_cell);
-
+  std::vector<Tensor<1,dim> >          elasticity_phi      (elasticity_dofs_per_cell);
+
+                                  // Then comes the main loop over all cells
+                                  // and, as in step-27, the initialization
+                                  // of the hp::FEValues object for the
+                                  // current cell and the extraction of a
+                                  // FEValues object that is appropriate for
+                                  // the current cell:
   typename hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
@@ -523,6 +625,29 @@ void FluidStructureProblem<dim>::assemble_system ()
                           cell->get_fe().dofs_per_cell);
       local_rhs.reinit (cell->get_fe().dofs_per_cell);
 
+                                      // With all of this done, we continue
+                                      // to assemble the cell terms for cells
+                                      // that are part of the Stokes and
+                                      // elastic regions. While we could in
+                                      // principle do this in one formula, in
+                                      // effect implementing the one bilinear
+                                      // form stated in the introduction, we
+                                      // realize that our finite element
+                                      // spaces are chosen in such a way that
+                                      // on each cell, one set of variables
+                                      // (either velocities and pressure, or
+                                      // displacements) are always zero, and
+                                      // consequently a more efficient way of
+                                      // computing local integrals is to do
+                                      // only what's necessary based on an
+                                      // <code>if</code> clause that tests
+                                      // which part of the domain we are in.
+                                      //
+                                      // The actual computation of the local
+                                      // matrix is the same as in step-22 as
+                                      // well as that given in the @ref
+                                      // vector_valued documentation module
+                                      // for the elasticity equations:
       if (cell_is_in_fluid_domain (cell))
        {
          const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
@@ -579,32 +704,100 @@ void FluidStructureProblem<dim>::assemble_system ()
            }
        }
 
+                                      // Once we have the contributions from
+                                      // cell integrals, we copy them into
+                                      // the global matrix (taking care of
+                                      // constraints right away, through the
+                                      // ConstraintMatrix::distribute_local_to_global
+                                      // function). Note that we have not
+                                      // written anything into the
+                                      // <code>local_rhs</code> variable,
+                                      // though we still need to pass it
+                                      // along since the elimination of
+                                      // nonzero boundary values requires the
+                                      // modification of local and
+                                      // consequently also global right hand
+                                      // side values:
       local_dof_indices.resize (cell->get_fe().dofs_per_cell);
       cell->get_dof_indices (local_dof_indices);
-
-                                      // local_rhs==0, but need to do
-                                      // this here because of
-                                      // boundary values
       constraints.distribute_local_to_global (local_matrix, local_rhs,
                                              local_dof_indices,
                                              system_matrix, system_rhs);
 
-                                      // see about face terms
+                                      // The more interesting part of this
+                                      // function is where we see about face
+                                      // terms along the interface between
+                                      // the two subdomains. To this end, we
+                                      // first have to make sure that we only
+                                      // assemble them once even though a
+                                      // loop over all faces of all cells
+                                      // would encounter each part of the
+                                      // interface twice. We arbitrarily make
+                                      // the decision that we will only
+                                      // evaluate interface terms if the
+                                      // current cell is part of the solid
+                                      // subdomain and if, consequently, a
+                                      // face is not at the boundary and the
+                                      // potential neighbor behind it is part
+                                      // of the fluid domain. Let's start
+                                      // with these conditions:
       if (cell_is_in_solid_domain (cell))
-                                        // we are on a solid cell
        for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
          if (cell->at_boundary(f) == false)
            {
+                                              // At this point we know that
+                                              // the current cell is a
+                                              // candidate for integration
+                                              // and that a neighbor behind
+                                              // face <code>f</code>
+                                              // exists. There are now three
+                                              // possibilities:
+                                              //
+                                              // - The neighbor is at the
+                                              //   same refinement level and
+                                              //   has no children.
+                                              // - The neighbor has children.
+                                              // - The neighbor is coarser.
+                                              //
+                                              // In all three cases, we are
+                                              // only interested in it if it
+                                              // is part of the fluid
+                                              // subdomain. So let us start
+                                              // with the first and simplest
+                                              // case: if the neighbor is at
+                                              // the same level, has no
+                                              // children, and is a fluid
+                                              // cell, then the two cells
+                                              // share a boundary that is
+                                              // part of the interface along
+                                              // which we want to integrate
+                                              // interface terms. All we have
+                                              // to do is initialize two
+                                              // FEFaceValues object with the
+                                              // current face and the face of
+                                              // the neighboring cell (note
+                                              // how we find out which face
+                                              // of the neighboring cell
+                                              // borders on the current cell)
+                                              // and pass things off to the
+                                              // function that evaluates the
+                                              // interface terms (the third
+                                              // through fifth arguments to
+                                              // this function provide it
+                                              // with scratch arrays). The
+                                              // result is then again copied
+                                              // into the global matrix,
+                                              // using a function that knows
+                                              // that the DoF indices of rows
+                                              // and columns of the local
+                                              // matrix result from different
+                                              // cells:
              if ((cell->neighbor(f)->level() == cell->level())
                  &&
                  (cell->neighbor(f)->has_children() == false)
                  &&
                  cell_is_in_fluid_domain (cell->neighbor(f)))
                {
-                                                  // same size
-                                                  // neighbors;
-                                                  // neighbor is
-                                                  // fluid cell
                  elasticity_fe_face_values.reinit (cell, f);
                  stokes_fe_face_values.reinit (cell->neighbor(f),
                                                cell->neighbor_of_neighbor(f));
@@ -619,16 +812,32 @@ void FluidStructureProblem<dim>::assemble_system ()
                                                         neighbor_dof_indices,
                                                         system_matrix);
                }
+
+                                              // The second case is if the
+                                              // neighbor has further
+                                              // children. In that case, we
+                                              // have to loop over all the
+                                              // children of the neighbor to
+                                              // see if they are part of the
+                                              // fluid subdomain. If they
+                                              // are, then we integrate over
+                                              // the common interface, which
+                                              // is a face for the neighbor
+                                              // and a subface of the current
+                                              // cell, requiring us to use an
+                                              // FEFaceValues for the
+                                              // neighbor and an
+                                              // FESubfaceValues for the
+                                              // current cell:
              else if ((cell->neighbor(f)->level() == cell->level())
                       &&
                       (cell->neighbor(f)->has_children() == true))
                {
-                                                  // neighbor has children. loop over
-                                                  // the cells adjacent to the commone
-                                                  // interface and see which subdomain
-                                                  // they belong to
-                 for (unsigned int subface=0; subface<cell->face(f)->n_children(); ++subface)
-                   if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface (f, subface)))
+                 for (unsigned int subface=0;
+                      subface<cell->face(f)->n_children();
+                      ++subface)
+                   if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
+                                                (f, subface)))
                      {
                        elasticity_fe_subface_values.reinit (cell,
                                                             f,
@@ -636,29 +845,42 @@ void FluidStructureProblem<dim>::assemble_system ()
                        stokes_fe_face_values.reinit (cell->neighbor_child_on_subface (f, subface),
                                                      cell->neighbor_of_neighbor(f));
 
-                       assemble_interface_term (elasticity_fe_subface_values, stokes_fe_face_values,
-                                                elasticity_phi, stokes_phi_grads_u, stokes_phi_p,
+                       assemble_interface_term (elasticity_fe_subface_values,
+                                                stokes_fe_face_values,
+                                                elasticity_phi,
+                                                stokes_phi_grads_u, stokes_phi_p,
                                                 local_interface_matrix);
 
-                       cell->neighbor_child_on_subface (f, subface)->get_dof_indices (neighbor_dof_indices);
+                       cell->neighbor_child_on_subface (f, subface)
+                         ->get_dof_indices (neighbor_dof_indices);
                        constraints.distribute_local_to_global(local_interface_matrix,
                                                               local_dof_indices,
                                                               neighbor_dof_indices,
                                                               system_matrix);
                      }
                }
+
+                                              // The last option is that the
+                                              // neighbor is coarser. In that
+                                              // case we have to use an
+                                              // FESubfaceValues object for
+                                              // the neighbor and a
+                                              // FEFaceValues for the current
+                                              // cell; the rest is the same
+                                              // as before:
              else if (cell->neighbor_is_coarser(f)
                       &&
                       cell_is_in_fluid_domain(cell->neighbor(f)))
                {
-                                                  // neighbor is coarser
                  elasticity_fe_face_values.reinit (cell, f);
                  stokes_fe_subface_values.reinit (cell->neighbor(f),
                                                   cell->neighbor_of_coarser_neighbor(f).first,
                                                   cell->neighbor_of_coarser_neighbor(f).second);
 
-                 assemble_interface_term (elasticity_fe_face_values, stokes_fe_subface_values,
-                                          elasticity_phi, stokes_phi_grads_u, stokes_phi_p,
+                 assemble_interface_term (elasticity_fe_face_values,
+                                          stokes_fe_subface_values,
+                                          elasticity_phi,
+                                          stokes_phi_grads_u, stokes_phi_p,
                                           local_interface_matrix);
 
                  cell->neighbor(f)->get_dof_indices (neighbor_dof_indices);
@@ -674,25 +896,47 @@ void FluidStructureProblem<dim>::assemble_system ()
 
 
 
+                                // In the function that assembles the global
+                                // system, we passed computing interface
+                                // terms to a separate function we discuss
+                                // here. The key is that even though we can't
+                                // predict the combination of FEFaceValues
+                                // and FESubfaceValues objects, they are both
+                                // derived from the FEFaceValuesBase class
+                                // and consequently we don't have to care:
+                                // the function is simply called with two
+                                // such objects denoting the values of the
+                                // shape functions on the quadrature points
+                                // of the two sides of the face. We then do
+                                // what we always do: we fill the scratch
+                                // arrays with the values of shape functions
+                                // and their derivatives, and then loop over
+                                // all entries of the matrix to compute the
+                                // local integrals. The details of the
+                                // bilinear form we evaluate here are given
+                                // in the introduction.
 template <int dim>
 void
-FluidStructureProblem<dim>::assemble_interface_term (const FEFaceValuesBase<dim>          &elasticity_fe_face_values,
-                                                    const FEFaceValuesBase<dim>          &stokes_fe_face_values,
-                                                    std::vector<Tensor<1,dim> >          &elasticity_phi,
-                                                    std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
-                                                    std::vector<double>                  &stokes_phi_p,
-                                                    FullMatrix<double>                   &local_interface_matrix) const
+FluidStructureProblem<dim>::
+assemble_interface_term (const FEFaceValuesBase<dim>          &elasticity_fe_face_values,
+                        const FEFaceValuesBase<dim>          &stokes_fe_face_values,
+                        std::vector<Tensor<1,dim> >          &elasticity_phi,
+                        std::vector<SymmetricTensor<2,dim> > &stokes_phi_grads_u,
+                        std::vector<double>                  &stokes_phi_p,
+                        FullMatrix<double>                   &local_interface_matrix) const
 {
   Assert (stokes_fe_face_values.n_quadrature_points ==
           elasticity_fe_face_values.n_quadrature_points,
          ExcInternalError());
-
+  const unsigned int n_face_quadrature_points
+    = elasticity_fe_face_values.n_quadrature_points;
+  
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
   const FEValuesExtractors::Vector displacements (dim+1);
 
   local_interface_matrix = 0;
-  for (unsigned int q=0; q<elasticity_fe_face_values.n_quadrature_points; ++q)
+  for (unsigned int q=0; q<n_face_quadrature_points; ++q)
     {
       const Tensor<1,dim> normal_vector = stokes_fe_face_values.normal_vector(q);
 
@@ -715,6 +959,16 @@ FluidStructureProblem<dim>::assemble_interface_term (const FEFaceValuesBase<dim>
 }
 
 
+                                 // @sect4{<code>FluidStructureProblem::solve</code>}
+
+                                // As discussed in the introduction, we use a
+                                // rather trivial solver here: we just pass
+                                // the linear system off to the
+                                // SparseDirectUMFPACK direct solver (see,
+                                // for example, step-29). The only thing we
+                                // have to do after solving is ensure that
+                                // hanging node and boundary value
+                                // constraints are correct.
 template <int dim>
 void
 FluidStructureProblem<dim>::solve ()
@@ -728,10 +982,21 @@ FluidStructureProblem<dim>::solve ()
 
 
 
+                                 // @sect4{<code>FluidStructureProblem::output_results</code>}
 
+                                // Generating graphical output is rather
+                                // trivial here: all we have to do is
+                                // identify which components of the solution
+                                // vector belong to scalars and/or vectors
+                                // (see, for example, step-22 for a previous
+                                // example), and then pass it all on to the
+                                // DataOut class (with the second template
+                                // argument equal to hp::DoFHandler instead
+                                // of the usual default DoFHandler):
 template <int dim>
 void
-FluidStructureProblem<dim>::output_results (const unsigned int refinement_cycle)  const
+FluidStructureProblem<dim>::
+output_results (const unsigned int refinement_cycle)  const
 {
   std::vector<std::string> solution_names (dim, "velocity");
   solution_names.push_back ("pressure");
@@ -765,21 +1030,37 @@ FluidStructureProblem<dim>::output_results (const unsigned int refinement_cycle)
 }
 
 
-
+                                 // @sect4{<code>FluidStructureProblem::refine_mesh</code>}
+
+                                // The next step is to refine the mesh. As
+                                // was discussed in the introduction, this is
+                                // a bit tricky primarily because the fluid
+                                // and the solid subdomains use variables
+                                // that have different physical dimensions
+                                // and for which the absolute magnitude of
+                                // error estimates is consequently not
+                                // directly comparable. We will therefore
+                                // have to scale them. At the top of the
+                                // function, we therefore first compute error
+                                // estimates for the different variables
+                                // separately (using the velocities but not
+                                // the pressure for the fluid domain, and the
+                                // displacements in the solid domain):
 template <int dim>
 void
 FluidStructureProblem<dim>::refine_mesh ()
 {
-  Vector<float> stokes_estimated_error_per_cell (triangulation.n_active_cells());
-  Vector<float> elasticity_estimated_error_per_cell (triangulation.n_active_cells());
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+  Vector<float>
+    stokes_estimated_error_per_cell (triangulation.n_active_cells());
+  Vector<float>
+    elasticity_estimated_error_per_cell (triangulation.n_active_cells());
 
-  const QGauss<dim-1> stokes_quadrature(stokes_degree+2);
-  const QGauss<dim-1> elasticity_quadrature(elasticity_degree+2);
+  const QGauss<dim-1> stokes_face_quadrature(stokes_degree+2);
+  const QGauss<dim-1> elasticity_face_quadrature(elasticity_degree+2);
 
   hp::QCollection<dim-1> face_q_collection;
-  face_q_collection.push_back (stokes_quadrature);
-  face_q_collection.push_back (elasticity_quadrature);
+  face_q_collection.push_back (stokes_face_quadrature);
+  face_q_collection.push_back (elasticity_face_quadrature);
 
   std::vector<bool> stokes_component_mask (dim+1+dim, false);
   for (unsigned int d=0; d<dim; ++d)
@@ -801,11 +1082,52 @@ FluidStructureProblem<dim>::refine_mesh ()
                                       elasticity_estimated_error_per_cell,
                                       elasticity_component_mask);
 
-  stokes_estimated_error_per_cell /= 0.25 * stokes_estimated_error_per_cell.l2_norm();
-  elasticity_estimated_error_per_cell /= elasticity_estimated_error_per_cell.l2_norm();
+                                  // We then normalize error estimates by
+                                  // dividing by their norm and scale the
+                                  // fluid error indicators by a factor of 4
+                                  // as discussed in the introduction. The
+                                  // results are then added together into a
+                                  // vector that contains error indicators
+                                  // for all cells:
+  stokes_estimated_error_per_cell
+    *= 4 . / stokes_estimated_error_per_cell.l2_norm();
+  elasticity_estimated_error_per_cell
+    *= 1. / elasticity_estimated_error_per_cell.l2_norm();
+
+  Vector<float>
+    estimated_error_per_cell (triangulation.n_active_cells());
+
   estimated_error_per_cell += stokes_estimated_error_per_cell;
   estimated_error_per_cell += elasticity_estimated_error_per_cell;
 
+                                  // The second to last part of the function,
+                                  // before actually refining the mesh,
+                                  // involves a heuristic that we have
+                                  // already mentioned in the introduction:
+                                  // because the solution is discontinuous,
+                                  // the KellyErrorEstimator class gets all
+                                  // confused about cells that sit at the
+                                  // boundary between subdomains: it believes
+                                  // that the error is large there because
+                                  // the jump in the gradient is large, even
+                                  // though this is entirely expected and a
+                                  // feature that is in fact present in the
+                                  // exact solution as well and therefore not
+                                  // indicative of any numerical error.
+                                  //
+                                  // Consequently, we set the error
+                                  // indicators to zero for all cells at the
+                                  // interface; the conditions determining
+                                  // which cells this affects are slightly
+                                  // awkward because we have to account for
+                                  // the possibility of adaptively refined
+                                  // meshes, meaning that the neighboring
+                                  // cell can be coarser than the current
+                                  // one, or could in fact be refined some
+                                  // more. The structure of these nested
+                                  // conditions is much the same as we
+                                  // encountered when assembling interface
+                                  // terms in <code>assemble_system</code>.
   {
     unsigned int cell_index = 0;
     for (typename hp::DoFHandler<dim>::active_cell_iterator
@@ -868,6 +1190,14 @@ FluidStructureProblem<dim>::refine_mesh ()
 
 
 
+                                 // @sect4{<code>FluidStructureProblem::run</code>}
+
+                                // This is, as usual, the function that
+                                // controls the overall flow of operation. If
+                                // you've read through tutorial programs
+                                // step-1 through step-6, for example, then
+                                // you are already quite familiar with the
+                                // following structure:
 template <int dim>
 void FluidStructureProblem<dim>::run ()
 {
@@ -898,6 +1228,11 @@ void FluidStructureProblem<dim>::run ()
 
 
 
+                                 // @sect4{The <code>main()</code> function}
+
+                                // This, final, function contains pretty much
+                                // exactly what most of the other tutorial
+                                // programs have:
 int main ()
 {
   try

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.