From ac85fa24a76a9a1c1f6a73c15ec1655bdc71ccfe Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 28 Aug 2015 07:29:34 -0500 Subject: [PATCH] Fix Issue 1128 and enable periodicity_03 --- source/dofs/dof_tools_constraints.cc | 22 ++++++--- tests/mpi/periodicity_03.cc | 73 +++++++++------------------- 2 files changed, 38 insertions(+), 57 deletions(-) diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 8b1021b092..e4fd5a47e1 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -2098,10 +2098,13 @@ namespace DoFTools const unsigned int n_dofs_per_face = face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face; - Assert(matrix.m() == 0 || matrix.m() == n_dofs_per_face || - matrix.m() == (int)spacedim, - ExcMessage ("matrix must have either size 0 or spacedim or the " - "size must be equal to the # of DoFs on the face")); + Assert(matrix.m() == 0 + || (first_vector_components.empty() && matrix.m() == n_dofs_per_face) + || (!first_vector_components.empty() && matrix.m() == (int)spacedim), + ExcMessage ("The matrix must have either size 0 or spacedim " + "(if first_vector_components is nonempty) " + "or the size must be equal to the # of DoFs on the face " + "(if first_vector_components is empty).")); } if (!face_2->has_children()) @@ -2110,10 +2113,13 @@ namespace DoFTools const unsigned int n_dofs_per_face = face_2->get_fe(face_2->nth_active_fe_index(0)).dofs_per_face; - Assert(matrix.m() == 0 || matrix.m() == n_dofs_per_face || - matrix.m() == (int)spacedim, - ExcMessage ("matrix must have either size 0 or spacedim or the " - "size must be equal to the # of DoFs on the face")); + Assert(matrix.m() == 0 + || (first_vector_components.empty() && matrix.m() == n_dofs_per_face) + || (!first_vector_components.empty() && matrix.m() == (int)spacedim), + ExcMessage ("The matrix must have either size 0 or spacedim " + "(if first_vector_components is nonempty) " + "or the size must be equal to the # of DoFs on the face " + "(if first_vector_components is empty).")); } #endif diff --git a/tests/mpi/periodicity_03.cc b/tests/mpi/periodicity_03.cc index 2d66b82dcc..e4c0d34f28 100644 --- a/tests/mpi/periodicity_03.cc +++ b/tests/mpi/periodicity_03.cc @@ -70,7 +70,7 @@ namespace Step22 void solve (); void get_point_value (const Point point, const int proc, Vector &value) const; - void check_periodicity(const unsigned int cycle) const; + void check_periodicity() const; void output_results (const unsigned int refinement_cycle) const; void refine_mesh (); @@ -78,7 +78,6 @@ namespace Step22 MPI_Comm mpi_communicator; - HyperShellBoundary boundary; parallel::distributed::Triangulation triangulation; FESystem fe; DoFHandler dof_handler; @@ -97,7 +96,6 @@ namespace Step22 FullMatrix rot_matrix; Tensor<1,dim> offset; - MappingQ mapping; }; @@ -142,9 +140,8 @@ namespace Step22 (TrilinosWrappers::MPI::Vector &dst, const TrilinosWrappers::MPI::Vector &src) const { - SolverControl solver_control (10*src.size(), 1e-6*src.l2_norm(), false, false); - TrilinosWrappers::SolverCG cg (solver_control, - TrilinosWrappers::SolverCG::AdditionalData()); + SolverControl solver_control (src.size(), 1e-6*src.l2_norm(), false, false); + SolverCG cg (solver_control); tmp = 0.; cg.solve (*matrix, tmp, src, *preconditioner); @@ -161,7 +158,6 @@ namespace Step22 const InverseMatrix &A_inverse, const IndexSet &owned_pres, - const IndexSet &relevant_pres, const MPI_Comm &mpi_communicator); void vmult (TrilinosWrappers::MPI::Vector &dst, @@ -182,7 +178,6 @@ namespace Step22 const InverseMatrix &A_inverse, const IndexSet &owned_vel, - const IndexSet &relevant_vel, const MPI_Comm &mpi_communicator) : system_matrix (&system_matrix), @@ -223,8 +218,7 @@ namespace Step22 std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), - rot_matrix(dim), - mapping(1) + rot_matrix(dim) { rot_matrix[0][dim-1] = -1.; rot_matrix[dim-1][0] = 1.; @@ -262,10 +256,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); @@ -281,34 +274,31 @@ namespace Step22 GridTools::collect_periodic_faces (dof_handler, 4, 0, 1, periodicity_vector, offset, - rot_matrix, first_vector_components); + rot_matrix); DoFTools::make_periodicity_constraints > - (periodicity_vector, constraints, fe.component_mask(velocities)); + (periodicity_vector, constraints, fe.component_mask(velocities), + first_vector_components); - VectorTools::interpolate_boundary_values (mapping, - dof_handler, + VectorTools::interpolate_boundary_values (dof_handler, 1, ZeroFunction(dim+1), constraints, fe.component_mask(velocities)); - VectorTools::interpolate_boundary_values (mapping, - dof_handler, + VectorTools::interpolate_boundary_values (dof_handler, 2, ZeroFunction(dim+1), constraints, fe.component_mask(velocities)); - VectorTools::interpolate_boundary_values (mapping, - dof_handler, + VectorTools::interpolate_boundary_values (dof_handler, 3, ZeroFunction(dim+1), constraints, fe.component_mask(velocities)); - VectorTools::interpolate_boundary_values (mapping, - dof_handler, + VectorTools::interpolate_boundary_values (dof_handler, 5, ZeroFunction(dim+1), constraints, @@ -346,8 +336,7 @@ namespace Step22 QGauss quadrature_formula(degree+2); - FEValues fe_values (mapping, - fe, quadrature_formula, + FEValues fe_values (fe, quadrature_formula, update_values | update_quadrature_points | update_JxW_values | @@ -392,7 +381,7 @@ namespace Step22 { for (unsigned int j=0; j<=i; ++j) { - local_matrix(i,j) += (symgrad_phi_u[i] * symgrad_phi_u[j] + local_matrix(i,j) += (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) - div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] + phi_p[i] * phi_p[j]) @@ -422,8 +411,6 @@ namespace Step22 system_matrix.compress (VectorOperation::add); system_rhs.compress (VectorOperation::add); - - pcout << " Computing preconditioner..." << std::endl << std::flush; } @@ -454,22 +441,16 @@ namespace Step22 SchurComplement schur_complement (system_matrix, A_inverse, - owned_partitioning[0], relevant_partitioning[0], + owned_partitioning[0], mpi_communicator); SolverControl solver_control (solution.block(1).size(), 1e-6*schur_rhs.l2_norm(), false, false); -// TrilinosWrappers::SolverCG cg (solver_control, -// TrilinosWrappers::SolverCG::AdditionalData(true)); SolverCG cg(solver_control); TrilinosWrappers::PreconditionAMG preconditioner; preconditioner.initialize (system_matrix.block(1,1)); - InverseMatrix - m_inverse (system_matrix.block(1,1), preconditioner, - owned_partitioning[1], mpi_communicator); cg.solve (schur_complement, tmp.block(1), @@ -499,10 +480,10 @@ namespace Step22 try { const typename DoFHandler::active_cell_iterator cell - = GridTools::find_active_cell_around_point (mapping, dof_handler, point).first; + = GridTools::find_active_cell_around_point (dof_handler, point); if (cell->is_locally_owned()) - VectorTools::point_value (mapping, dof_handler, solution, + VectorTools::point_value (dof_handler, solution, point, value); } catch (GridTools::ExcPointNotFound &p) @@ -524,17 +505,17 @@ namespace Step22 } template - void StokesProblem::check_periodicity (const unsigned int cycle) const + void StokesProblem::check_periodicity () const { AssertThrow(false, ExcNotImplemented()); } template <> - void StokesProblem<3>::check_periodicity (const unsigned int cycle) const + void StokesProblem<3>::check_periodicity () const { const unsigned int dim = 3; Quadrature q_dummy (fe.base_element(0).get_unit_face_support_points()); - FEFaceValues fe_face_values (mapping,fe, q_dummy, update_quadrature_points); + FEFaceValues fe_face_values (fe, q_dummy, update_quadrature_points); DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); DoFHandler::active_cell_iterator endc= dof_handler.end(); @@ -771,12 +752,10 @@ namespace Step22 void StokesProblem::refine_mesh () { - Vector estimated_error_per_cell (triangulation.n_active_cells()); FEValuesExtractors::Scalar pressure(dim); - KellyErrorEstimator::estimate (mapping, - dof_handler, + KellyErrorEstimator::estimate (dof_handler, QGauss(degree+1), typename FunctionMap::type(), solution, @@ -800,12 +779,8 @@ namespace Step22 std::vector::cell_iterator> > periodicity_vector; - std::vector first_vector_components; - first_vector_components.push_back(0); - GridTools::collect_periodic_faces - (triangulation, 4, 0, 1, periodicity_vector, offset, - rot_matrix, first_vector_components); + (triangulation, 4, 0, 1, periodicity_vector, offset, rot_matrix); triangulation.add_periodicity(periodicity_vector); triangulation.refine_global (1); @@ -827,7 +802,7 @@ namespace Step22 solve (); // output_results (refinement_cycle); - check_periodicity(refinement_cycle); + check_periodicity(); pcout << std::endl; } -- 2.39.5