From: Daniel Arndt Date: Mon, 6 Jul 2015 09:18:34 +0000 (+0200) Subject: Fix smaller errors for rotated periodicity constraints and allow for constraints... X-Git-Tag: v8.3.0-rc1~12^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1ee1c50e7784243bfd1c0b64ed114d5065bf36b2;p=dealii.git Fix smaller errors for rotated periodicity constraints and allow for constraints on the same DoF in this case --- diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 68adca2558..5c580c7e0e 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -1781,7 +1781,7 @@ namespace DoFTools cell_to_rotated_face_index[cell_index] = i; } - // loop over all dofs on face 2 and constrain them again the ones on face 1 + // loop over all dofs on face 2 and constrain them against the ones on face 1 for (unsigned int i=0; i > *constraint_entries + = constraint_matrix.get_constraint_entries(dofs_1[j]); + if (constraint_entries->size()==1 && (*constraint_entries)[0].first == dofs_2[i]) + { + if ((is_identity_constrained && std::abs((*constraint_entries)[0].second-1) > eps) || + (is_inverse_constrained && std::abs((*constraint_entries)[0].second+1) > eps)) + { + //this pair of constraints means that both dofs have to be constrained to 0. + constraint_matrix.add_line(dofs_2[i]); + } + } + else + enter_constraint = true; + } + + if (enter_constraint) + { + constraint_matrix.add_line(dofs_2[i]); + constraint_matrix.add_entry(dofs_2[i], dofs_1[j], is_identity_constrained?1.0:-1.0); + } + } + else + { + // this is just a regular constraint. enter it piece by piece + constraint_matrix.add_line(dofs_2[i]); + for (unsigned int jj=0; jj quadrature (fe.get_unit_face_support_points()); - FEFaceValues fe_face_values (fe, quadrature, update_q_points); // have an array that stores the location of each vector-dof tuple // we want to rotate. @@ -1927,6 +2004,7 @@ namespace DoFTools // find corresponding other components of vector DoFTuple vector_dofs; vector_dofs[0] = i; + unsigned int n_found = 1; Assert(*comp_it + spacedim <= fe.n_components(), ExcMessage("Error: the finite element does not have enough components " @@ -1943,7 +2021,9 @@ namespace DoFTools vector_dofs[fe.face_system_to_component_index(k).first - first_vector_component] = k; - break; + n_found++; + if (n_found==dim) + break; } // ... and rotate all dofs belonging to vector valued diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index b157fcae93..9d289634d0 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2530,7 +2530,7 @@ next_cell: if (matrix.m() == spacedim) for (int i = 0; i < spacedim; ++i) for (int j = 0; j < spacedim; ++j) - distance(i) = matrix(i,j) * point1(j); + distance(i) += matrix(i,j) * point1(j); else distance = point1; diff --git a/tests/mpi/periodicity_02.cc b/tests/mpi/periodicity_02.cc index ef207ba153..1750ee28a3 100644 --- a/tests/mpi/periodicity_02.cc +++ b/tests/mpi/periodicity_02.cc @@ -325,10 +325,9 @@ namespace Step22 DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs); relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u)); relevant_partitioning.push_back(locally_relevant_dofs.get_view(n_u, n_u+n_p)); - } - - { + constraints.clear (); + constraints.reinit(locally_relevant_dofs); FEValuesExtractors::Vector velocities(0); FEValuesExtractors::Scalar pressure(dim);