From db0af0507954c86947f2fbff8de08db4e61d27f8 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 1 Sep 2018 18:08:10 -0400 Subject: [PATCH] step-37: Redo the second Dirichlet boundary algorithm. Doing this with blocks will take a lot more work than the original version implied, so this example now presents an alternative that constructs a second MatrixFree object instead. --- examples/step-37/doc/results.dox | 106 ++++++++++++++----------------- 1 file changed, 46 insertions(+), 60 deletions(-) diff --git a/examples/step-37/doc/results.dox b/examples/step-37/doc/results.dox index baf5b54ed5..db5fb60b7b 100644 --- a/examples/step-37/doc/results.dox +++ b/examples/step-37/doc/results.dox @@ -621,72 +621,58 @@ residual with respect to the solution variable.
Use LaplaceOperator with a second AffineConstraints object without Dirichlet conditions
-A second alternative to get the right hand side that re-uses the -@p LaplaceOperator::apply_add() function is to instead add a second constraint -matrix that skips Dirichlet constraints on the read operation. To do this, we -initialize a MatrixFree object in a more extended way with two different -DoFHandler - AffineConstraints combinations. The zeroth component includes -Dirichlet conditions for solving the linear system, whereas first component -also reads from Dirichlet-constrained degrees of freedom for the right hand -side assembly: - -@code - constraints.clear(); - constraints.reinit(locally_relevant_dofs); - DoFTools::make_hanging_node_constraints(dof_handler, - constraints); - VectorTools::interpolate_boundary_values (dof_handler, - 0, - BoundaryValues(), - constraints); - constraints.close(); - constraints_without_dirichlet.clear(); - constraints_without_dirichlet.reinit(locally_relevant_dofs); - DoFTools::make_hanging_node_constraints(dof_handler, - constraints_without_dirichlet); - constraints_without_dirichlet.close(); - - std::vector *> dof_handlers(2, &dof_handler); - { - std::vector constraint(2); - constraint[0] = &constraints; - constraint[1] = &constraints_without_dirichlet; - typename MatrixFree::AdditionalData additional_data; - additional_data.mapping_update_flags = (update_gradients | update_JxW_values | - update_quadrature_points); - - std::shared_ptr > matrix_free; - matrix_free->reinit (mapping, dof_handlers, constraint, - QGauss<1>(fe.degree+1), additional_data); - - // select zeroth block in matrix_free for the main matrix - std::vector selected_block {0}; - laplace_operator.initialize(matrix_free, selected_block); - } -@endcode - -This @p matrix_free object is then passed to a @p LaplaceOperator class -instance @p laplace_operator that gets used in the linear solver. Alongside, -we create a second @p LaplaceOperator object that fills the right hand side: +A second alternative to get the right hand side that re-uses the @p +LaplaceOperator::apply_add() function is to instead add a second LaplaceOperator +that skips Dirichlet constraints. To do this, we initialize a second MatrixFree +object which does not have any boundary value constraints. This @p matrix_free +object is then passed to a @p LaplaceOperator class instance @p +inhomogeneous_operator that is only used to create the right hand side: @code template void LaplaceProblem::assemble_rhs() { - LaplaceOperator laplace_operator_inhomogenous; - - // select first block in matrix_free to use constraints_without_dirichlet - std::vector selected_block{1}; - laplace_operator_inhomogeneous.initialize(matrix_free, selected_block); - solution = 0; + system_rhs = 0; + AffineConstraints no_constraints; + no_constraints.close(); + LaplaceOperator inhomogeneous_operator; + + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + std::shared_ptr> matrix_free( + new MatrixFree()); + matrix_free->reinit(dof_handler, + no_constraints, + QGauss<1>(fe.degree + 1), + additional_data); + inhomogeneous_operator.initialize(matrix_free); + + solution = 0.0; constraints.distribute(solution); - laplace_operator_inhomogeneous.vmult(system_rhs, solution); - system_rhs *= -1.; + inhomogeneous_operator.evaluate_coefficient(Coefficient()); + inhomogeneous_operator.vmult(system_rhs, solution); + system_rhs *= -1.0; - // proceed as usual with integration of right hand side function... + FEEvaluation phi( + *inhomogeneous_operator.get_matrix_free()); + for (unsigned int cell = 0; + cell < inhomogeneous_operator.get_matrix_free()->n_macro_cells(); + ++cell) + { + phi.reinit(cell); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_value(make_vectorized_array(1.0), q); + phi.integrate(true, false); + phi.distribute_local_to_global(system_rhs); + } + system_rhs.compress(VectorOperation::add); } @endcode -Instead of adding a second DoFHandler - AffineConstraints pair to -the same MatrixFree::reinit() call, one could of course also construct an -independent MatrixFree object that feeds the second @p LaplaceOperator instance, -see also the discussion in MatrixFreeOperators::Base. +A more sophisticated implementation of this technique could reuse the original +MatrixFree object. This can be done by initializing the MatrixFree object with +multiple blocks, where each block corresponds to a different AffineConstraints +object. Doing this would require making substantial modifications to the +LaplaceOperator class, but the MatrixFreeOperators::LaplaceOperator class that +comes with the library can do this. See the discussion on blocks in +MatrixFreeOperators::Base for more information on how to set up blocks. -- 2.39.5