From b9f30454d70382c129421b688611f4be29a0fc31 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 12 Nov 2017 16:30:29 -0700 Subject: [PATCH] Add a test. --- tests/hp/step-3a_mapping_collection.cc | 293 +++++++++++++++++++++ tests/hp/step-3a_mapping_collection.output | 150 +++++++++++ 2 files changed, 443 insertions(+) create mode 100644 tests/hp/step-3a_mapping_collection.cc create mode 100644 tests/hp/step-3a_mapping_collection.output diff --git a/tests/hp/step-3a_mapping_collection.cc b/tests/hp/step-3a_mapping_collection.cc new file mode 100644 index 0000000000..953c53df53 --- /dev/null +++ b/tests/hp/step-3a_mapping_collection.cc @@ -0,0 +1,293 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// like the _3a test, but using a hp::MappingCollection for +// VectorTools::interpolate_boundary_values(). this function existed +// before, but was not instantiated + + +#include "../tests.h" +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +std::ofstream logfile("output"); + + + +class LaplaceProblem +{ +public: + LaplaceProblem (); + + void run (); + +private: + void make_grid_and_dofs (); + void assemble_system (); + void solve (); + void output_results () const; + + Triangulation<2> triangulation; + hp::FECollection<2> fe; + hp::MappingCollection<2> mappings; + hp::DoFHandler<2> dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + // Although we do not have h-refinement, + // hanging nodes will inevitably appear + // due to different polynomial degrees. + ConstraintMatrix hanging_node_constraints; + + Vector solution; + Vector system_rhs; +}; + + +LaplaceProblem::LaplaceProblem () : + dof_handler (triangulation) +{} + + + +void LaplaceProblem::make_grid_and_dofs () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (2); + deallog << "Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + deallog << "Total number of cells: " + << triangulation.n_cells() + << std::endl; + + hp::DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active (), + endc = dof_handler.end (); + + unsigned int cell_no = 0; + for (; cell != endc; ++cell) + { + if (cell_no < triangulation.n_active_cells () / 2) + cell->set_active_fe_index (1); + else + cell->set_active_fe_index (0); + deallog << "Cell " << cell << " has fe_index=" << cell->active_fe_index() + << std::endl; + ++cell_no; + } + + dof_handler.distribute_dofs (fe); + deallog << "Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + // Create sparsity pattern. + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + + // Create constraints which stem from + // the different polynomial degrees on + // the different elements. + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + + hanging_node_constraints.print (deallog.get_file_stream()); + + hanging_node_constraints.close (); + hanging_node_constraints.condense (sparsity_pattern); + + sparsity_pattern.compress(); + system_matrix.reinit (sparsity_pattern); +} + + + +void LaplaceProblem::assemble_system () +{ + hp::QCollection<2> quadrature_formula(QGauss<2>(6)); + hp::FEValues<2> x_fe_values (fe, quadrature_formula, + update_values | update_gradients | update_JxW_values); + + const unsigned int max_dofs_per_cell = fe.max_dofs_per_cell (); + const unsigned int n_q_points = quadrature_formula[0].size(); + + FullMatrix cell_matrix (max_dofs_per_cell, max_dofs_per_cell); + Vector cell_rhs (max_dofs_per_cell); + + std::vector local_dof_indices (max_dofs_per_cell); + + hp::DoFHandler<2>::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + x_fe_values.reinit (cell); + + const FEValues<2> &fe_values = x_fe_values.get_present_fe_values(); + + cell_matrix = 0; + cell_rhs = 0; + + const unsigned int dofs_per_cell = cell->get_fe ().dofs_per_cell; + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + for (unsigned int i=0; i boundary_values; + Functions::ZeroFunction<2> zero; + std::map*> + b_v_functions {{ + types::boundary_id(0), + &zero + } + }; + + VectorTools::interpolate_boundary_values (mappings, + dof_handler, + b_v_functions, + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +} + + + +void LaplaceProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + + cg.solve (system_matrix, solution, system_rhs, + PreconditionIdentity()); + + solution.print (deallog.get_file_stream()); + + hanging_node_constraints.distribute (solution); +} + + + +void LaplaceProblem::output_results () const +{ + DataOut<2,hp::DoFHandler<2> > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + + data_out.write_gnuplot (deallog.get_file_stream()); +} + + + +void LaplaceProblem::run () +{ + FE_Q<2> fe_1 (1), + fe_2 (2), + fe_3 (3), + fe_4 (4); + + fe.push_back (fe_1); + fe.push_back (fe_2); + fe.push_back (fe_3); + fe.push_back (fe_4); + + mappings.push_back (MappingQ1<2>()); + mappings.push_back (MappingQ1<2>()); + mappings.push_back (MappingQ1<2>()); + mappings.push_back (MappingQ1<2>()); + + make_grid_and_dofs (); + assemble_system (); + solve (); + output_results (); +} + + + +int main () +{ + logfile.precision(6); + + deallog.attach(logfile); + + LaplaceProblem laplace_problem; + laplace_problem.run (); + + return 0; +} diff --git a/tests/hp/step-3a_mapping_collection.output b/tests/hp/step-3a_mapping_collection.output new file mode 100644 index 0000000000..6429376647 --- /dev/null +++ b/tests/hp/step-3a_mapping_collection.output @@ -0,0 +1,150 @@ + +DEAL::Number of active cells: 16 +DEAL::Total number of cells: 21 +DEAL::Cell 2.0 has fe_index=1 +DEAL::Cell 2.1 has fe_index=1 +DEAL::Cell 2.2 has fe_index=1 +DEAL::Cell 2.3 has fe_index=1 +DEAL::Cell 2.4 has fe_index=1 +DEAL::Cell 2.5 has fe_index=1 +DEAL::Cell 2.6 has fe_index=1 +DEAL::Cell 2.7 has fe_index=1 +DEAL::Cell 2.8 has fe_index=0 +DEAL::Cell 2.9 has fe_index=0 +DEAL::Cell 2.10 has fe_index=0 +DEAL::Cell 2.11 has fe_index=0 +DEAL::Cell 2.12 has fe_index=0 +DEAL::Cell 2.13 has fe_index=0 +DEAL::Cell 2.14 has fe_index=0 +DEAL::Cell 2.15 has fe_index=0 +DEAL::Number of degrees of freedom: 55 + 17 40: 0.500000 + 17 41: 0.500000 + 20 41: 0.500000 + 20 44: 0.500000 + 35 44: 0.500000 + 35 49: 0.500000 + 38 49: 0.500000 + 38 51: 0.500000 +DEAL:cg::Starting value 0.634648 +DEAL:cg::Convergence step 16 value 0 +0.000e+00 0.000e+00 0.000e+00 1.814e-01 0.000e+00 1.122e-01 0.000e+00 1.109e-01 7.270e-02 0.000e+00 2.303e-01 1.395e-01 0.000e+00 2.170e-01 1.331e-01 0.000e+00 2.157e-01 0.000e+00 1.299e-01 2.790e-01 0.000e+00 2.646e-01 0.000e+00 1.814e-01 1.122e-01 0.000e+00 2.170e-01 1.331e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 1.109e-01 7.270e-02 2.157e-01 0.000e+00 2.646e-01 0.000e+00 0.000e+00 1.299e-01 0.000e+00 2.421e-01 0.000e+00 1.926e-01 3.080e-01 2.409e-01 0.000e+00 0.000e+00 0.000e+00 2.421e-01 1.926e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +-1.00000 -1.00000 0.00000 +-0.500000 -1.00000 0.00000 + +-1.00000 -0.500000 0.00000 +-0.500000 -0.500000 0.181364 + + +-0.500000 -1.00000 0.00000 +0.00000 -1.00000 0.00000 + +-0.500000 -0.500000 0.181364 +0.00000 -0.500000 0.230304 + + +-1.00000 -0.500000 0.00000 +-0.500000 -0.500000 0.181364 + +-1.00000 0.00000 0.00000 +-0.500000 0.00000 0.242082 + + +-0.500000 -0.500000 0.181364 +0.00000 -0.500000 0.230304 + +-0.500000 0.00000 0.242082 +0.00000 0.00000 0.307982 + + +0.00000 -1.00000 0.00000 +0.500000 -1.00000 0.00000 + +0.00000 -0.500000 0.230304 +0.500000 -0.500000 0.181364 + + +0.500000 -1.00000 0.00000 +1.00000 -1.00000 0.00000 + +0.500000 -0.500000 0.181364 +1.00000 -0.500000 0.00000 + + +0.00000 -0.500000 0.230304 +0.500000 -0.500000 0.181364 + +0.00000 0.00000 0.307982 +0.500000 0.00000 0.242082 + + +0.500000 -0.500000 0.181364 +1.00000 -0.500000 0.00000 + +0.500000 0.00000 0.242082 +1.00000 0.00000 0.00000 + + +-1.00000 0.00000 0.00000 +-0.500000 0.00000 0.242082 + +-1.00000 0.500000 0.00000 +-0.500000 0.500000 0.192624 + + +-0.500000 0.00000 0.242082 +0.00000 0.00000 0.307982 + +-0.500000 0.500000 0.192624 +0.00000 0.500000 0.240924 + + +-1.00000 0.500000 0.00000 +-0.500000 0.500000 0.192624 + +-1.00000 1.00000 0.00000 +-0.500000 1.00000 0.00000 + + +-0.500000 0.500000 0.192624 +0.00000 0.500000 0.240924 + +-0.500000 1.00000 0.00000 +0.00000 1.00000 0.00000 + + +0.00000 0.00000 0.307982 +0.500000 0.00000 0.242082 + +0.00000 0.500000 0.240924 +0.500000 0.500000 0.192624 + + +0.500000 0.00000 0.242082 +1.00000 0.00000 0.00000 + +0.500000 0.500000 0.192624 +1.00000 0.500000 0.00000 + + +0.00000 0.500000 0.240924 +0.500000 0.500000 0.192624 + +0.00000 1.00000 0.00000 +0.500000 1.00000 0.00000 + + +0.500000 0.500000 0.192624 +1.00000 0.500000 0.00000 + +0.500000 1.00000 0.00000 +1.00000 1.00000 0.00000 + + -- 2.39.5