diff --git a/include/bem_problem.h b/include/bem_problem.h index 0545a8e7..961fa9c5 100644 --- a/include/bem_problem.h +++ b/include/bem_problem.h @@ -128,6 +128,8 @@ class BEMProblem : public ParameterAcceptor /// in our computations (assemble system and compute_normals-gradients). void reinit(); + const Quadrature & get_singular_quadrature(const unsigned int index) const; + /// This function compute a very specific case, a double node that has a /// dirichlet-dirichlet condition. In this case there is a constraint for /// the normal derivative since we want a conitnuos velocity thus a conitnuos diff --git a/source/bem_problem.cc b/source/bem_problem.cc index 7ecd8fd2..d3b472d0 100644 --- a/source/bem_problem.cc +++ b/source/bem_problem.cc @@ -316,6 +316,44 @@ void BEMProblem::reinit() } + +template<> +const Quadrature<2> &BEMProblem<3>::get_singular_quadrature(const unsigned int index) const +{ + Assert(index < fe->dofs_per_cell, + ExcIndexRange(0, fe->dofs_per_cell, index)); + + + + static std::vector > quadratures; + { + if (quadratures.size() == 0) + for (unsigned int i=0; idofs_per_cell; ++i) + { + quadratures.push_back(QSplit<2> (QDuffy (singular_quadrature_order,1.),fe->get_unit_support_points()[i])); + } + } + + return quadratures[index]; + +} + +template<> +const Quadrature<1> &BEMProblem<2>::get_singular_quadrature(const unsigned int index) const +{ + Assert(index < fe->dofs_per_cell, + ExcIndexRange(0, fe->dofs_per_cell, index)); + + static std::vector > quadratures; + if (quadratures.size() == 0) + for (unsigned int i=0; idofs_per_cell; ++i) + { + quadratures.push_back(QTelles<1>(singular_quadrature_order, + fe->get_unit_support_points()[i])); + } + return quadratures[index]; +} + template void BEMProblem::declare_parameters (ParameterHandler &prm) { @@ -566,25 +604,6 @@ void BEMProblem::assemble_system() dirichlet_matrix = 0; - std::vector > sing_quadratures; - for (unsigned int i=0; idofs_per_cell; ++i) - { - if (fe->degree > 1) - { - sing_quadratures.push_back(QIterated(QGauss<1> (singular_quadrature_order),fe->degree)); - } - else - { - sing_quadratures.push_back - (QTelles(singular_quadrature_order, - fe->get_unit_support_points()[i])); - // - // Usage of alternative singular quadrature formula - // (QGaussOneOverR(singular_quadrature_order, - // fe->get_unit_support_points()[i],true)); - } - } - // Next, we initialize an FEValues // object with the quadrature @@ -945,8 +964,7 @@ void BEMProblem::assemble_system() const Quadrature * singular_quadrature - = dynamic_cast*>( - &sing_quadratures[singular_index]); + = &(get_singular_quadrature(singular_index)); Assert(singular_quadrature, ExcInternalError()); FEValues fe_v_singular (*mapping, *fe, *singular_quadrature, diff --git a/tests/tests.h b/tests/tests.h index a187bebe..f4bb2438 100644 --- a/tests/tests.h +++ b/tests/tests.h @@ -1,7 +1,6 @@ // --------------------------------------------------------------------- -// $Id$ // -// Copyright (C) 2004 - 2014 by the deal.II authors +// Copyright (C) 2004 - 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -14,28 +13,54 @@ // // --------------------------------------------------------------------- -#ifndef __tests_tests_h -#define __tests_tests_h +#ifndef dealii_tests_h +#define dealii_tests_h // common definitions used in all the tests #include #include #include +#include #include #include #include #include -#include +#include +#include #include #include #include #include #include #include -#include +#include + +#if defined(DEBUG) && defined(DEAL_II_HAVE_FP_EXCEPTIONS) +# include +#endif + + +// silence extra diagnostics in the testsuite +#ifdef DEAL_II_DISABLE_EXTRA_DIAGNOSTICS +DEAL_II_DISABLE_EXTRA_DIAGNOSTICS +#endif +#ifdef DEAL_II_MSVC +// Under windows tests will hang and show a debugging dialog box from the +// debug CRT if an exception is encountered. Disable this: +#include + +struct DisableWindowsDebugRuntimeDialog +{ + DisableWindowsDebugRuntimeDialog () + { + _set_abort_behavior( 0, _WRITE_ABORT_MSG); + } +} deal_II_windows_crt_dialog; +#endif + // implicitly use the deal.II namespace everywhere, without us having to say // so in each and every testcase using namespace dealii; @@ -43,27 +68,102 @@ using namespace dealii; // ------------------------------ Utility functions used in tests ----------------------- +/** + * Go through the input stream @p in and filter out binary data for the key @p key . + * The filtered stream is returned in @p out. + */ +void filter_out_xml_key(std::istream &in, const std::string &key, std::ostream &out) +{ + std::string line; + bool found = false; + const std::string opening = "<" + key; + const std::string closing = "" but keep things after "", 0)+1) + line.substr(pos); + found = false; + } + else + line = line.substr(0, line.find(">", 0)+1); + out << line << std::endl; + } + else if (line.find(closing) != std::string::npos) + { + found = false; + // remove everything before "<" + line = line.substr(line.find("<",0)); + out << line << std::endl; + } + else if (!found) + out << line << std::endl; + } +} + +/** + * A function to return real part of the number and check that + * its imaginary part is zero. + */ +#ifdef DEAL_II_WITH_PETSC +#include +PetscReal get_real_assert_zero_imag(const PETScWrappers::internal::VectorReference &a) +{ + Assert (a.imag() == 0.0, ExcInternalError()); + return a.real(); +} +#endif + +template +number get_real_assert_zero_imag(const std::complex &a) +{ + Assert (a.imag() == 0.0, ExcInternalError()); + return a.real(); +} + +template +number get_real_assert_zero_imag(const number &a) +{ + return a; +} + + // Cygwin has a different implementation for rand() which causes many tests to fail. // This here is a reimplementation that gives the same sequence of numbers as a program // that uses rand() on a typical linux machine. // we put this into a namespace to not conflict with stdlib namespace Testing { - int rand(bool reseed=false, int seed=1) throw() + int rand(const bool reseed=false, + const int seed=1) { static int r[32]; static int k; static bool inited=false; + if (!inited || reseed) { //srand treats a seed 0 as 1 for some reason r[0]=(seed==0)?1:seed; + long int word = r[0]; for (int i=1; i<31; i++) { - r[i] = (16807LL * r[i-1]) % 2147483647; - if (r[i] < 0) - r[i] += 2147483647; + // This does: + // r[i] = (16807 * r[i-1]) % 2147483647; + // buit avoids overflowing 31 bits. + const long int hi = word / 127773; + const long int lo = word % 127773; + word = 16807 * lo - 2836 * hi; + if (word < 0) + word += 2147483647; + r[i] = word; } k=31; for (int i=31; i<34; i++) @@ -88,8 +188,8 @@ namespace Testing return (unsigned int)ret >> 1; } -// reseed our random number generator - void srand(int seed) throw() + // reseed our random number generator + void srand(const int seed) { rand(true, seed); } @@ -97,6 +197,31 @@ namespace Testing +// Get a uniformly distributed random value between min and max +template +T random_value(const T &min=static_cast(0), + const T &max=static_cast(1)) +{ + return min+(max-min)*(static_cast(Testing::rand())/static_cast(RAND_MAX)); +} + + + +// Construct a uniformly distributed random point, with each coordinate +// between min and max +template +inline Point random_point(const double &min=0.0, + const double &max=1.0) +{ + Assert(max >= min, ExcMessage("Make sure max>=min")); + Point p; + for (unsigned int i=0; i +unsigned int checksum(const IT &begin, const IT &end) +{ + AssertThrow(sizeof(unsigned int)==4, ExcInternalError()); + AssertThrow(sizeof(*begin)==1, ExcInternalError()); + + unsigned int a = 1; + unsigned int b = 0; + + IT it = begin; + + while (it != end) + { + a = (a + (unsigned char)*it) % 65521; + b = (a + b) % 65521; + ++it; + } + + return (b << 16) | a; } @@ -139,14 +289,65 @@ void sort_file_contents (const std::string &filename) * Also, while GCC prepends the name by "virtual " if the function is virtual, * Intel's ICC does not do that, so filter that out as well. */ -void unify_pretty_function (const std::string &filename) +std::string unify_pretty_function (const std::string &text) { - int error = std::system ((std::string ("sed -i -e 's/ \\&/ \\& /g' -e 's/ & ,/\\&,/g' -e 's/ \\& )/\\&)/g' -e 's/ \\& /\\& /g' -e 's/^DEAL::virtual /DEAL::/g' ") + filename).c_str()); - - Assert (error == 0, ExcInternalError()); + std::string t=text; + t=Utilities::replace_in_string(t, " &", " & "); + t=Utilities::replace_in_string(t, " & ,", "&,"); + t=Utilities::replace_in_string(t, " & )", "&)"); + t=Utilities::replace_in_string(t, " & ", "& "); + t=Utilities::replace_in_string(t, "virtual ", ""); + return t; } +/* + * Test that a solver converged within a certain range of iteration steps. + * + * SolverType_COMMAND is the command to issue, CONTROL_COMMAND a function call + * that returns the number of iterations (castable to unsigned int), and + * MIN_ALLOWED, MAX_ALLOWED is the inclusive range of allowed iteration + * steps. + */ + +#define check_solver_within_range(SolverType_COMMAND, CONTROL_COMMAND, MIN_ALLOWED, MAX_ALLOWED) \ + { \ + const unsigned int previous_depth = deallog.depth_file(0); \ + try \ + { \ + SolverType_COMMAND; \ + } \ + catch (SolverControl::NoConvergence &exc) \ + {} \ + deallog.depth_file(previous_depth); \ + const unsigned int steps = CONTROL_COMMAND; \ + if (steps >= MIN_ALLOWED && steps <= MAX_ALLOWED) \ + { \ + deallog << "Solver stopped within " << MIN_ALLOWED << " - " \ + << MAX_ALLOWED << " iterations" << std::endl; \ + } \ + else \ + { \ + deallog << "Solver stopped after " << steps << " iterations" \ + << std::endl; \ + } \ + } + +/* + * Allow a test program to define a number that is very small to a given + * tolerance to be output as zero. This is used e.g. for the output of float + * numbers where roundoff difference can make the error larger than what we + * have set for numdiff (that is appropriate for double variables). + */ +template +Number filter_out_small_numbers (const Number number, const double tolerance) +{ + if (std::abs(number) < tolerance) + return Number(); + else + return number; +} + // ------------------------------ Functions used in initializing subsystems ------------------- @@ -155,21 +356,19 @@ void unify_pretty_function (const std::string &filename) * each of them runs 64 threads, then we get astronomical loads. * Limit concurrency to a fixed (small) number of threads, independent * of the core count. - * - * Note that we can't do this if we run in MPI mode because then - * MPI_InitFinalize already calls this function. Since every test - * calls MPI_InitFinalize itself, we can't adjust the thread count - * for this here. */ -#ifndef DEAL_II_WITH_MPI +inline unsigned int testing_max_num_threads() +{ + return 5; +} + struct LimitConcurrency { LimitConcurrency () { - multithread_info.set_thread_limit (5); + MultithreadInfo::set_thread_limit (testing_max_num_threads()); } } limit_concurrency; -#endif @@ -180,6 +379,7 @@ namespace { void check_petsc_allocations() { +#if DEAL_II_PETSC_VERSION_GTE(3, 2, 0) PetscStageLog stageLog; PetscLogGetStageLog(&stageLog); @@ -208,6 +408,7 @@ namespace if (errors) throw dealii::ExcMessage("PETSc memory leak"); +#endif } } #endif @@ -230,17 +431,13 @@ initlog(bool console=false) deallogname = "output"; deallogfile.open(deallogname.c_str()); deallog.attach(deallogfile); - if (!console) - deallog.depth_console(0); - -//TODO: Remove this line and replace by test_mode() - deallog.threshold_float(1.e-8); + deallog.depth_console(console?10:0); } inline void -mpi_initlog(bool console=false) +mpi_initlog(const bool console=false) { #ifdef DEAL_II_WITH_MPI unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); @@ -249,13 +446,10 @@ mpi_initlog(bool console=false) deallogname = "output"; deallogfile.open(deallogname.c_str()); deallog.attach(deallogfile); - if (!console) - deallog.depth_console(0); - -//TODO: Remove this line and replace by test_mode() - deallog.threshold_float(1.e-8); + deallog.depth_console(console?10:0); } #else + (void)console; // can't use this function if not using MPI Assert (false, ExcInternalError()); #endif @@ -263,26 +457,38 @@ mpi_initlog(bool console=false) -/* helper class to include the deallogs of all processors - on proc 0 */ +/** + * A helper class that gives each MPI process its own output file + * for the `deallog` stream, and at the end of the program (or, + * more correctly, the end of the current object), concatenates them + * all into the output file used on processor 0. + */ struct MPILogInitAll { - MPILogInitAll(bool console=false) + MPILogInitAll(const bool console=false) { #ifdef DEAL_II_WITH_MPI - unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - deallogname = "output"; - if (myid != 0) - deallogname = deallogname + Utilities::int_to_string(myid); - deallogfile.open(deallogname.c_str()); - deallog.attach(deallogfile); - if (!console) - deallog.depth_console(0); - -//TODO: Remove this line and replace by test_mode() - deallog.threshold_float(1.e-8); + const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + if (myid == 0) + { + if (!deallog.has_file()) + { + deallogfile.open("output"); + deallog.attach(deallogfile); + } + } + else + { + deallogname = "output" + Utilities::int_to_string(myid); + deallogfile.open(deallogname.c_str()); + deallog.attach(deallogfile); + } + + deallog.depth_console(console ? 10 : 0); + deallog.push(Utilities::int_to_string(myid)); #else + (void)console; // can't use this function if not using MPI Assert (false, ExcInternalError()); #endif @@ -291,9 +497,10 @@ struct MPILogInitAll ~MPILogInitAll() { #ifdef DEAL_II_WITH_MPI - unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); - unsigned int nproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + const unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD); + const unsigned int nproc = Utilities::MPI::n_mpi_processes (MPI_COMM_WORLD); + // pop the prefix for the MPI rank of the current process deallog.pop(); if (myid!=0) @@ -314,19 +521,11 @@ struct MPILogInitAll for (unsigned int i=1; i &v) +{ + for (unsigned int i=0; i &v) +{ + for (unsigned int i=0; i &v) +{ + for (unsigned int i=0; i