]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Ignore constraints computed for degrees of freedom that are already constrained.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 17 Sep 2010 12:43:00 +0000 (12:43 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 17 Sep 2010 12:43:00 +0000 (12:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@22011 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/include/numerics/vectors.templates.h

index 3105748402a576b33ef6b56bcd0f2724067041d9..aa5103da74227137cb7064e568edb7808ea1e39b 100644 (file)
@@ -707,23 +707,22 @@ class VectorTools
                                      * constrained to, inhomogeneities) is
                                      * kept and nothing happens.
                                      *
-                                     * Please note that when combining
-                                     * adaptively refined meshes with hanging
-                                     * node constraints and inhomogeneous
-                                     * Dirichlet boundary conditions within
+                                     * @note When combining adaptively
+                                     * refined meshes with hanging node
+                                     * constraints and boundary conditions
+                                     * like from the current function within
                                      * one ConstraintMatrix object, the
                                      * hanging node constraints should always
-                                     * be set first, and then Dirichlet
-                                     * boundary conditions should be
-                                     * interpolated. This makes sure that the
-                                     * discretization remains H<sup>1</sup>
-                                     * conforming as is needed e.g. for the
-                                     * Laplace equation in 3D, as hanging
-                                     * nodes on the boundary should be still
-                                     * set to the weighted average of
-                                     * neighbors, and not the actual
-                                     * Dirichlet value.
- *                                   *
+                                     * be set first, and then the boundary
+                                     * conditions since boundary conditions
+                                     * are not set in the second operation on
+                                     * degrees of freedom that are already
+                                     * constrained. This makes sure that the
+                                     * discretization remains conforming as
+                                     * is needed. See the discussion on
+                                     * conflicting constraints in the module
+                                     * on @ref constraints .
+                                     *
                                      * The parameter @p boundary_component
                                      * corresponds to the number @p
                                      * boundary_indicator of the face.  255
@@ -919,22 +918,21 @@ class VectorTools
                                      * inhomogeneities) is kept and nothing
                                      * happens.
                                      *
-                                     * Please note that when combining
-                                     * adaptively refined meshes with hanging
-                                     * node constraints and inhomogeneous
-                                     * Dirichlet boundary conditions within
+                                     * @note When combining adaptively
+                                     * refined meshes with hanging node
+                                     * constraints and boundary conditions
+                                     * like from the current function within
                                      * one ConstraintMatrix object, the
                                      * hanging node constraints should always
-                                     * be set first, and then Dirichlet
-                                     * boundary conditions should be
-                                     * interpolated. This makes sure that the
-                                     * discretization remains H<sup>1</sup>
-                                     * conforming as is needed e.g. for the
-                                     * Laplace equation in 3D, as hanging
-                                     * nodes on the boundary should be still
-                                     * set to the weighted average of
-                                     * neighbors, and not the actual
-                                     * Dirichlet value.
+                                     * be set first, and then the boundary
+                                     * conditions since boundary conditions
+                                     * are not set in the second operation on
+                                     * degrees of freedom that are already
+                                     * constrained. This makes sure that the
+                                     * discretization remains conforming as
+                                     * is needed. See the discussion on
+                                     * conflicting constraints in the module
+                                     * on @ref constraints .
                                      *
                                      * If @p component_mapping is empty, it
                                      * is assumed that the number of
@@ -1006,7 +1004,9 @@ class VectorTools
                                      * first, and then completed by hanging
                                      * node constraints, in order to make sure
                                      * that the discretization remains
-                                     * consistent.
+                                     * consistent. See the discussion on
+                                     * conflicting constraints in the
+                                     * module on @ref constraints .
                                      *
                                      * This function is explecitly written to
                                      * use with the FE_Nedelec elements. Thus
@@ -1177,6 +1177,22 @@ class VectorTools
                                      * compute the normal vector $\vec n$ at
                                      * the boundary points.
                                      *
+                                     * @note When combining adaptively
+                                     * refined meshes with hanging node
+                                     * constraints and boundary conditions
+                                     * like from the current function within
+                                     * one ConstraintMatrix object, the
+                                     * hanging node constraints should always
+                                     * be set first, and then the boundary
+                                     * conditions since boundary conditions
+                                     * are not set in the second operation on
+                                     * degrees of freedom that are already
+                                     * constrained. This makes sure that the
+                                     * discretization remains conforming as
+                                     * is needed. See the discussion on
+                                     * conflicting constraints in the module
+                                     * on @ref constraints .
+                                     *
                                      *
                                      * <h4>Computing constraints in 2d</h4>
                                      *
