]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Finish writing intro and results sections. Re-scale pictures.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 20 Apr 2011 14:21:36 +0000 (14:21 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 20 Apr 2011 14:21:36 +0000 (14:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@23616 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-46/doc/intro.dox
deal.II/examples/step-46/doc/results.dox
deal.II/examples/step-46/doc/step-46.displacement.png
deal.II/examples/step-46/doc/step-46.pressure.png
deal.II/examples/step-46/doc/step-46.velocity-magnitude.png
deal.II/examples/step-46/doc/step-46.velocity.png
deal.II/examples/step-46/step-46.cc

index 604242143e7d8fe425f3ff591ee7dc783c24e46b..7d7fec1019d9722961f70fc8247204d4d3fc3c62 100644 (file)
@@ -437,9 +437,94 @@ flow, yielding the implicit no-stress condition $(2\eta
 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 &mdash; 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.
+
index 6796b94fffafe5f6a66fc24394b530c50bb32815..d37bf00fe3d633086153159bdb359906f0e187e3 100644 (file)
 <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 &mdash; in particular one that scales well and
 will also work for realistic 3d problems. This shouldn't actually be
@@ -70,3 +160,14 @@ separately. These are well known, however: for Stokes, we can use the
 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.
index 84596915ef2785fe5213443b5587195890814d21..a41c7039c099a1948f9f81bfba27c592d8d95bfb 100644 (file)
Binary files a/deal.II/examples/step-46/doc/step-46.displacement.png and b/deal.II/examples/step-46/doc/step-46.displacement.png differ
index 353a9b90ba9153eeb3160265743bc29bfbb767e1..bc117c72b64988114401dad9a1f563009d6c8631 100644 (file)
Binary files a/deal.II/examples/step-46/doc/step-46.pressure.png and b/deal.II/examples/step-46/doc/step-46.pressure.png differ
index 5ef29b25447d3cc99af69998d9ddb1a6045d4ff8..ef3d5a5eaff6578379d16643205620473fd9ce47 100644 (file)
Binary files a/deal.II/examples/step-46/doc/step-46.velocity-magnitude.png and b/deal.II/examples/step-46/doc/step-46.velocity-magnitude.png differ
index 057a60ad93d706a4ead56e002d577482f9461b64..3a432e16ec85e32440c4719e953aa7c6cf11031e 100644 (file)
Binary files a/deal.II/examples/step-46/doc/step-46.velocity.png and b/deal.II/examples/step-46/doc/step-46.velocity.png differ
index da3e82cf31f1b51cedeb61ecb337aedc5a721418..6edf976d59e57ab6225a5e2c0dbc1e20a494869c 100644 (file)
@@ -188,8 +188,9 @@ RightHandSide<dim>::vector_value (const Point<dim> &p,
 
 
 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),
@@ -329,6 +330,7 @@ void FluidStructureProblem<dim>::setup_dofs ()
 }
 
 
+
 template <int dim>
 void
 FluidStructureProblem<dim>::setup_subdomains ()
@@ -349,6 +351,7 @@ FluidStructureProblem<dim>::setup_subdomains ()
 }
 
 
+
 template <int dim>
 void FluidStructureProblem<dim>::assemble_system ()
 {
@@ -371,17 +374,21 @@ 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;
@@ -395,7 +402,7 @@ void FluidStructureProblem<dim>::assemble_system ()
   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);
@@ -700,10 +707,65 @@ FluidStructureProblem<dim>::refine_mesh ()
                                       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);

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.