From be723b7edeebe4430c5df701ec467c058e573334 Mon Sep 17 00:00:00 2001 From: Zhuoran Wang Date: Mon, 16 Sep 2019 10:41:18 -0600 Subject: [PATCH] Add Darcy velocity. --- examples/step-61/step-61.cc | 60 +++++++++++++++++++++++++++++++++++-- 1 file changed, 57 insertions(+), 3 deletions(-) diff --git a/examples/step-61/step-61.cc b/examples/step-61/step-61.cc index 2ee5436a25..0ea6c6dc19 100644 --- a/examples/step-61/step-61.cc +++ b/examples/step-61/step-61.cc @@ -43,6 +43,7 @@ #include #include #include +#include #include #include #include @@ -64,7 +65,7 @@ // program into its own namespace: namespace Step61 { - using namespace dealii; + using namespace dealii; // @sect3{The WGDarcyEquation class template} @@ -96,6 +97,9 @@ namespace Step61 FESystem fe; DoFHandler dof_handler; + + FE_DGRaviartThomas fe_dgrt; + DoFHandler dof_handler_dgrt; AffineConstraints constraints; @@ -104,6 +108,8 @@ namespace Step61 Vector solution; Vector system_rhs; + + Vector darcy_velocity; }; @@ -254,7 +260,10 @@ namespace Step61 template WGDarcyEquation::WGDarcyEquation(const unsigned int degree) : fe(FE_DGQ(degree), 1, FE_FaceQ(degree), 1) - , dof_handler(triangulation) + , dof_handler(triangulation), + + fe_dgrt(0), + dof_handler_dgrt (triangulation) {} @@ -293,12 +302,15 @@ namespace Step61 void WGDarcyEquation::setup_system() { dof_handler.distribute_dofs(fe); + dof_handler_dgrt.distribute_dofs (fe_dgrt); std::cout << " Number of pressure degrees of freedom: " << dof_handler.n_dofs() << std::endl; solution.reinit(dof_handler.n_dofs()); system_rhs.reinit(dof_handler.n_dofs()); + + darcy_velocity.reinit (dof_handler_dgrt.n_dofs()); { constraints.clear(); @@ -683,9 +695,19 @@ namespace Step61 update_values | update_normal_vectors | update_quadrature_points | update_JxW_values); + + FEValues fe_values_dgrt (fe_dgrt, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEFaceValues fe_face_values_dgrt (fe_dgrt, face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int dofs_per_cell_rt = fe_rt.dofs_per_cell; + const unsigned int dofs_per_cell_dgrt = fe_dgrt.dofs_per_cell; const unsigned int n_q_points = fe_values.get_quadrature().size(); const unsigned int n_q_points_rt = fe_values_rt.get_quadrature().size(); @@ -696,6 +718,7 @@ namespace Step61 std::vector local_dof_indices(dofs_per_cell); + std::vector local_dof_indices_dgrt(dofs_per_cell_dgrt); FullMatrix cell_matrix_M(dofs_per_cell_rt, dofs_per_cell_rt); FullMatrix cell_matrix_G(dofs_per_cell_rt, dofs_per_cell); @@ -717,6 +740,7 @@ namespace Step61 const FEValuesExtractors::Scalar pressure(dim); const FEValuesExtractors::Scalar interior(0); const FEValuesExtractors::Scalar face(1); + const FEValuesExtractors::Vector velocities_dgrt(0); const ExactVelocity exact_velocity; @@ -736,9 +760,14 @@ namespace Step61 // components. Then, we multiply all these coefficients and call them beta. // The numerical velocity is the product of beta and the basis functions of // the Raviart-Thomas space. - for (const auto &cell : dof_handler.active_cell_iterators()) + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(), + cell_dgrt = dof_handler_dgrt.begin_active(); + for (; cell!=endc; ++cell,++cell_dgrt) { fe_values.reinit(cell); + fe_values_dgrt.reinit(cell_dgrt); const typename Triangulation::active_cell_iterator cell_rt = cell; fe_values_rt.reinit(cell_rt); @@ -832,6 +861,15 @@ namespace Step61 for (unsigned int i = 0; i < dofs_per_cell; ++i) cell_velocity(k) += -(cell_solution(i) * cell_matrix_C(i, j) * cell_matrix_D(k, j)); + + // We compute Darcy velocity. + // This is same as cell_velocity but used to graph Darcy velocity. + cell_dgrt->get_dof_indices(local_dof_indices_dgrt); + for (unsigned int k = 0; k solution_names (dim, "velocity"); + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + + DataOut data_out_dgrt; + data_out_dgrt.attach_dof_handler (dof_handler_dgrt); + data_out_dgrt.add_data_vector (darcy_velocity, solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out_dgrt.build_patches (fe_dgrt.degree); + std::ofstream dgrt_output ("Darcy_velocity.vtk"); + data_out_dgrt.write_vtk (dgrt_output); + } } -- 2.39.5