From 52f1bd1e5b8240ca4beb20680164a86882b4205b Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 24 Jun 2010 13:29:15 +0000 Subject: [PATCH] Talk about an alternative way of doing things. git-svn-id: https://svn.dealii.org/trunk@21316 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-45/doc/results.dox | 113 +++++++++++++++++++++++ 1 file changed, 113 insertions(+) diff --git a/deal.II/examples/step-45/doc/results.dox b/deal.II/examples/step-45/doc/results.dox index b4f55e3127..c3acba4c1b 100644 --- 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. + + + +

Possibilities for extensions

+ +The function LaplaceProblem::make_periodicity_constraints 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 +make_periodicity_constraint_recursively() function — an +implementation detail, not an external interface — put into an anonymous +namespace): +@code +namespace +{ + template + void + make_periodicity_constraints_recursively + (const typename DoFHandler::face_iterator &face_1, + const typename DoFHandler::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; cn_children(); ++c) + make_periodicity_constraints_recursively (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 local_dof_indices_1 (dofs_per_face); + face_1->get_dof_indices (local_dof_indices_1); + + std::vector local_dof_indices_2 (dofs_per_face); + face_2->get_dof_indices (local_dof_indices_2); + + for (unsigned int i=0; i (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 +LaplaceProblem::make_periodicity_constraints() 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. -- 2.39.5