2020-10-29-deal.II-step-26

2020.10.29 Step-26 Study

The step-26 is mainly talking about how deal.ii is doing heat transfer.

We can run the default code and derive the results.

The result:

The mesh:

Now I will go through some the possible extensions thought out by myself.

1. Save time-dependent results for Paraview
   By default, the output results are a bunch of vtk files, which are very difficult to be visualized in Paraview because a user needs to open them one by one and export results one by one. The reason for that is that the vtk files do not contain time information. And they are separated files for Paraview. However, I found the following code in tutorial step-52, where the output files are vtu files and a pvd file that contains the time information. So basically, when one wants to postprocess the results, he/she only needs to open the pvd file in Paraview and, therefore, all the time steps are read automatically.

   template <int dim>// zzd
     void HeatEquation<dim>::output_results(const double time,
         const unsigned int time_step) const
     {
       DataOut<dim> data_out;
   
       data_out.attach_dof_handler(dof_handler);
       data_out.add_data_vector(solution, "U");
   
       data_out.build_patches();
   
       data_out.set_flags(DataOutBase::VtkFlags(time, time_step)); // zzd
   //    data_out.set_flags(DataOutBase::VtkFlags(time, timestep_number));
   
       const std::string filename =
       "solution-" + Utilities::int_to_string(time_step, 3) + ".vtu"; //zzd
   //    "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu";
       std::ofstream output(filename);
       data_out.write_vtu(output); // zzd
       //    data_out.write_vtk(output);
   
       static std::vector<std::pair<double, std::string>> times_and_names; // zzd
       static std::string pvd_filename = "solution.pvd"; // zzd
       times_and_names.emplace_back(time, filename);
       std::ofstream pvd_output(pvd_filename);
       DataOutBase::write_pvd_record(pvd_output,times_and_names);
     }

A screenshot of the generated files is depicts as follows:

2. Moving circular heat source
   The initial problem of step-26 is a L-shape simulation domain under two alternatively flashing heat sources. So I think it is possible to make multiple heat sources alternatively falsh and, accordingly, we can produce a moving heat source in this fashion. The following is an example of this idea.

   template <int dim>
     double RightHandSide<dim>::value (const Point<dim> &p,
         const unsigned int component) const
         {
       (void) component;
       AssertIndexRange(component, 1);
       Assert(dim == 2, ExcNotImplemented());
   
       const double time = this->get_time ();
       const double point_within_period =
                                          (time / period - std::floor (
                                               time / period));
   
       if ((point_within_period >= 0.0) && (point_within_period <= 0.2))
           {
         if (std::pow (p[0] + 0.9, 2.0) + std::pow (p[1] - 0.9, 2.0) <= 0.01)
           return 1;
         else
           return 0;
       }
       else if ((point_within_period >= 0.2) && (point_within_period <= 0.4))
                {
         if (std::pow (p[0] + 0.5, 2.0) + std::pow (p[1] - 0.9, 2.0) <= 0.01)
           return 1;
         else
           return 0;
       }
       else if ((point_within_period >= 0.4) && (point_within_period <= 0.6))
                {
         if (std::pow (p[0] + 0.1, 2.0) + std::pow (p[1] - 0.9, 2.0) <= 0.01)
           return 1;
         else
           return 0;
       }
       else if ((point_within_period >= 0.6) && (point_within_period <= 0.8))
                {
         if (std::pow (p[0] - 0.5, 2.0) + std::pow (p[1] - 0.9, 2.0) <= 0.01)
           return 1;
         else
           return 0;
       }
       else if ((point_within_period >= 0.8) && (point_within_period <= 1.0))
                {
         if (std::pow (p[0] - 0.9, 2.0) + std::pow (p[1] - 0.9, 2.0) <= 0.01)
           return 1;
         else
           return 0;
       }
       else
         return 0;
     }

3. Improving adaptive mesh
   As mentioned in the original tutorial document https://www.dealii.org/current/doxygen/deal.II/step_26.html, the adaptively refined meshes look rather random. There are some islands where cells have been refined but that are surrounded by coarser cells, as shown in the following picture:

These can be fixed by considering putting an argument to the Triangulation class, which is as the following:

HeatEquation<dim>::HeatEquation ()
      : triangulation (Triangulation<dim>::maximum_smoothing) // newly added
            , fe (1)
            , dof_handler (triangulation)
            , time (0.0)
            , time_step (1. / 500)
                         , timestep_number (0)
            , theta (0.5)
  {
  }

The maximum_smoothing argument includes almost all smoothing algorithms, as stated in the file of tria.h. I quote it here for reference:

