From: bangerth Date: Wed, 27 Apr 2011 15:28:10 +0000 (+0000) Subject: Revamp computing the sparsity pattern. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ad1f43effa6ac55b9b1f392aa1e8ae5b00979002;p=dealii-svn.git Revamp computing the sparsity pattern. git-svn-id: https://svn.dealii.org/trunk@23655 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-46/doc/intro.dox b/deal.II/examples/step-46/doc/intro.dox index 51bb22cc12..4eb097985b 100644 --- a/deal.II/examples/step-46/doc/intro.dox +++ b/deal.II/examples/step-46/doc/intro.dox @@ -279,11 +279,14 @@ mapping the problem on to the hp framework makes sense: More specifically, in the program we have to address the following points: - Implementing the bilinear form, and in particular dealing with the - interface term. + interface term, both in the matrix and the sparsity pattern. - Implementing Dirichlet boundary conditions on the external and internal parts of the boundaryies $\partial\Omega_f,\partial\Omega_s$. + +

Dealing with the interface terms

+ Let us first discuss implementing the bilinear form, which at the discrete level we recall to be @f{multline*} @@ -338,15 +341,102 @@ integrate over parts of a face. Take a look at the implementation below on how to deal with this. As an additional complication, the matrix entries that result from this term -need to be added to the sparsity pattern of the matrix somehow. This, however, -is not too difficult: the term couples degrees of freedom from two adjacent -cells along a face, which is exactly the kind of thing one would do in -discontinuous Galerkin schemes for which the function -DoFTools::make_flux_sparsity_pattern was written. In the current context it -computes a superset of matrix entries compared to the usual -DoFTools::make_sparsity_pattern: it will also add the ones for faces where the -solution is discontinuous, which is exactlty the case at the interface between -the two subdomains. +need to be added to the sparsity pattern of the matrix somehow. This is the +realm of various functions in the DoFTools namespace like +DoFTools::make_sparsity_pattern and +DoFTools::make_flux_sparsity_pattern. Essentially, what these functions do is +simulate what happens during assembly of the system matrix: whenever assembly +would write a nonzero entry into the global matrix, the functions in DoFTools +would add an entry to the sparsity pattern. We could therefore do the +following: let DoFTools::make_sparsity_pattern add all those entries to the +sparsity pattern that arise from the regular cell-by-cell integration, and +then do the same by hand that arise from the interface terms. If you look at +the implementation of the interface integrals in the program below, it should +be obvious how to do that and would require no more than maybe 100 lines of +code at most. + +But we're lazy people: the interface term couples degrees of freedom from two +adjacent cells along a face, which is exactly the kind of thing one would do +in discontinuous Galerkin schemes for which the function +DoFTools::make_flux_sparsity_pattern was written. This is a superset of matrix +entries compared to the usual DoFTools::make_sparsity_pattern: it will also +add all entries that result from computing terms coupling the degrees of +freedom from both sides of all faces. Unfortunately, for the simplest version +of this function, this is a pretty big superset. Consider for example the +following mesh with two cells and a $Q_1$ finite element: +@code + 2---3---5 + | | | + 0---1---4 +@endcode +Here, the sparsity pattern produced by DoFTools::make_sparsity_pattern will +only have entries for degrees of freedom that couple on a cell. However, it +will not have sparsity pattern entries $(0,4),(0,5),(2,4),(2,5)$. The sparsity +pattern generated by DoFTools::make_flux_sparsity_pattern will have these +entries, however: it assumes that you want to build a sparsity pattern for a +bilinear form that couples all degrees of freedom from adjacent +cells. This is not what we want: our interface term acts only on a small +subset of cells, and we certainly don't need all the extra couplings between +two adjacent fluid cells, or two adjacent solid cells. Furthermore, the fact that we +use higher order elements means that we would really generate many many more +entries than we actually need: on the coarsest mesh, in 2d, 44,207 nonzero +entries instead of 16,635 for DoFTools::make_sparsity_pattern, leading to +plenty of zeros in the matrix we later build (of course, the 16,635 are not +enough since they don't include the interface entries). This ratio would be +even worse in 3d. + +So being extremely lazy comes with a cost: too many entries in the matrix. But +we can get away with being moderately lazy: there is a variant of +DoFTools::make_flux_sparsity_pattern that allows us +to specify which vector components of the finite element couple with which +other components, both in cell terms as well as in face terms. For cells that +are in the solid subdomain, we couple all displacements with each other; for +fluid cells, all velocities with all velocities and the pressure, but not the +pressure with itself. Since no cell has both sets of +variables, there is no need to distinguish between the two kinds of cells, so +we can write the mask like this: +@code + Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(), + fe_collection.n_components()); + + for (unsigned int c=0; c=dim+1) && (d>=dim+1))) + cell_coupling[c][d] = DoFTools::Coupling::always; +@endcode +Here, we have used the fact that the first dim components of the +finite element are the velocities, then the pressure, and then the +dim displacements. (We could as well have stated that the +velocities/pressure also couple with the displacements since no cell ever has +both sets of variables.) On the other hand, the interface terms require a mask +like this: +@code + Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(), + fe_collection.n_components()); + + for (unsigned int c=0; c=dim+1) && (dc@>=dim+1) couple with all velocity and pressure shape functions +on the other side of an interface. This is not entirely true, though close: in +fact, the exact form of the interface term only those pressure displacement +shape functions that are indeed nonzero on the common interface, which is not +true for all shape functions; on the other hand, it really couples all +velocities (since the integral involves gradients of the velocity shape +functions, which are all nonzero on all faces of the cell). However, the mask we +build above, is not capable of these subtleties. Nevertheless, through these +masks we manage to get the number of sparsity pattern entries down to 21,028 +— good enough for now. + + + +

Velocity boundary conditions on the interface

The second difficulty is that while we know how to enforce a zero velocity or stress on the external boundary (using diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index bea373cf3a..58a934a62e 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -501,15 +501,36 @@ FluidStructureProblem::setup_dofs () << dof_handler.n_dofs() << std::endl; - // The rest of this function is standard: - // Create a sparsity pattern and use it to - // initialize the matrix; then also set - // vectors to their correct sizes. + // In the rest of this function we create a + // sparsity pattern as discussed + // extensively in the introduction, and use + // it to initialize the matrix; then also + // set vectors to their correct sizes: { CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(), dof_handler.n_dofs()); - DoFTools::make_flux_sparsity_pattern (dof_handler, csp, constraints, false); + Table<2,DoFTools::Coupling> cell_coupling (fe_collection.n_components(), + fe_collection.n_components()); + Table<2,DoFTools::Coupling> face_coupling (fe_collection.n_components(), + fe_collection.n_components()); + + for (unsigned int c=0; c=dim+1) && (d>=dim+1))) + cell_coupling[c][d] = DoFTools::Coupling::always; + + if ((c>=dim+1) && (d