From 5fc8cf2cea45eff1913408a7be3eb1254dca6840 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 24 Aug 2011 13:52:08 +0000 Subject: [PATCH] Check in a fix provided by Martin Kronbichler for the problem that different processors disagree as to which DoF to constrain, leading to an incredibly difficult to find bug. This happened with step-32, so is a common enough case trying to work around. git-svn-id: https://svn.dealii.org/trunk@24176 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/numerics/vectors.templates.h | 68 +++++++++++++++++-- 1 file changed, 63 insertions(+), 5 deletions(-) diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h index 62292846b4..b75d91c0e7 100644 --- a/deal.II/include/deal.II/numerics/vectors.templates.h +++ b/deal.II/include/deal.II/numerics/vectors.templates.h @@ -2213,11 +2213,55 @@ namespace internal // can also deal with cases like // x0 = 0 // if necessary + // + // there is a problem if we have a + // normal vector of the form + // (a,a,small) or (a,a,a). Depending on + // round-off we may choose the first or + // second component (or third, in the + // latter case) as the largest one, and + // depending on our choice one or + // another degree of freedom will be + // constrained. On a single processor + // this is not much of a problem, but + // it's a nightmare when we run in + // parallel and two processors disagree + // on which DoF should be + // constrained. This led to an + // incredibly difficult to find bug in + // step-32 when running in parallel + // with 9 or more processors. + // + // in practice, such normal vectors of + // the form (a,a,small) or (a,a,a) + // happen not infrequently since they + // lie on the diagonals where vertices + // frequently happen to land upon mesh + // refinement if one starts from a + // symmetric and regular body. we work + // around this problem in the following + // way: if we have a normal vector of + // the form (a,b) (similarly algorithm + // in 3d), we choose 'a' as the largest + // coefficient not if a>b but if + // a>b+1e-10. this shifts the problem + // away from the frequently visited + // diagonal to a line that is off the + // diagonal. there will of course be + // problems where the exact values of a + // and b differ by exactly 1e-10 and we + // get into the same instability, but + // from a practical viewpoint such + // problems should be much rarer. in + // particular, meshes have to be very + // very fine for a vertex to land on + // this line if the original body had a + // vertex on the diagonal as well switch (dim) { case 2: { - if (std::fabs(constraining_vector[0]) > std::fabs(constraining_vector[1])) + if (std::fabs(constraining_vector[0]) > std::fabs(constraining_vector[1]) + 1e-10) { if (!constraints.is_constrained(dof_indices.dof_indices[0]) && @@ -2252,9 +2296,9 @@ namespace internal case 3: { - if ((std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[1])) + if ((std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[1])+1e-10) && - (std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[2]))) + (std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[2])+2e-10)) { if (!constraints.is_constrained(dof_indices.dof_indices[0]) && @@ -2276,9 +2320,9 @@ namespace internal } } else - if ((std::fabs(constraining_vector[1]) >= std::fabs(constraining_vector[0])) + if ((std::fabs(constraining_vector[1])+1e-10 >= std::fabs(constraining_vector[0])) && - (std::fabs(constraining_vector[1]) >= std::fabs(constraining_vector[2]))) + (std::fabs(constraining_vector[1]) >= std::fabs(constraining_vector[2])+1e-10)) { if (!constraints.is_constrained(dof_indices.dof_indices[1]) && @@ -4419,6 +4463,12 @@ compute_no_normal_flux_constraints (const DH &dof_handler, normal += x->second.first; normal /= normal.norm(); + // normalize again + for (unsigned int d=0; d & @@ -4719,6 +4769,14 @@ compute_no_normal_flux_constraints (const DH &dof_handler, internal::VectorTools:: compute_orthonormal_vectors (average_tangent, constraining_normals); + // normalize again + for (unsigned int e=0; e