From 3353f0b79071fcd177ce3822414defef5b0fb12d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 17 Sep 2010 12:43:00 +0000 Subject: [PATCH] Ignore constraints computed for degrees of freedom that are already constrained. git-svn-id: https://svn.dealii.org/trunk@22011 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/numerics/vectors.h | 76 ++++++---- .../include/numerics/vectors.templates.h | 135 +++++++++++------- 2 files changed, 126 insertions(+), 85 deletions(-) diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 3105748402..aa5103da74 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -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 H1 - * 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 H1 - * 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 . + * * *

Computing constraints in 2d

* diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index c975d44e56..54fc7c1213 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -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 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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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::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 &dof_handler, // constrained. enter // this into the // constraint matrix + // + // ignore dofs already + // constrained for (unsigned int i=0; ifirst.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; } -- 2.39.5