From 17c4960b77501ecb188b698e386f55e87e1ef801 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 6 Jun 2012 16:12:44 +0000 Subject: [PATCH] Add a section. git-svn-id: https://svn.dealii.org/trunk@25608 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/doxygen/headers/constraints.h | 70 +++++++++++++++++++++-- 1 file changed, 64 insertions(+), 6 deletions(-) diff --git a/deal.II/doc/doxygen/headers/constraints.h b/deal.II/doc/doxygen/headers/constraints.h index ae120a6138..b4493e1564 100644 --- a/deal.II/doc/doxygen/headers/constraints.h +++ b/deal.II/doc/doxygen/headers/constraints.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2010 by the deal.II authors +// Copyright (C) 2010, 2012 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -304,13 +304,13 @@ * system matrix and the right-hand-side with the parameter * use_inhomogeneities_for_rhs = false (default) * - Set the solution to zero in the inhomogeneous constrained components - * using the ConstraintMatrix::set_zero() function (or start with a solution + * using the ConstraintMatrix::set_zero() function (or start with a solution * vector equal to zero) * - solve() the linear system * - Apply ConstraintMatrix::distribute() to the solution - * + * * Second approach: - * - Use the ConstraintMatrix::distribute_local_to_global() function with the parameter + * - Use the ConstraintMatrix::distribute_local_to_global() function with the parameter * use_inhomogeneities_for_rhs = true and apply it to * the system matrix and the right-hand-side * - Set the concerning components of the solution to the inhomogeneous @@ -324,6 +324,64 @@ * and you may want to call ConstraintMatrix::distribute() anyway if you have * additional constraints such as from hanging nodes. * + * Of course, both approaches lead to the same final answer but in different + * ways. Using approach (i.e., when using use_inhomogeneities_for_rhs = false + * in ConstraintMatrix::distribute_local_to_global()), the linear system we + * build has zero entries in the right hand side in all those places where a + * degree of freedom is constrained, and some positive value on the matrix + * diagonal of these lines. Consequently, the solution vector of the linear + * system will have a zero value for inhomogeneously constrained degrees of + * freedom and we need to call ConstraintMatrix::distribute() to give these + * degrees of freedom their correct nonzero values. + * + * On the other hand, in the second approach, the matrix diagonal element and + * corresponding right hand side entry for inhomogeneously constrained degrees + * of freedom are so that the solution of the linear system already has the + * correct value (e.g., if the constraint is that $x_{13}=42$ then row $13$ if + * the matrix is empty with the exception of the diagonal entry, and + * $b_{13}/A_{13,13}=42$ so that the solution of $Ax=b$ must satisfy + * $x_{13}=42$ as desired). As a consequence, we do not need to call + * ConstraintMatrix::distribute() after solving to fix up inhomogeneously + * constrained components of the solution, though there is also no harm in + * doing so. + * + * There remains the question of which of the approaches to take and why we + * need to set to zero the values of the solution vector in the first + * approach. The answer to both questions has to do with how iterative solvers + * solve the linear system. To this end, consider that we typically stop + * iterations when the residual has dropped below a certain fraction of the + * norm of the right hand side, or, alternatively, a certain fraction of the + * norm of the initial residual. Now consider this: + * + * - In the first approach, the right hand side entries for constrained + * degrees of freedom are zero, i.e., the norm of the right hand side + * really only consists of those parts that we care about. On the other + * hand, if we start with a solution vector that is not zero in + * constrained entries, then the initial residual is very large because + * the value that is currently in the solution vector does not match the + * solution of the linear system (which is zero in these components). + * Thus, if we stop iterations once we have reduced the initial residual + * by a certain factor, we may reach the threshold after a single + * iteration because constrained degrees of freedom are resolved by + * iterative solvers in just one iteration. If the initial residual + * was dominated by these degrees of freedom, then we see a steep + * reduction in the first step although we did not really make much + * progress on the remainder of the linear system in this just one + * iteration. We can avoid this problem by either stopping iterations + * once the norm of the residual reaches a certain fraction of the + * norm of the right hand side, or we can set the solution + * components to zero (thus reducing the initial residual) and iterating + * until we hit a certain fraction of the norm of the initial + * residual. + * - In the second approach, we get the same problem if the starting vector + * in the iteration is zero, since then then the residual may be + * dominated by constrained degrees of freedom having values that do not + * match the values we want for them at the solution. We can again + * circumvent this problem by setting the corresponding elements of the + * solution vector to their correct values, by calling + * ConstraintMatrix::distribute() before solving the linear system + * (and then, as necessary, a second time after solving). + * * *

Dealing with conflicting constraints

* @@ -348,12 +406,12 @@ * hanging node constraints at this point and consequently would not * satisfy the regularity properties of the element chosen (e.g. would not * be continuous despite using a $Q_1$ element). - * - The situation becomes completely hopeless if you consider + * - The situation becomes completely hopeless if you consider * curved boundaries since then the edge midpoint (i.e. the hanging node) * does in general not lie on the mother edge. Consequently, the solution * will not be $H^1$ conforming anyway, regardless of the priority of * the two competing constraints. If the hanging node constraint wins, then - * the solution will be neither conforming, nor have the right boundary + * the solution will be neither conforming, nor have the right boundary * values. * In other words, it is not entirely clear what the "correct" solution would * be. In most cases, it will not matter much: in either case, the error -- 2.39.5