enum MeshSmoothing
  {
    /**
     * No mesh smoothing at all, except that meshes have to remain one-
     * irregular.
     */
    none = 0x0,
    /**
     * It can be shown, that degradation of approximation occurs if the
     * triangulation contains vertices which are member of cells with levels
     * differing by more than one. One such example is the following:
     *
     * @image html limit_level_difference_at_vertices.png ""
     *
     * It would seem that in two space dimensions, the maximum jump in levels
     * between cells sharing a common vertex is two (as in the example above).
     * However, this is not true if more than four cells meet at a vertex. It
     * is not uncommon that a
     * @ref GlossCoarseMesh "coarse (initial) mesh" contains vertices at which
     * six or even eight cells meet, when small features of the domain have to
     * be resolved even on the coarsest mesh. In that case, the maximum
     * difference in levels is three or four, respectively. The problem gets
     * even worse in three space dimensions.
     *
     * Looking at an interpolation of the second derivative of the finite
     * element solution (assuming bilinear finite elements), one sees that the
     * numerical solution is almost totally wrong, compared with the true
     * second derivative. Indeed, on regular meshes, there exist sharp
     * estimations that the H<sup>2</sup>-error is only of order one, so we
     * should not be surprised; however, the numerical solution may show a
     * value for the second derivative which may be a factor of ten away from
     * the true value. These problems are located on the small cell adjacent
     * to the center vertex, where cells of non-subsequent levels meet, as
     * well as on the upper and right neighbor of this cell (but with a less
     * degree of deviation from the true value).
     *
     * If the smoothing indicator given to the constructor contains the bit
     * for #limit_level_difference_at_vertices, situations as the above one
     * are eliminated by also marking the upper right cell for refinement.
     *
     * In case of anisotropic refinement, the level of a cell is not linked to
     * the refinement of a cell as directly as in case of isotropic
     * refinement. Furthermore, a cell can be strongly refined in one
     * direction and not or at least much less refined in another. Therefore,
     * it is very difficult to decide, which cases should be excluded from the
     * refinement process. As a consequence, when using anisotropic
     * refinement, the #limit_level_difference_at_vertices flag must not be
     * set. On the other hand, the implementation of multigrid methods in
     * deal.II requires that this bit be set.
     */
    limit_level_difference_at_vertices = 0x1,
    /**
     * Single cells which are not refined and are surrounded by cells which
     * are refined usually also lead to a sharp decline in approximation
     * properties locally. The reason is that the nodes on the faces between
     * unrefined and refined cells are not real degrees of freedom but carry
     * constraints. The patch without additional degrees of freedom is thus
     * significantly larger then the unrefined cell itself. If in the
     * parameter passed to the constructor the bit for
     * #eliminate_unrefined_islands is set, all cells which are not flagged
     * for refinement but which are surrounded by more refined cells than
     * unrefined cells are flagged for refinement. Cells which are not yet
     * refined but flagged for that are accounted for the number of refined
     * neighbors. Cells on the boundary are not accounted for at all. An
     * unrefined island is, by this definition also a cell which (in 2D) is
     * surrounded by three refined cells and one unrefined one, or one
     * surrounded by two refined cells, one unrefined one and is at the
     * boundary on one side. It is thus not a true island, as the name of the
     * flag may indicate. However, no better name came to mind to the author
     * by now.
     */
    eliminate_unrefined_islands = 0x2,
    /**
     * A triangulation of patch level 1 consists of patches, i.e. of cells
     * that are refined once. This flag ensures that a mesh of patch level 1
     * is still of patch level 1 after coarsening and refinement. It is,
     * however, the user's responsibility to ensure that the mesh is of patch
     * level 1 before calling
     * Triangulation::execute_coarsening_and_refinement() the first time. The
     * easiest way to achieve this is by calling global_refine(1) straight
     * after creation of the triangulation.  It follows that if at least one
     * of the children of a cell is or will be refined than all children need
     * to be refined. If the #patch_level_1 flag is set, than the flags
     * #eliminate_unrefined_islands, #eliminate_refined_inner_islands and
     * #eliminate_refined_boundary_islands will be ignored as they will be
     * fulfilled automatically.
     */
    patch_level_1 = 0x4,
    /**
     * Each
     * @ref GlossCoarseMesh "coarse grid"
     * cell is refined at least once,
     * i.e., the triangulation
     * might have active cells on level 1 but not on level 0. This flag
     * ensures that a mesh which has coarsest_level_1 has still
     * coarsest_level_1 after coarsening and refinement. It is, however, the
     * user's responsibility to ensure that the mesh has coarsest_level_1
     * before calling execute_coarsening_and_refinement the first time. The
     * easiest way to achieve this is by calling global_refine(1) straight
     * after creation of the triangulation. It follows that active cells on
     * level 1 may not be coarsened.
     *
     * The main use of this flag is to ensure that each cell has at least one
     * neighbor in each coordinate direction (i.e. each cell has at least a
     * left or right, and at least an upper or lower neighbor in 2d). This is
     * a necessary precondition for some algorithms that compute finite
     * differences between cells. The DerivativeApproximation class is one of
     * these algorithms that require that a triangulation is coarsest_level_1
     * unless all cells already have at least one neighbor in each coordinate
     * direction on the coarsest level.
     */
    coarsest_level_1 = 0x8,
    /**
     * This flag is not included in @p maximum_smoothing. The flag is
     * concerned with the following case: consider the case that an unrefined
     * and a refined cell share a common face and that one of the children of
     * the refined cell along the common face is flagged for further
     * refinement. In that case, the resulting mesh would have more than one
     * hanging node along one or more of the edges of the triangulation, a
     * situation that is not allowed. Consequently, in order to perform the
     * refinement, the coarser of the two original cells is also going to be
     * refined.
     *
     * However, in many cases it is sufficient to refine the coarser of the
     * two original cells in an anisotropic way to avoid the case of multiple
     * hanging vertices on a single edge. Doing only the minimal anisotropic
     * refinement can save cells and degrees of freedom. By specifying this
     * flag, the library can produce these anisotropic refinements.
     *
     * The flag is not included by default since it may lead to
     * anisotropically refined meshes even though no cell has ever been
     * refined anisotropically explicitly by a user command. This surprising
     * fact may lead to programs that do the wrong thing since they are not
     * written for the additional cases that can happen with anisotropic
     * meshes, see the discussion in the introduction to step-30.
     */
    allow_anisotropic_smoothing = 0x10,
    /**
     * This algorithm seeks for isolated cells which are refined or flagged
     * for refinement. This definition is unlike that for
     * #eliminate_unrefined_islands, which would mean that an island is
     * defined as a cell which is refined but more of its neighbors are not
     * refined than are refined. For example, in 2D, a cell's refinement would
     * be reverted if at most one of its neighbors is also refined (or refined
     * but flagged for coarsening).
     *
     * The reason for the change in definition of an island is, that this
     * option would be a bit dangerous, since if you consider a chain of
     * refined cells (e.g. along a kink in the solution), the cells at the two
     * ends would be coarsened, after which the next outermost cells would
     * need to be coarsened. Therefore, only one loop of flagging cells like
     * this could be done to avoid eating up the whole chain of refined cells
     * (`chain reaction'...).
     *
     * This algorithm also takes into account cells which are not actually
     * refined but are flagged for refinement. If necessary, it takes away the
     * refinement flag.
     *
     * Actually there are two versions of this flag,
     * #eliminate_refined_inner_islands and
     * #eliminate_refined_boundary_islands. There first eliminates islands
     * defined by the definition above which are in the interior of the
     * domain, while the second eliminates only those islands if the cell is
     * at the boundary. The reason for this split of flags is that one often
     * wants to eliminate such islands in the interior while those at the
     * boundary may well be wanted, for example if one refines the mesh
     * according to a criterion associated with a boundary integral or if one
     * has rough boundary data.
     */
    eliminate_refined_inner_islands = 0x100,
    /**
     * The result of this flag is very similar to
     * #eliminate_refined_inner_islands. See the documentation there.
     */
    eliminate_refined_boundary_islands = 0x200,
    /**
     * This flag prevents the occurrence of unrefined islands. In more detail:
     * It prohibits the coarsening of a cell if 'most of the neighbors' will
     * be refined after the step.
     */
    do_not_produce_unrefined_islands = 0x400,

    /**
     * This flag sums up all smoothing algorithms which may be performed upon
     * refinement by flagging some more cells for refinement.
     */
    smoothing_on_refinement =
      (limit_level_difference_at_vertices | eliminate_unrefined_islands),
    /**
     * This flag sums up all smoothing algorithms which may be performed upon
     * coarsening by flagging some more cells for coarsening.
     */
    smoothing_on_coarsening =
      (eliminate_refined_inner_islands | eliminate_refined_boundary_islands |
       do_not_produce_unrefined_islands),

    /**
     * This flag includes all the above ones (therefore combines all
     * smoothing algorithms implemented), with the exception of anisotropic
     * smoothing.
     */
    maximum_smoothing = 0xffff ^ allow_anisotropic_smoothing
  };

