From: bangerth Date: Wed, 13 Dec 2006 00:26:41 +0000 (+0000) Subject: Reimplement the algorithm by which we decide which will be the master dofs. This... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6883ba9f0ddcb45516013b7c9c153fd6be9b466b;p=dealii-svn.git Reimplement the algorithm by which we decide which will be the master dofs. This should now work, and indeed hp/crash_12 now passes. git-svn-id: https://svn.dealii.org/trunk@14228 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 5d447075bd..3aeaad1958 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -1564,6 +1564,117 @@ namespace internal { namespace { + inline + bool + check_master_dof_list (const FullMatrix &face_interpolation_matrix, + const std::vector &master_dof_list) + { + const unsigned int N = master_dof_list.size(); + + FullMatrix tmp (N,N); + for (unsigned int i=0; i p(N); + for (unsigned int i=0; i max) + { + max = std::fabs(tmp(i,j)); + r = i; + } + } + // check whether the + // pivot is too small. if + // that is the case, then + // the matrix is singular + // and we shouldn't use + // this set of master + // dofs + if (max < 1.e-12*typical_diagonal_element) + return false; + + // row interchange + if (r>j) + { + for (unsigned int k=0; k void select_master_dofs_for_face_restriction (const FiniteElement &fe1, const FiniteElement &fe2, + const FullMatrix &face_interpolation_matrix, std::vector &master_dof_mask) { Assert (fe1.dofs_per_face >= fe2.dofs_per_face, @@ -1601,9 +1733,6 @@ namespace internal Assert (fe2.dofs_per_quad <= fe1.dofs_per_quad, ExcInternalError()); - // initialize all with false - std::fill (master_dof_mask.begin(), master_dof_mask.end(), false); - // the idea here is to designate as // many DoFs in fe1 per object // (vertex, line, quad) as master as @@ -1612,33 +1741,125 @@ namespace internal // to avoid the 'unsigned int < 0 is // always false warning for the cases // at the bottom in 1d and 2d) + // + // as mentioned in the paper, + // it is not always easy to + // find a set of master dofs + // that produces an + // invertible matrix. to this + // end, we check in each step + // whether the matrix is + // still invertible. we have + // never had a problem with + // vertex and line dofs, so + // we only check inside an + // Assert statement, without + // providing any code that + // can cope with the + // situation if this goes + // wrong. + // + // on the other hand, we did + // have trouble with adding + // more quad dofs when Q3 and + // Q4 elements meet at a + // refined face in 3d (see + // the hp/crash_12 test that + // tests that we can do + // exactly this, and failed + // before we had code to + // compensate for this + // case). therefore, for quad + // dofs (which appear only in + // 3d anyway), we have a + // piece of code that can add + // dofs until we have finally + // found a suitable set + std::vector master_dof_list; unsigned int index = 0; - for (int v=0; v(GeometryInfo::vertices_per_face); ++v) + for (int v=0; + v(GeometryInfo::vertices_per_face); + ++v) { for (unsigned int i=0; i(GeometryInfo::lines_per_face); ++l) + for (int l=0; + l(GeometryInfo::lines_per_face); + ++l) { for (unsigned int i=0; i(GeometryInfo::quads_per_face); ++q) + for (int q=0; + q(GeometryInfo::quads_per_face); + ++q) { - for (unsigned int i=0; i(fe2.dofs_per_face), + Assert (master_dof_list.size() == fe2.dofs_per_face, ExcInternalError()); + + // finally copy the list into the + // mask + std::fill (master_dof_mask.begin(), master_dof_mask.end(), false); + for (std::vector::const_iterator i=master_dof_list.begin(); + i!=master_dof_list.end(); ++i) + master_dof_mask[*i] = true; } @@ -2720,9 +2941,13 @@ namespace internal // master and slave // components. invert the // master component +//TODO[WB]: cache this mask as well std::vector master_dof_mask (cell->get_fe().dofs_per_face); select_master_dofs_for_face_restriction (cell->get_fe(), dominating_fe, + (*face_interpolation_matrices + [dominating_fe_index] + [cell->active_fe_index()]), master_dof_mask); ensure_existence_of_split_face_matrix