From 83b0395f383e6b606e3a1b94d00587e6bd34d0fa Mon Sep 17 00:00:00 2001 From: bangerth Date: Sun, 21 Oct 2007 06:16:36 +0000 Subject: [PATCH] Use a better refinement algorithm. Re-structure the code that refines the mesh up front. Also, add the term to the bilinear form for the T-equation. This term is zero on the continuous level, but isn't on the discrete level. git-svn-id: https://svn.dealii.org/trunk@15359 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-22/step-22.cc | 74 ++++++++++++++++------------- 1 file changed, 41 insertions(+), 33 deletions(-) diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index c55e7377d9..492ca1e63e 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -685,6 +685,11 @@ void BoussinesqFlowProblem::assemble_rhs_T () std::vector > present_solution_values(n_q_points, Vector(dim+2)); std::vector > present_solution_values_face(n_face_q_points, Vector(dim+2)); + std::vector > > + present_solution_grads(n_q_points, + std::vector >(dim+2)); + + std::vector neighbor_temperature (n_face_q_points); std::vector local_dof_indices (dofs_per_cell); @@ -701,6 +706,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () fe_values.get_function_values (old_solution, old_solution_values); fe_values.get_function_values (solution, present_solution_values); + fe_values.get_function_gradients (solution, present_solution_grads); for (unsigned int q=0; q::assemble_rhs_T () for (unsigned int d=0; d grad_phi_i_T = extract_grad_T(fe_values, i, q); local_rhs(i) += (time_step * old_T * - present_u * - grad_phi_i_T + (present_u * + grad_phi_i_T + + + present_div_u * + phi_i_T) + old_T * phi_i_T) * fe_values.JxW(q); } + //TODO: unify the code that actually does the assembly down below for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) @@ -764,7 +778,6 @@ void BoussinesqFlowProblem::assemble_rhs_T () for (unsigned int q=0; q::assemble_rhs_T () for (unsigned int q=0; q present_u_face; @@ -1023,7 +1036,7 @@ void BoussinesqFlowProblem::output_results () const data_out.add_data_vector (solution, solution_names, DataOut::type_dof_data, data_component_interpretation); - + data_out.build_patches (); std::ostringstream filename; @@ -1046,11 +1059,18 @@ BoussinesqFlowProblem::refine_mesh () old_solution, estimated_error_per_cell, dim+1); - + + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index) + estimated_error_per_cell(cell_index) *= cell->diameter(); + GridRefinement::refine_and_coarsen_fixed_fraction (triangulation, estimated_error_per_cell, 0.3, 0.03, - triangulation.n_active_cells()); + static_cast + (triangulation.n_active_cells()*1.1)); SolutionTransfer soltrans(dof_handler); @@ -1125,7 +1145,7 @@ void BoussinesqFlowProblem::run () static HalfHyperShellBoundary boundary; triangulation.set_boundary (0, boundary); - triangulation.refine_global (3); + triangulation.refine_global (4); break; } @@ -1148,38 +1168,24 @@ void BoussinesqFlowProblem::run () } - for (unsigned int pre_refinement=0; pre_refinement<3-dim; ++pre_refinement) + setup_dofs(false); + + VectorTools::project (dof_handler, + hanging_node_constraints, + QGauss(degree+2), + InitialValues(), + old_solution); + + for (unsigned int pre_refinement=0; pre_refinement<4-dim; ++pre_refinement) { - setup_dofs(false); - + refine_mesh (); + VectorTools::project (dof_handler, hanging_node_constraints, QGauss(degree+2), InitialValues(), old_solution); - - Vector estimated_error_per_cell (triangulation.n_active_cells()); - -//TODO do this better - DerivativeApproximation::approximate_gradient (dof_handler, - old_solution, - estimated_error_per_cell, - dim+1); - - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.03); - - triangulation.execute_coarsening_and_refinement (); } - - setup_dofs(true); - VectorTools::project (dof_handler, - hanging_node_constraints, - QGauss(degree+2), - InitialValues(), - old_solution); - timestep_number = 0; double time = 0; @@ -1206,6 +1212,8 @@ void BoussinesqFlowProblem::run () std::cout << std::endl; + break; + if (timestep_number % 10 == 0) refine_mesh (); } -- 2.39.5