From db0105bab6205bdc084931793f257162223654c1 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 5 Oct 2016 10:44:28 +0200 Subject: [PATCH] minor cleanup of fe_enriched_step-36 unit test --- tests/fe/fe_enriched_step-36.cc | 154 ++++++-------------------------- 1 file changed, 25 insertions(+), 129 deletions(-) diff --git a/tests/fe/fe_enriched_step-36.cc b/tests/fe/fe_enriched_step-36.cc index 22702ae43e..9a20201fe7 100644 --- a/tests/fe/fe_enriched_step-36.cc +++ b/tests/fe/fe_enriched_step-36.cc @@ -62,19 +62,19 @@ using namespace dealii; * Coulomb potential */ template -class PotentialFunction : public dealii::Function +class PotentialFunction : public Function { public: PotentialFunction() - : dealii::Function(1) + : Function(1) {} - virtual double value(const dealii::Point &point, + virtual double value(const Point &point, const unsigned int component = 0 ) const; }; template -double PotentialFunction::value(const dealii::Point &p, +double PotentialFunction::value(const Point &p, const unsigned int ) const { Assert (p.square() > 0., @@ -159,7 +159,6 @@ namespace Step36 std::pair setup_system (); void assemble_system (); std::pair solve (); - void constrain_pou_dofs(); void estimate_error (); void refine_grid (); void output_results (const unsigned int cycle) const; @@ -192,11 +191,9 @@ namespace Step36 EnrichmentFunction enrichment; const FEValuesExtractors::Scalar fe_extractor; - const unsigned int fe_group; const unsigned int fe_fe_index; const unsigned int fe_material_id; const FEValuesExtractors::Scalar pou_extractor; - const unsigned int pou_group; const unsigned int pou_fe_index; const unsigned int pou_material_id; @@ -209,8 +206,8 @@ namespace Step36 : dof_handler (triangulation), mpi_communicator(MPI_COMM_WORLD), - n_mpi_processes(dealii::Utilities::MPI::n_mpi_processes(mpi_communicator)), - this_mpi_process(dealii::Utilities::MPI::this_mpi_process(mpi_communicator)), + n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator)), + this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator)), pcout (std::cout, (this_mpi_process == 0)), number_of_eigenvalues(1), @@ -218,15 +215,13 @@ namespace Step36 /*Z*/1.0, /*radius*/2.5), // radius is set such that 8 cells are marked as enriched fe_extractor(/*dofs start at...*/0), - fe_group(/*in FE*/0), fe_fe_index(0), fe_material_id(0), pou_extractor(/*dofs start at (scalar fields!)*/1), - pou_group(/*FE*/1), pou_fe_index(1), pou_material_id(1) { - dealii::GridGenerator::hyper_cube (triangulation, -10, 10); + GridGenerator::hyper_cube (triangulation, -10, 10); triangulation.refine_global (2); // 64 cells for (auto cell= triangulation.begin_active(); cell != triangulation.end(); ++cell) @@ -363,49 +358,6 @@ namespace Step36 return cell->material_id() == pou_material_id; } - template - void - EigenvalueProblem::constrain_pou_dofs() - { - std::vector local_face_dof_indices(fe_collection[pou_fe_index].dofs_per_face); - for (typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (cell_is_pou(cell)) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (!cell->at_boundary(f)) - { - bool face_is_on_interface = false; - if ((cell->neighbor(f)->has_children() == false /* => it is active */) - && - (! cell_is_pou(cell->neighbor(f)))) - face_is_on_interface = true; - else if (cell->neighbor(f)->has_children() == true) - { - // The neighbor does - // have - // children. See if - // any of the cells - // on the other - // side are vanilla FEM - for (unsigned int sf=0; sfface(f)->n_children(); ++sf) - if (!cell_is_pou (cell->neighbor_child_on_subface(f, sf))) - { - face_is_on_interface = true; - break; - } - } - // add constraints - if (face_is_on_interface) - { - cell->face(f)->get_dof_indices (local_face_dof_indices, pou_fe_index); - for (unsigned int i=0; i 0) - constraints.add_line (local_face_dof_indices[i]); - } - } - } template void @@ -414,13 +366,10 @@ namespace Step36 stiffness_matrix = 0; mass_matrix = 0; - dealii::FullMatrix cell_stiffness_matrix; - dealii::FullMatrix cell_mass_matrix; + FullMatrix cell_stiffness_matrix; + FullMatrix cell_mass_matrix; std::vector local_dof_indices; std::vector potential_values; - std::vector enrichment_values; - std::vector> enrichment_gradients; - hp::FEValues fe_values_hp(fe_collection, q_collection, update_values | update_gradients | update_quadrature_points | update_JxW_values); @@ -439,8 +388,6 @@ namespace Step36 const unsigned int &n_q_points = fe_values.n_quadrature_points; potential_values.resize(n_q_points); - enrichment_values.resize(n_q_points); - enrichment_gradients.resize(n_q_points); local_dof_indices.resize (dofs_per_cell); cell_stiffness_matrix.reinit (dofs_per_cell,dofs_per_cell); @@ -492,17 +439,17 @@ namespace Step36 mass_matrix); } - stiffness_matrix.compress (dealii::VectorOperation::add); - mass_matrix.compress (dealii::VectorOperation::add); + stiffness_matrix.compress (VectorOperation::add); + mass_matrix.compress (VectorOperation::add); } template std::pair EigenvalueProblem::solve() { - dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-9,false,false); - dealii::SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control, - mpi_communicator); + SolverControl solver_control (dof_handler.n_dofs(), 1e-9,false,false); + SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control, + mpi_communicator); eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL); eigensolver.set_problem_type (EPS_GHEP); @@ -533,7 +480,7 @@ namespace Step36 { { std::vector sol (number_of_eigenvalues); - std::vector *> error (number_of_eigenvalues); + std::vector *> error (number_of_eigenvalues); for (unsigned int i = 0; i < number_of_eigenvalues; i++) { @@ -542,8 +489,8 @@ namespace Step36 } hp::QCollection face_quadrature_formula; - face_quadrature_formula.push_back (dealii::QGauss(3)); - face_quadrature_formula.push_back (dealii::QGauss(3)); + face_quadrature_formula.push_back (QGauss(3)); + face_quadrature_formula.push_back (QGauss(3)); KellyErrorEstimator::estimate (dof_handler, face_quadrature_formula, @@ -575,62 +522,12 @@ namespace Step36 triangulation.execute_coarsening_and_refinement (); } - template - class Postprocessor : public DataPostprocessorScalar - { - public: - Postprocessor(const EnrichmentFunction &enrichment); - - virtual - void - compute_derived_quantities_vector (const std::vector > &uh, - const std::vector > > &duh, - const std::vector > > &dduh, - const std::vector > &normals, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const; - - private: - const EnrichmentFunction &enrichment; - }; - - template - Postprocessor::Postprocessor(const EnrichmentFunction &enrichment) - : - DataPostprocessorScalar ("total_solution", - update_values | update_q_points), - enrichment(enrichment) - {} - - template - void - Postprocessor:: - compute_derived_quantities_vector (const std::vector > &uh, - const std::vector > > &/*duh*/, - const std::vector > > &/*dduh*/, - const std::vector > &/*normals*/, - const std::vector > &evaluation_points, - std::vector > &computed_quantities) const - { - const unsigned int n_quadrature_points = uh.size(); - Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); - for (unsigned int q=0; q void EigenvalueProblem::output_results (const unsigned int cycle) const { - dealii::Vector fe_index(triangulation.n_active_cells()); + Vector fe_index(triangulation.n_active_cells()); { - typename dealii::hp::DoFHandler::active_cell_iterator + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (unsigned int index=0; cell!=endc; ++cell, ++index) @@ -645,7 +542,6 @@ namespace Step36 filename += ".vtk"; std::ofstream output (filename.c_str()); - Postprocessor postprocessor(enrichment); // has to live until the DataOut object is destroyed; objects are destroyed in reverse order of declaration DataOut > data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (eigenfunctions_locally_relevant[0], "solution"); @@ -744,14 +640,14 @@ int main (int argc,char **argv) try { - dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); { Step36::EigenvalueProblem step36; - dealii::PETScWrappers::set_option_value("-eps_target","-1.0"); - dealii::PETScWrappers::set_option_value("-st_type","sinvert"); - dealii::PETScWrappers::set_option_value("-st_ksp_type","cg"); - dealii::PETScWrappers::set_option_value("-st_pc_type", "gamg"); - dealii::PETScWrappers::set_option_value("-st_ksp_tol", "1e-11"); + PETScWrappers::set_option_value("-eps_target","-1.0"); + PETScWrappers::set_option_value("-st_type","sinvert"); + PETScWrappers::set_option_value("-st_ksp_type","cg"); + PETScWrappers::set_option_value("-st_pc_type", "gamg"); + PETScWrappers::set_option_value("-st_ksp_tol", "1e-11"); step36.run(); } -- 2.39.5