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*}
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
+— 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
<< 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);
}
projects / dealii.git / commitdiff