From: Wolfgang Bangerth Date: Thu, 24 May 2007 01:28:41 +0000 (+0000) Subject: New test X-Git-Tag: v8.0.0~10329 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=02a0f4075c978e4f67fe4ded466de2a3afa03083;p=dealii.git New test git-svn-id: https://svn.dealii.org/trunk@14694 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/fail/hp-crash_18.cc b/tests/fail/hp-crash_18.cc new file mode 100644 index 0000000000..8282ce3266 --- /dev/null +++ b/tests/fail/hp-crash_18.cc @@ -0,0 +1,735 @@ +//---------------------------- crash_18.cc --------------------------- +// $Id: crash_18.cc 12732 2006-03-28 23:15:45Z wolf $ +// Version: $Name$ +// +// Copyright (C) 2006, 2007 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. +// +//---------------------------- crash_18.cc --------------------------- + + +// a modified version of step-27 that crashes because of circular +// constraints. like crash_17, but only check the second cycle. there, +// we have dofs that are constrained against itself, and we don't +// currently detect this... + +char logname[] = "hp-crash_18/output"; + + +#include "../tests.h" + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include + + // Finally, this is as in previous + // programs: +using namespace dealii; + +template +class LaplaceProblem +{ + public: + LaplaceProblem (); + ~LaplaceProblem (); + + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void create_coarse_grid (); + void refine_grid (); + void estimate_smoothness (Vector &smoothness_indicators) const; + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + hp::QCollection face_quadrature_collection; + + ConstraintMatrix hanging_node_constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + Timer distr, condense, hang; +}; + + + +template +class RightHandSide : public Function +{ + public: + RightHandSide () : Function () {}; + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double product = 1; + for (unsigned int d=0; d +LaplaceProblem::LaplaceProblem () : + dof_handler (triangulation) +{ + for (unsigned int degree=2; degree<5; ++degree) + { + fe_collection.push_back (FE_Q(degree)); + quadrature_collection.push_back (QGauss(degree+2)); + face_quadrature_collection.push_back (QGauss(degree+2)); + } +} + + +template +LaplaceProblem::~LaplaceProblem () +{ + dof_handler.clear (); +} + +template +void LaplaceProblem::setup_system () +{ + distr.reset(); + distr.start(); + dof_handler.distribute_dofs (fe_collection); + distr.stop(); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + hanging_node_constraints.clear (); + hang.reset(); + hang.start(); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + + // check for self-referential constraints + + hanging_node_constraints.close (); + hang.stop(); + + if (dim < 3) + { + 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); + + condense.reset(); + condense.start(); + hanging_node_constraints.condense (sparsity_pattern); + condense.stop(); + sparsity_pattern.compress(); + } + else + { + CompressedSparsityPattern csp (dof_handler.n_dofs(), + dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, csp); + + condense.reset(); + condense.start(); + hanging_node_constraints.condense (csp); + condense.stop(); + + sparsity_pattern.copy_from (csp); + } + + system_matrix.reinit (sparsity_pattern); +} + +template +void LaplaceProblem::assemble_system () +{ + hp::FEValues hp_fe_values (fe_collection, + quadrature_collection, + update_values | update_gradients | + update_q_points | update_JxW_values); + + const RightHandSide rhs_function; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + cell_matrix = 0; + cell_rhs = 0; + + hp_fe_values.reinit (cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values (); + + std::vector rhs_values (fe_values.n_quadrature_points); + rhs_function.value_list (fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction(), + boundary_values); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); +} + +template +void LaplaceProblem::solve () +{ + SolverControl solver_control (system_rhs.size(), + 1e-8*system_rhs.l2_norm()); + SolverCG<> cg (solver_control); + + PreconditionSSOR<> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + + hanging_node_constraints.distribute (solution); +} + + +unsigned int +int_pow (const unsigned int x, + const unsigned int n) +{ + unsigned int p=1; + for (unsigned int i=0; i +void +LaplaceProblem:: +estimate_smoothness (Vector &smoothness_indicators) const +{ + const unsigned int N = 4; + + // form all the Fourier vectors + // that we want to + // consider. exclude k=0 to avoid + // problems with |k|^{-mu} and also + // logarithms of |k| + std::vector > k_vectors; + std::vector k_vectors_magnitude; + switch (dim) + { + case 2: + { + for (unsigned int i=0; i(deal_II_numbers::PI * i, + deal_II_numbers::PI * j)); + k_vectors_magnitude.push_back (i*i+j*j); + } + + break; + } + + case 3: + { + for (unsigned int i=0; i(deal_II_numbers::PI * i, + deal_II_numbers::PI * j, + deal_II_numbers::PI * k)); + k_vectors_magnitude.push_back (i*i+j*j+k*k); + } + + break; + } + + default: + Assert (false, ExcNotImplemented()); + } + + const unsigned n_fourier_modes = k_vectors.size(); + std::vector ln_k (n_fourier_modes); + for (unsigned int i=0; i base_quadrature (2); + QIterated quadrature (base_quadrature, N); + + std::vector > > + fourier_transform_matrices (fe_collection.size()); + for (unsigned int fe=0; fe sum = 0; + for (unsigned int q=0; q x_q = quadrature.point(q); + sum += std::exp(std::complex(0,1) * + (k_vectors[k] * x_q)) * + fe_collection[fe].shape_value(i,x_q) * + quadrature.weight(q); + } + fourier_transform_matrices[fe](k,i) + = sum / std::pow(2*deal_II_numbers::PI, 1.*dim/2); + } + } + + // the next thing is to loop over + // all cells and do our work there, + // i.e. to locally do the Fourier + // transform and estimate the decay + // coefficient + std::vector > fourier_coefficients (n_fourier_modes); + Vector local_dof_values; + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + { + local_dof_values.reinit (cell->get_fe().dofs_per_cell); + cell->get_dof_values (solution, local_dof_values); + + // first compute the Fourier + // transform of the local + // solution + std::fill (fourier_coefficients.begin(), fourier_coefficients.end(), 0); + for (unsigned int f=0; fget_fe().dofs_per_cell; ++i) + fourier_coefficients[f] += + fourier_transform_matrices[cell->active_fe_index()](f,i) + * + local_dof_values(i); + + // enter the Fourier + // coefficients into a map with + // the k-magnitudes, to make + // sure that we get only the + // largest magnitude for each + // value of |k| + std::map k_to_max_U_map; + for (unsigned int f=0; f +void LaplaceProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + KellyErrorEstimator::estimate (dof_handler, + face_quadrature_collection, + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + + Vector smoothness_indicators (triangulation.n_active_cells()); + estimate_smoothness (smoothness_indicators); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + + float max_smoothness = 0, + min_smoothness = smoothness_indicators.linfty_norm(); + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + if (cell->refine_flag_set()) + { + max_smoothness = std::max (max_smoothness, + smoothness_indicators(index)); + min_smoothness = std::min (min_smoothness, + smoothness_indicators(index)); + } + } + const float cutoff_smoothness = (max_smoothness + min_smoothness) / 2; + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + if (cell->refine_flag_set() + && + (smoothness_indicators(index) > cutoff_smoothness) + && + !(cell->active_fe_index() == fe_collection.size() - 1)) + { + cell->clear_refine_flag(); + cell->set_active_fe_index (std::min (cell->active_fe_index() + 1, + fe_collection.size() - 1)); + } + } + + triangulation.execute_coarsening_and_refinement (); +} + + + +template +void LaplaceProblem::output_results (const unsigned int cycle) const +{ + Assert (cycle < 10, ExcNotImplemented()); + + { + const std::string filename = "grid-" + + Utilities::int_to_string (cycle, 2) + + ".eps"; + std::ofstream output (filename.c_str()); + + GridOut grid_out; + grid_out.write_eps (triangulation, output); + } + + { + Vector smoothness_indicators (triangulation.n_active_cells()); + estimate_smoothness (smoothness_indicators); + + Vector smoothness_field (dof_handler.n_dofs()); + DoFTools::distribute_cell_to_dof_vector (dof_handler, + smoothness_indicators, + smoothness_field); + + Vector fe_indices (triangulation.n_active_cells()); + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + { + + fe_indices(index) = cell->active_fe_index(); +// smoothness_indicators(index) *= std::sqrt(cell->diameter()); + } + } + + const std::string filename = "solution-" + + Utilities::int_to_string (cycle, 2) + + ".vtk"; + DataOut > data_out; + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.add_data_vector (smoothness_indicators, "smoothness1"); + data_out.add_data_vector (smoothness_field, "smoothness2"); + data_out.add_data_vector (fe_indices, "fe_index"); + data_out.build_patches (); + + std::ofstream output (filename.c_str()); + data_out.write_vtk (output); + } +} + + +template <> +void LaplaceProblem<2>::create_coarse_grid () +{ + const unsigned int dim = 2; + + static const Point<2> vertices_1[] + = { Point<2> (-1., -1.), + Point<2> (-1./2, -1.), + Point<2> (0., -1.), + Point<2> (+1./2, -1.), + Point<2> (+1, -1.), + + Point<2> (-1., -1./2.), + Point<2> (-1./2, -1./2.), + Point<2> (0., -1./2.), + Point<2> (+1./2, -1./2.), + Point<2> (+1, -1./2.), + + Point<2> (-1., 0.), + Point<2> (-1./2, 0.), + Point<2> (+1./2, 0.), + Point<2> (+1, 0.), + + Point<2> (-1., 1./2.), + Point<2> (-1./2, 1./2.), + Point<2> (0., 1./2.), + Point<2> (+1./2, 1./2.), + Point<2> (+1, 1./2.), + + Point<2> (-1., 1.), + Point<2> (-1./2, 1.), + Point<2> (0., 1.), + Point<2> (+1./2, 1.), + Point<2> (+1, 1.) }; + const unsigned int + n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]); + const std::vector > vertices (&vertices_1[0], + &vertices_1[n_vertices]); + static const int cell_vertices[][GeometryInfo::vertices_per_cell] + = {{0, 1, 5, 6}, + {1, 2, 6, 7}, + {2, 3, 7, 8}, + {3, 4, 8, 9}, + {5, 6, 10, 11}, + {8, 9, 12, 13}, + {10, 11, 14, 15}, + {12, 13, 17, 18}, + {14, 15, 19, 20}, + {15, 16, 20, 21}, + {16, 17, 21, 22}, + {17, 18, 22, 23}}; + const unsigned int + n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]); + + std::vector > cells (n_cells, CellData()); + for (unsigned int i=0; i::vertices_per_cell; + ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + triangulation.create_triangulation (vertices, + cells, + SubCellData()); + triangulation.refine_global (3); +} + + + +template <> +void LaplaceProblem<3>::create_coarse_grid () +{ + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (1); +} + + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<2; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + create_coarse_grid (); + else + refine_grid (); + + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + Timer all; + all.start(); + setup_system (); + + deallog << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + deallog << " Number of constraints : " + << hanging_node_constraints.n_constraints() + << std::endl; + + assemble_system (); + solve (); + all.stop(); + + } +} + +int main () +{ + std::ofstream logfile(logname); + logfile.precision (3); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + try + { + LaplaceProblem<3> laplace_problem_2d; + laplace_problem_2d.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +}