From 2f68e10c84da825e55dbd005e3e5b14ab657518c Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Sat, 11 May 2024 10:29:37 -0400 Subject: [PATCH] step-72: use AffineConstraints only Same as done in step-15 in #16981 --- examples/step-72/step-72.cc | 171 +++++++++--------------------------- 1 file changed, 41 insertions(+), 130 deletions(-) diff --git a/examples/step-72/step-72.cc b/examples/step-72/step-72.cc index b0479b77e5..b1fe79105c 100644 --- a/examples/step-72/step-72.cc +++ b/examples/step-72/step-72.cc @@ -146,13 +146,12 @@ namespace Step72 void run(const int formulation, const double tolerance); private: - void setup_system(const bool initial_step); + void setup_system(); void assemble_system_unassisted(); void assemble_system_with_residual_linearization(); void assemble_system_using_energy_functional(); void solve(); void refine_mesh(); - void set_boundary_values(); double compute_residual(const double alpha) const; double determine_step_length() const; void output_results(const unsigned int refinement_cycle) const; @@ -163,7 +162,8 @@ namespace Step72 const FE_Q fe; const QGauss quadrature_formula; - AffineConstraints hanging_node_constraints; + AffineConstraints nonzero_constraints; + AffineConstraints zero_constraints; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; @@ -212,26 +212,30 @@ namespace Step72 // data structures, namely the DoFHandler, the hanging node constraints // applied to the problem, and the linear system. template - void MinimalSurfaceProblem::setup_system(const bool initial_step) + void MinimalSurfaceProblem::setup_system() { - if (initial_step) - { - dof_handler.distribute_dofs(fe); - current_solution.reinit(dof_handler.n_dofs()); - - hanging_node_constraints.clear(); - DoFTools::make_hanging_node_constraints(dof_handler, - hanging_node_constraints); - hanging_node_constraints.close(); - } - + dof_handler.distribute_dofs(fe); + current_solution.reinit(dof_handler.n_dofs()); newton_update.reinit(dof_handler.n_dofs()); system_rhs.reinit(dof_handler.n_dofs()); - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, dsp); + zero_constraints.clear(); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + zero_constraints); + DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints); + zero_constraints.close(); - hanging_node_constraints.condense(dsp); + nonzero_constraints.clear(); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + BoundaryValues(), + nonzero_constraints); + nonzero_constraints.close(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, zero_constraints); sparsity_pattern.copy_from(dsp); system_matrix.reinit(sparsity_pattern); @@ -385,15 +389,8 @@ namespace Step72 const std::vector &local_dof_indices = copy_data.local_dof_indices[0]; - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int j = 0; j < dofs_per_cell; ++j) - system_matrix.add(local_dof_indices[i], - local_dof_indices[j], - cell_matrix(i, j)); - - system_rhs(local_dof_indices[i]) += cell_rhs(i); - } + zero_constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); }; // We have all of the required functions definitions in place, so @@ -410,22 +407,6 @@ namespace Step72 sample_scratch_data, sample_copy_data, MeshWorker::assemble_own_cells); - - // And finally, as is done in step-15, we remove hanging nodes from the - // system and apply zero boundary values to the linear system that defines - // the Newton updates $\delta u^n$. - hanging_node_constraints.condense(system_matrix); - hanging_node_constraints.condense(system_rhs); - - std::map boundary_values; - VectorTools::interpolate_boundary_values(dof_handler, - 0, - Functions::ZeroFunction(), - boundary_values); - MatrixTools::apply_boundary_values(boundary_values, - system_matrix, - newton_update, - system_rhs); } // @sect5{Assembly via differentiation of the residual vector} @@ -617,15 +598,8 @@ namespace Step72 const std::vector &local_dof_indices = copy_data.local_dof_indices[0]; - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int j = 0; j < dofs_per_cell; ++j) - system_matrix.add(local_dof_indices[i], - local_dof_indices[j], - cell_matrix(i, j)); - - system_rhs(local_dof_indices[i]) += cell_rhs(i); - } + zero_constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); }; MeshWorker::mesh_loop(dof_handler.active_cell_iterators(), @@ -634,19 +608,6 @@ namespace Step72 sample_scratch_data, sample_copy_data, MeshWorker::assemble_own_cells); - - hanging_node_constraints.condense(system_matrix); - hanging_node_constraints.condense(system_rhs); - - std::map boundary_values; - VectorTools::interpolate_boundary_values(dof_handler, - 0, - Functions::ZeroFunction(), - boundary_values); - MatrixTools::apply_boundary_values(boundary_values, - system_matrix, - newton_update, - system_rhs); } // @sect5{Assembly via differentiation of the energy functional} @@ -785,15 +746,8 @@ namespace Step72 const std::vector &local_dof_indices = copy_data.local_dof_indices[0]; - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int j = 0; j < dofs_per_cell; ++j) - system_matrix.add(local_dof_indices[i], - local_dof_indices[j], - cell_matrix(i, j)); - - system_rhs(local_dof_indices[i]) += cell_rhs(i); - } + zero_constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); }; MeshWorker::mesh_loop(dof_handler.active_cell_iterators(), @@ -802,19 +756,6 @@ namespace Step72 sample_scratch_data, sample_copy_data, MeshWorker::assemble_own_cells); - - hanging_node_constraints.condense(system_matrix); - hanging_node_constraints.condense(system_rhs); - - std::map boundary_values; - VectorTools::interpolate_boundary_values(dof_handler, - 0, - Functions::ZeroFunction(), - boundary_values); - MatrixTools::apply_boundary_values(boundary_values, - system_matrix, - newton_update, - system_rhs); } @@ -833,7 +774,7 @@ namespace Step72 solver.solve(system_matrix, newton_update, system_rhs, preconditioner); - hanging_node_constraints.distribute(newton_update); + zero_constraints.distribute(newton_update); const double alpha = determine_step_length(); current_solution.add(alpha, newton_update); @@ -862,46 +803,21 @@ namespace Step72 0.03); triangulation.prepare_coarsening_and_refinement(); - SolutionTransfer solution_transfer(dof_handler); - solution_transfer.prepare_for_coarsening_and_refinement(current_solution); - triangulation.execute_coarsening_and_refinement(); - - dof_handler.distribute_dofs(fe); - Vector tmp(dof_handler.n_dofs()); - solution_transfer.interpolate(current_solution, tmp); - current_solution = tmp; + SolutionTransfer solution_transfer(dof_handler); + const Vector coarse_solution = current_solution; + solution_transfer.prepare_for_coarsening_and_refinement(coarse_solution); - hanging_node_constraints.clear(); - DoFTools::make_hanging_node_constraints(dof_handler, - hanging_node_constraints); - hanging_node_constraints.close(); + triangulation.execute_coarsening_and_refinement(); - set_boundary_values(); + setup_system(); - setup_system(false); + solution_transfer.interpolate(coarse_solution, current_solution); + nonzero_constraints.distribute(current_solution); } - // @sect4{MinimalSurfaceProblem::set_boundary_values} - - // The choice of boundary conditions remains identical to step-15... - template - void MinimalSurfaceProblem::set_boundary_values() - { - std::map boundary_values; - VectorTools::interpolate_boundary_values(dof_handler, - 0, - BoundaryValues(), - boundary_values); - for (auto &boundary_value : boundary_values) - current_solution(boundary_value.first) = boundary_value.second; - - hanging_node_constraints.distribute(current_solution); - } - - // @sect4{MinimalSurfaceProblem::compute_residual} // ... as does the function used to compute the residual during the @@ -951,16 +867,11 @@ namespace Step72 } cell->get_dof_indices(local_dof_indices); - for (unsigned int i = 0; i < dofs_per_cell; ++i) - residual(local_dof_indices[i]) += cell_residual(i); + zero_constraints.distribute_local_to_global(cell_residual, + local_dof_indices, + residual); } - hanging_node_constraints.condense(residual); - - for (const types::global_dof_index i : - DoFTools::extract_boundary_dofs(dof_handler)) - residual(i) = 0; - return residual.l2_norm(); } @@ -1033,8 +944,8 @@ namespace Step72 GridGenerator::hyper_ball(triangulation); triangulation.refine_global(2); - setup_system(/*first time=*/true); - set_boundary_values(); + setup_system(); + nonzero_constraints.distribute(current_solution); double last_residual_norm = std::numeric_limits::max(); unsigned int refinement_cycle = 0; -- 2.39.5