From 5fdc999995dbccce05586b2c906eb7516b263aea Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 20 Apr 2011 14:21:36 +0000 Subject: [PATCH] Finish writing intro and results sections. Re-scale pictures. git-svn-id: https://svn.dealii.org/trunk@23616 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-46/doc/intro.dox | 85 +++++++++++++++++++ deal.II/examples/step-46/doc/results.dox | 101 +++++++++++++++++++++++ deal.II/examples/step-46/step-46.cc | 88 +++++++++++++++++--- 3 files changed, 261 insertions(+), 13 deletions(-) diff --git a/deal.II/examples/step-46/doc/intro.dox b/deal.II/examples/step-46/doc/intro.dox index 604242143e..7d7fec1019 100644 --- a/deal.II/examples/step-46/doc/intro.dox +++ b/deal.II/examples/step-46/doc/intro.dox @@ -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$. + + +

Linear solvers

+ 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 results section. + +

Mesh refinement

+ +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 stokes_estimated_error_per_cell (triangulation.n_active_cells()); + Vector elasticity_estimated_error_per_cell (triangulation.n_active_cells()); + + std::vector stokes_component_mask (dim+1+dim, false); + for (unsigned int d=0; d::estimate (dof_handler, + face_q_collection, + typename FunctionMap::type(), + solution, + stokes_estimated_error_per_cell, + stokes_component_mask); + + std::vector elasticity_component_mask (dim+1+dim, false); + for (unsigned int d=0; d::estimate (dof_handler, + face_q_collection, + typename FunctionMap::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 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. + diff --git a/deal.II/examples/step-46/doc/results.dox b/deal.II/examples/step-46/doc/results.dox index 6796b94fff..d37bf00fe3 100644 --- a/deal.II/examples/step-46/doc/results.dox +++ b/deal.II/examples/step-46/doc/results.dox @@ -1,10 +1,100 @@

Results

+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: + + + + + + + + + + + + +
+ @image html step-46.velocity-magnitude.png + +

+ Magnitude of the fluid velocity +

+ 		+ @image html step-46.pressure.png + +

+ 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. +

+
+ @image html step-46.velocity.png + +

+ Fluid velocity. +

+ 		+ @image html step-46.displacement.png + +

+ Solid displacement. +

+
+ +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.

Possibilities for extensions

+

Linear solvers and preconditioners

+ 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 @@ -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. + + +

Refinement indicators

+ +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. diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index da3e82cf31..6edf976d59 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -188,8 +188,9 @@ RightHandSide::vector_value (const Point &p, template -FluidStructureProblem::FluidStructureProblem (const unsigned int stokes_degree, - const unsigned int elasticity_degree) +FluidStructureProblem:: +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::setup_dofs () } + template void FluidStructureProblem::setup_subdomains () @@ -349,6 +351,7 @@ FluidStructureProblem::setup_subdomains () } + template void FluidStructureProblem::assemble_system () { @@ -371,17 +374,21 @@ void FluidStructureProblem::assemble_system () const QGauss face_quadrature(std::max (stokes_degree+2, elasticity_degree+2)); - FEFaceValues stokes_fe_face_values (stokes_fe, face_quadrature, - update_JxW_values | - update_normal_vectors | - update_gradients); - FEFaceValues elasticity_fe_face_values (elasticity_fe, face_quadrature, - update_values); - FESubfaceValues stokes_fe_subface_values (stokes_fe, face_quadrature, + FEFaceValues stokes_fe_face_values (stokes_fe, + face_quadrature, + update_JxW_values | + update_normal_vectors | + update_gradients); + FEFaceValues elasticity_fe_face_values (elasticity_fe, + face_quadrature, + update_values); + FESubfaceValues stokes_fe_subface_values (stokes_fe, + face_quadrature, update_JxW_values | update_normal_vectors | update_gradients); - FESubfaceValues elasticity_fe_subface_values (elasticity_fe, face_quadrature, + FESubfaceValues 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::assemble_system () std::vector local_dof_indices; std::vector neighbor_dof_indices (stokes_dofs_per_cell); - const RightHandSide right_hand_side; + const RightHandSide right_hand_side; const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); @@ -700,10 +707,65 @@ FluidStructureProblem::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::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell, ++cell_index) + for (unsigned int f=0; f::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); -- 2.39.5