From: bangerth Date: Sat, 23 Apr 2011 03:18:46 +0000 (+0000) Subject: Some more cleanups, a few more comments. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=123cf4d3e97483d1a68ee43f07cb1724a430340d;p=dealii-svn.git Some more cleanups, a few more comments. git-svn-id: https://svn.dealii.org/trunk@23633 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 ebaf9f1937..5b7032f770 100644 --- a/deal.II/examples/step-46/step-46.cc +++ b/deal.II/examples/step-46/step-46.cc @@ -13,6 +13,12 @@ // @sect3{Include files} + // The include files for this program are the + // same as for many others before. The only + // new one is the one that declares + // FE_Nothing as discussed in the + // introduction. The ones in the hp directory + // have already been discussed in step-27. #include #include @@ -52,7 +58,31 @@ using namespace dealii; - + // @sect3{The FluidStructureProblem class template} + + // This is the main class. It is, if you + // want, a combination of step-8 and step-22 + // in that it has member variables that + // either address the global problem (the + // Triangulation and hp::DoFHandler objects, + // as well as the hp::FECollection and + // various linear algebra objects) or that + // pertain to either the elasticity or Stokes + // sub-problems. The general structure of the + // class, however, is like that of most of + // the other programs implementing stationary + // problems. + // + // There are a few helper functions + // (cell_is_in_fluid_domain, + // cell_is_in_solid_domain) of + // self-explanatory nature and a few + // functions (make_grid, + // setup_subdomains, + // assemble_interface_terms) that have + // been broken out of other functions and + // will be discussed as we get to their + // implementation. template class FluidStructureProblem { @@ -69,6 +99,7 @@ class FluidStructureProblem cell_is_in_solid_domain (const typename hp::DoFHandler::cell_iterator &cell); + void make_grid (); void setup_subdomains (); void setup_dofs (); void assemble_system (); @@ -96,21 +127,33 @@ class FluidStructureProblem SparsityPattern sparsity_pattern; SparseMatrix system_matrix; - Vector solution; - Vector system_rhs; + Vector solution; + Vector system_rhs; - const double viscosity; - const double lambda; - const double mu; + const double viscosity; + const double lambda; + const double mu; }; + // @sect3{Boundary values and right hand side} + // The following classes do as their names + // suggest. The boundary values for the + // velocity are $\mathbf u=(0, \sin(\pi + // x))^T$ in 2d and $\mathbf u=(0, 0, + // \sin(\pi x)\sin(\pi y))^T$ in 3d, + // respectively. The remaining boundary + // conditions for this problem are all + // homogenous and have been discussed in the + // introduction. The right hand side forcing + // term is zero for both the fluid and the + // solid. template -class BoundaryValues : public Function +class StokesBoundaryValues : public Function { public: - BoundaryValues () : Function(dim+1+dim) {} + StokesBoundaryValues () : Function(dim+1+dim) {} virtual double value (const Point &p, const unsigned int component = 0) const; @@ -122,7 +165,7 @@ class BoundaryValues : public Function template double -BoundaryValues::value (const Point &p, +StokesBoundaryValues::value (const Point &p, const unsigned int component) const { Assert (component < this->n_components, @@ -145,11 +188,11 @@ BoundaryValues::value (const Point &p, template void -BoundaryValues::vector_value (const Point &p, +StokesBoundaryValues::vector_value (const Point &p, Vector &values) const { for (unsigned int c=0; cn_components; ++c) - values(c) = BoundaryValues::value (p, c); + values(c) = StokesBoundaryValues::value (p, c); } @@ -189,8 +232,11 @@ RightHandSide::vector_value (const Point &p, + // @sect3{The FluidStructureProblem implementation} - + // Let's now get to the implementation of the + // primary class of this program. The first + // few functions are the constructor and template @@ -238,13 +284,52 @@ cell_is_in_solid_domain (const typename hp::DoFHandler::cell_iterator &cell +template +void +FluidStructureProblem::make_grid () +{ + GridGenerator::subdivided_hyper_cube (triangulation, 1, -1, 1); + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary() + && + (cell->face(f)->center()[dim-1] == 1)) + cell->face(f)->set_all_boundary_indicators(1); + + triangulation.refine_global (5-dim); +} + + +template +void +FluidStructureProblem::setup_subdomains () +{ + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (((std::fabs(cell->center()[0]) < 0.25) + && + (cell->center()[dim-1] > 0.5)) + || + ((std::fabs(cell->center()[0]) >= 0.25) + && + (cell->center()[dim-1] > -0.5))) + cell->set_active_fe_index (0); + else + cell->set_active_fe_index (1); +} + + + template void FluidStructureProblem::setup_dofs () { system_matrix.clear (); dof_handler.distribute_dofs (fe_collection); - + { constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, @@ -255,7 +340,7 @@ void FluidStructureProblem::setup_dofs () velocity_mask[d] = true; VectorTools::interpolate_boundary_values (dof_handler, 1, - BoundaryValues(), + StokesBoundaryValues(), constraints, velocity_mask); std::vector elasticity_mask (dim+1+dim, false); @@ -336,27 +421,6 @@ void FluidStructureProblem::setup_dofs () -template -void -FluidStructureProblem::setup_subdomains () -{ - for (typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (((std::fabs(cell->center()[0]) < 0.25) - && - (cell->center()[dim-1] > 0.5)) - || - ((std::fabs(cell->center()[0]) >= 0.25) - && - (cell->center()[dim-1] > -0.5))) - cell->set_active_fe_index (0); - else - cell->set_active_fe_index (1); -} - - - template void FluidStructureProblem::assemble_system () { @@ -782,16 +846,7 @@ FluidStructureProblem::refine_mesh () template void FluidStructureProblem::run () { - GridGenerator::hyper_cube (triangulation, -1, 1); - for (typename Triangulation::active_cell_iterator - cell = triangulation.begin_active(); - cell != triangulation.end(); ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->at_boundary() - && - (cell->face(f)->center()[dim-1] == 1)) - cell->face(f)->set_all_boundary_indicators(1); - triangulation.refine_global (5-dim); + make_grid (); for (unsigned int refinement_cycle = 0; refinement_cycle<10-2*dim; ++refinement_cycle) @@ -804,16 +859,16 @@ void FluidStructureProblem::run () setup_subdomains (); setup_dofs (); - std::cout << " Assembling..." << std::endl << std::flush; + std::cout << " Assembling..." << std::endl; assemble_system (); - std::cout << " Solving..." << std::flush; + std::cout << " Solving..." << std::endl; solve (); - std::cout << " Writing output..." << std::flush; + std::cout << " Writing output..." << std::endl; output_results (refinement_cycle); - std::cout << std::endl << std::endl; + std::cout << std::endl; } }