From: kronbichler Date: Wed, 12 Nov 2008 09:03:32 +0000 (+0000) Subject: Added implementation of ConstraintMatrix::condense() for CompressedSimpleSparsityPatt... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=55443f24e2023572c6a25c1543cb7df5c097ad21;p=dealii-svn.git Added implementation of ConstraintMatrix::condense() for CompressedSimpleSparsityPattern by copying the code. Changed step-22, step-31 and step-32 to that sparsity pattern format, since it is faster than CompressedSetSparsityPattern and uses less memory. git-svn-id: https://svn.dealii.org/trunk@17560 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index 137e941665..06ec28b472 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -31,9 +31,11 @@ template class FullMatrix; class SparsityPattern; class CompressedSparsityPattern; class CompressedSetSparsityPattern; +class CompressedSimpleSparsityPattern; class BlockSparsityPattern; class BlockCompressedSparsityPattern; class BlockCompressedSetSparsityPattern; +class BlockCompressedSimpleSparsityPattern; template class SparseMatrix; template class BlockSparseMatrix; class BlockIndices; @@ -719,6 +721,14 @@ class ConstraintMatrix : public Subscriptor */ void condense (CompressedSetSparsityPattern &sparsity) const; + /** + * Same function as above, but + * condenses compressed + * sparsity patterns, which are + * based on the ''simple'' aproach. + */ + void condense (CompressedSimpleSparsityPattern &sparsity) const; + /** * Same function as above, but * condenses square compressed @@ -752,6 +762,13 @@ class ConstraintMatrix : public Subscriptor * sparsity patterns. */ void condense (BlockCompressedSetSparsityPattern &sparsity) const; + + /** + * Same function as above, but + * condenses square compressed + * sparsity patterns. + */ + void condense (BlockCompressedSimpleSparsityPattern &sparsity) const; /** diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 5dbd6097a9..2421304424 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -1095,6 +1095,7 @@ void ConstraintMatrix::condense (CompressedSparsityPattern &sparsity) const } + void ConstraintMatrix::condense (CompressedSetSparsityPattern &sparsity) const { Assert (sorted == true, ExcMatrixNotClosed()); @@ -1185,6 +1186,141 @@ void ConstraintMatrix::condense (CompressedSetSparsityPattern &sparsity) const +void ConstraintMatrix::condense (CompressedSimpleSparsityPattern &sparsity) const +{ + Assert (sorted == true, ExcMatrixNotClosed()); + Assert (sparsity.n_rows() == sparsity.n_cols(), + ExcNotQuadratic()); + + // store for each index whether it must be + // distributed or not. If entry is + // numbers::invalid_unsigned_int, + // no distribution is necessary. + // otherwise, the number states which line + // in the constraint matrix handles this + // index + std::vector distribute(sparsity.n_rows(), + numbers::invalid_unsigned_int); + + for (unsigned int c=0; c distribute (sparsity.n_rows(), + numbers::invalid_unsigned_int); + + for (unsigned int c=0; c(c); + + const unsigned int n_rows = sparsity.n_rows(); + for (unsigned int row=0; row + block_index = index_mapping.global_to_local(row); + const unsigned int block_row = block_index.first; + const unsigned int local_row = block_index.second; + + if (distribute[row] == numbers::invalid_unsigned_int) + // regular line. loop over + // all columns and see + // whether this column must + // be distributed. note that + // as we proceed to + // distribute cols, the loop + // over cols may get longer. + // + // don't try to be clever + // here as in the algorithm + // for the + // CompressedSparsityPattern, + // as that would be much more + // complicated here. after + // all, we know that + // compressed patterns are + // inefficient... + { + + // to loop over all entries + // in this row, we have to + // loop over all blocks in + // this blockrow and the + // corresponding row + // therein + for (unsigned int block_col=0; block_col2D calculations -Running the program with the space dimension set to 2 in -main() yields the following output: +Running the program with the space dimension set to 2 in main() +yields the following output (when the flag is set to optimized in the +Makefile): @code examples/step-22> make run ============================ Remaking Makefile.dep @@ -54,16 +55,9 @@ Refinement cycle 5 Assembling... Computing preconditioner... Solving... 11 outer CG Schur complement iterations for pressure - -Refinement cycle 6 - Number of active cells: 8896 - Number of degrees of freedom: 83885 (74474+9411) - Assembling... - Computing preconditioner... - Solving... 11 outer CG Schur complement iterations for pressure @endcode -The entire computation above takes about 30 seconds on a reasonably +The entire computation above takes about 20 seconds on a reasonably quick (for 2007 standards) machine. What we see immediately from this is that the number of (outer) @@ -181,12 +175,7 @@ Refinement cycle 5 Again, we see that the number of outer iterations does not increase as we refine the mesh. Nevertheless, the compute time increases significantly: for each of the iterations above separately, it takes a -few seconds, a few seconds, 1min, 5min, 21min, and 1h35. (One level -more can be run on machines with more memory, though one may have to -replace the BlockCompressedSetSparsityPattern class by the -BlockCompressedSparsityPattern class as the former needs significantly -more memory and becomes the bottleneck; the latter class, however, has -a superlinear runtime complexity.) This overall +few seconds, a few seconds, 1min, 5min, 21min, and 1h35. This overall superlinear (in the number of unknowns) increase in runtime is due to the fact that our inner solver is not ${\cal O}(N)$: a simple experiment shows that as we keep refining the mesh, the average number of diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index a6e73ce1e3..205b8cc518 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -336,22 +336,20 @@ InverseMatrix::InverseMatrix (const Matrix &m, // This is the implementation of the // vmult function. - // Note that we use a rather large - // tolerance for the solver - // control. The reason for this is - // that the function is used very - // frequently, and hence, any - // additional effort to make the - // residual in the CG solve smaller - // makes the solution more - // expensive. Note that we do not - // only use this class as a - // preconditioner for the Schur - // complement, but also when forming - // the inverse of the Laplace matrix - - // which needs to be accurate in - // order to obtain a solution to the - // right problem. + // In this class we use a rather large + // tolerance for the solver control. The + // reason for this is that the function is + // used very frequently, and hence, any + // additional effort to make the residual + // in the CG solve smaller makes the + // solution more expensive. Note that we do + // not only use this class as a + // preconditioner for the Schur complement, + // but also when forming the inverse of the + // Laplace matrix – which is hence + // directly responsible for the accuracy of + // the solution itself, so we can't choose + // a too large tolerance, either. template void InverseMatrix::vmult (Vector &dst, const Vector &src) const @@ -599,95 +597,88 @@ void StokesProblem::setup_dofs () << std::endl; // The next task is to allocate a - // sparsity pattern for the system - // matrix we will create. We could - // do this in the same way as in - // step-20, i.e. directly build an - // object of type SparsityPattern - // through + // sparsity pattern for the system matrix + // we will create. We could do this in + // the same way as in step-20, + // i.e. directly build an object of type + // SparsityPattern through // DoFTools::make_sparsity_pattern. However, - // there is a major reason not to - // do so: In 3D, the function + // there is a major reason not to do so: + // In 3D, the function // DoFTools::max_couplings_between_dofs - // yields a conservative but rather - // large number for the coupling - // between the individual dofs, so - // that the memory initially - // provided for the creation of the - // sparsity pattern of the matrix - // is far too much -- so much - // actually that the initial - // sparsity pattern won't even fit - // into the physical memory of most - // systems already for - // moderately-sized 3D problems, - // see also the discussion in - // step-18. Instead, we first - // build a temporary object that - // uses a different data structure - // that doesn't require allocating - // more memory than necessary but - // isn't suitable for use as a - // basis of SparseMatrix or - // BlockSparseMatrix objects; in a - // second step we then copy this - // object into an object of - // BlockSparsityPattern. This is - // entirely analgous to what we - // already did in step-11 and - // step-18. + // yields a conservative but rather large + // number for the coupling between the + // individual dofs, so that the memory + // initially provided for the creation of + // the sparsity pattern of the matrix is + // far too much -- so much actually that + // the initial sparsity pattern won't + // even fit into the physical memory of + // most systems already for + // moderately-sized 3D problems, see also + // the discussion in step-18. Instead, + // we first build a temporary object that + // uses a different data structure that + // doesn't require allocating more memory + // than necessary but isn't suitable for + // use as a basis of SparseMatrix or + // BlockSparseMatrix objects; in a second + // step we then copy this object into an + // object of BlockSparsityPattern. This + // is entirely analgous to what we + // already did in step-11 and step-18. // - // There is one snag again here, - // though: just as in step-27, it - // turns out that using the - // CompressedSparsityPattern (or - // the block version - // BlockCompressedSparsityPattern - // we would use here) - // has a bottleneck that makes the - // algorithm to build the sparsity - // pattern be quadratic in the - // number of degrees of - // freedom. This doesn't become - // noticable until we get well into - // the range of several 100,000 - // degrees of freedom, but - // eventually dominates the setup - // of the linear system when we get - // to more than a million degrees - // of freedom. This is due to the - // data structures used in the + // There is one snag again here, though: + // it turns out that using the + // CompressedSparsityPattern (or the + // block version + // BlockCompressedSparsityPattern we + // would use here) has a bottleneck that + // makes the algorithm to build the + // sparsity pattern be quadratic in the + // number of degrees of freedom. This + // doesn't become noticable until we get + // well into the range of several 100,000 + // degrees of freedom, but eventually + // dominates the setup of the linear + // system when we get to more than a + // million degrees of freedom. This is + // due to the data structures used in the // CompressedSparsityPattern class, // nothing that can easily be - // changed. Fortunately, there is - // an easy solution, as already - // pointed out in step-27: the - // CompressedSetSparsityPattern - // class (and its block variant - // BlockCompressedSetSparsityPattern) - // has exactly the same interface, - // uses a different internal data - // structure, is slightly slower - // for smaller numbers of degrees - // of freedom (but there we don't - // care that much anyway) but is - // linear in the number of degrees - // of freedom and therefore much - // more efficient for large - // problems. + // changed. Fortunately, there is an easy + // solution: the + // CompressedSimpleSparsityPattern class + // (and its block variant + // BlockCompressedSimpleSparsityPattern) + // has exactly the same interface, uses a + // different internal data structure, is + // slightly slower for smaller numbers of + // degrees of freedom (but there we don't + // care that much anyway) but is linear + // in the number of degrees of freedom + // and therefore much more efficient for + // large problems. As another + // alternative, we could also have chosen + // the class + // BlockCompressedSetSparsityPattern that + // uses yet another strategy for internal + // memory management. Though, that class + // turns out to be more memory-demanding + // than + // BlockCompressedSimpleSparsityPattern + // for this example. // - // Consequently, this is the class - // that we will use for our - // intermediate sparsity - // representation. All this is done - // inside a new scope, which means - // that the memory of - // csp will be - // released once the information + // Consequently, this is the class that + // we will use for our intermediate + // sparsity representation. All this is + // done inside a new scope, which means + // that the memory of csp + // will be released once the information // has been copied to // sparsity_pattern. { - BlockCompressedSetSparsityPattern csp (2,2); + BlockCompressedSimpleSparsityPattern csp (2,2); csp.block(0,0).reinit (n_u, n_u); csp.block(1,0).reinit (n_p, n_u); @@ -858,23 +849,20 @@ void StokesProblem::assemble_system () } } - // Note that in the above - // computation of the local - // matrix contribution we added - // the term phi_p[i] * - // phi_p[j] , yielding a + // Note that in the above computation + // of the local matrix contribution + // we added the term phi_p[i] + // * phi_p[j] , yielding a // pressure mass matrix in the - // $(1,1)$ block of the matrix - // as discussed in the - // introduction. That this term - // only ends up in the $(1,1)$ - // block stems from the fact - // that both of the factors in - // phi_p[i] * - // phi_p[j] are only - // non-zero when all the other - // terms vanish (and the other - // way around). + // $(1,1)$ block of the matrix as + // discussed in the + // introduction. That this term only + // ends up in the $(1,1)$ block stems + // from the fact that both of the + // factors in phi_p[i] * + // phi_p[j] are only non-zero + // when all the other terms vanish + // (and the other way around). // // Note also that operator* is // overloaded for symmetric @@ -923,33 +911,28 @@ void StokesProblem::assemble_system () } // After the addition of the local - // contributions, we have to - // condense the hanging node - // constraints and interpolate - // Dirichlet boundary conditions. - // Further down below where we set - // up the mesh, we will associate + // contributions, we have to condense the + // hanging node constraints and + // interpolate Dirichlet boundary + // conditions. Further down below where + // we set up the mesh, we will associate // the top boundary where we impose - // Dirichlet boundary conditions - // with boundary indicator 1. We - // will have to pass this boundary - // indicator as second argument to - // the function below interpolating - // boundary values. There is one - // more thing, though. The - // function describing the - // Dirichlet conditions was defined - // for all components, both - // velocity and pressure. However, - // the Dirichlet conditions are to - // be set for the velocity only. - // To this end, we use a + // Dirichlet boundary conditions with + // boundary indicator 1. We will have to + // pass this boundary indicator as second + // argument to the function below + // interpolating boundary values. There + // is one more thing, though. The + // function describing the Dirichlet + // conditions was defined for all + // components, both velocity and + // pressure. However, the Dirichlet + // conditions are to be set for the + // velocity only. To this end, we use a // component_mask that - // filters out the pressure - // component, so that the - // condensation is performed on - // velocity degrees of freedom - // only: + // filters out the pressure component, so + // that the condensation is performed on + // velocity degrees of freedom only: hanging_node_constraints.condense (system_matrix); hanging_node_constraints.condense (system_rhs); @@ -973,17 +956,16 @@ void StokesProblem::assemble_system () // linear system, we generate a // preconditioner for the // velocity-velocity matrix, i.e., - // block(0,0) in the - // system matrix. As mentioned - // above, this depends on the - // spatial dimension. Since the two - // classes described by the - // InnerPreconditioner@ :: type - // typedef have the same interface, - // we do not have to do anything - // different whether we want to use - // a sparse direct solver or an - // ILU: + // block(0,0) in the system + // matrix. As mentioned above, this + // depends on the spatial + // dimension. Since the two classes + // described by the + // InnerPreconditioner@ :: + // type typedef have the same + // interface, we do not have to do + // anything different whether we want to + // use a sparse direct solver or an ILU: std::cout << " Computing preconditioner..." << std::endl << std::flush; A_preconditioner @@ -1041,49 +1023,50 @@ void StokesProblem::solve () 1e-6*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); - // Now to the preconditioner to - // the Schur complement. As - // explained in the introduction, - // the preconditioning is done by - // a mass matrix in the pressure - // variable. It is stored in the - // $(1,1)$ block of the system - // matrix (that is not used + // Now to the preconditioner to the + // Schur complement. As explained in + // the introduction, the + // preconditioning is done by a mass + // matrix in the pressure variable. It + // is stored in the $(1,1)$ block of + // the system matrix (that is not used // anywhere else but in // preconditioning). // - // Actually, the solver needs to have the - // preconditioner in the form $P^{-1}$, so - // we need to create an inverse - // operation. Once again, we use an - // object of the class + // Actually, the solver needs to have + // the preconditioner in the form + // $P^{-1}$, so we need to create an + // inverse operation. Once again, we + // use an object of the class // InverseMatrix, which // implements the vmult // operation that is needed by the // solver. In this case, we have to - // invert the pressure mass matrix. As it - // already turned out in earlier tutorial - // programs, the inversion of a mass - // matrix is a rather cheap and + // invert the pressure mass matrix. As + // it already turned out in earlier + // tutorial programs, the inversion of + // a mass matrix is a rather cheap and // straight-forward operation (compared // to, e.g., a Laplace matrix). The CG // method with ILU preconditioning // converges in 5-10 steps, - // independently on the mesh size. This - // is precisely what we do here: We - // choose another ILU preconditioner + // independently on the mesh size. + // This is precisely what we do here: + // We choose another ILU preconditioner // and take it along to the // InverseMatrix object via the // corresponding template parameter. A // CG solver is then called within the - // vmult operation of the inverse matrix. + // vmult operation of the inverse + // matrix. // - // An alternative that is cheaper to build, - // but needs more iterations afterwards, - // would be to choose a SSOR preconditioner - // with factor 1.2. It needs about twice - // the number of iterations, but the costs - // for its generation are almost neglible. + // An alternative that is cheaper to + // build, but needs more iterations + // afterwards, would be to choose a + // SSOR preconditioner with factor + // 1.2. It needs about twice the number + // of iterations, but the costs for its + // generation are almost neglible. SparseILU preconditioner; preconditioner.initialize (system_matrix.block(1,1), SparseILU::AdditionalData()); @@ -1091,21 +1074,20 @@ void StokesProblem::solve () InverseMatrix,SparseILU > m_inverse (system_matrix.block(1,1), preconditioner); - // With the Schur complement and - // an efficient preconditioner at - // hand, we can solve the - // respective equation for the - // pressure (i.e. block 0 in the - // solution vector) in the usual + // With the Schur complement and an + // efficient preconditioner at hand, we + // can solve the respective equation + // for the pressure (i.e. block 0 in + // the solution vector) in the usual // way: cg.solve (schur_complement, solution.block(1), schur_rhs, m_inverse); - // After this first solution step, - // the hanging node constraints have - // to be distributed to the solution - // in order to achieve a consistent - // pressure field. + // After this first solution step, the + // hanging node constraints have to be + // distributed to the solution in order + // to achieve a consistent pressure + // field. hanging_node_constraints.distribute (solution); std::cout << " " @@ -1115,19 +1097,17 @@ void StokesProblem::solve () << std::endl; } - // As in step-20, we finally need - // to solve for the velocity - // equation where we plug in the - // solution to the pressure - // equation. This involves only - // objects we already know - so we - // simply multiply $p$ by $B^T$, - // subtract the right hand side and - // multiply by the inverse of - // $A$. At the end, we need to + // As in step-20, we finally need to + // solve for the velocity equation where + // we plug in the solution to the + // pressure equation. This involves only + // objects we already know - so we simply + // multiply $p$ by $B^T$, subtract the + // right hand side and multiply by the + // inverse of $A$. At the end, we need to // distribute the constraints from - // hanging nodes in order to obtain - // a constistent flow field: + // hanging nodes in order to obtain a + // constistent flow field: { system_matrix.block(0,1).vmult (tmp, solution.block(1)); tmp *= -1; @@ -1144,43 +1124,41 @@ void StokesProblem::solve () // The next function generates graphical // output. In this example, we are going to - // use the VTK file format. We attach names - // to the individual variables in the problem: - // velocity to the dim - // components of velocity and pressure - // to the pressure. + // use the VTK file format. We attach + // names to the individual variables in the + // problem: velocity to the + // dim components of velocity + // and pressure to the + // pressure. // - // Not all visualization programs - // have the ability to group - // individual vector components into - // a vector to provide vector plots; - // in particular, this holds for some - // VTK-based visualization - // programs. In this case, the - // logical grouping of components - // into vectors should already be - // described in the file containing - // the data. In other words, what we - // need to do is provide our output - // writers with a way to know which - // of the components of the finite - // element logically form a vector - // (with $d$ components in $d$ space - // dimensions) rather than letting - // them assume that we simply have a - // bunch of scalar fields. This is - // achieved using the members of the + // Not all visualization programs have the + // ability to group individual vector + // components into a vector to provide + // vector plots; in particular, this holds + // for some VTK-based visualization + // programs. In this case, the logical + // grouping of components into vectors + // should already be described in the file + // containing the data. In other words, + // what we need to do is provide our output + // writers with a way to know which of the + // components of the finite element + // logically form a vector (with $d$ + // components in $d$ space dimensions) + // rather than letting them assume that we + // simply have a bunch of scalar fields. + // This is achieved using the members of + // the // DataComponentInterpretation - // namespace: as with the filename, - // we create a vector in which the - // first dim components - // refer to the velocities and are - // given the tag + // namespace: as with the filename, we + // create a vector in which the first + // dim components refer to the + // velocities and are given the tag // DataComponentInterpretation::component_is_part_of_vector; // we finally push one tag // DataComponentInterpretation::component_is_scalar - // to describe the grouping of the - // pressure variable. + // to describe the grouping of the pressure + // variable. // The rest of the function is then // the same as in step-20. @@ -1216,19 +1194,18 @@ StokesProblem::output_results (const unsigned int refinement_cycle) const // @sect4{StokesProblem::refine_mesh} - // This is the last interesting function - // of the StokesProblem class. + // This is the last interesting function of + // the StokesProblem class. // As indicated by its name, it takes the - // solution to the problem and - // refines the mesh where this is - // needed. The procedure is the same - // as in the respective step in - // step-6, with the exception that - // we base the refinement only on the - // change in pressure, i.e., we call - // the Kelly error estimator with a - // mask object. Additionally, we do - // not coarsen the grid again: + // solution to the problem and refines the + // mesh where this is needed. The procedure + // is the same as in the respective step in + // step-6, with the exception that we base + // the refinement only on the change in + // pressure, i.e., we call the Kelly error + // estimator with a mask + // object. Additionally, we do not coarsen + // the grid again: template void StokesProblem::refine_mesh () @@ -1253,28 +1230,25 @@ StokesProblem::refine_mesh () // @sect4{StokesProblem::run} - // The last step in the Stokes class - // is, as usual, the function that generates - // the initial grid and calls the other + // The last step in the Stokes class is, as + // usual, the function that generates the + // initial grid and calls the other // functions in the respective order. // - // We start off with a rectangle of - // size $4 \times 1$ (in 2d) or $4 - // \times 1 \times 1$ (in 3d), placed - // in $R^2/R^3$ as + // We start off with a rectangle of size $4 + // \times 1$ (in 2d) or $4 \times 1 \times + // 1$ (in 3d), placed in $R^2/R^3$ as // $(-2,2)\times(-1,0)$ or // $(-2,2)\times(0,1)\times(-1,1)$, - // respectively. It is natural to - // start with equal mesh size in each - // direction, so we subdivide the - // initial rectangle four times in - // the first coordinate direction. To - // limit the scope of the variables - // involved in the creation of the - // mesh to the range where we - // actually need them, we put the - // entire block between a pair of - // braces: + // respectively. It is natural to start + // with equal mesh size in each direction, + // so we subdivide the initial rectangle + // four times in the first coordinate + // direction. To limit the scope of the + // variables involved in the creation of + // the mesh to the range where we actually + // need them, we put the entire block + // between a pair of braces: template void StokesProblem::run () { @@ -1309,22 +1283,19 @@ void StokesProblem::run () cell->face(f)->set_all_boundary_indicators(1); - // We then apply an initial - // refinement before solving for - // the first time. In 3D, there are - // going to be more degrees of - // freedom, so we refine less - // there: + // We then apply an initial refinement + // before solving for the first time. In + // 3D, there are going to be more degrees + // of freedom, so we refine less there: triangulation.refine_global (4-dim); - // As first seen in step-6, we - // cycle over the different - // refinement levels and refine - // (except for the first cycle), + // As first seen in step-6, we cycle over + // the different refinement levels and + // refine (except for the first cycle), // setup the degrees of freedom and - // matrices, assemble, solve and - // create output: - for (unsigned int refinement_cycle = 0; refinement_cycle<4; + // matrices, assemble, solve and create + // output: + for (unsigned int refinement_cycle = 0; refinement_cycle<6; ++refinement_cycle) { std::cout << "Refinement cycle " << refinement_cycle << std::endl; @@ -1359,7 +1330,7 @@ int main () { deallog.depth_console (0); - StokesProblem<3> flow_problem(1); + StokesProblem<2> flow_problem(1); flow_problem.run (); } catch (std::exception &exc) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 369abc7b6f..a224297f41 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -406,12 +406,11 @@ namespace LinearSolvers // SIAM J. Numer. Anal., 31 (1994), // pp. 1352-1367). // - // Replacing P by $\tilde{P}$ does - // not change the situation - // dramatically. The product $P^{-1} A$ - // will still be close to a matrix with - // eigenvalues 0 and 1, which lets us - // hope to be able to get a number of + // Replacing P by $\tilde{P}$ + // keeps that spirit alive: the product + // $P^{-1} A$ will still be close to a + // matrix with eigenvalues 1, which lets + // us hope to be able to get a number of // GMRES iterations that does not depend // on the problem size. // @@ -1129,7 +1128,7 @@ void BoussinesqFlowProblem::setup_dofs () // step-22, we choose to create the pattern // not as in the first few tutorial // programs, but by using the blocked - // version of CompressedSetSparsityPattern. + // version of CompressedSimpleSparsityPattern. // The reason for doing this is mainly // memory, that is, the SparsityPattern // class would consume too much memory when @@ -1139,56 +1138,59 @@ void BoussinesqFlowProblem::setup_dofs () // So, we first release the memory stored // in the matrices, then set up an object // of type - // BlockCompressedSetSparsityPattern + // BlockCompressedSimpleSparsityPattern // consisting of $2\times 2$ blocks (for // the Stokes system matrix and // preconditioner) or - // CompressedSparsityPattern (for the - // temperature part). We then fill these - // sparsity patterns with the nonzero + // CompressedSimpleSparsityPattern (for + // the temperature part). We then fill + // these objects with the nonzero // pattern, taking into account that for // the Stokes system matrix, there are no // entries in the pressure-pressure block // (but all velocity vector components // couple with each other and with the // pressure). Similarly, in the Stokes - // preconditioner matrix, only the diagonal - // blocks are nonzero, since we use the - // vector Laplacian as discussed in the - // introduction. This operator only couples - // each vector component of the Laplacian - // with itself, but not with the other - // vector components. (Application of the + // preconditioner matrix, only the + // diagonal blocks are nonzero, since we + // use the vector Laplacian as discussed + // in the introduction. This operator + // only couples each vector component of + // the Laplacian with itself, but not + // with the other vector + // components. (Application of the // constraints resulting from the no-flux // boundary conditions will couple vector // components at the boundary again, // however.) // - // When generating the sparsity pattern, we - // directly apply the constraints from + // When generating the sparsity pattern, + // we directly apply the constraints from // hanging nodes and no-flux boundary // conditions. This approach was already // used in step-27, but is different from - // the one in early tutorial programs where - // we first built the original sparsity - // pattern and only then added the entries - // resulting from constraints. The reason - // for doing so is that later during - // assembly we are going to distribute the - // constraints immediately when - // transferring local to global - // dofs. Consequently, there will be no - // data written at positions of constrained - // degrees of freedom, so we can let the - // DoFTools::make_sparsity_pattern function - // omit these entries by setting the last - // boolean flag to false. Once - // the sparsity pattern is ready, we can - // use it to initialize the Trilinos - // matrices. Note that the Trilinos - // matrices store the sparsity pattern - // internally, so there is no need to keep - // the sparsity pattern around after the + // the one in early tutorial programs + // where we first built the original + // sparsity pattern and only then added + // the entries resulting from + // constraints. The reason for doing so + // is that later during assembly we are + // going to distribute the constraints + // immediately when transferring local to + // global dofs. Consequently, there will + // be no data written at positions of + // constrained degrees of freedom, so we + // can let the + // DoFTools::make_sparsity_pattern + // function omit these entries by setting + // the last boolean flag to + // false. Once the sparsity + // pattern is ready, we can use it to + // initialize the Trilinos + // matrices. Since the Trilinos matrices + // store the sparsity pattern internally, + // there is no need to keep the sparsity + // pattern around after the // initialization of the matrix. stokes_block_sizes.resize (2); stokes_block_sizes[0] = n_u; @@ -1196,7 +1198,7 @@ void BoussinesqFlowProblem::setup_dofs () { stokes_matrix.clear (); - BlockCompressedSetSparsityPattern csp (2,2); + BlockCompressedSimpleSparsityPattern csp (2,2); csp.block(0,0).reinit (n_u, n_u); csp.block(0,1).reinit (n_u, n_p); @@ -1225,7 +1227,7 @@ void BoussinesqFlowProblem::setup_dofs () Mp_preconditioner.reset (); stokes_preconditioner_matrix.clear (); - BlockCompressedSetSparsityPattern csp (2,2); + BlockCompressedSimpleSparsityPattern csp (2,2); csp.block(0,0).reinit (n_u, n_u); csp.block(0,1).reinit (n_u, n_p); @@ -1256,14 +1258,14 @@ void BoussinesqFlowProblem::setup_dofs () // discretization) follows the generation // of the Stokes matrix – except // that it is much easier here since we - // do not need to take care of any - // blocks or coupling between components: + // do not need to take care of any blocks + // or coupling between components: { temperature_mass_matrix.clear (); temperature_stiffness_matrix.clear (); temperature_matrix.clear (); - CompressedSetSparsityPattern csp (n_T, n_T); + CompressedSimpleSparsityPattern csp (n_T, n_T); DoFTools::make_sparsity_pattern (temperature_dof_handler, csp, temperature_constraints, false); diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 20871d259a..61d9629c11 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -633,7 +633,7 @@ void BoussinesqFlowProblem::setup_dofs () { stokes_matrix.clear (); - BlockCompressedSetSparsityPattern csp (2,2); + BlockCompressedSimpleSparsityPattern csp (2,2); csp.block(0,0).reinit (n_u, n_u); csp.block(0,1).reinit (n_u, n_p); @@ -662,7 +662,7 @@ void BoussinesqFlowProblem::setup_dofs () Mp_preconditioner.reset (); stokes_preconditioner_matrix.clear (); - BlockCompressedSetSparsityPattern csp (2,2); + BlockCompressedSimpleSparsityPattern csp (2,2); csp.block(0,0).reinit (n_u, n_u); csp.block(0,1).reinit (n_u, n_p); @@ -697,7 +697,7 @@ void BoussinesqFlowProblem::setup_dofs () temperature_stiffness_matrix.clear (); temperature_matrix.clear (); - CompressedSetSparsityPattern csp (n_T, n_T); + CompressedSimpleSparsityPattern csp (n_T, n_T); DoFTools::make_sparsity_pattern (temperature_dof_handler, csp, temperature_constraints, false);