From 37874497544af5b0eb4dc7897e8cabfbf955a6a6 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 11 Aug 2008 22:27:21 +0000 Subject: [PATCH] Provide the infrastructure for building the preconditioner as a system of three decoupled blocks each of which consists of a scalar Laplace operator. git-svn-id: https://svn.dealii.org/trunk@16510 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 128 ++++++++++++++++++++++++---- 1 file changed, 111 insertions(+), 17 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index addc5e237f..e186f05786 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -33,6 +33,7 @@ #include #include +#include #include #include #include @@ -123,14 +124,24 @@ class BoussinesqFlowProblem const unsigned int degree; Triangulation triangulation; + FESystem fe; DoFHandler dof_handler; + FE_Q scalar_velocity_fe; + DoFHandler scalar_velocity_dof_handler; + ConstraintMatrix hanging_node_constraints; + ConstraintMatrix scalar_velocity_hanging_node_constraints; BlockSparsityPattern sparsity_pattern; BlockSparseMatrix system_matrix; + SparsityPattern scalar_velocity_sparsity_pattern; + SparseMatrix scalar_velocity_system_matrix; + std::vector scalar_to_vector_velocities[dim]; + + double time_step; unsigned int timestep_number; @@ -582,6 +593,8 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) FE_Q(degree), 1, FE_DGQ(degree-1), 1), dof_handler (triangulation), + scalar_velocity_fe (degree+1), + scalar_velocity_dof_handler (triangulation), time_step (0), rebuild_matrices (true), rebuild_preconditioner (true) @@ -679,6 +692,18 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) << std::endl << std::endl; + + // also initialize the DoFHandler + // for the scalar version of the + // Laplace equation on the velocity + // variables + scalar_velocity_dof_handler.distribute_dofs (scalar_velocity_fe); + DoFRenumbering::Cuthill_McKee (scalar_velocity_dof_handler); + scalar_velocity_hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (scalar_velocity_dof_handler, + scalar_velocity_hanging_node_constraints); + scalar_velocity_hanging_node_constraints.close (); + // The next step is to // create the sparsity // pattern for the system matrix @@ -732,28 +757,79 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) // the sparsity pattern. if (setup_matrices == true) { - system_matrix.clear (); + { + system_matrix.clear (); - BlockCompressedSetSparsityPattern csp (3,3); + BlockCompressedSetSparsityPattern csp (3,3); - csp.block(0,0).reinit (n_u, n_u); - csp.block(0,1).reinit (n_u, n_p); - csp.block(0,2).reinit (n_u, n_T); - csp.block(1,0).reinit (n_p, n_u); - csp.block(1,1).reinit (n_p, n_p); - csp.block(1,2).reinit (n_p, n_T); - csp.block(2,0).reinit (n_T, n_u); - csp.block(2,1).reinit (n_T, n_p); - csp.block(2,2).reinit (n_T, n_T); + csp.block(0,0).reinit (n_u, n_u); + csp.block(0,1).reinit (n_u, n_p); + csp.block(0,2).reinit (n_u, n_T); + csp.block(1,0).reinit (n_p, n_u); + csp.block(1,1).reinit (n_p, n_p); + csp.block(1,2).reinit (n_p, n_T); + csp.block(2,0).reinit (n_T, n_u); + csp.block(2,1).reinit (n_T, n_p); + csp.block(2,2).reinit (n_T, n_T); + + csp.collect_sizes (); - csp.collect_sizes (); + DoFTools::make_sparsity_pattern (dof_handler, csp); + hanging_node_constraints.condense (csp); + sparsity_pattern.copy_from (csp); + + system_matrix.reinit (sparsity_pattern); + } - DoFTools::make_sparsity_pattern (dof_handler, csp); - hanging_node_constraints.condense (csp); - sparsity_pattern.copy_from (csp); + { + scalar_velocity_system_matrix.clear (); - system_matrix.reinit (sparsity_pattern); + CompressedSetSparsityPattern csp (scalar_velocity_dof_handler.n_dofs(), + scalar_velocity_dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (scalar_velocity_dof_handler, csp); + scalar_velocity_hanging_node_constraints.condense (csp); + scalar_velocity_sparsity_pattern.copy_from (csp); + + scalar_velocity_system_matrix.reinit (scalar_velocity_sparsity_pattern); + } } + + { + InterGridMap > scalar_to_vector; + scalar_to_vector.make_mapping (scalar_velocity_dof_handler, + dof_handler); + std::vector + cell_dof_indices (fe.dofs_per_cell), + scalar_velocity_cell_dof_indices (scalar_velocity_fe.dofs_per_cell); + + for (unsigned int component=0; component::active_cell_iterator + cell = scalar_velocity_dof_handler.begin_active(); + cell != scalar_velocity_dof_handler.end(); + ++cell) + { + cell->get_dof_indices (scalar_velocity_cell_dof_indices); + scalar_to_vector[cell]->get_dof_indices (cell_dof_indices); + + for (unsigned int i=0; i::assemble_system () std::cout << " Rebuilding preconditioner..." << std::flush; - A_preconditioner + MatrixTools::create_laplace_matrix (scalar_velocity_dof_handler, + QGauss (degree+1), + scalar_velocity_system_matrix); + scalar_velocity_hanging_node_constraints.condense (scalar_velocity_system_matrix); + { + std::map boundary_values; + VectorTools::interpolate_boundary_values (scalar_velocity_dof_handler, + 0, + ZeroFunction(), + boundary_values); + Vector dummy (scalar_velocity_dof_handler.n_dofs()); + MatrixTools::apply_boundary_values (boundary_values, + scalar_velocity_system_matrix, + dummy, + dummy); + } + + + A_preconditioner = boost::shared_ptr::type> (new typename InnerPreconditioner::type()); A_preconditioner->initialize (system_matrix.block(0,0), -- 2.39.5