]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Revamp computing the sparsity pattern.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 27 Apr 2011 15:28:10 +0000 (15:28 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 27 Apr 2011 15:28:10 +0000 (15:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@23655 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-46/doc/intro.dox
deal.II/examples/step-46/step-46.cc

index 51bb22cc12f84b7c31b9c646ecdaa8c7315638d7..4eb097985b84c9a713a0ec695d8372b8f926be70 100644 (file)
@@ -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$.
 
+
+<h4>Dealing with the interface terms</h4>
+
 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 <i>all</i> 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<fe_collection.n_components(); ++c)
+      for (unsigned int d=0; d<fe_collection.n_components(); ++d)
+       if (((c<dim+1) && (d<dim+1)
+            && !((c==dim) && (d==dim)))
+           ||
+           ((c>=dim+1) && (d>=dim+1)))
+         cell_coupling[c][d] = DoFTools::Coupling::always;
+@endcode
+Here, we have used the fact that the first <code>dim</code> components of the
+finite element are the velocities, then the pressure, and then the
+<code>dim</code> 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<fe_collection.n_components(); ++c)
+      for (unsigned int d=0; d<fe_collection.n_components(); ++d)
+       if ((c>=dim+1) && (d<dim+1))
+         face_coupling[c][d] = DoFTools::Coupling::always;
+@endcode
+In other words, all displacement test functions (components
+<code>c@>=dim+1</code>) 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
+&mdash; good enough for now.
+
+
+
+<h4>Velocity boundary conditions on the interface</h4>
 
 The second difficulty is that while we know how to enforce a zero
 velocity or stress on the external boundary (using
index bea373cf3a41b6369566cbf95ad27db182c42c17..58a934a62ef44eeda63d884e1c2a3e2a445fcbcb 100644 (file)
@@ -501,15 +501,36 @@ FluidStructureProblem<dim>::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<fe_collection.n_components(); ++c)
+      for (unsigned int d=0; d<fe_collection.n_components(); ++d)
+       {
+         if (((c<dim+1) && (d<dim+1)
+              && !((c==dim) && (d==dim)))
+             ||
+             ((c>=dim+1) && (d>=dim+1)))
+           cell_coupling[c][d] = DoFTools::Coupling::always;
+
+         if ((c>=dim+1) && (d<dim+1))
+           face_coupling[c][d] = DoFTools::Coupling::always;
+       }
+    
+    DoFTools::make_flux_sparsity_pattern (dof_handler, csp,
+                                         cell_coupling, face_coupling);
+    constraints.condense (csp);
     sparsity_pattern.copy_from (csp);
   }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.