From: bangerth Date: Tue, 22 Jan 2008 23:34:03 +0000 (+0000) Subject: Check in the benchmark that mimics three candles. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2b81c00e5a8c03e7079537602d266f44bde36c41;p=dealii-svn.git Check in the benchmark that mimics three candles. git-svn-id: https://svn.dealii.org/trunk@15642 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index cfd44f6b53..9527fef181 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -163,7 +163,6 @@ class InitialValues : public Function virtual void vector_value (const Point &p, Vector &value) const; - }; @@ -171,9 +170,51 @@ template double InitialValues::value (const Point &p, const unsigned int component) const +{ + return 0; +} + + +template +void +InitialValues::vector_value (const Point &p, + Vector &values) const +{ + for (unsigned int c=0; cn_components; ++c) + values(c) = InitialValues::value (p, c); +} + + + +template +class RightHandSide : public Function +{ + public: + RightHandSide () : Function(dim+2) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int component) const { if (component == dim+1) - return (p.distance (Point()) < .6 ? 1 : 0); + return ((p.distance (Point(.3,.1)) < 1./32) + || + (p.distance (Point(.45,.1)) < 1./32) + || + (p.distance (Point(.75,.1)) < 1./32) + ? + 1 + : + 0); else return 0; } @@ -181,11 +222,11 @@ InitialValues::value (const Point &p, template void -InitialValues::vector_value (const Point &p, +RightHandSide::vector_value (const Point &p, Vector &values) const { for (unsigned int c=0; cn_components; ++c) - values(c) = InitialValues::value (p, c); + values(c) = RightHandSide::value (p, c); } @@ -763,10 +804,10 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints); - std::set no_normal_flux_boundaries; - no_normal_flux_boundaries.insert (0); - compute_no_normal_flux_constraints (dof_handler, 0, no_normal_flux_boundaries, - hanging_node_constraints); +// std::set no_normal_flux_boundaries; +// no_normal_flux_boundaries.insert (0); +// compute_no_normal_flux_constraints (dof_handler, 0, no_normal_flux_boundaries, +// hanging_node_constraints); hanging_node_constraints.close (); std::vector dofs_per_component (dim+2); @@ -934,8 +975,7 @@ void BoussinesqFlowProblem::assemble_system () } } - const Point gravity = fe_values.quadrature_point(q) / - fe_values.quadrature_point(q).norm(); + const Point gravity (0,1); local_rhs(i) += (Raleigh_number * gravity * phi_i_u * old_temperature)* @@ -1010,13 +1050,13 @@ void BoussinesqFlowProblem::assemble_system () break; } -// std::vector component_mask (dim+2, true); -// component_mask[dim] = component_mask[dim+1] = false; -// VectorTools::interpolate_boundary_values (dof_handler, -// 0, -// ZeroFunction(dim+2), -// boundary_values, -// component_mask); + std::vector component_mask (dim+2, true); + component_mask[dim] = component_mask[dim+1] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(dim+2), + boundary_values, + component_mask); MatrixTools::apply_boundary_values (boundary_values, system_matrix, @@ -1052,8 +1092,8 @@ void BoussinesqFlowProblem::assemble_system () template void BoussinesqFlowProblem::assemble_rhs_T () { - QGauss quadrature_formula(degree+1); - QGauss face_quadrature_formula(degree+1); + QGauss quadrature_formula(degree+2); + QGauss face_quadrature_formula(degree+2); FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); @@ -1118,7 +1158,9 @@ void BoussinesqFlowProblem::assemble_rhs_T () const double phi_i_T = extract_T(fe_values, i, q); const Tensor<1,dim> grad_phi_i_T = extract_grad_T(fe_values, i, q); - + + const Point p = fe_values.quadrature_point(q); + local_rhs(i) += (time_step * old_T * (present_u * @@ -1127,7 +1169,11 @@ void BoussinesqFlowProblem::assemble_rhs_T () present_div_u * phi_i_T) + - old_T * phi_i_T) + old_T * phi_i_T + + + time_step * + RightHandSide().value (p, dim+1) + * phi_i_T) * fe_values.JxW(q); } @@ -1381,8 +1427,9 @@ void BoussinesqFlowProblem::solve () hanging_node_constraints.distribute (solution); } + // for DGQ1 needs to be /15 time_step = GridTools::minimal_cell_diameter(triangulation) / - get_maximal_velocity() / 2; + std::max (get_maximal_velocity(), .05) / 2; assemble_rhs_T (); { @@ -1409,7 +1456,11 @@ void BoussinesqFlowProblem::solve () std::cout << " " << solver_control.last_step() << " CG iterations for temperature." - << std::endl; + << std::endl; + std::cout << " Max temperature: " + << *std::max_element (solution.block(2).begin(), + solution.block(2).end()) + << std::endl; } } @@ -1543,13 +1594,9 @@ void BoussinesqFlowProblem::run () { case 2: { - GridGenerator::hyper_shell (triangulation, - Point(), 0.5, 1.0); - - static HyperShellBoundary boundary; - triangulation.set_boundary (0, boundary); + GridGenerator::hyper_cube (triangulation); - triangulation.refine_global (4); + triangulation.refine_global (5); break; }