The conditions on the interface between the two domains has been
discussed above already.
+For simplicity, we choose the material parameters to be $\eta=\lambda=\mu=1$.
+
+
+<h4>Linear solvers</h4>
+
This program is primarily intended to show how to deal with different
physics in different parts of the domain, and how to implement such
models in deal.II. As a consequence, we won't bother coming up with a
good solver: we'll just use the SparseDirectUMFPACK class which always
works, even if not with optimal complexity. We will, however, comment
on possible other solvers in the <a href="#Results">results</a> section.
+
+<h4>Mesh refinement</h4>
+
+One of the trickier aspects of this program is how to estimate the
+error. Because it works on almost any program, we'd like to use the
+KellyErrorEstimator, and we can relatively easily do that here as well using
+code like the following:
+@code
+ Vector<float> stokes_estimated_error_per_cell (triangulation.n_active_cells());
+ Vector<float> elasticity_estimated_error_per_cell (triangulation.n_active_cells());
+
+ std::vector<bool> stokes_component_mask (dim+1+dim, false);
+ for (unsigned int d=0; d<dim; ++d)
+ stokes_component_mask[d] = true;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ face_q_collection,
+ typename FunctionMap<dim>::type(),
+ solution,
+ stokes_estimated_error_per_cell,
+ stokes_component_mask);
+
+ std::vector<bool> elasticity_component_mask (dim+1+dim, false);
+ for (unsigned int d=0; d<dim; ++d)
+ elasticity_component_mask[dim+1+d] = true;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ face_q_collection,
+ typename FunctionMap<dim>::type(),
+ solution,
+ elasticity_estimated_error_per_cell,
+ elasticity_component_mask);
+@endcode
+This gives us two sets of error indicators for each cell. We would then
+somehow combine them into one for mesh refinement, for example using something
+like the following (note that we normalize the squared error indicator in the
+two vectors because error quantities have physical units that do not match in
+the current situation, leading to error indicators that may differ by orders
+of magnitude between the two subdomains):
+@code
+ stokes_estimated_error_per_cell /= stokes_estimated_error_per_cell.l2_norm();
+ elasticity_estimated_error_per_cell /= elasticity_estimated_error_per_cell.l2_norm();
+
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ estimated_error_per_cell += stokes_estimated_error_per_cell;
+ estimated_error_per_cell += elasticity_estimated_error_per_cell;
+@endcode
+(In the code, we actually weigh the error indicators 4:1 in favor of the ones
+computed on the Stokes subdomain since refinement is otherwise heavily biased
+towards the elastic subdomain, but this is just a technicality.)
+
+While this principle is sound, it doesn't quite work as expected. The reason
+is that the KellyErrorEstimator class computes error indicators by integrating
+the jump in the solution's gradient around the faces of each cell. This jump
+is likely to be very large at the locations where the solution is
+discontinuous and extended by zero; it also doesn't become smaller as the mesh
+is refined. The KellyErrorEstimator class can't just ignore the interface
+because it essentially only sees an hp::DoFHandler where the element type
+changes from one cell to another — precisely the thing that the
+hp::DoFHandler was designed for, the interface in the current program looks no
+different than the interfaces in step-27, for example, and certainly no less
+legitimate. Be that as it may, the end results is that there is a layer of
+cells on both sides of the interface between the two subdomains where error
+indicators are irrationally large. Consequently, most of the mesh refinement
+is focused on the interface.
+
+This clearly wouldn't happen if we had a refinement indicator that actually
+understood something about the problem and simply ignore the interface between
+subdomains when integrating jump terms.
+On the other hand, this program is
+about showing how to represent problems where we have different physics in
+different subdomains, not about the peculiarities of the KellyErrorEstimator,
+and so we resort to the big hammer called "heuristics": we simply set the
+error indicators of cells at the interface to zero. This cuts off the spikes
+in the error indicators. At first sight one would also think that it prevents
+the mesh from being refined at the interface, but the requirement that
+neighboring cells may only differ by one level of refinement will still lead
+to a reasonably refined mesh.
+
+While this is clearly a suboptimal solution, it works for now and leaves room
+for future improvement.
+
<a name="Results"></a>
<h1>Results</h1>
+When running the program, you should get output like the following:
+@code
+Refinement cycle 0
+ Number of active cells: 64
+ Number of degrees of freedom: 531
+ Assembling...
+ Solving...
+Refinement cycle 1
+ Number of active cells: 184
+ Number of degrees of freedom: 1565
+ Assembling...
+ Solving...
+Refinement cycle 2
+ Number of active cells: 568
+ Number of degrees of freedom: 4849
+ Assembling...
+ Solving...
+Refinement cycle 3
+ Number of active cells: 1528
+ Number of degrees of freedom: 11605
+ Assembling...
+ Solving...
+Refinement cycle 4
+ Number of active cells: 3628
+ Number of degrees of freedom: 22309
+ Assembling...
+ Solving...
+Refinement cycle 5
+ Number of active cells: 8440
+ Number of degrees of freedom: 45139
+ Assembling...
+ Solving...
+@endcode
+
+The results are easily visualized:
+
+<TABLE WIDTH="60%" ALIGN="center">
+ <tr valign="top">
+ <td valign="top" align="center">
+ @image html step-46.velocity-magnitude.png
+
+ <p align="center">
+ Magnitude of the fluid velocity
+ </p>
+ </td>
+
+ <td valign="top" align="center">
+ @image html step-46.pressure.png
+
+ <p align="center">
+ Fluid pressure. The dynamic range has been truncated to cut off the
+ pressure singularities at the top left and right corners of the domain
+ as well as the top corners of the solid that forms re-entrant corners
+ into the fluid domain.
+ </p>
+ </td>
+ </tr>
+ <tr valign="top">
+ <td valign="top" align="center">
+ @image html step-46.velocity.png
+
+ <p align="center">
+ Fluid velocity.
+ </p>
+ </td>
+
+ <td valign="top" align="center">
+ @image html step-46.displacement.png
+
+ <p align="center">
+ Solid displacement.
+ </p>
+ </td>
+ </tr>
+</table>
+
+In all figures, we have applied a mask to only show the original field, not
+the one extended by zero: for example, to plot the pressure, we have selected
+that part of the domain where the magnitude of the velocity is greater than
+$10^{-7}$.
+
+The plots are easily interpreted: as the flow drives down on the left side and
+up on the right side of the upright part of the solid, it produces a shear
+force that pulls the left side down and the right side up. An additional part
+force comes from the pressure, which bears down on the left side of the top
+and pulls up on the right side. Both forces yield a net torque on the solid
+that bends it to the left, as confirmed by the plot of the displacement
+vectors.
<a name="extensions"></a>
<h3>Possibilities for extensions</h3>
+<h4>Linear solvers and preconditioners</h4>
+
An obvious place to improve the program would be to use a more
sophisticated solver — in particular one that scales well and
will also work for realistic 3d problems. This shouldn't actually be
preconditioner discussed in the results section of step-22; for
elasticity, a good preconditioner would be a single V-cycle of a
geometric or algebraic multigrid.
+
+
+<h4>Refinement indicators</h4>
+
+As mentioned in the introduction, the refinement indicator we use for this
+program is rather ad hoc. A better one would understand that the jump in the
+gradient of the solution across the interface is not indicative of the error
+but to be expected and ignore the interface when integrating the jump
+terms. Nevertheless, this is not what the KellyErrorEstimator class
+does. Consequently, an obvious possibility for improving the program would be
+to implement a better refinement criterion.
template <int dim>
-FluidStructureProblem<dim>::FluidStructureProblem (const unsigned int stokes_degree,
- const unsigned int elasticity_degree)
+FluidStructureProblem<dim>::
+FluidStructureProblem (const unsigned int stokes_degree,
+ const unsigned int elasticity_degree)
:
stokes_degree (stokes_degree),
elasticity_degree (elasticity_degree),
}
+
template <int dim>
void
FluidStructureProblem<dim>::setup_subdomains ()
}
+
template <int dim>
void FluidStructureProblem<dim>::assemble_system ()
{
const QGauss<dim-1> face_quadrature(std::max (stokes_degree+2,
elasticity_degree+2));
- FEFaceValues<dim> stokes_fe_face_values (stokes_fe, face_quadrature,
- update_JxW_values |
- update_normal_vectors |
- update_gradients);
- FEFaceValues<dim> elasticity_fe_face_values (elasticity_fe, face_quadrature,
- update_values);
- FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe, face_quadrature,
+ FEFaceValues<dim> stokes_fe_face_values (stokes_fe,
+ face_quadrature,
+ update_JxW_values |
+ update_normal_vectors |
+ update_gradients);
+ FEFaceValues<dim> elasticity_fe_face_values (elasticity_fe,
+ face_quadrature,
+ update_values);
+ FESubfaceValues<dim> stokes_fe_subface_values (stokes_fe,
+ face_quadrature,
update_JxW_values |
update_normal_vectors |
update_gradients);
- FESubfaceValues<dim> elasticity_fe_subface_values (elasticity_fe, face_quadrature,
+ FESubfaceValues<dim> elasticity_fe_subface_values (elasticity_fe,
+ face_quadrature,
update_values);
const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell;
std::vector<unsigned int> local_dof_indices;
std::vector<unsigned int> neighbor_dof_indices (stokes_dofs_per_cell);
- const RightHandSide<dim> right_hand_side;
+ const RightHandSide<dim> right_hand_side;
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
elasticity_estimated_error_per_cell,
elasticity_component_mask);
- stokes_estimated_error_per_cell /= stokes_estimated_error_per_cell.linfty_norm();
- elasticity_estimated_error_per_cell /= elasticity_estimated_error_per_cell.linfty_norm();
+ stokes_estimated_error_per_cell /= 0.25 * stokes_estimated_error_per_cell.l2_norm();
+ elasticity_estimated_error_per_cell /= elasticity_estimated_error_per_cell.l2_norm();
estimated_error_per_cell += stokes_estimated_error_per_cell;
estimated_error_per_cell += elasticity_estimated_error_per_cell;
+
+ {
+ unsigned int cell_index = 0;
+ for (typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active();
+ cell != dof_handler.end(); ++cell, ++cell_index)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell_is_in_solid_domain (cell))
+ {
+ if ((cell->at_boundary(f) == false)
+ &&
+ (((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == false)
+ &&
+ cell_is_in_fluid_domain (cell->neighbor(f)))
+ ||
+ ((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == true)
+ &&
+ (cell_is_in_fluid_domain (cell->neighbor_child_on_subface
+ (f, 0))))
+ ||
+ (cell->neighbor_is_coarser(f)
+ &&
+ cell_is_in_fluid_domain(cell->neighbor(f)))
+ ))
+ estimated_error_per_cell(cell_index) = 0;
+ }
+ else
+ {
+ if ((cell->at_boundary(f) == false)
+ &&
+ (((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == false)
+ &&
+ cell_is_in_solid_domain (cell->neighbor(f)))
+ ||
+ ((cell->neighbor(f)->level() == cell->level())
+ &&
+ (cell->neighbor(f)->has_children() == true)
+ &&
+ (cell_is_in_solid_domain (cell->neighbor_child_on_subface
+ (f, 0))))
+ ||
+ (cell->neighbor_is_coarser(f)
+ &&
+ cell_is_in_solid_domain(cell->neighbor(f)))
+ ))
+ estimated_error_per_cell(cell_index) = 0;
+ }
+ }
+
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
0.3, 0.0);
projects / dealii.git / commitdiff