From: Wolfgang Bangerth Date: Mon, 25 Apr 2011 00:00:59 +0000 (+0000) Subject: Finish writing documentation. X-Git-Tag: v8.0.0~4113 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4f37cbed1bb66e3337e33f50741d5b0e6c8f0c0d;p=dealii.git Finish writing documentation. git-svn-id: https://svn.dealii.org/trunk@23644 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-46/step-46.cc b/deal.II/examples/step-46/step-46.cc index bf2e6c8202..dd9e9e2f6c 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -249,11 +249,25 @@ RightHandSide::vector_value (const Point &p, // @sect3{The FluidStructureProblem implementation} + // @sect4{Constructors and helper functions} + // Let's now get to the implementation of the // primary class of this program. The first - // few functions are the constructor and - - + // few functions are the constructor and the + // helper functions that can be used to + // determine which part of the domain a cell + // is in. Given the discussion of these + // topics in the introduction, their + // implementation is rather obvious. In the + // constructor, note that we have to + // construct the hp::FECollection object from + // the base elements for Stokes and + // elasticity; using the + // hp::FECollection::push_back function + // assigns them spots zero and one in this + // collection, an order that we have to + // remember and use consistently in the rest + // of the program. template FluidStructureProblem:: FluidStructureProblem (const unsigned int stokes_degree, @@ -298,13 +312,32 @@ cell_is_in_solid_domain (const typename hp::DoFHandler::cell_iterator &cell } - + // @sect4{Meshes and assigning subdomains} + + // The next pair of functions deals with + // generating a mesh and making sure all + // flags that denote subdomains are + // correct. make_grid, as + // discussed in the introduction, generates + // an $8\times 8$ mesh (or an $8\times + // 8\times 8$ mesh in 3d) to make sure that + // each coarse mesh cell is completely within + // one of the subdomains. After generating + // this mesh, we loop over its boundary and + // set the boundary indicator to one at the + // top boundary, the only place where we set + // nonzero Dirichlet boundary + // conditions. After this, we loop again over + // all cells to set the material indicator + // — used to denote which part of the + // domain we are in, to either the fluid or + // solid indicator. template void FluidStructureProblem::make_grid () { -// not quite what we want... GridGenerator::subdivided_hyper_cube (triangulation, 8, -1, 1); + for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); cell != triangulation.end(); ++cell) @@ -331,7 +364,24 @@ FluidStructureProblem::make_grid () } - + // The second part of this pair of functions + // determines which finite element to use on + // each cell. Above we have set the material + // indicator for each coarse mesh cell, and + // as mentioned in the introduction, this + // information is inherited from mother to + // child cell upon mesh refinement. + // + // In other words, whenever we have refined + // (or created) the mesh, we can rely on the + // material indicators to be a correct + // description of which part of the domain a + // cell is in. We then use this to set the + // active FE index of the cell to the + // corresponding element of the + // hp::FECollection member variable of this + // class: zero for fluid cells, one for solid + // cells. template void FluidStructureProblem::set_active_fe_indices () @@ -339,18 +389,34 @@ FluidStructureProblem::set_active_fe_indices () for (typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell) - if (cell_is_in_fluid_domain(cell)) - cell->set_active_fe_index (0); - else if (cell_is_in_solid_domain(cell)) - cell->set_active_fe_index (1); - else - Assert (false, ExcNotImplemented()); + { + if (cell_is_in_fluid_domain(cell)) + cell->set_active_fe_index (0); + else if (cell_is_in_solid_domain(cell)) + cell->set_active_fe_index (1); + else + Assert (false, ExcNotImplemented()); + } } - + // @sect4{FluidStructureProblem::setup_dofs} + + // The next step is to setup the data + // structures for the linear system. To this + // end, we first have to set the active FE + // indices with the function immediately + // above, then distribute degrees of freedom, + // and then determine constraints on the + // linear system. The latter includes hanging + // node constraints as usual, but also the + // inhomogenous boundary values at the top + // fluid boundary, and zero boundary values + // along the perimeter of the solid + // subdomain. template -void FluidStructureProblem::setup_dofs () +void +FluidStructureProblem::setup_dofs () { set_active_fe_indices (); dof_handler.distribute_dofs (fe_collection); @@ -368,6 +434,7 @@ void FluidStructureProblem::setup_dofs () StokesBoundaryValues(), constraints, velocity_mask); + std::vector elasticity_mask (dim+1+dim, false); for (unsigned int d=dim+1; d::setup_dofs () elasticity_mask); } - // make sure velocity is zero at - // the interface + // There are more constraints we have to + // handle, though: we have to make sure + // that the velocity is zero at the + // interface between fluid and solid. The + // following piece of code was already + // presented in the introduction: { std::vector local_face_dof_indices (stokes_fe.dofs_per_face); for (typename hp::DoFHandler::active_cell_iterator @@ -397,14 +468,9 @@ void FluidStructureProblem::setup_dofs () face_is_on_interface = true; else if (cell->neighbor(f)->has_children() == true) { - // neighbor does - // have - // children. see if - // any of the cells - // on the other - // side are elastic for (unsigned int sf=0; sfface(f)->n_children(); ++sf) - if (cell_is_in_solid_domain (cell->neighbor_child_on_subface(f, sf))) + if (cell_is_in_solid_domain (cell->neighbor_child_on_subface + (f, sf))) { face_is_on_interface = true; break; @@ -420,7 +486,11 @@ void FluidStructureProblem::setup_dofs () } } - + // At the end of all this, we can declare + // to the constraints object that we now + // have all constraints ready to go and + // that the object can rebuild its internal + // data structures for better efficiency: constraints.close (); std::cout << " Number of active cells: " @@ -430,6 +500,10 @@ void FluidStructureProblem::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. { CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(), dof_handler.n_dofs()); @@ -446,6 +520,20 @@ void FluidStructureProblem::setup_dofs () + // @sect4{FluidStructureProblem::assemble_system} + + // Following is the central function of this + // program: the one that assembles the linear + // system. It has a long section of setting + // up auxiliary functions at the beginning: + // from creating the quadrature formulas and + // setting up the FEValues, FEFaceValues and + // FESubfaceValues objects necessary to + // integrate the cell terms as well as the + // interface terms for the case where cells + // along the interface come together at same + // size or with differing levels of + // refinement... template void FluidStructureProblem::assemble_system () { @@ -465,42 +553,50 @@ void FluidStructureProblem::assemble_system () update_JxW_values | update_gradients); - const QGauss face_quadrature(std::max (stokes_degree+2, - elasticity_degree+2)); + const QGauss common_face_quadrature(std::max (stokes_degree+2, + elasticity_degree+2)); FEFaceValues stokes_fe_face_values (stokes_fe, - face_quadrature, + common_face_quadrature, update_JxW_values | update_normal_vectors | update_gradients); FEFaceValues elasticity_fe_face_values (elasticity_fe, - face_quadrature, + common_face_quadrature, update_values); FESubfaceValues stokes_fe_subface_values (stokes_fe, - face_quadrature, + common_face_quadrature, update_JxW_values | update_normal_vectors | update_gradients); FESubfaceValues elasticity_fe_subface_values (elasticity_fe, - face_quadrature, + common_face_quadrature, update_values); - const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell; - const unsigned int elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell; + // ...to objects that are needed to + // describe the local contributions to the + // global linear system... + const unsigned int stokes_dofs_per_cell = stokes_fe.dofs_per_cell; + const unsigned int elasticity_dofs_per_cell = elasticity_fe.dofs_per_cell; - FullMatrix local_matrix; - FullMatrix local_interface_matrix (elasticity_dofs_per_cell, - stokes_dofs_per_cell); - Vector local_rhs; + FullMatrix local_matrix; + FullMatrix local_interface_matrix (elasticity_dofs_per_cell, + stokes_dofs_per_cell); + Vector local_rhs; 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); - const FEValuesExtractors::Vector displacements (dim+1); + // ...to variables that allow us to extract + // certain components of the shape + // functions and cache their values rather + // than having to recompute them at every + // quadrature point: + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + const FEValuesExtractors::Vector displacements (dim+1); std::vector > stokes_phi_grads_u (stokes_dofs_per_cell); std::vector stokes_div_phi_u (stokes_dofs_per_cell); @@ -508,8 +604,14 @@ void FluidStructureProblem::assemble_system () std::vector > elasticity_phi_grad (elasticity_dofs_per_cell); std::vector elasticity_phi_div (elasticity_dofs_per_cell); - std::vector > elasticity_phi (elasticity_dofs_per_cell); - + std::vector > elasticity_phi (elasticity_dofs_per_cell); + + // Then comes the main loop over all cells + // and, as in step-27, the initialization + // of the hp::FEValues object for the + // current cell and the extraction of a + // FEValues object that is appropriate for + // the current cell: typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -523,6 +625,29 @@ void FluidStructureProblem::assemble_system () cell->get_fe().dofs_per_cell); local_rhs.reinit (cell->get_fe().dofs_per_cell); + // With all of this done, we continue + // to assemble the cell terms for cells + // that are part of the Stokes and + // elastic regions. While we could in + // principle do this in one formula, in + // effect implementing the one bilinear + // form stated in the introduction, we + // realize that our finite element + // spaces are chosen in such a way that + // on each cell, one set of variables + // (either velocities and pressure, or + // displacements) are always zero, and + // consequently a more efficient way of + // computing local integrals is to do + // only what's necessary based on an + // if clause that tests + // which part of the domain we are in. + // + // The actual computation of the local + // matrix is the same as in step-22 as + // well as that given in the @ref + // vector_valued documentation module + // for the elasticity equations: if (cell_is_in_fluid_domain (cell)) { const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; @@ -579,32 +704,100 @@ void FluidStructureProblem::assemble_system () } } + // Once we have the contributions from + // cell integrals, we copy them into + // the global matrix (taking care of + // constraints right away, through the + // ConstraintMatrix::distribute_local_to_global + // function). Note that we have not + // written anything into the + // local_rhs variable, + // though we still need to pass it + // along since the elimination of + // nonzero boundary values requires the + // modification of local and + // consequently also global right hand + // side values: local_dof_indices.resize (cell->get_fe().dofs_per_cell); cell->get_dof_indices (local_dof_indices); - - // local_rhs==0, but need to do - // this here because of - // boundary values constraints.distribute_local_to_global (local_matrix, local_rhs, local_dof_indices, system_matrix, system_rhs); - // see about face terms + // The more interesting part of this + // function is where we see about face + // terms along the interface between + // the two subdomains. To this end, we + // first have to make sure that we only + // assemble them once even though a + // loop over all faces of all cells + // would encounter each part of the + // interface twice. We arbitrarily make + // the decision that we will only + // evaluate interface terms if the + // current cell is part of the solid + // subdomain and if, consequently, a + // face is not at the boundary and the + // potential neighbor behind it is part + // of the fluid domain. Let's start + // with these conditions: if (cell_is_in_solid_domain (cell)) - // we are on a solid cell for (unsigned int f=0; f::faces_per_cell; ++f) if (cell->at_boundary(f) == false) { + // At this point we know that + // the current cell is a + // candidate for integration + // and that a neighbor behind + // face f + // exists. There are now three + // possibilities: + // + // - The neighbor is at the + // same refinement level and + // has no children. + // - The neighbor has children. + // - The neighbor is coarser. + // + // In all three cases, we are + // only interested in it if it + // is part of the fluid + // subdomain. So let us start + // with the first and simplest + // case: if the neighbor is at + // the same level, has no + // children, and is a fluid + // cell, then the two cells + // share a boundary that is + // part of the interface along + // which we want to integrate + // interface terms. All we have + // to do is initialize two + // FEFaceValues object with the + // current face and the face of + // the neighboring cell (note + // how we find out which face + // of the neighboring cell + // borders on the current cell) + // and pass things off to the + // function that evaluates the + // interface terms (the third + // through fifth arguments to + // this function provide it + // with scratch arrays). The + // result is then again copied + // into the global matrix, + // using a function that knows + // that the DoF indices of rows + // and columns of the local + // matrix result from different + // cells: if ((cell->neighbor(f)->level() == cell->level()) && (cell->neighbor(f)->has_children() == false) && cell_is_in_fluid_domain (cell->neighbor(f))) { - // same size - // neighbors; - // neighbor is - // fluid cell elasticity_fe_face_values.reinit (cell, f); stokes_fe_face_values.reinit (cell->neighbor(f), cell->neighbor_of_neighbor(f)); @@ -619,16 +812,32 @@ void FluidStructureProblem::assemble_system () neighbor_dof_indices, system_matrix); } + + // The second case is if the + // neighbor has further + // children. In that case, we + // have to loop over all the + // children of the neighbor to + // see if they are part of the + // fluid subdomain. If they + // are, then we integrate over + // the common interface, which + // is a face for the neighbor + // and a subface of the current + // cell, requiring us to use an + // FEFaceValues for the + // neighbor and an + // FESubfaceValues for the + // current cell: else if ((cell->neighbor(f)->level() == cell->level()) && (cell->neighbor(f)->has_children() == true)) { - // neighbor has children. loop over - // the cells adjacent to the commone - // interface and see which subdomain - // they belong to - for (unsigned int subface=0; subfaceface(f)->n_children(); ++subface) - if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface (f, subface))) + for (unsigned int subface=0; + subfaceface(f)->n_children(); + ++subface) + if (cell_is_in_fluid_domain (cell->neighbor_child_on_subface + (f, subface))) { elasticity_fe_subface_values.reinit (cell, f, @@ -636,29 +845,42 @@ void FluidStructureProblem::assemble_system () stokes_fe_face_values.reinit (cell->neighbor_child_on_subface (f, subface), cell->neighbor_of_neighbor(f)); - assemble_interface_term (elasticity_fe_subface_values, stokes_fe_face_values, - elasticity_phi, stokes_phi_grads_u, stokes_phi_p, + assemble_interface_term (elasticity_fe_subface_values, + stokes_fe_face_values, + elasticity_phi, + stokes_phi_grads_u, stokes_phi_p, local_interface_matrix); - cell->neighbor_child_on_subface (f, subface)->get_dof_indices (neighbor_dof_indices); + cell->neighbor_child_on_subface (f, subface) + ->get_dof_indices (neighbor_dof_indices); constraints.distribute_local_to_global(local_interface_matrix, local_dof_indices, neighbor_dof_indices, system_matrix); } } + + // The last option is that the + // neighbor is coarser. In that + // case we have to use an + // FESubfaceValues object for + // the neighbor and a + // FEFaceValues for the current + // cell; the rest is the same + // as before: else if (cell->neighbor_is_coarser(f) && cell_is_in_fluid_domain(cell->neighbor(f))) { - // neighbor is coarser elasticity_fe_face_values.reinit (cell, f); stokes_fe_subface_values.reinit (cell->neighbor(f), cell->neighbor_of_coarser_neighbor(f).first, cell->neighbor_of_coarser_neighbor(f).second); - assemble_interface_term (elasticity_fe_face_values, stokes_fe_subface_values, - elasticity_phi, stokes_phi_grads_u, stokes_phi_p, + assemble_interface_term (elasticity_fe_face_values, + stokes_fe_subface_values, + elasticity_phi, + stokes_phi_grads_u, stokes_phi_p, local_interface_matrix); cell->neighbor(f)->get_dof_indices (neighbor_dof_indices); @@ -674,25 +896,47 @@ void FluidStructureProblem::assemble_system () + // In the function that assembles the global + // system, we passed computing interface + // terms to a separate function we discuss + // here. The key is that even though we can't + // predict the combination of FEFaceValues + // and FESubfaceValues objects, they are both + // derived from the FEFaceValuesBase class + // and consequently we don't have to care: + // the function is simply called with two + // such objects denoting the values of the + // shape functions on the quadrature points + // of the two sides of the face. We then do + // what we always do: we fill the scratch + // arrays with the values of shape functions + // and their derivatives, and then loop over + // all entries of the matrix to compute the + // local integrals. The details of the + // bilinear form we evaluate here are given + // in the introduction. template void -FluidStructureProblem::assemble_interface_term (const FEFaceValuesBase &elasticity_fe_face_values, - const FEFaceValuesBase &stokes_fe_face_values, - std::vector > &elasticity_phi, - std::vector > &stokes_phi_grads_u, - std::vector &stokes_phi_p, - FullMatrix &local_interface_matrix) const +FluidStructureProblem:: +assemble_interface_term (const FEFaceValuesBase &elasticity_fe_face_values, + const FEFaceValuesBase &stokes_fe_face_values, + std::vector > &elasticity_phi, + std::vector > &stokes_phi_grads_u, + std::vector &stokes_phi_p, + FullMatrix &local_interface_matrix) const { Assert (stokes_fe_face_values.n_quadrature_points == elasticity_fe_face_values.n_quadrature_points, ExcInternalError()); - + const unsigned int n_face_quadrature_points + = elasticity_fe_face_values.n_quadrature_points; + const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); const FEValuesExtractors::Vector displacements (dim+1); local_interface_matrix = 0; - for (unsigned int q=0; q normal_vector = stokes_fe_face_values.normal_vector(q); @@ -715,6 +959,16 @@ FluidStructureProblem::assemble_interface_term (const FEFaceValuesBase } + // @sect4{FluidStructureProblem::solve} + + // As discussed in the introduction, we use a + // rather trivial solver here: we just pass + // the linear system off to the + // SparseDirectUMFPACK direct solver (see, + // for example, step-29). The only thing we + // have to do after solving is ensure that + // hanging node and boundary value + // constraints are correct. template void FluidStructureProblem::solve () @@ -728,10 +982,21 @@ FluidStructureProblem::solve () + // @sect4{FluidStructureProblem::output_results} + // Generating graphical output is rather + // trivial here: all we have to do is + // identify which components of the solution + // vector belong to scalars and/or vectors + // (see, for example, step-22 for a previous + // example), and then pass it all on to the + // DataOut class (with the second template + // argument equal to hp::DoFHandler instead + // of the usual default DoFHandler): template void -FluidStructureProblem::output_results (const unsigned int refinement_cycle) const +FluidStructureProblem:: +output_results (const unsigned int refinement_cycle) const { std::vector solution_names (dim, "velocity"); solution_names.push_back ("pressure"); @@ -765,21 +1030,37 @@ FluidStructureProblem::output_results (const unsigned int refinement_cycle) } - + // @sect4{FluidStructureProblem::refine_mesh} + + // The next step is to refine the mesh. As + // was discussed in the introduction, this is + // a bit tricky primarily because the fluid + // and the solid subdomains use variables + // that have different physical dimensions + // and for which the absolute magnitude of + // error estimates is consequently not + // directly comparable. We will therefore + // have to scale them. At the top of the + // function, we therefore first compute error + // estimates for the different variables + // separately (using the velocities but not + // the pressure for the fluid domain, and the + // displacements in the solid domain): template void FluidStructureProblem::refine_mesh () { - Vector stokes_estimated_error_per_cell (triangulation.n_active_cells()); - Vector elasticity_estimated_error_per_cell (triangulation.n_active_cells()); - Vector estimated_error_per_cell (triangulation.n_active_cells()); + Vector + stokes_estimated_error_per_cell (triangulation.n_active_cells()); + Vector + elasticity_estimated_error_per_cell (triangulation.n_active_cells()); - const QGauss stokes_quadrature(stokes_degree+2); - const QGauss elasticity_quadrature(elasticity_degree+2); + const QGauss stokes_face_quadrature(stokes_degree+2); + const QGauss elasticity_face_quadrature(elasticity_degree+2); hp::QCollection face_q_collection; - face_q_collection.push_back (stokes_quadrature); - face_q_collection.push_back (elasticity_quadrature); + face_q_collection.push_back (stokes_face_quadrature); + face_q_collection.push_back (elasticity_face_quadrature); std::vector stokes_component_mask (dim+1+dim, false); for (unsigned int d=0; d::refine_mesh () elasticity_estimated_error_per_cell, elasticity_component_mask); - 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(); + // We then normalize error estimates by + // dividing by their norm and scale the + // fluid error indicators by a factor of 4 + // as discussed in the introduction. The + // results are then added together into a + // vector that contains error indicators + // for all cells: + stokes_estimated_error_per_cell + *= 4 . / stokes_estimated_error_per_cell.l2_norm(); + elasticity_estimated_error_per_cell + *= 1. / 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; + // The second to last part of the function, + // before actually refining the mesh, + // involves a heuristic that we have + // already mentioned in the introduction: + // because the solution is discontinuous, + // the KellyErrorEstimator class gets all + // confused about cells that sit at the + // boundary between subdomains: it believes + // that the error is large there because + // the jump in the gradient is large, even + // though this is entirely expected and a + // feature that is in fact present in the + // exact solution as well and therefore not + // indicative of any numerical error. + // + // Consequently, we set the error + // indicators to zero for all cells at the + // interface; the conditions determining + // which cells this affects are slightly + // awkward because we have to account for + // the possibility of adaptively refined + // meshes, meaning that the neighboring + // cell can be coarser than the current + // one, or could in fact be refined some + // more. The structure of these nested + // conditions is much the same as we + // encountered when assembling interface + // terms in assemble_system. { unsigned int cell_index = 0; for (typename hp::DoFHandler::active_cell_iterator @@ -868,6 +1190,14 @@ FluidStructureProblem::refine_mesh () + // @sect4{FluidStructureProblem::run} + + // This is, as usual, the function that + // controls the overall flow of operation. If + // you've read through tutorial programs + // step-1 through step-6, for example, then + // you are already quite familiar with the + // following structure: template void FluidStructureProblem::run () { @@ -898,6 +1228,11 @@ void FluidStructureProblem::run () + // @sect4{The main() function} + + // This, final, function contains pretty much + // exactly what most of the other tutorial + // programs have: int main () { try