From: Martin Kronbichler Date: Tue, 30 Apr 2013 14:56:37 +0000 (+0000) Subject: Make partial initialization work properly, including a minimum functionality in case... X-Git-Tag: v8.0.0~602 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5a9de4f1449e9bed5db322ba9ca33ab5d971c9fa;p=dealii.git Make partial initialization work properly, including a minimum functionality in case of no index initializations. git-svn-id: https://svn.dealii.org/trunk@29408 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/matrix_free/dof_info.h b/deal.II/include/deal.II/matrix_free/dof_info.h index 90c7fccdab..b7a791dd2b 100644 --- a/deal.II/include/deal.II/matrix_free/dof_info.h +++ b/deal.II/include/deal.II/matrix_free/dof_info.h @@ -354,6 +354,13 @@ namespace internal */ std::vector plain_dof_indices; + /** + * Stores the dimension of the underlying DoFHandler. Since the indices + * are not templated, this is the variable that makes the dimension + * accessible in the (rare) cases it is needed inside this class. + */ + unsigned int dimension; + /** * Stores the number of components in the DoFHandler where the indices * have been read from. diff --git a/deal.II/include/deal.II/matrix_free/dof_info.templates.h b/deal.II/include/deal.II/matrix_free/dof_info.templates.h index 2430dcecdb..64623405d5 100644 --- a/deal.II/include/deal.II/matrix_free/dof_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/dof_info.templates.h @@ -131,6 +131,7 @@ namespace internal constrained_dofs (dof_info_in.constrained_dofs), row_starts_plain_indices (dof_info_in.row_starts_plain_indices), plain_dof_indices (dof_info_in.plain_dof_indices), + dimension (dof_info_in.dimension), n_components (dof_info_in.n_components), dofs_per_cell (dof_info_in.dofs_per_cell), dofs_per_face (dof_info_in.dofs_per_face), @@ -153,6 +154,7 @@ namespace internal ghost_dofs.clear(); dofs_per_cell.clear(); dofs_per_face.clear(); + dimension = 2; n_components = 0; row_starts_plain_indices.clear(); plain_dof_indices.clear(); diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h index 9acb29a336..187716af64 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -687,22 +687,6 @@ protected: */ unsigned int cell_data_number; - /** - * If the present cell chunk for vectorization is not completely filled up - * with data, this field stores how many physical cells are underlying. Is - * between 1 and VectorizedArray::n_array_elements-1 (inclusive). - */ - unsigned int n_irreg_components_filled; - - /** - * Stores whether the present cell chunk used in vectorization is not - * completely filled up with physical cells. E.g. if vectorization dictates - * that four cells should be worked with but only three physical cells are - * left, this flag will be set to true, otherwise to false. Mainly used for - * internal checking when reading from vectors or writing to vectors. - */ - bool at_irregular_cell; - /** * Debug information to track whether dof values have been initialized * before accessed. Used to control exceptions when uninitialized data is @@ -1465,17 +1449,14 @@ FEEvaluationBase jacobian_grad_upper(0), cell (numbers::invalid_unsigned_int), cell_type (internal::MatrixFreeFunctions::undefined), - cell_data_number (0), - n_irreg_components_filled (0), - at_irregular_cell (false) + cell_data_number (0) { - Assert (matrix_info.indices_initialized() == true, - ExcNotInitialized()); Assert (matrix_info.mapping_initialized() == true, ExcNotInitialized()); AssertDimension (matrix_info.get_size_info().vectorization_length, VectorizedArray::n_array_elements); Assert (n_fe_components == 1 || + n_components == 1 || n_components == n_fe_components, ExcMessage ("The underlying FE is vector-valued. In this case, the " "template argument n_components must be a the same " @@ -1553,9 +1534,6 @@ FEEvaluationBase } } - n_irreg_components_filled = - std_cxx1x::get<2>(dof_info.row_starts[cell_in]); - at_irregular_cell = n_irreg_components_filled > 0; #ifdef DEBUG dof_values_initialized = false; values_quad_initialized = false; @@ -1832,6 +1810,8 @@ FEEvaluationBase // into the local data field or write local data into the vector. Certain // operations are no-ops for the given use case. + Assert (matrix_info.indices_initialized() == true, + ExcNotInitialized()); Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); // loop over all local dofs. ind_local holds local number on cell, index @@ -1844,6 +1824,10 @@ FEEvaluationBase dof_info.end_indicators(cell); unsigned int ind_local = 0; + const unsigned int n_irreg_components_filled = + std_cxx1x::get<2>(dof_info.row_starts[cell]); + const bool at_irregular_cell = n_irreg_components_filled > 0; + // scalar case (or case when all components have the same degrees of freedom // and sit on a different vector each) if (n_fe_components == 1) @@ -2447,6 +2431,8 @@ FEEvaluationBase { // this is different from the other three operations because we do not use // constraints here, so this is a separate function. + Assert (matrix_info.indices_initialized() == true, + ExcNotInitialized()); Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized()); Assert (dof_info.store_plain_indices == true, ExcNotInitialized()); @@ -2455,6 +2441,10 @@ FEEvaluationBase // points to the global indices stored in index_local_to_global const unsigned int *dof_indices = dof_info.begin_indices_plain(cell); + const unsigned int n_irreg_components_filled = + std_cxx1x::get<2>(dof_info.row_starts[cell]); + const bool at_irregular_cell = n_irreg_components_filled > 0; + // scalar case (or case when all components have the same degrees of freedom // and sit on a different vector each) if (n_fe_components == 1) diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h index ccb397894f..2014ab8435 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h @@ -159,6 +159,32 @@ internal_reinit(const Mapping &mapping, initialize_indices (constraint, locally_owned_set); } + // initialize bare structures + else if (dof_info.size() != dof_handler.size()) + { + initialize_dof_handlers(dof_handler, additional_data.level_mg_handler); + std::vector dummy; + size_info.make_layout (cell_level_index.size(), + VectorizedArray::n_array_elements, + dummy, dummy); + for (unsigned int i=0; iget_fe().element_multiplicity(0); + dof_info[i].dofs_per_cell.push_back(dof_handler[i]->get_fe().dofs_per_cell); + dof_info[i].row_starts.resize(size_info.n_macro_cells+1); + std_cxx1x::get<2>(dof_info[i].row_starts.back()) = + cell_level_index.size() % VectorizedArray::n_array_elements; + + // if indices are not initialized, the cell_level_index might not be + // divisible by the vectorization length. But it must be for + // mapping_info... + while (cell_level_index.size() % VectorizedArray::n_array_elements + != 0) + cell_level_index.push_back(cell_level_index.back()); + } + } + // Reads out the FE information and stores the shape function values, // gradients and Hessians for quadrature points. const unsigned int n_fe = dof_handler.size(); @@ -254,6 +280,33 @@ internal_reinit(const Mapping &mapping, initialize_indices (constraint, locally_owned_set); } + // initialize bare structures + else if (dof_info.size() != dof_handler.size()) + { + initialize_dof_handlers(dof_handler, additional_data.level_mg_handler); + std::vector dummy; + size_info.make_layout (cell_level_index.size(), + VectorizedArray::n_array_elements, + dummy, dummy); + for (unsigned int i=0; iget_fe().size() == 1, ExcNotImplemented()); + dof_info[i].dimension = dim; + dof_info[i].n_components = dof_handler[i]->get_fe()[0].element_multiplicity(0); + dof_info[i].dofs_per_cell.push_back(dof_handler[i]->get_fe()[0].dofs_per_cell); + dof_info[i].row_starts.resize(size_info.n_macro_cells+1); + std_cxx1x::get<2>(dof_info[i].row_starts.back()) = + cell_level_index.size() % VectorizedArray::n_array_elements; + + // if indices are not initialized, the cell_level_index might not be + // divisible by the vectorization length. But it must be for + // mapping_info... + while (cell_level_index.size() % VectorizedArray::n_array_elements + != 0) + cell_level_index.push_back(cell_level_index.back()); + } + } + // Reads out the FE information and stores the shape function values, // gradients and Hessians for quadrature points. const unsigned int n_components = dof_handler.size(); @@ -460,8 +513,8 @@ void MatrixFree::initialize_indices // cache number of finite elements and dofs_per_cell dof_info[no].dofs_per_cell.push_back (fe.dofs_per_cell); dof_info[no].dofs_per_face.push_back (fe.dofs_per_face); - dof_info[no].n_components = n_fe_components; - + dof_info[no].dimension = dim; + dof_info[no].n_components = n_fe_components; // get permutation that gives lexicographic renumbering of the cell // dofs renumber (this is necessary for FE_Q, for example, since diff --git a/tests/matrix_free/get_functions_common.h b/tests/matrix_free/get_functions_common.h index b859b07386..f3e66e8c25 100644 --- a/tests/matrix_free/get_functions_common.h +++ b/tests/matrix_free/get_functions_common.h @@ -59,7 +59,7 @@ class MatrixFreeTest virtual ~MatrixFreeTest () {} - + // make function virtual to allow derived // classes to define a different function virtual void @@ -269,4 +269,3 @@ int main () deallog.pop(); } } - diff --git a/tests/matrix_free/no_index_initialize.cc b/tests/matrix_free/no_index_initialize.cc new file mode 100644 index 0000000000..f3736b8894 --- /dev/null +++ b/tests/matrix_free/no_index_initialize.cc @@ -0,0 +1,183 @@ +//------------------ no_index_initialize.cc ------------------------ +// $Id$ +// Version: $Name$ +// +//------------------ no_index_initialize.cc ------------------------ + + +// check that FEEvaluation can be evaluated without indices initialized (and +// throws an exception when trying to read/write from/to vectors) + +#include "../tests.h" + + +std::ofstream logfile("no_index_initialize/output"); + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +// forward declare this function +template +void test (); + + + +template +class MatrixFreeTest +{ + public: + bool read_vector; + + MatrixFreeTest(const MatrixFree &data_in): + read_vector(false), + data (data_in) + {}; + + void + operator () (const MatrixFree &data, + Vector &, + const Vector &src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval (data); + + for(unsigned int cell=cell_range.first;cell(1.), e); + fe_eval.evaluate (true,true,true); + + // values should evaluate to one, derivatives to zero + for (unsigned int q=0; q::n_array_elements; ++d) + { + Assert(fe_eval.get_value(q)[d] == 1., ExcInternalError()); + for (unsigned int e=0; e &src) const + { + Vector dst_dummy; + data.cell_loop (&MatrixFreeTest::operator(), this, dst_dummy, src); + deallog << "OK" << std::endl; + }; + +protected: + const MatrixFree &data; +}; + + + + +template +void do_test (const DoFHandler &dof, + const ConstraintMatrix&constraints) +{ + // use this for info on problem + //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() + // << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + Vector solution (dof.n_dofs()); + + // create vector with random entries + for (unsigned int i=0; i mf_data; + { + const QGauss<1> quad (fe_degree+1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_gradients | update_second_derivatives; + data.initialize_indices = false; + mf_data.reinit (dof, constraints, quad, data); + } + + deallog << "Testing " << dof.get_fe().get_name() << " without read" << std::endl; + MatrixFreeTest mf (mf_data); + mf.test_functions(solution); + + deallog << "Testing " << dof.get_fe().get_name() << " with read" << std::endl; + try + { + mf.read_vector = true; + mf.test_functions(solution); + } + catch (ExceptionBase &e) + { + deallog << e.get_exc_name() << std::endl; + } +} + + +int main () +{ + deal_II_exceptions::disable_abort_on_exception(); + deallog.attach(logfile); + deallog.depth_console(0); + + deallog << std::setprecision (3); + test<2,1>(); +} + + +template +void test () +{ + Triangulation tria; + GridGenerator::hyper_ball (tria); + static const HyperBallBoundary boundary; + tria.set_boundary (0, boundary); + // refine first and last cell + tria.begin(tria.n_levels()-1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.refine_global (4-dim); + + FE_Q fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + constraints.close(); + + do_test (dof, constraints); +} diff --git a/tests/matrix_free/no_index_initialize/cmp/generic b/tests/matrix_free/no_index_initialize/cmp/generic new file mode 100644 index 0000000000..1a9723e8a1 --- /dev/null +++ b/tests/matrix_free/no_index_initialize/cmp/generic @@ -0,0 +1,5 @@ + +DEAL::Testing FE_Q<2>(1) without read +DEAL::OK +DEAL::Testing FE_Q<2>(1) with read +DEAL::ExcNotInitialized() diff --git a/tests/matrix_free/update_mapping_only.cc b/tests/matrix_free/update_mapping_only.cc new file mode 100644 index 0000000000..d0a619fc27 --- /dev/null +++ b/tests/matrix_free/update_mapping_only.cc @@ -0,0 +1,84 @@ +//------------------ update_mapping_only.cc ------------------------ +// $Id$ +// Version: $Name$ +// +//------------------ update_mapping_only.cc ------------------------ + + +// correctness of matrix-free initialization when the mapping changes and the +// indices are not updated, only the mapping info + +#include "../tests.h" + +#include +#include + + +std::ofstream logfile("update_mapping_only/output"); + +#include "get_functions_common.h" + + +template +void test () +{ + Triangulation tria; + GridGenerator::hyper_cube (tria); + tria.refine_global(1); + + FE_Q fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + + ConstraintMatrix constraints; + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + // use this for info on problem + //std::cout << "Number of cells: " << dof.get_tria().n_active_cells() + // << std::endl; + //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl; + //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl; + + Vector solution (dof.n_dofs()); + + // create vector with random entries + for (unsigned int i=0; i fe_sys(dof.get_fe(), dim); + DoFHandler dofh_eulerian(dof.get_tria()); + dofh_eulerian.distribute_dofs(fe_sys); + + MatrixFree mf_data; + Vector shift(dofh_eulerian.n_dofs()); + for (unsigned int i=0; i<2; ++i) + { + if (i == 1) + { + shift(0) = 0.121; + shift(1) = -0.005; + } + MappingQEulerian mapping(2,shift,dofh_eulerian); + + { + const QGauss<1> quad (fe_degree+1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_gradients | update_second_derivatives; + if (i==1) + data.initialize_indices = false; + mf_data.reinit (mapping, dof, constraints, quad, data); + } + + MatrixFreeTest mf (mf_data, mapping); + mf.test_functions(solution); + } +} diff --git a/tests/matrix_free/update_mapping_only/cmp/generic b/tests/matrix_free/update_mapping_only/cmp/generic new file mode 100644 index 0000000000..8586990170 --- /dev/null +++ b/tests/matrix_free/update_mapping_only/cmp/generic @@ -0,0 +1,79 @@ + +DEAL:2d::Testing FE_Q<2>(1) +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(2) +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(3) +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Testing FE_Q<2>(4) +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:2d::Error function values: 0 +DEAL:2d::Error function gradients: 0 +DEAL:2d::Error function Laplacians: 0 +DEAL:2d::Error function diagonal of Hessian: 0 +DEAL:2d::Error function Hessians: 0 +DEAL:2d:: +DEAL:3d::Testing FE_Q<3>(1) +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d:: +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d:: +DEAL:3d::Testing FE_Q<3>(2) +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d:: +DEAL:3d::Error function values: 0 +DEAL:3d::Error function gradients: 0 +DEAL:3d::Error function Laplacians: 0 +DEAL:3d::Error function diagonal of Hessian: 0 +DEAL:3d::Error function Hessians: 0 +DEAL:3d::