From 8a8f14745a02f437d8de92a30b85549716dcfb43 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 6 Jan 2010 10:34:24 +0000 Subject: [PATCH] Use homogenous boundary values, to avoid an assertion in the constraint matrix. git-svn-id: https://svn.dealii.org/trunk@20309 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-42/step-42.cc | 337 ++++++++++++++-------------- 1 file changed, 168 insertions(+), 169 deletions(-) diff --git a/deal.II/examples/step-42/step-42.cc b/deal.II/examples/step-42/step-42.cc index 794db35448..3d97a248d6 100644 --- a/deal.II/examples/step-42/step-42.cc +++ b/deal.II/examples/step-42/step-42.cc @@ -12,8 +12,8 @@ // @sect3{Include files} - - // As usual, we start by including + + // As usual, we start by including // some well-known files: #include #include @@ -51,7 +51,7 @@ #include #include - + #include #include @@ -64,39 +64,39 @@ #include #include #include - + #include #include - + using namespace dealii; - + template struct InnerPreconditioner; - + template <> -struct InnerPreconditioner<2> +struct InnerPreconditioner<2> { typedef SparseDirectUMFPACK type; }; - + template <> -struct InnerPreconditioner<3> +struct InnerPreconditioner<3> { typedef SparseILU type; }; template -class StokesProblem +class StokesProblem { public: StokesProblem (const unsigned int degree); void run (); - + private: void setup_dofs (); void assemble_system (); @@ -105,19 +105,19 @@ class StokesProblem void solve_block (); void solve_mg (); - void find_dofs_on_lower_level (std::vector > &lower_dofs, + void find_dofs_on_lower_level (std::vector > &lower_dofs, std::vector > &boundary_dofs); void output_results (const unsigned int refinement_cycle) const; void refine_mesh (); - + const unsigned int degree; - + Triangulation triangulation; FESystem fe; MGDoFHandler dof_handler; ConstraintMatrix constraints; - + BlockSparsityPattern sparsity_pattern; BlockSparseMatrix system_matrix; @@ -130,22 +130,22 @@ class StokesProblem MGLevelObject > mg_interface_matrices; std::vector > mg_dofs_per_component; - + std_cxx1x::shared_ptr::type> A_preconditioner; std::vector::type> > mg_A_preconditioner; }; - + template -class BoundaryValues : public Function +class BoundaryValues : public Function { public: BoundaryValues () : Function(dim+1) {} - + virtual double value (const Point &p, const unsigned int component = 0) const; - virtual void vector_value (const Point &p, + virtual void vector_value (const Point &p, Vector &value) const; }; @@ -153,11 +153,11 @@ class BoundaryValues : public Function template double BoundaryValues::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { Assert (component < this->n_components, ExcIndexRange (component, 0, this->n_components)); - + if (component == 0 && p[0] == 0) return (dim == 2 ? - p[1]*(p[1]-1.) : p[1]*(p[1]-1.) * p[2]*(p[2]-1.)); return 0; @@ -167,7 +167,7 @@ BoundaryValues::value (const Point &p, template void BoundaryValues::vector_value (const Point &p, - Vector &values) const + Vector &values) const { for (unsigned int c=0; cn_components; ++c) values(c) = BoundaryValues::value (p, c); @@ -175,43 +175,43 @@ BoundaryValues::vector_value (const Point &p, - + template -class RightHandSide : public Function +class RightHandSide : public Function { public: RightHandSide () : Function(dim+1) {} - + virtual double value (const Point &p, const unsigned int component = 0) const; - virtual void vector_value (const Point &p, + virtual void vector_value (const Point &p, Vector &value) const; - + }; template double RightHandSide::value (const Point &/*p*/, - const unsigned int /*component*/) const + const unsigned int component) const { - return 0; + return (component == 0 ? 1 : 0); } template void RightHandSide::vector_value (const Point &p, - Vector &values) const + Vector &values) const { for (unsigned int c=0; cn_components; ++c) values(c) = RightHandSide::value (p, c); } - - + + template class InverseMatrix : public Subscriptor { @@ -264,10 +264,10 @@ class BlockSchurPreconditioner : public Subscriptor private: const SmartPointer > system_matrix; - const SmartPointer, + const SmartPointer, PreconditionerMp > > m_inverse; const PreconditionerA &a_preconditioner; - + mutable Vector tmp; }; @@ -286,7 +286,7 @@ BlockSchurPreconditioner::BlockSchurPrecondit {} // Now the interesting function, the multiplication of - // the preconditioner with a BlockVector. + // the preconditioner with a BlockVector. template void BlockSchurPreconditioner::vmult ( BlockVector &dst, @@ -294,13 +294,13 @@ void BlockSchurPreconditioner::vmult ( { // Form u_new = A^{-1} u a_preconditioner.vmult (dst.block(0), src.block(0)); - // Form tmp = - B u_new + p + // Form tmp = - B u_new + p // (SparseMatrix::residual // does precisely this) system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); // Change sign in tmp tmp *= -1; - // Multiply by approximate Schur complement + // Multiply by approximate Schur complement // (i.e. a pressure mass matrix) m_inverse->vmult (dst.block(1), tmp); } @@ -319,7 +319,7 @@ class SchurComplement : public Subscriptor private: const SmartPointer > system_matrix; const SmartPointer, Preconditioner> > A_inverse; - + mutable Vector tmp1, tmp2; }; @@ -347,7 +347,7 @@ void SchurComplement::vmult (Vector &dst, } - + template StokesProblem::StokesProblem (const unsigned int degree) : @@ -359,31 +359,31 @@ StokesProblem::StokesProblem (const unsigned int degree) {} - - - + + + template void StokesProblem::setup_dofs () { A_preconditioner.reset (); mg_A_preconditioner.resize (0); system_matrix.clear (); - - dof_handler.distribute_dofs (fe); + + dof_handler.distribute_dofs (fe); // DoFRenumbering::Cuthill_McKee (dof_handler); std::vector block_component (dim+1,0); block_component[dim] = 1; DoFRenumbering::component_wise (dof_handler, block_component); - + { constraints.clear (); std::vector component_mask (dim+1, true); component_mask[dim] = false; VectorTools::interpolate_boundary_values (dof_handler, 0, - BoundaryValues(), + ZeroFunction(dim+1), //TODO: go back to BoundaryValues constraints, component_mask); DoFTools::make_hanging_node_constraints (dof_handler, @@ -392,10 +392,10 @@ void StokesProblem::setup_dofs () constraints.close (); - + std::vector dofs_per_block (2); DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, - block_component); + block_component); const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1]; @@ -406,10 +406,10 @@ void StokesProblem::setup_dofs () << dof_handler.n_dofs() << " (" << n_u << '+' << n_p << ')' << std::endl; - - - - + + + + { BlockCompressedSimpleSparsityPattern csp (2,2); @@ -417,23 +417,23 @@ void StokesProblem::setup_dofs () csp.block(1,0).reinit (n_p, n_u); csp.block(0,1).reinit (n_u, n_p); csp.block(1,1).reinit (n_p, n_p); - - csp.collect_sizes(); - + + csp.collect_sizes(); + DoFTools::make_sparsity_pattern ( - static_cast&>(dof_handler), + static_cast&>(dof_handler), csp, constraints, false); sparsity_pattern.copy_from (csp); } - - + + system_matrix.reinit (sparsity_pattern); - + solution.reinit (2); solution.block(0).reinit (n_u); solution.block(1).reinit (n_p); solution.collect_sizes (); - + system_rhs.reinit (2); system_rhs.block(0).reinit (n_u); system_rhs.block(1).reinit (n_p); @@ -460,7 +460,7 @@ void StokesProblem::setup_dofs () for (unsigned int level=0; level::setup_dofs () } - + template -void StokesProblem::assemble_system () +void StokesProblem::assemble_system () { system_matrix=0; system_rhs=0; - + QGauss quadrature_formula(degree+2); FEValues fe_values (fe, quadrature_formula, @@ -484,42 +484,42 @@ void StokesProblem::assemble_system () update_quadrature_points | update_JxW_values | update_gradients); - + const unsigned int dofs_per_cell = fe.dofs_per_cell; - + const unsigned int n_q_points = quadrature_formula.size(); FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); Vector local_rhs (dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); - + const RightHandSide right_hand_side; std::vector > rhs_values (n_q_points, Vector(dim+1)); - + const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); - - + + std::vector > phi_grads_u (dofs_per_cell); std::vector div_phi_u (dofs_per_cell); std::vector phi_p (dofs_per_cell); - + typename MGDoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) - { + { fe_values.reinit (cell); local_matrix = 0; local_rhs = 0; - + right_hand_side.vector_value_list(fe_values.get_quadrature_points(), rhs_values); - + for (unsigned int q=0; q::assemble_system () - div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] + phi_p[i] * phi_p[j]) - * fe_values.JxW(q); - + * fe_values.JxW(q); } const unsigned int component_i = fe.system_to_component_index(i).first; - local_rhs(i) += fe_values.shape_value(i,q) * + local_rhs(i) += fe_values.shape_value(i,q) * rhs_values[q](component_i) * fe_values.JxW(q); } } - + cell->get_dof_indices (local_dof_indices); constraints.distribute_local_to_global (local_matrix, local_rhs, - local_dof_indices, + local_dof_indices, system_matrix, system_rhs); } - - + + std::cout << " Computing preconditioner..." << std::endl << std::flush; - + A_preconditioner = std_cxx1x::shared_ptr::type>(new typename InnerPreconditioner::type()); A_preconditioner->initialize (system_matrix.block(0,0), @@ -570,7 +569,7 @@ void StokesProblem::assemble_system () template -void StokesProblem::assemble_multigrid () +void StokesProblem::assemble_multigrid () { QGauss quadrature_formula(degree+2); FEValues fe_values (fe, quadrature_formula, @@ -578,22 +577,22 @@ void StokesProblem::assemble_multigrid () update_quadrature_points | update_JxW_values | update_gradients); - + const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); - + const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); - - + + std::vector > phi_grads_u (dofs_per_cell); std::vector div_phi_u (dofs_per_cell); std::vector phi_p (dofs_per_cell); - + std::vector > interface_dofs; std::vector > boundary_interface_dofs; for (unsigned int level = 0; level::assemble_multigrid () cell = dof_handler.begin(), endc = dof_handler.end(); for (; cell!=endc; ++cell) - { + { // Remember the level of the // current cell. const unsigned int level = cell->level(); @@ -640,7 +639,7 @@ void StokesProblem::assemble_multigrid () // by update flags above. fe_values.reinit (cell); local_matrix = 0; - + for (unsigned int q=0; q::assemble_multigrid () - div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] + phi_p[i] * phi_p[j]) - * fe_values.JxW(q); + * fe_values.JxW(q); } - + cell->get_mg_dof_indices (local_dof_indices); // for (unsigned int i=0; i::assemble_multigrid () for (unsigned int i=0; i::assemble_multigrid () } } - + template -void StokesProblem::solve_block () +void StokesProblem::solve_block () { SparseMatrix pressure_mass_matrix; pressure_mass_matrix.reinit(sparsity_pattern.block(1,1)); @@ -706,14 +705,14 @@ void StokesProblem::solve_block () system_matrix.block(1,1) = 0; SparseILU pmass_preconditioner; - pmass_preconditioner.initialize (pressure_mass_matrix, + pmass_preconditioner.initialize (pressure_mass_matrix, SparseILU::AdditionalData()); InverseMatrix,SparseILU > m_inverse (pressure_mass_matrix, pmass_preconditioner); BlockSchurPreconditioner::type, - SparseILU > + SparseILU > preconditioner (system_matrix, m_inverse, *A_preconditioner); SolverControl solver_control (system_matrix.m(), @@ -740,14 +739,14 @@ template class SchurComplementSmoother { public: - struct AdditionalData + struct AdditionalData { const InnerPreconditioner *A_preconditioner; }; void initialize (const BlockSparseMatrix &system_matrix, const AdditionalData &data); - + void vmult (BlockVector &dst, const BlockVector &src) const; @@ -755,7 +754,7 @@ class SchurComplementSmoother const BlockVector &src) const; void clear (); - + private: SmartPointer > system_matrix; SmartPointer A_preconditioner; @@ -784,51 +783,51 @@ vmult (BlockVector &dst, const InverseMatrix,InnerPreconditioner> A_inverse (system_matrix->block(0,0), *A_preconditioner); Vector tmp (dst.block(0).size()); - - + + { Vector schur_rhs (dst.block(1).size()); A_inverse.vmult (tmp, src.block(0)); system_matrix->block(1,0).vmult (schur_rhs, tmp); schur_rhs -= src.block(1); - + SchurComplement schur_complement (*system_matrix, A_inverse); - + // The usual control structures for // the solver call are created... SolverControl solver_control (dst.block(1).size(), 1e-6*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); - - - + + + SparseILU preconditioner; - preconditioner.initialize (system_matrix->block(1,1), + preconditioner.initialize (system_matrix->block(1,1), SparseILU::AdditionalData()); - + InverseMatrix,SparseILU > m_inverse (system_matrix->block(1,1), preconditioner); - - + + cg.solve (schur_complement, dst.block(1), schur_rhs, m_inverse); - + // no constraints to be taken care of here - + std::cout << " " << solver_control.last_step() << " CG Schur complement iterations in smoother" << std::flush - << std::endl; + << std::endl; } - - + + { system_matrix->block(0,1).vmult (tmp, dst.block(1)); tmp *= -1; tmp += src.block(0); - + A_inverse.vmult (dst.block(0), tmp); // no constraints here either @@ -856,7 +855,7 @@ Tvmult (BlockVector &, template -void StokesProblem::solve_mg () +void StokesProblem::solve_mg () { assemble_multigrid (); typedef PreconditionMG, MGTransferPrebuilt > > @@ -881,7 +880,7 @@ void StokesProblem::solve_mg () typedef SchurComplementSmoother::type> Smoother; - + MGSmootherPrecondition, Smoother, BlockVector > @@ -930,58 +929,58 @@ void StokesProblem::solve_mg () template -void StokesProblem::solve_schur () +void StokesProblem::solve_schur () { const InverseMatrix, typename InnerPreconditioner::type> A_inverse (system_matrix.block(0,0), *A_preconditioner); Vector tmp (solution.block(0).size()); - - + + { Vector schur_rhs (solution.block(1).size()); A_inverse.vmult (tmp, system_rhs.block(0)); system_matrix.block(1,0).vmult (schur_rhs, tmp); schur_rhs -= system_rhs.block(1); - + SchurComplement::type> schur_complement (system_matrix, A_inverse); - + // The usual control structures for // the solver call are created... SolverControl solver_control (solution.block(1).size(), 1e-6*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); - - - + + + SparseILU preconditioner; - preconditioner.initialize (system_matrix.block(1,1), + preconditioner.initialize (system_matrix.block(1,1), SparseILU::AdditionalData()); - + InverseMatrix,SparseILU > m_inverse (system_matrix.block(1,1), preconditioner); - - + + cg.solve (schur_complement, solution.block(1), schur_rhs, m_inverse); - - + + constraints.distribute (solution); - + std::cout << " " << solver_control.last_step() << " outer CG Schur complement iterations for pressure" << std::flush - << std::endl; + << std::endl; } - - + + { system_matrix.block(0,1).vmult (tmp, solution.block(1)); tmp *= -1; tmp += system_rhs.block(0); - + A_inverse.vmult (solution.block(0), tmp); constraints.distribute (solution); @@ -1010,14 +1009,14 @@ void StokesProblem::find_dofs_on_lower_level (std::vector std::vector boundary_cell_dofs(dofs_per_cell); const unsigned int level = cell->level(); cell->get_mg_dof_indices (local_dof_indices); - for (unsigned int face_nr=0; + for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) { typename DoFHandler::face_iterator face = cell->face(face_nr); if(!cell->at_boundary(face_nr)) { //interior face - typename MGDoFHandler::cell_iterator neighbor + typename MGDoFHandler::cell_iterator neighbor = cell->neighbor(face_nr); // Do refinement face // from the coarse side @@ -1038,7 +1037,7 @@ void StokesProblem::find_dofs_on_lower_level (std::vector if(!cell->at_boundary(neighbor_face_nr)) { //other face is interior - typename MGDoFHandler::cell_iterator neighbor + typename MGDoFHandler::cell_iterator neighbor = cell->neighbor(neighbor_face_nr); // Do refinement face // from the coarse side @@ -1063,27 +1062,27 @@ void StokesProblem::find_dofs_on_lower_level (std::vector } } } - + template void StokesProblem::output_results (const unsigned int refinement_cycle) const { std::vector solution_names (dim, "velocity"); solution_names.push_back ("pressure"); - + std::vector data_component_interpretation (dim, DataComponentInterpretation::component_is_part_of_vector); data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); - + DataOut data_out; - data_out.attach_dof_handler (dof_handler); + data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, solution_names, DataOut::type_dof_data, data_component_interpretation); data_out.build_patches (); - + std::ostringstream filename; filename << "solution-" << Utilities::int_to_string (refinement_cycle, 2) @@ -1094,10 +1093,10 @@ StokesProblem::output_results (const unsigned int refinement_cycle) const } - + template void -StokesProblem::refine_mesh () +StokesProblem::refine_mesh () { Vector estimated_error_per_cell (triangulation.n_active_cells()); @@ -1117,9 +1116,9 @@ StokesProblem::refine_mesh () } - + template -void StokesProblem::run () +void StokesProblem::run () { { std::vector subdivisions (dim, 1); @@ -1131,42 +1130,42 @@ void StokesProblem::run () const Point top_right = (dim == 2 ? Point(1,1) : Point(1,1,1)); - + GridGenerator::subdivided_hyper_rectangle (triangulation, subdivisions, bottom_left, top_right); } - - + + for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) for (unsigned int f=0; f::faces_per_cell; ++f) if (cell->face(f)->center()[0] == 1) cell->face(f)->set_all_boundary_indicators(1); - - - + + + triangulation.refine_global (4-dim); - + for (unsigned int refinement_cycle = 0; refinement_cycle<6; ++refinement_cycle) { std::cout << "Refinement cycle " << refinement_cycle << std::endl; - + if (refinement_cycle > 0) refine_mesh (); - + setup_dofs (); std::cout << " Assembling..." << std::endl << std::flush; - assemble_system (); + assemble_system (); std::cout << " Solving..." << std::flush; solve_mg (); - + output_results (refinement_cycle); std::cout << std::endl; @@ -1174,12 +1173,12 @@ void StokesProblem::run () } - -int main () + +int main () { try { -// deallog.depth_console (0); + deallog.depth_console (0); StokesProblem<2> flow_problem(1); flow_problem.run (); @@ -1194,10 +1193,10 @@ int main () << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; - + return 1; } - catch (...) + catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" -- 2.39.5