Talk about an alternative way of doing things.

git-svn-id: https://svn.dealii.org/trunk@21316 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/examples/step-45/doc/results.dox b/deal.II/examples/step-45/doc/results.dox
index b4f55e3127fb93dca9200e6e1b3c53032742cf83..c3acba4c1b1e22aeb46c878e12ce9884df658e7c 100644 (file)
--- a/deal.II/examples/step-45/doc/results.dox
+++ b/deal.II/examples/step-45/doc/results.dox
@@ -15,3 +15,116 @@ required by the homogeneous Dirichlet boundary conditions. On the left and
 right parts the values coincide with each other, just as we wanted:
 
 @image html step-45.solution.png
+
+Note also that the solution is clearly not left-right symmetric and so would
+not likely have been periodic had we prescribed, for example, homogeneous
+Neumann boundary condition. However, it is periodic thanks to the constraints
+imposed.
+
+
+<a name="extensions"></a>
+<h3>Possibilities for extensions</h3>
+
+The function <code>LaplaceProblem::make_periodicity_constraints</code> is
+relatively simple in that it just matches the location of degrees of
+freedom. This makes it flexible when the periodicity boundary conditions are
+posed not just on opposite faces of the unit rectangle but on separate parts
+of a possibly more complicated domain. Or, if one wanted to "twist" the
+boundary condition by prescribing, for example,
+@f{align*}
+   u(0,y) &= u(1,1-y) \qquad &\text{for }y\in(0,1).
+@f}
+
+On the other hand, the function is somewhat limited by the assumption that the
+domain is two-dimensional and that we only use $Q_1$ elements. The former
+assumption is easily lifted by looping over all four vertices of a face in 3d,
+but the latter is somewhat more complicated to lift because we have assumed
+that degrees of freedom are only located in vertices. In the following,
+therefore, let us describe a function that computes the same constraints but
+in a dimension-independent way and for any finite element one may want to
+consider.
+
+The idea is to work recursively on pairs of faces. For example, let us start
+with the left and right face of the (single) coarse mesh cell. They need to
+match, but they are not active (i.e. they are further refined) and so there
+are no degrees of freedom associated with these faces. However, if the two
+current faces are periodic, then so are the zeroth children of the two as well
+as the respective first children, etc. We can then in turn work on each of
+these pairs of faces. If they are not active, we may recurse further into this
+refinement tree until we find a pair of active faces. In that case, we enter
+equivalences between matching degrees of freedom on the two active faces.
+
+An implementation of this idea would look like follows (with the
+<code>make_periodicity_constraint_recursively()</code> function &mdash; an
+implementation detail, not an external interface &mdash; put into an anonymous
+namespace):
+@code
+namespace 
+{
+  template <int dim>
+  void
+  make_periodicity_constraints_recursively
+  (const typename DoFHandler<dim>::face_iterator &face_1,
+   const typename DoFHandler<dim>::face_iterator &face_2,
+   ConstraintMatrix &constraints)
+  {
+    Assert (face_1->n_children() == face_2->n_children(),
+           ExcNotImplemented());
+    if (face_1->has_children())
+      {
+       for (unsigned int c=0; c<face_1->n_children(); ++c)
+         make_periodicity_constraints_recursively<dim> (face_1->child(c),
+                                                        face_2->child(c),
+                                                        constraints);
+      }
+    else
+      {
+       const unsigned int dofs_per_face
+         = face_1->get_dof_handler().get_fe().dofs_per_face;
+       
+       std::vector<unsigned int> local_dof_indices_1 (dofs_per_face);
+       face_1->get_dof_indices (local_dof_indices_1);
+       
+       std::vector<unsigned int> local_dof_indices_2 (dofs_per_face);
+       face_2->get_dof_indices (local_dof_indices_2);
+
+       for (unsigned int i=0; i<dofs_per_face; ++i)
+         {
+           constraints.add_line (local_dof_indices_1[i]);
+           constraints.add_entry (local_dof_indices_1[i],
+                                  local_dof_indices_2[i],
+                                  1.0);
+         }
+      }
+  }
+}
+
+
+void LaplaceProblem::make_periodicity_constraints ()
+{
+  make_periodicity_constraints_recursively<2> (dof_handler.begin(0)->face(0),
+                                               dof_handler.begin(0)->face(1),
+                                               constraints);
+}
+@endcode
+
+The implementation of the recursive function should be mostly self explanatory
+given the discussion above. The
+<code>LaplaceProblem::make_periodicity_constraints()</code> function simply
+calls the former with matching faces of the first (and only) coarse mesh cell
+on refinement level 0. Note that when calling the recursive function we have
+to explicitly specify the template argument since the compiler can not deduce
+it (the template argument is only used in a "non-deducible context").
+
+This function is now dimension and finite element independent, but it still
+has the restriction that it assumes that the mesh is uniformly refined (or, in
+fact, only that matching periodic faces are refined equally). We check this at
+the beginning by asserting that both faces have the same number of children
+(that includes that neither have any children, i.e. that both are active).
+On the other hand, the function above can be extended to also allow this sort
+of thing. In that case, if we encounter the situation that only one cell is
+refined, we would have to recurse into its children and interpolate their
+degrees of freedom with respect to the degrees of freedom to the coarser
+matching face. This can use the same facilities the finite element classes
+already provide for computing constraints based on hanging nodes. We leave
+implementing this as an exercise, however.