pattern. It tells the sparsity pattern that constrained entries should not
be inserted into the sparsity pattern. This makes sense for the way we're
going to do assembly, which does the elimination of constraints already
-locally. Hence, we will never write touch the constrained entries. This is
-the second key to increase this program's performance, since the system
-matrix will have less entries — this saves both memory and computing
-time when doing matrix-vector multiplications or applying a preconditioner.
+locally. Hence, we will never write the constrained entries. This is the
+second key to increase this program's performance, since the system matrix
+will have less entries — this saves both memory and computing time
+when doing matrix-vector multiplications or applying a preconditioner.
In a similar vein, we have to slightly modify the way we copy local
contributions into global matrices and vectors. In previous tutorial programs,
hanging_node_constraints.condense (system_rhs);
@endcode
We now replace copying and later condensing into one step, as already shown in
-@ref step_17 "step-17" and @ref step_18 "step-18":
+@ref step_17 "step-17", @ref step_18 "step-18", and @ref step_22 "step-22":
@code
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
cell->get_dof_indices (local_dof_indices);
hanging_node_constraints
- .distribute_local_to_global (cell_matrix,
+ .distribute_local_to_global (cell_matrix, cell_rhs,
local_dof_indices,
- system_matrix);
-
- hanging_node_constraints
- .distribute_local_to_global (cell_rhs,
- local_dof_indices,
- system_rhs);
+ system_matrix, system_rhs);
}
@endcode
Essentially, what the
-<code>hanging_node_constraints.distribute_local_to_global</code> calls do is
-to implement the same loop as before, but whenever we hit a constrained degree
+<code>hanging_node_constraints.distribute_local_to_global</code> call does is
+to implement the same loops as before, but whenever we hit a constrained degree
of freedom, the function does the right thing and already condenses it away.
Using this technique of eliminating constrained nodes already when
of having to go back later and change these objects. Timing these operations
shows that this makes the overall algorithms faster.
+In our program, we will also treat the boundary conditions as (possibly
+inhomogeneous) constraints and eliminate the matrix rows and columns to
+those as well. All we have to do for this is to call the function that
+interpolates the Dirichlet boundary conditions already in the setup phase in
+order to tell the ConstraintMatrix object about them, and then do the
+transfer from local to global data on matrix and vector simultaneously. This
+is exactly what we've shown in the @ref step_22 "step-22" tutorial program.
+
<h3>The test case</h3>
hp::QCollection<dim> quadrature_collection;
hp::QCollection<dim-1> face_quadrature_collection;
- ConstraintMatrix hanging_node_constraints;
+ ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
//
// This function is again an almost
// verbatim copy of what we already did in
- // step-6, with the main difference that we
- // don't directly build the sparsity
- // pattern, but first create an
- // intermediate object that we later copy
- // into the right data structure. In
+ // step-6. The first change is that we
+ // append the Dirichlet boundary conditions
+ // to the ConstraintMatrix object, which we
+ // consequently call just
+ // <code>constraints</code> instead of
+ // <code>hanging_node_constraints</code>. The
+ // second difference is that we don't
+ // directly build the sparsity pattern, but
+ // first create an intermediate object that
+ // we later copy into the usual
+ // SparsityPattern data structure, since
+ // this is more efficient for the problem
+ // with many entries per row (and different
+ // number of entries in different rows). In
// another slight deviation, we do not
- // first build the sparsity pattern then
- // condense away constrained degrees of
- // freedom, but pass the constraint matrix
- // object directly to the function that
- // builds the sparsity pattern. We disable
- // the insertion of constrained entries
- // with <tt>false</tt> as fourth argument
- // in the DoFTools::make_sparsity_pattern
- // function. Both of these changes are
+ // first build the sparsity pattern and
+ // then condense away constrained degrees
+ // of freedom, but pass the constraint
+ // matrix object directly to the function
+ // that builds the sparsity pattern. We
+ // disable the insertion of constrained
+ // entries with <tt>false</tt> as fourth
+ // argument in the
+ // DoFTools::make_sparsity_pattern
+ // function. All of these changes are
// explained in the introduction of this
// program.
//
- // The second change, maybe hidden in plain
+ // The last change, maybe hidden in plain
// sight, is that the dof_handler variable
// here is an hp object -- nevertheless all
// the function calls we had before still
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
- hanging_node_constraints.clear ();
+ constraints.clear ();
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(),
+ constraints);
DoFTools::make_hanging_node_constraints (dof_handler,
- hanging_node_constraints);
- hanging_node_constraints.close ();
+ constraints);
+ constraints.close ();
CompressedSetSparsityPattern csp (dof_handler.n_dofs(),
dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (dof_handler, csp,
- hanging_node_constraints, false);
+ DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
sparsity_pattern.copy_from (csp);
system_matrix.reinit (sparsity_pattern);
Vector<double> cell_rhs;
std::vector<unsigned int> local_dof_indices;
-
+
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
local_dof_indices.resize (dofs_per_cell);
cell->get_dof_indices (local_dof_indices);
- hanging_node_constraints
- .distribute_local_to_global (cell_matrix,
- local_dof_indices,
- system_matrix);
-
- hanging_node_constraints
- .distribute_local_to_global (cell_rhs,
- local_dof_indices,
- system_rhs);
+ constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+ local_dof_indices,
+ system_matrix, system_rhs);
}
- // After the steps already discussed above,
- // all we still have to do is to treat
- // Dirichlet boundary values
- // correctly. This works in exactly the
- // same way as for non-<i>hp</i> objects:
- std::map<unsigned int,double> boundary_values;
- VectorTools::interpolate_boundary_values (dof_handler,
- 0,
- ZeroFunction<dim>(),
- boundary_values);
- MatrixTools::apply_boundary_values (boundary_values,
- system_matrix,
- solution,
- system_rhs);
+ // Now with the loop over all cells
+ // finished, we are done for this
+ // function. The steps we still had to do
+ // at this point in earlier tutorial
+ // programs, namely condensing hanging
+ // node constraints and applying
+ // Dirichlet boundary conditions, have
+ // been taken care of by the
+ // ConstraintMatrix object
+ // <code>constraints</code> on the fly.
}
cg.solve (system_matrix, solution, system_rhs,
preconditioner);
- hanging_node_constraints.distribute (solution);
+ constraints.distribute (solution);
}
// VTK format):
const std::string filename = "solution-" +
Utilities::int_to_string (cycle, 2) +
- ".gmv";
+ ".vtk";
std::ofstream output (filename.c_str());
- data_out.write_gmv (output);
+ data_out.write_vtk (output);
}
// After this, we would like to actually
<< dof_handler.n_dofs()
<< std::endl
<< " Number of constraints : "
- << hanging_node_constraints.n_constraints()
+ << constraints.n_constraints()
<< std::endl;
assemble_system ();
projects / dealii-svn.git / commitdiff