From: Martin Kronbichler Date: Wed, 9 Dec 2009 18:11:22 +0000 (+0000) Subject: Correct counting of dofs per component. Improve speed of counting dofs per component... X-Git-Tag: v8.0.0~6752 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=241d7332245ac937e666038dc03369451d017db8;p=dealii.git Correct counting of dofs per component. Improve speed of counting dofs per component. Correct bug in extract_constant_modes. git-svn-id: https://svn.dealii.org/trunk@20219 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index 82919b65be..4752f92904 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -1106,55 +1106,44 @@ class DoFTools std::vector &selected_dofs); /** - * Extract a vector that - * represents the constant - * modes of the DoFHandler - * for the components - * chosen by component_select. - * The constant modes on a - * discretization are the - * null space of a Laplace - * operator on the selected - * components with Neumann - * boundary conditions - * applied. The null space is - * a necessary ingredient for - * obtaining a good AMG - * preconditioner when using - * the class + * Extract a vector that represents the + * constant modes of the DoFHandler for + * the components chosen by + * component_select. The + * constant modes on a discretization are + * the null space of a Laplace operator + * on the selected components with + * Neumann boundary conditions + * applied. The null space is a necessary + * ingredient for obtaining a good AMG + * preconditioner when using the class * TrilinosWrappers::PreconditionAMG. - * Since the ML AMG package - * only works on algebraic - * properties of the respective - * matrix, it has no chance - * to detect whether the matrix - * comes from a scalar or a - * vector valued problem. However, - * a near null space supplies - * exactly the needed information - * about these components. - * The null space will consist - * of as many vectors as there + * Since the ML AMG package only works on + * algebraic properties of the respective + * matrix, it has no chance to detect + * whether the matrix comes from a scalar + * or a vector valued problem. However, a + * near null space supplies exactly the + * needed information about these + * components. The null space will + * consist of as many vectors as there * are true arguments in - * component_select, - * each of which will be one - * in one component and - * zero in all others. We store - * this object in a vector of - * vectors, where the outer - * vector is of the size of - * the number of selected - * components, and each inner - * vector has as many components - * as there are degrees of + * component_select, each of + * which will be one in one component and + * zero in all others. We store this + * object in a vector of vectors, where + * the outer vector is of the size of the + * number of selected components, and + * each inner vector has as many + * components as there are degrees of * freedom in the selected - * components. Note that any - * matrix associated with this - * null space must have been - * constructed using the - * same - * component_select - * argument. + * components. Note that any matrix + * associated with this null space must + * have been constructed using the same + * component_select argument, + * since the numbering of DoFs is done + * relative to the selected dofs, not to + * all dofs. * * The main reason for this * program is the use of the diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 9d0ad85840..645de5e5dc 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -2987,6 +2987,116 @@ void DoFTools::distribute_cell_to_dof_vector ( +namespace internal +{ + // this internal function assigns to each dof + // the respective component of the vector + // system. if use_blocks is set, then the + // assignment is done by blocks, not by + // components, as specified by + // component_select. The additional argument + // component_select only is used for + // non-primitive FEs, where we need it since + // more components couple, and no unique + // component can be assigned. Then, we sort + // them to the first selected component of the + // vector system. + template + inline + void + extract_dofs_by_component (const dealii::DoFHandler &dof, + const std::vector &component_select, + const bool sort_by_blocks, + std::vector &dofs_by_component) + { + const FiniteElement &fe = dof.get_fe(); + Assert (fe.n_components() < 256, ExcNotImplemented()); + Assert(dofs_by_component.size() == dof.n_dofs(), + ExcDimensionMismatch(dofs_by_component.size(), dof.n_dofs())); + + + // next set up a table for the + // degrees of freedom on each of + // the cells whether it is + // something interesting or not + std::vector local_component_association (fe.dofs_per_cell); + if (sort_by_blocks == true) + { + for (unsigned int i=0; i indices(fe.dofs_per_cell); + typename dealii::DoFHandler::active_cell_iterator c; + for (c=dof.begin_active(); c!=dof.end(); ++ c) + { + c->get_dof_indices(indices); + for (unsigned int i=0; i + inline + void + extract_dofs_by_component (const dealii::hp::DoFHandler &, + const std::vector &, + const bool , + std::vector&) + { + Assert (false, ExcNotImplemented()); + } +} + + + template void DoFTools::extract_dofs ( @@ -3032,71 +3142,20 @@ DoFTools::extract_dofs ( // preset all values by false std::fill_n (selected_dofs.begin(), dof.n_dofs(), false); - // next set up a table for the - // degrees of freedom on each of - // the cells whether it is - // something interesting or not - std::vector local_selected_dofs (fe.dofs_per_cell, false); - for (unsigned int i=0; i indices(fe.dofs_per_cell); - typename DoFHandler::active_cell_iterator c; - for (c=dof.begin_active(); c!=dof.end(); ++ c) - { - c->get_dof_indices(indices); - for (unsigned int i=0; i dofs_by_component (dof.n_dofs()); + internal::extract_dofs_by_component (dof, component_select, count_by_blocks, + dofs_by_component); + + for (unsigned int i=0; i void -DoFTools::extract_boundary_dofs (const DH &dof_handler, - const std::vector &component_select, - std::vector &selected_dofs, +DoFTools::extract_boundary_dofs (const DH &dof_handler, + const std::vector &component_select, + std::vector &selected_dofs, const std::set &boundary_indicators) { Assert (component_select.size() == n_components(dof_handler), @@ -3756,48 +3815,36 @@ DoFTools::extract_constant_modes (const DH &dof_handler, Assert (n_components == component_select.size(), ExcDimensionMismatch(n_components, component_select.size())); + std::vector localized_component (n_components, + numbers::invalid_unsigned_int); + unsigned int n_components_selected = 0; + for (unsigned int i=0; i dofs_by_component (dof_handler.n_dofs()); + internal::extract_dofs_by_component (dof_handler, component_select, false, + dofs_by_component); + unsigned int n_selected_dofs = 0; + for (unsigned int i=0; i(n_selected_dofs, + false)); std::vector component_list (n_components, 0); for (unsigned int d=0; d dofs_per_block(2); - count_dofs_per_block(dof_handler, dofs_per_block, component_list); - - const unsigned int n_u = dofs_per_block[1]; - std::vector selection_dof_list (dof_handler.n_dofs(), false); - std::vector temporary_dof_list (dof_handler.n_dofs(), false); - extract_dofs (dof_handler, component_select, selection_dof_list); + component_list[d] = component_select[d]; - constant_modes.resize (n_components_selected, std::vector(n_u, false)); - - for (unsigned int component=0, component_used=0; - component < n_components; ++component, ++component_used) - if (component_select[component]) + unsigned int counter = 0; + for (unsigned int i=0; i selection_mask (n_components, false); - selection_mask[component] = true; - extract_dofs (dof_handler, selection_mask, temporary_dof_list); - - unsigned int counter = 0; - for (unsigned int i=0; i subdomain_association (dof_handler.n_dofs()); get_subdomain_association (dof_handler, subdomain_association); - std::vector component_association (dof_handler.n_dofs()); + std::vector component_association (dof_handler.n_dofs()); + internal::extract_dofs_by_component (dof_handler, std::vector(), false, + component_association); + for (unsigned int c=0; c component_mask (dof_handler.get_fe().n_components(), false); - component_mask[c] = true; - extract_dofs (dof_handler, component_mask, component_association); - for (unsigned int i=0; i(c))) ++n_dofs_on_subdomain[c]; } } +namespace internal +{ + template + void + resolve_components (const FiniteElement&fe, + const std::vector &dofs_by_component, + const std::vector &target_component, + const bool only_once, + std::vector &dofs_per_component, + unsigned int &component) + { + for (unsigned int b=0;b& base = fe.base_element(b); + // Dimension of base element + unsigned int d = base.n_components(); + + for (unsigned int m=0;m 1) + resolve_components(base, dofs_by_component, target_component, + only_once, dofs_per_component, component); + else + { + for (unsigned int dd=0;dd void DoFTools::count_dofs_per_component ( @@ -4018,7 +4108,6 @@ DoFTools::count_dofs_per_component ( fe.n_components())); - const unsigned int max_component = *std::max_element (target_component.begin(), target_component.end()); @@ -4041,48 +4130,15 @@ DoFTools::count_dofs_per_component ( // otherwise determine the number // of dofs in each component // separately. do so in parallel - std::vector > - dofs_in_component (n_components, - std::vector(dof_handler.n_dofs(), false)); - std::vector > - component_select (n_components, - std::vector(n_components, false)); - Threads::TaskGroup<> tasks; - for (unsigned int i=0; i &, - const std::vector &, - std::vector &, - bool) - = &DoFTools::template extract_dofs; - component_select[i][i] = true; - tasks += Threads::new_task (fun_ptr, - dof_handler, component_select[i], - dofs_in_component[i], false); - } - tasks.join_all (); + std::vector dofs_by_component (dof_handler.n_dofs()); + internal::extract_dofs_by_component (dof_handler, std::vector(), false, + dofs_by_component); // next count what we got unsigned int component = 0; - for (unsigned int b=0;b& base = fe.base_element(b); - // Dimension of base element - unsigned int d = base.n_components(); - - for (unsigned int m=0;m void DoFTools::count_dofs_per_block ( - const DoFHandler& dof_handler, + const DoFHandler& dof_handler, std::vector& dofs_per_block, std::vector target_block) { @@ -4144,31 +4200,15 @@ DoFTools::count_dofs_per_block ( // otherwise determine the number // of dofs in each block // separately. do so in parallel - std::vector > - dofs_in_block (n_blocks, std::vector(dof_handler.n_dofs(), false)); - std::vector > - block_select (n_blocks, std::vector(n_blocks, false)); - Threads::TaskGroup<> tasks; - for (unsigned int i=0; i &, - const std::vector &, - std::vector &, - bool) - = &DoFTools::template extract_dofs; - block_select[i][i] = true; - tasks += Threads::new_task (fun_ptr, - dof_handler, block_select[i], - dofs_in_block[i], true); - }; - tasks.join_all (); + std::vector dofs_by_block (dof_handler.n_dofs()); + internal::extract_dofs_by_component (dof_handler, std::vector(), true, + dofs_by_block); // next count what we got for (unsigned int block=0;block