From: Guido Kanschat Date: Mon, 24 Jun 2013 16:08:05 +0000 (+0000) Subject: reduce example X-Git-Tag: v8.0.0~235 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f20a0e0b7b3ca16ac22c6dcc3d16307c4c252d66;p=dealii.git reduce example git-svn-id: https://svn.dealii.org/trunk@29883 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/mpi/flux_edge_01.cc b/tests/mpi/flux_edge_01.cc index 070faebecf..5b68d0e1f3 100644 --- a/tests/mpi/flux_edge_01.cc +++ b/tests/mpi/flux_edge_01.cc @@ -1,13 +1,16 @@ -/* Author: Guido Kanschat, Texas A&M University, 2009 */ - -/* $Id: step-39.cc 28552 2013-02-25 15:18:15Z kanschat $ */ -/* */ -/* Copyright (C) 2010-2013 by the deal.II authors */ -/* */ -/* This file is subject to QPL and may not be distributed */ -/* without copyright and license information. Please refer */ -/* to the file deal.II/doc/license.html for the text and */ -/* further information on this license. */ +//---------------------------------------------------------------------- +// $Id: mesh_worker_01.cc 28614 2013-02-27 15:20:34Z kanschat $ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + +// Test problem in MGTools::make_flux_sparsity_pattern_edge #include "../tests.h" #include @@ -55,283 +58,6 @@ namespace Step39 { using namespace dealii; - Functions::SlitSingularityFunction<2> exact_solution; - - - - - template - class MatrixIntegrator : public MeshWorker::LocalIntegrator - { - public: - void cell(MeshWorker::DoFInfo &dinfo, - typename MeshWorker::IntegrationInfo &info) const; - void boundary(MeshWorker::DoFInfo &dinfo, - typename MeshWorker::IntegrationInfo &info) const; - void face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - typename MeshWorker::IntegrationInfo &info1, - typename MeshWorker::IntegrationInfo &info2) const; - }; - - - template - void MatrixIntegrator::cell( - MeshWorker::DoFInfo &dinfo, - typename MeshWorker::IntegrationInfo &info) const - { - LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values()); - } - - - template - void MatrixIntegrator::boundary( - MeshWorker::DoFInfo &dinfo, - typename MeshWorker::IntegrationInfo &info) const - { - const unsigned int deg = info.fe_values(0).get_fe().tensor_degree(); - LocalIntegrators::Laplace::nitsche_matrix( - dinfo.matrix(0,false).matrix, info.fe_values(0), - LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg)); - } - - template - void MatrixIntegrator::face( - MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - typename MeshWorker::IntegrationInfo &info1, - typename MeshWorker::IntegrationInfo &info2) const - { - const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree(); - LocalIntegrators::Laplace::ip_matrix( - dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, - dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix, - info1.fe_values(0), info2.fe_values(0), - LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); - } - - template - class RHSIntegrator : public MeshWorker::LocalIntegrator - { - public: - void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; - void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; - void face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - typename MeshWorker::IntegrationInfo &info1, - typename MeshWorker::IntegrationInfo &info2) const; - }; - - - template - void RHSIntegrator::cell(MeshWorker::DoFInfo &, typename MeshWorker::IntegrationInfo &) const - {} - - - template - void RHSIntegrator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const - { - const FEValuesBase &fe = info.fe_values(); - Vector &local_vector = dinfo.vector(0).block(0); - - std::vector boundary_values(fe.n_quadrature_points); - exact_solution.value_list(fe.get_quadrature_points(), boundary_values); - - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); - - for (unsigned k=0; k - void RHSIntegrator::face(MeshWorker::DoFInfo &, - MeshWorker::DoFInfo &, - typename MeshWorker::IntegrationInfo &, - typename MeshWorker::IntegrationInfo &) const - {} - - - template - class Estimator : public MeshWorker::LocalIntegrator - { - public: - void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; - void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; - void face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - typename MeshWorker::IntegrationInfo &info1, - typename MeshWorker::IntegrationInfo &info2) const; - }; - - - template - void Estimator::cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const - { - const FEValuesBase &fe = info.fe_values(); - - const std::vector > &DDuh = info.hessians[0][0]; - for (unsigned k=0; kdiameter() * trace(DDuh[k]); - dinfo.value(0) += t*t * fe.JxW(k); - } - dinfo.value(0) = std::sqrt(dinfo.value(0)); - } - - template - void Estimator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const - { - const FEValuesBase &fe = info.fe_values(); - - std::vector boundary_values(fe.n_quadrature_points); - exact_solution.value_list(fe.get_quadrature_points(), boundary_values); - - const std::vector &uh = info.values[0][0]; - - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); - - for (unsigned k=0; k - void Estimator::face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - typename MeshWorker::IntegrationInfo &info1, - typename MeshWorker::IntegrationInfo &info2) const - { - const FEValuesBase &fe = info1.fe_values(); - const std::vector &uh1 = info1.values[0][0]; - const std::vector &uh2 = info2.values[0][0]; - const std::vector > &Duh1 = info1.gradients[0][0]; - const std::vector > &Duh2 = info2.gradients[0][0]; - - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); - const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); - const double penalty = penalty1 + penalty2; - const double h = dinfo1.face->measure(); - - for (unsigned k=0; k - class ErrorIntegrator : public MeshWorker::LocalIntegrator - { - public: - void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; - void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; - void face(MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - typename MeshWorker::IntegrationInfo &info1, - typename MeshWorker::IntegrationInfo &info2) const; - }; - - - template - void ErrorIntegrator::cell( - MeshWorker::DoFInfo &dinfo, - typename MeshWorker::IntegrationInfo &info) const - { - const FEValuesBase &fe = info.fe_values(); - std::vector > exact_gradients(fe.n_quadrature_points); - std::vector exact_values(fe.n_quadrature_points); - - exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients); - exact_solution.value_list(fe.get_quadrature_points(), exact_values); - - const std::vector > &Duh = info.gradients[0][0]; - const std::vector &uh = info.values[0][0]; - - for (unsigned k=0; k - void ErrorIntegrator::boundary( - MeshWorker::DoFInfo &dinfo, - typename MeshWorker::IntegrationInfo &info) const - { - const FEValuesBase &fe = info.fe_values(); - - std::vector exact_values(fe.n_quadrature_points); - exact_solution.value_list(fe.get_quadrature_points(), exact_values); - - const std::vector &uh = info.values[0][0]; - - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); - - for (unsigned k=0; k - void ErrorIntegrator::face( - MeshWorker::DoFInfo &dinfo1, - MeshWorker::DoFInfo &dinfo2, - typename MeshWorker::IntegrationInfo &info1, - typename MeshWorker::IntegrationInfo &info2) const - { - const FEValuesBase &fe = info1.fe_values(); - const std::vector &uh1 = info1.values[0][0]; - const std::vector &uh2 = info2.values[0][0]; - - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); - const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); - const double penalty = penalty1 + penalty2; - - for (unsigned k=0; k class InteriorPenaltyProblem { @@ -344,13 +70,6 @@ namespace Step39 private: void setup_system (); - void assemble_matrix (); - void assemble_mg_matrix (); - void assemble_right_hand_side (); - void error (); - double estimate (); - void solve (); - void output_results (const unsigned int cycle) const; parallel::distributed::Triangulation triangulation; const MappingQ1 mapping; @@ -360,9 +79,6 @@ namespace Step39 IndexSet locally_relevant_set; TrilinosWrappers::SparseMatrix matrix; - TrilinosWrappers::MPI::Vector solution; - TrilinosWrappers::MPI::Vector right_hand_side; - BlockVector estimates; MGLevelObject mg_matrix; @@ -394,9 +110,7 @@ namespace Step39 dof_handler.distribute_mg_dofs (fe); DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_set); - solution.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); - right_hand_side.reinit(dof_handler.locally_owned_dofs(), MPI_COMM_WORLD); - + CompressedSimpleSparsityPattern c_sparsity(dof_handler.n_dofs(), dof_handler.n_dofs()); DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity); matrix.reinit(dof_handler.locally_owned_dofs(), c_sparsity, MPI_COMM_WORLD, true); @@ -438,268 +152,6 @@ namespace Step39 } - template - void - InteriorPenaltyProblem::assemble_matrix() - { - MeshWorker::IntegrationInfoBox info_box; - UpdateFlags update_flags = update_values | update_gradients; - info_box.add_update_flags_all(update_flags); - info_box.initialize(fe, mapping); - - MeshWorker::DoFInfo dof_info(dof_handler); - - MeshWorker::Assembler::MatrixSimple assembler; - assembler.initialize(matrix); - - FilteredIterator::active_cell_iterator> - begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active()); - FilteredIterator::active_cell_iterator> - end(IteratorFilters::LocallyOwnedCell(), dof_handler.end()); - - MatrixIntegrator integrator; - MeshWorker::integration_loop( - begin, end, - dof_info, info_box, - integrator, assembler); - - matrix.compress(VectorOperation::add); - } - - - template - void - InteriorPenaltyProblem::assemble_mg_matrix() - { - MeshWorker::IntegrationInfoBox info_box; - UpdateFlags update_flags = update_values | update_gradients; - info_box.add_update_flags_all(update_flags); - info_box.initialize(fe, mapping); - - MeshWorker::DoFInfo dof_info(dof_handler); - - MeshWorker::Assembler::MGMatrixSimple assembler; - assembler.initialize(mg_matrix); - assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); - - FilteredIterator::level_cell_iterator> - begin(IteratorFilters::LocallyOwnedLevelCell(), dof_handler.begin()); - FilteredIterator::level_cell_iterator> - end(IteratorFilters::LocallyOwnedLevelCell(), dof_handler.end()); - - MatrixIntegrator integrator; - MeshWorker::integration_loop ( - begin, end, - dof_info, info_box, - integrator, assembler); - - for (unsigned int level=mg_matrix.min_level(); level <= mg_matrix.max_level();++level) - { - mg_matrix[level].compress(VectorOperation::add); - if (level > mg_matrix.min_level()) - { - mg_matrix_dg_up[level].compress(VectorOperation::add); - mg_matrix_dg_down[level].compress(VectorOperation::add); - } - } - } - - - template - void - InteriorPenaltyProblem::assemble_right_hand_side() - { - MeshWorker::IntegrationInfoBox info_box; - UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; - info_box.add_update_flags_all(update_flags); - info_box.initialize(fe, mapping); - - MeshWorker::DoFInfo dof_info(dof_handler); - - MeshWorker::Assembler::ResidualSimple assembler; - NamedData data; - TrilinosWrappers::MPI::Vector *rhs = &right_hand_side; - data.add(rhs, "RHS"); - assembler.initialize(data); - - FilteredIterator::active_cell_iterator> - begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active()); - FilteredIterator::active_cell_iterator> - end(IteratorFilters::LocallyOwnedCell(), dof_handler.end()); - - - RHSIntegrator integrator; - MeshWorker::integration_loop(begin, end, - dof_info, info_box, - integrator, assembler); - - right_hand_side.compress(VectorOperation::add); - right_hand_side *= -1.; - } - - - template - void - InteriorPenaltyProblem::solve() - { - SolverControl control(1000, 1.e-12); - SolverCG solver(control); - - MGTransferPrebuilt mg_transfer; - mg_transfer.build_matrices(dof_handler); - - SolverControl coarse_solver_control (1000, 1e-10, false, false); - SolverCG coarse_solver(coarse_solver_control); - PreconditionIdentity identity; - TrilinosWrappers::SparseMatrix coarse_matrix = mg_matrix[0]; - MGCoarseGridLACIteration,TrilinosWrappers::MPI::Vector> - coarse_grid_solver(coarse_solver, coarse_matrix, identity); - - typedef TrilinosWrappers::PreconditionJacobi Smoother; - MGSmootherPrecondition mg_smoother; - mg_smoother.initialize(mg_matrix, .7); - mg_smoother.set_steps(2); - - MGMatrix mgmatrix(&mg_matrix); - MGMatrix mgdown(&mg_matrix_dg_down); - MGMatrix mgup(&mg_matrix_dg_up); - - Multigrid mg(dof_handler, mgmatrix, - coarse_grid_solver, mg_transfer, - mg_smoother, mg_smoother); - mg.set_edge_flux_matrices(mgdown, mgup); - - PreconditionMG > - preconditioner(dof_handler, mg, mg_transfer); - solver.solve(matrix, solution, right_hand_side, preconditioner); - } - - - template - double - InteriorPenaltyProblem::estimate() - { - TrilinosWrappers::MPI::Vector ghost; - ghost.reinit(locally_relevant_set, MPI_COMM_WORLD); - ghost = solution; - - std::vector old_user_indices; - triangulation.save_user_indices(old_user_indices); - - estimates.block(0).reinit(triangulation.n_active_cells()); - - unsigned int i=0; - for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); - cell != triangulation.end(); ++cell,++i) - cell->set_user_index(i); - - MeshWorker::IntegrationInfoBox info_box; - const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; - info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); - - NamedData solution_data; - solution_data.add(&ghost, "solution"); - - info_box.cell_selector.add("solution", false, false, true); - info_box.boundary_selector.add("solution", true, true, false); - info_box.face_selector.add("solution", true, true, false); - - info_box.add_update_flags_boundary(update_quadrature_points); - info_box.initialize(fe, mapping, solution_data); - - MeshWorker::DoFInfo dof_info(dof_handler); - - MeshWorker::Assembler::CellsAndFaces assembler; - NamedData* > out_data; - BlockVector *est = &estimates; - out_data.add(est, "cells"); - assembler.initialize(out_data, false); - - FilteredIterator::active_cell_iterator> - begin(IteratorFilters::LocallyOwnedCell(), dof_handler.begin_active()); - FilteredIterator::active_cell_iterator> - end(IteratorFilters::LocallyOwnedCell(), dof_handler.end()); - - Estimator integrator; - MeshWorker::integration_loop ( - begin, end, - dof_info, info_box, - integrator, assembler); - - triangulation.load_user_indices(old_user_indices); - return estimates.block(0).l2_norm(); - } - - - // template - // void - // InteriorPenaltyProblem::error() - // { - // BlockVector errors(2); - // errors.block(0).reinit(triangulation.n_active_cells()); - // errors.block(1).reinit(triangulation.n_active_cells()); - // unsigned int i=0; - // for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(); - // cell != triangulation.end(); ++cell,++i) - // cell->set_user_index(i); - - // MeshWorker::IntegrationInfoBox info_box; - // const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; - // info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); - - // NamedData solution_data; - // solution_data.add(&solution, "solution"); - - // info_box.cell_selector.add("solution", true, true, false); - // info_box.boundary_selector.add("solution", true, false, false); - // info_box.face_selector.add("solution", true, false, false); - - // info_box.add_update_flags_cell(update_quadrature_points); - // info_box.add_update_flags_boundary(update_quadrature_points); - // info_box.initialize(fe, mapping, solution_data); - - // MeshWorker::DoFInfo dof_info(dof_handler); - - // MeshWorker::Assembler::CellsAndFaces assembler; - // NamedData* > out_data; - // BlockVector *est = &errors; - // out_data.add(est, "cells"); - // assembler.initialize(out_data, false); - - // ErrorIntegrator integrator; - // MeshWorker::integration_loop ( - // dof_handler.begin_active(), dof_handler.end(), - // dof_info, info_box, - // integrator, assembler); - - // deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl; - // deallog << "L2-error: " << errors.block(1).l2_norm() << std::endl; - // } - - - // template - // void InteriorPenaltyProblem::output_results (const unsigned int cycle) const - // { - // char *fn = new char[100]; - // sprintf(fn, "step-39/sol-%02d", cycle); - - // std::string filename(fn); - // filename += ".gnuplot"; - // deallog << "Writing solution to <" << filename << ">..." - // << std::endl << std::endl; - // std::ofstream gnuplot_output (filename.c_str()); - - // DataOut data_out; - // data_out.attach_dof_handler (dof_handler); - // data_out.add_data_vector (solution, "u"); - // data_out.add_data_vector (estimates.block(0), "est"); - - // data_out.build_patches (); - - // data_out.write_gnuplot(gnuplot_output); - // } - template void InteriorPenaltyProblem::run(unsigned int n_steps) @@ -708,16 +160,13 @@ namespace Step39 for (unsigned int s=0; sset_refine_flag(); - // parallel::distributed::GridRefinement::refine_and_coarsen_fixed_fraction (triangulation, - // estimates.block(0), - // 0.5, 0.0); triangulation.execute_coarsening_and_refinement (); } @@ -730,18 +179,6 @@ namespace Step39 for (unsigned int l=0; l fe1(2); + FE_DGQ<2> fe1(0); InteriorPenaltyProblem<2> test1(fe1); test1.run(20); }