After considering the smoothing algorithm, the mesh becomes much better as shown in the below picture:

4. A user can convert to 3D very easily by just assigning

   HeatEquation<3> heat_equation_solver;

It should be noted that the following line should be commented off. This is just an assertion that the developers are encouraging the user to contribute to the library. If one is interested, he/she can submit an pull request for the 3D model.

//    Assert(dim == 2, ExcNotImplemented());  // commented off intentionally

5. 3D moving heat source
   

6. Make top surface a Neuuman boundary

   for (const auto &cell : triangulation.active_cell_iterators ())
         for (const auto &face : cell->face_iterators ())
         {
           if (face->at_boundary () == true &&
               face->center ()[2] == 1.0)
               {
             face->set_boundary_id (6);
           }
         }

7. Multi-track scanning
   

8. A part-level simulation for the first time
   

Future work:

1. The FEM information flow control should recieve more attention, such as how exactly the boundary_values are applied on the final linear equations system.

2. The “friend” and “virtual” in C++ is still not crystal clear and should be further studied.

3. I don’t understand how the error estimation inequality is derived.
   

4. Manifold unclear either.

5. 3D moving heat source.

6. 3D moving heat source with nonlinear material properties.

7. Mechanical simulation, plasticity

The official tutorial document is as follows,
https://www.dealii.org/current/doxygen/deal.II/step_26.html