From: Wolfgang Bangerth Date: Mon, 8 Jul 2013 14:01:00 +0000 (+0000) Subject: Minimal form of the test. This fails with my patches to work_stream.h because we... X-Git-Tag: v8.0.0~167 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ed371a8332418e5a4b01edadbe3e4d9932a1ecbd;p=dealii.git Minimal form of the test. This fails with my patches to work_stream.h because we step on some other array's memory... git-svn-id: https://svn.dealii.org/trunk@29953 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/base/work_stream_03.cc b/tests/base/work_stream_03.cc index 61c094d4ef..cea42e56fd 100644 --- a/tests/base/work_stream_03.cc +++ b/tests/base/work_stream_03.cc @@ -20,7 +20,6 @@ // verifying that it indeed works #include "../tests.h" -#include #include #include #include @@ -28,29 +27,12 @@ #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 @@ -60,39 +42,29 @@ char logname[] = "work_stream_03/output"; template -class F : public Function +double value (const Point &p) { - public: - virtual double value (const Point &p, - const unsigned int component) const - { - Assert ((component == 0) && (this->n_components == 1), - ExcInternalError()); - double val = 0; - for (unsigned int d=0; d + template struct Scratch { - Scratch (const ::dealii::hp::FECollection &fe, - const UpdateFlags update_flags, - const ::dealii::hp::QCollection &quadrature) + Scratch (const FiniteElement &fe, + const Quadrature &quadrature) : fe_collection (fe), quadrature_collection (quadrature), x_fe_values (fe_collection, quadrature_collection, - update_flags), - rhs_values(quadrature_collection.max_n_quadrature_points()), - update_flags (update_flags) + update_q_points), + rhs_values(quadrature_collection.size()) {} Scratch (const Scratch &data) @@ -101,19 +73,16 @@ namespace quadrature_collection (data.quadrature_collection), x_fe_values (fe_collection, quadrature_collection, - data.update_flags), - rhs_values (data.rhs_values), - update_flags (data.update_flags) + update_q_points), + rhs_values (data.rhs_values) {} - const ::dealii::hp::FECollection &fe_collection; - const ::dealii::hp::QCollection &quadrature_collection; + const FiniteElement &fe_collection; + const Quadrature &quadrature_collection; - ::dealii::hp::FEValues x_fe_values; + FEValues x_fe_values; std::vector rhs_values; - - const UpdateFlags update_flags; }; @@ -123,108 +92,67 @@ namespace }; } -template +template void - mass_assembler(const CellIterator &cell, Scratch &data, + mass_assembler(const CellIterator &cell, Scratch &data, CopyData ©_data) { data.x_fe_values.reinit(cell); - const FEValues &fe_values = - data.x_fe_values.get_present_fe_values(); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = - fe_values.n_quadrature_points; - const FiniteElement &fe = fe_values.get_fe(); - const unsigned int n_components = fe.n_components(); - - { - copy_data.cell_rhs.reinit(dofs_per_cell); - } - - data.rhs_values.resize(n_q_points); - F().value_list(fe_values.get_quadrature_points(), - data.rhs_values); - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - const unsigned int component_i = fe.system_to_component_index(i).first; - const double *phi_i = &fe_values.shape_value(i, 0); - - double add_data = 0; - for (unsigned int point = 0; point < n_q_points; ++point) - add_data += phi_i[point] * fe_values.JxW(point) - * data.rhs_values[point]; - copy_data.cell_rhs(i) = add_data; - } + + // this appears to be the key: the following line overwrites some of the memory + // in which we store the quadrature point location. if the line is moved down, + // the comparison in the if() always succeeds... + copy_data.cell_rhs = 0; + if (cell->center().distance (data.x_fe_values.quadrature_point(0)) >= 1e-6 *cell->diameter()) + std::cout << '.' << std::flush; + else + std::cout << '*' << std::flush; + + copy_data.cell_rhs[0] = value(data.x_fe_values.quadrature_point(0)); } -template void copy_local_to_global (const CopyData &data, - VectorType *right_hand_side) + double *sum) { -// std::cout << data.cell_rhs.l2_norm() << ' '; - right_hand_side->push_back (data.cell_rhs.l2_norm()); + *sum += data.cell_rhs[0]; } - -template -void create_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - std::vector &rhs_vector) +void do_project () { - hp::FECollection fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - Scratch - assembler_data (fe_collection, - update_values | - update_JxW_values | update_quadrature_points, - q_collection); - CopyData copy_data; - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - - dealii::WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &mass_assembler::active_cell_iterator>, - std_cxx1x::bind(& - copy_local_to_global >, - std_cxx1x::_1, &rhs_vector), - assembler_data, - copy_data, - 2*multithread_info.n_default_threads, - 1); -} + static const int dim = 3; - -template -void do_project (const FiniteElement &fe) -{ Triangulation triangulation; GridGenerator::hyper_cube (triangulation); triangulation.refine_global (2); + FE_Nothing<3> fe; DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs (fe); - std::vector tmp; - - create_mass_matrix (dof_handler, QGauss(3), - tmp); - - double sum=0; - for (unsigned int i=0; i q; -template -void test_no_hanging_nodes (const FiniteElement &fe) -{ for (unsigned int i=0; i<12; ++i) - do_project (fe); + { + std::vector tmp; + + double sum = 0; + Scratch assembler_data (dof_handler.get_fe(), q); + CopyData copy_data; + copy_data.cell_rhs.reinit (8); + dealii::WorkStream::run (dof_handler.begin_active(), + dof_handler.end(), + &mass_assembler::active_cell_iterator>, + std_cxx1x::bind(& + copy_local_to_global, + std_cxx1x::_1, &sum), + assembler_data, + copy_data, + 8, + 1); + printf ("\nCheck: %5.3f\n", sum); + } } @@ -237,5 +165,5 @@ int main () deallog.depth_console(0); deallog.threshold_double(1.e-10); - test_no_hanging_nodes (FE_Q<3>(1)); + do_project (); }