From b3fcb5015920d94cbcb639d3f2fb3257562510bc Mon Sep 17 00:00:00 2001 From: heister Date: Sat, 3 Nov 2012 02:51:05 +0000 Subject: [PATCH] step-6 now uses ConstraintMatrix::distribute_local_to_global() git-svn-id: https://svn.dealii.org/trunk@27325 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 + deal.II/examples/step-6/doc/intro.dox | 15 +++ deal.II/examples/step-6/step-6.cc | 179 +++++++++++--------------- 3 files changed, 93 insertions(+), 107 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 94f3c94277..9ff00210fb 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -108,6 +108,12 @@ never working correctly and it is not used.
    +
  1. step-6 now uses ConstraintMatrix::distribute_local_to_global() +instead of condense(), which is the preferred way to use a ConstraintMatrix + (and the only sensible way in parallel). +
    +(Timo Heister, 2012/11/02) +
  2. Fixed: DoFTools::make_flux_sparsity_pattern wasn't prepared to deal with adaptively refined meshes in 1d.
    diff --git a/deal.II/examples/step-6/doc/intro.dox b/deal.II/examples/step-6/doc/intro.dox index 1d20efef62..0fc9752982 100644 --- a/deal.II/examples/step-6/doc/intro.dox +++ b/deal.II/examples/step-6/doc/intro.dox @@ -54,3 +54,18 @@ The only other new thing is a method to catch exceptions in the main function in order to output some information in case the program crashes for some reason. + +

    The ConstraintMatrix

    + +As explained above, we are going to use an object called ConstraintMatrix +that will store the constraints at the hanging nodes to insure the solution +is continuous at these nodes. We could first assemble the system as normal +and then condense out the degrees of freedom that need to be constrained. +This is also explained @ref constraints documentation +module. Instead we will go the more efficient route and eliminate the +constrained entries while we are copying from the local to the global system. +Because boundary conditions can be treated in the same way, we will +incorporate the them as constraints in the same ConstraintMatrix object. +This way, we don't need to apply the boundary conditions after assembly +(like we did in the earlier steps). + \ No newline at end of file diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc index caf09197cd..90c31101f1 100644 --- a/deal.II/examples/step-6/step-6.cc +++ b/deal.II/examples/step-6/step-6.cc @@ -72,7 +72,8 @@ // sure that the degrees of freedom // on hanging nodes conform to some // constraints such that the global - // solution is continuous. The + // solution is continuous. We are +// also going to store the boundary conditions in this object. The // following file contains a class // which is used to handle these // constraints: @@ -114,8 +115,7 @@ using namespace dealii; // (instead of the global refinement // in the previous examples), and a // variable which will hold the - // constraints associated to the - // hanging nodes. In addition, we + // constraints. In addition, we // have added a destructor to the // class for reasons that will become // clear when we discuss its @@ -144,9 +144,10 @@ class Step6 // This is the new variable in // the main class. We need an // object which holds a list of - // constraints originating from - // the hanging nodes: - ConstraintMatrix hanging_node_constraints; + // constraints to hold the + // hanging nodes and the + // boundary conditions. + ConstraintMatrix constraints; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; @@ -421,7 +422,7 @@ void Step6::setup_system () // hanging nodes. In the class // desclaration, we have already // allocated space for an object - // hanging_node_constraints + // constraints // that will hold a list of these // constraints (they form a matrix, // which is reflected in the name @@ -435,23 +436,38 @@ void Step6::setup_system () // over from computations on the // previous mesh before the last // adaptive refinement): - hanging_node_constraints.clear (); + constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, - hanging_node_constraints); + constraints); + + + // Now we are ready to interpolate the ZeroFunction to our + // boundary with indicator 0 (the whole boundary) and store + // the resulting constraints in our constraints + // object. Note that we do not to apply the boundary + // conditions after assembly, like we did in earlier steps. + // As almost all the stuff, + // the interpolation of boundary + // values works also for higher + // order elements without the need + // to change your code for that. We + // note that for proper results, it + // is important that the + // elimination of boundary nodes + // from the system of equations + // happens *after* the elimination + // of hanging nodes. For that reason + // we are filling the boundary values into + // the ContraintMatrix after the hanging node + // constraints. + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + constraints); + // The next step is closing - // this object. For this note that, - // in principle, the - // ConstraintMatrix class can - // hold other constraints as well, - // i.e. constraints that do not - // stem from hanging - // nodes. Sometimes, it is useful - // to use such constraints, in - // which case they may be added to - // the ConstraintMatrix object - // after the hanging node - // constraints were computed. After + // this object. After // all constraints have been added, // they need to be sorted and // rearranged to perform some @@ -460,7 +476,7 @@ void Step6::setup_system () // close() function, after which // no further constraints may be // added any more: - hanging_node_constraints.close (); + constraints.close (); // Now we first build our // compressed sparsity pattern like @@ -468,26 +484,18 @@ void Step6::setup_system () // examples. Nevertheless, we do // not copy it to the final // sparsity pattern immediately. + // Note that we call a variant of make_sparsity_pattern that takes + // the ConstraintMatrix as the third argument. We are letting + // the routine know, the we will never write into the locations + // given by constraints by setting the argument + // keep_constrained_dofs to false. If we were to + // condense the constraints after assembling, we would have + // to pass true instead. CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, c_sparsity); - - // The constrained hanging nodes - // will later be eliminated from - // the linear system of - // equations. When doing so, some - // additional entries in the global - // matrix will be set to non-zero - // values, so we have to reserve - // some space for them here. Since - // the process of elimination of - // these constrained nodes is - // called condensation, the - // functions that eliminate them - // are called condense for both - // the system matrix and right hand - // side, as well as for the - // sparsity pattern. - hanging_node_constraints.condense (c_sparsity); + DoFTools::make_sparsity_pattern(dof_handler, + c_sparsity, + constraints, + false /*keep_constrained_dofs*/); // Now all non-zero entries of the // matrix are known (i.e. those @@ -509,9 +517,10 @@ void Step6::setup_system () // @sect4{Step6::assemble_system} // Next, we have to assemble the - // matrix again. There are no code - // changes compared to step-5 except - // for a single place: We have to use + // matrix again. There are two code + // changes compared to step-5: +// +// First, we have to use // a higher-order quadrature formula // to account for the higher // polynomial degree in the finite @@ -523,11 +532,15 @@ void Step6::setup_system () // we had two points for bilinear // elements. Now we should use three // points for biquadratic elements. - // + +// Second, to copy the local matrix and vector on each cell +// into the global system, we are no longer using a hand-written +// loop. Instead, we use ConstraintMatrix::distribute_local_to_global +// that internally executes this loop and eliminates all the constraints +// at the same time. +// // The rest of the code that forms - // the local contributions and - // transfers them into the global - // objects remains unchanged. It is + // the local contributions remains unchanged. It is // worth noting, however, that under // the hood several things are // different than before. First, the @@ -596,67 +609,19 @@ void Step6::assemble_system () } cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; icondense - // function modifies the system so - // that the values in the solution - // corresponding to constrained - // nodes are invalid, but that the - // system still has a well-defined - // solution; we compute the correct + // are invalid. We compute the correct // values for these nodes at the - // end of the solve function). - - // As almost all the stuff before, - // the interpolation of boundary - // values works also for higher - // order elements without the need - // to change your code for that. We - // note that for proper results, it - // is important that the - // elimination of boundary nodes - // from the system of equations - // happens *after* the elimination - // of hanging nodes. - std::map boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(), - boundary_values); - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - system_rhs); + // end of the solve function. } @@ -671,8 +636,8 @@ void Step6::assemble_system () // have to incorporate hanging node // constraints. As mentioned above, // the degrees of freedom - // corresponding to hanging node - // constraints have been removed from + // from the ConstraintMatrix corresponding to hanging node + // constraints and boundary values have been removed from // the linear system by giving the // rows and columns of the matrix a // special treatment. This way, the @@ -684,7 +649,7 @@ void Step6::assemble_system () // constraints to assign to them the // values that they should have. This // process, called distributing - // hanging nodes, computes the values + // constraints, computes the values // of constrained nodes from the // values of the unconstrained ones, // and requires only a single @@ -703,7 +668,7 @@ void Step6::solve () solver.solve (system_matrix, solution, system_rhs, preconditioner); - hanging_node_constraints.distribute (solution); + constraints.distribute (solution); } @@ -1068,7 +1033,7 @@ void Step6::run () } // After we have finished computing - // the solution on the finesh mesh, + // the solution on the finest mesh, // and writing all the grids to // disk, we want to also write the // actual solution on this final -- 2.39.5