index c975d44e566c97f28435dc8edd5bee0e9e4ef9da..54fc7c1213091b87e5f7a5beebecfe2a2590a1af 100644 (file)
@@ -2122,6 +2122,10 @@ namespace internal
                                      * indices, and $\vec n$ by the
                                      * vector specified as the second
                                      * argument.
+                                     *
+                                     * The function does not add constraints
+                                     * if a degree of freedom is already
+                                     * constrained in the constraints object.
                                      */
     template <int dim>
     void
@@ -2129,6 +2133,7 @@ namespace internal
                    const Tensor<1,dim>       &constraining_vector,
                    ConstraintMatrix          &constraints)
     {
+      
                                       // choose the DoF that has the
                                       // largest component in the
                                       // constraining_vector as the
@@ -2152,23 +2157,29 @@ namespace internal
          {
            if (std::fabs(constraining_vector[0]) > std::fabs(constraining_vector[1]))
              {
-               constraints.add_line (dof_indices.dof_indices[0]);
+               if (!constraints.is_constrained(dof_indices.dof_indices[0]))
+                 {
+                   constraints.add_line (dof_indices.dof_indices[0]);
 
-               if (std::fabs (constraining_vector[1]/constraining_vector[0])
-                   > std::numeric_limits<double>::epsilon())
-                 constraints.add_entry (dof_indices.dof_indices[0],
-                                        dof_indices.dof_indices[1],
-                                        -constraining_vector[1]/constraining_vector[0]);
+                   if (std::fabs (constraining_vector[1]/constraining_vector[0])
+                       > std::numeric_limits<double>::epsilon())
+                     constraints.add_entry (dof_indices.dof_indices[0],
+                                            dof_indices.dof_indices[1],
+                                            -constraining_vector[1]/constraining_vector[0]);
+                 }
              }
            else
              {
-               constraints.add_line (dof_indices.dof_indices[1]);
+               if (!constraints.is_constrained(dof_indices.dof_indices[1]))
+                 {
+                   constraints.add_line (dof_indices.dof_indices[1]);
 
-               if (std::fabs (constraining_vector[0]/constraining_vector[1])
-                   > std::numeric_limits<double>::epsilon())
-                 constraints.add_entry (dof_indices.dof_indices[1],
-                                        dof_indices.dof_indices[0],
-                                        -constraining_vector[0]/constraining_vector[1]);
+                   if (std::fabs (constraining_vector[0]/constraining_vector[1])
+                       > std::numeric_limits<double>::epsilon())
+                     constraints.add_entry (dof_indices.dof_indices[1],
+                                            dof_indices.dof_indices[0],
+                                            -constraining_vector[0]/constraining_vector[1]);
+                 }
              }
            break;
          }
@@ -2179,54 +2190,63 @@ namespace internal
                &&
                (std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[2])))
              {
-               constraints.add_line (dof_indices.dof_indices[0]);
-
-               if (std::fabs (constraining_vector[1]/constraining_vector[0])
-                   > std::numeric_limits<double>::epsilon())
-                 constraints.add_entry (dof_indices.dof_indices[0],
-                                        dof_indices.dof_indices[1],
-                                        -constraining_vector[1]/constraining_vector[0]);
-
-               if (std::fabs (constraining_vector[2]/constraining_vector[0])
-                   > std::numeric_limits<double>::epsilon())
-                 constraints.add_entry (dof_indices.dof_indices[0],
-                                        dof_indices.dof_indices[2],
-                                        -constraining_vector[2]/constraining_vector[0]);
+               if (!constraints.is_constrained(dof_indices.dof_indices[0]))
+                 {
+                   constraints.add_line (dof_indices.dof_indices[0]);
+
+                   if (std::fabs (constraining_vector[1]/constraining_vector[0])
+                       > std::numeric_limits<double>::epsilon())
+                     constraints.add_entry (dof_indices.dof_indices[0],
+                                            dof_indices.dof_indices[1],
+                                            -constraining_vector[1]/constraining_vector[0]);
+
+                   if (std::fabs (constraining_vector[2]/constraining_vector[0])
+                       > std::numeric_limits<double>::epsilon())
+                     constraints.add_entry (dof_indices.dof_indices[0],
+                                            dof_indices.dof_indices[2],
+                                            -constraining_vector[2]/constraining_vector[0]);
+                 }
              }
            else
              if ((std::fabs(constraining_vector[1]) >= std::fabs(constraining_vector[0]))
                  &&
                  (std::fabs(constraining_vector[1]) >= std::fabs(constraining_vector[2])))
                {
-                 constraints.add_line (dof_indices.dof_indices[1]);
-
-                 if (std::fabs (constraining_vector[0]/constraining_vector[1])
-                     > std::numeric_limits<double>::epsilon())
-                   constraints.add_entry (dof_indices.dof_indices[1],
-                                          dof_indices.dof_indices[0],
-                                          -constraining_vector[0]/constraining_vector[1]);
-
-                 if (std::fabs (constraining_vector[2]/constraining_vector[1])
-                     > std::numeric_limits<double>::epsilon())
-                   constraints.add_entry (dof_indices.dof_indices[1],
-                                          dof_indices.dof_indices[2],
-                                          -constraining_vector[2]/constraining_vector[1]);
+                 if (!constraints.is_constrained(dof_indices.dof_indices[1]))
+                   {
+                     constraints.add_line (dof_indices.dof_indices[1]);
+
+                     if (std::fabs (constraining_vector[0]/constraining_vector[1])
+                         > std::numeric_limits<double>::epsilon())
+                       constraints.add_entry (dof_indices.dof_indices[1],
+                                              dof_indices.dof_indices[0],
+                                              -constraining_vector[0]/constraining_vector[1]);
+
+                     if (std::fabs (constraining_vector[2]/constraining_vector[1])
+                         > std::numeric_limits<double>::epsilon())
+                       constraints.add_entry (dof_indices.dof_indices[1],
+                                              dof_indices.dof_indices[2],
+                                              -constraining_vector[2]/constraining_vector[1]);
+                   }
                }
              else
                {
-                 constraints.add_line (dof_indices.dof_indices[2]);
-
-                 if (std::fabs (constraining_vector[0]/constraining_vector[2])
-                     > std::numeric_limits<double>::epsilon())
-                   constraints.add_entry (dof_indices.dof_indices[2],
-                                          dof_indices.dof_indices[0],
-                                          -constraining_vector[0]/constraining_vector[2]);
-
-                 if (std::fabs (constraining_vector[1]/constraining_vector[2])
-                     > std::numeric_limits<double>::epsilon())
-                   constraints.add_entry (dof_indices.dof_indices[2],
-                                          dof_indices.dof_indices[1],
-                                          -constraining_vector[1]/constraining_vector[2]);
+                 if (!constraints.is_constrained(dof_indices.dof_indices[2]))
+                   {
+                     constraints.add_line (dof_indices.dof_indices[2]);
+
+                     if (std::fabs (constraining_vector[0]/constraining_vector[2])
+                         > std::numeric_limits<double>::epsilon())
+                       constraints.add_entry (dof_indices.dof_indices[2],
+                                              dof_indices.dof_indices[0],
+                                              -constraining_vector[0]/constraining_vector[2]);
+
+                     if (std::fabs (constraining_vector[1]/constraining_vector[2])
+                         > std::numeric_limits<double>::epsilon())
+                       constraints.add_entry (dof_indices.dof_indices[2],
+                                              dof_indices.dof_indices[1],
+                                              -constraining_vector[1]/constraining_vector[2]);
+                   }
                }
 
            break;
@@ -3788,11 +3808,16 @@ compute_no_normal_flux_constraints (const DH<dim,spacedim>         &dof_handler,
                                             // constrained. enter
                                             // this into the
                                             // constraint matrix
+                                            //
+                                            // ignore dofs already
+                                            // constrained
            for (unsigned int i=0; i<dim; ++i)
-             {
-               constraints.add_line (same_dof_range[0]->first.dof_indices[i]);
-                                                // no add_entries here
-             }
+             if (!constraints.is_constrained (same_dof_range[0]
+                                              ->first.dof_indices[i]))
+               {
+                 constraints.add_line (same_dof_range[0]->first.dof_indices[i]);
+                                                  // no add_entries here
+               }
 
            break;
          }

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.