From: bangerth Date: Mon, 8 Oct 2012 18:36:08 +0000 (+0000) Subject: Merge from mainline. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2cdb83c17f7f5f57d72c35e48713b945c6392fb3;p=dealii-svn.git Merge from mainline. git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@27015 0785d39b-7218-0410-832d-ea1e28bc413d --- 2cdb83c17f7f5f57d72c35e48713b945c6392fb3 diff --cc deal.II/include/deal.II/base/index_set.h index 46710d65b3,2baa112708..b1de672aed --- a/deal.II/include/deal.II/base/index_set.h +++ b/deal.II/include/deal.II/base/index_set.h @@@ -64,7 -64,7 +64,7 @@@ class IndexSe * overall size of the index * range. */ - IndexSet (const types::global_dof_index size); - explicit IndexSet (const unsigned int size); ++ explicit IndexSet (const types::global_dof_index size); /** * Remove all indices from this diff --cc deal.II/include/deal.II/dofs/dof_renumbering.h index 294c13ea88,341d1da45e..16328740bb --- a/deal.II/include/deal.II/dofs/dof_renumbering.h +++ b/deal.II/include/deal.II/dofs/dof_renumbering.h @@@ -667,8 -688,13 +688,13 @@@ namespace DoFRenumberin Cuthill_McKee (MGDoFHandler &dof_handler, const unsigned int level, const bool reversed_numbering = false, - const std::vector &starting_indices = std::vector ()); + const std::vector &starting_indices = std::vector ()); + /** + * @name Component-wise numberings + * @{ + */ + /** * Sort the degrees of freedom by * vector component. The @@@ -1245,9 -1398,18 +1398,18 @@@ */ template void - compute_random (std::vector &new_dof_indices, + compute_random (std::vector &new_dof_indices, const DH& dof_handler); + /** + * @} + */ + + /** + * @name Numberings based on cell attributes + * @{ + */ + /** * Renumber the degrees of * freedom such that they are @@@ -1298,9 -1460,13 +1460,13 @@@ */ template void - compute_subdomain_wise (std::vector &new_dof_indices, + compute_subdomain_wise (std::vector &new_dof_indices, const DH &dof_handler); + /** + * @} + */ + /** * Exception * diff --cc deal.II/include/deal.II/dofs/dof_tools.h index b6f7b2a531,3565e81af0..b3db9a13d6 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@@ -737,8 -738,8 +738,8 @@@ namespace DoFTool template void make_boundary_sparsity_pattern (const DH &dof, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &sparsity_pattern); - const std::vector &dof_to_boundary_mapping, - SparsityPattern &sparsity_pattern); ++ const std::vector &dof_to_boundary_mapping, ++ SparsityPattern &sparsity_pattern); /** * Write the sparsity structure of the @@@ -773,9 -774,9 +774,9 @@@ template void make_boundary_sparsity_pattern (const DH &dof, - const typename FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &sparsity); - const typename FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &sparsity); ++ const typename FunctionMap::type &boundary_indicators, ++ const std::vector &dof_to_boundary_mapping, ++ SparsityPattern &sparsity); /** * Generate sparsity pattern for @@@ -1234,22 -1250,61 +1250,61 @@@ template void extract_dofs (const hp::DoFHandler &dof_handler, - const std::vector &select, - std::vector &selected_dofs, - const bool blocks = false); + const ComponentMask &component_mask, + std::vector &selected_dofs); - /** - * Do the same thing as - * extract_dofs() for one level - * of a multi-grid DoF numbering. - */ + /** + * This function is the equivalent to the DoFTools::extract_dofs() functions above + * except that the selection of which degrees of freedom to extract is not done + * based on components (see @ref GlossComponent) but instead based on whether they + * are part of a particular block (see @ref GlossBlock). Consequently, the second + * argument is not a ComponentMask but a BlockMask object. + * + * @param dof_handler The DoFHandler object from which to extract degrees of freedom + * @param block_mask The block mask that describes which blocks to consider (see + * @ref GlossBlockMask) + * @param selected_dofs A vector of length DoFHandler::n_dofs() in which those + * entries are true that correspond to the selected blocks. + */ + template + void + extract_dofs (const DoFHandler &dof_handler, + const BlockMask &block_mask, + std::vector &selected_dofs); + + /** + * The same function as above, + * but for a hp::DoFHandler. + */ + template + void + extract_dofs (const hp::DoFHandler &dof_handler, + const BlockMask &block_mask, + std::vector &selected_dofs); + + /** + * Do the same thing as the corresponding + * extract_dofs() functino for one level + * of a multi-grid DoF numbering. + */ - template ++ template + void + extract_level_dofs (const unsigned int level, - const MGDoFHandler &dof, ++ const DH &dof, + const ComponentMask &component_mask, + std::vector &selected_dofs); + + /** + * Do the same thing as the corresponding + * extract_dofs() functino for one level + * of a multi-grid DoF numbering. + */ - template + template void extract_level_dofs (const unsigned int level, - const DH &dof, - const std::vector &select, - std::vector &selected_dofs, - const bool blocks = false); - const MGDoFHandler &dof, ++ const DH &dof, + const BlockMask &component_mask, + std::vector &selected_dofs); /** * Extract all degrees of freedom @@@ -1748,7 -1767,8 +1767,7 @@@ template IndexSet dof_indices_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id subdomain); + const types::subdomain_id subdomain); - // @} /** * @name Dof indices for patches @@@ -2095,10 -2115,10 +2114,10 @@@ template void count_dofs_per_component (const DH & dof_handler, - std::vector& dofs_per_component, - const bool vector_valued_once = false, - std::vector target_component - = std::vector()); - std::vector& dofs_per_component, - const bool vector_valued_once = false, - std::vector target_component - = std::vector()); ++ std::vector& dofs_per_component, ++ const bool vector_valued_once = false, ++ std::vector target_component ++ = std::vector()); /** * Count the degrees of freedom @@@ -2149,8 -2169,8 +2168,8 @@@ template void count_dofs_per_component (const DoFHandler& dof_handler, - std::vector& dofs_per_component, - std::vector target_component); - std::vector& dofs_per_component, - std::vector target_component); ++ std::vector& dofs_per_component, ++ std::vector target_component); /** * This function can be used when @@@ -2438,7 -2458,7 +2457,7 @@@ template void map_dof_to_boundary_indices (const DH &dof_handler, - std::vector &mapping); - std::vector &mapping); ++ std::vector &mapping); /** * Same as the previous function, @@@ -2454,8 -2474,8 +2473,8 @@@ template void map_dof_to_boundary_indices (const DH &dof_handler, - const std::set &boundary_indicators, - std::vector &mapping); - const std::set &boundary_indicators, - std::vector &mapping); ++ const std::set &boundary_indicators, ++ std::vector &mapping); /** * Return a list of support @@@ -2479,13 -2509,12 +2498,13 @@@ template void map_dofs_to_support_points (const Mapping &mapping, - const DoFHandler &dof_handler, - std::vector > &support_points); + const DoFHandler &dof_handler, + std::vector > &support_points); /** - * Same as the previous function but for the hp case. + * Same as above for the hp case. */ + template void map_dofs_to_support_points (const dealii::hp::MappingCollection &mapping, @@@ -2524,8 -2595,8 +2543,8 @@@ template void map_support_points_to_dofs (const Mapping &mapping, - const DH &dof_handler, - std::map, types::global_dof_index, Comp> &point_to_index_map); - const DH &dof_handler, - std::map, unsigned int, Comp> &point_to_index_map); ++ const DH &dof_handler, ++ std::map, types::global_dof_index, Comp> &point_to_index_map); /** * Map a coupling table from the diff --cc deal.II/include/deal.II/multigrid/mg_constrained_dofs.h index 4288c6d92c,e011893a8e..5b90025982 --- a/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h @@@ -59,14 -58,9 +59,14 @@@ class MGConstrainedDoFs : public Subscr */ template void initialize(const MGDoFHandler& dof, - const typename FunctionMap::type& function_map, - const std::vector& component_mask = std::vector()); + const typename FunctionMap::type& function_map, + const ComponentMask &component_mask = ComponentMask()); + template + void initialize(const DoFHandler& dof, + const typename FunctionMap::type& function_map, + const std::vector& component_mask = std::vector()); + /** * Reset the data structures. */ diff --cc deal.II/include/deal.II/multigrid/mg_tools.h index b9b2b94cd8,3dcf1b751e..c724e72280 --- a/deal.II/include/deal.II/multigrid/mg_tools.h +++ b/deal.II/include/deal.II/multigrid/mg_tools.h @@@ -173,11 -173,11 +173,11 @@@ namespace MGTool * for each block (access is * result[level][block]). */ - template + template void - count_dofs_per_block (const MGDoFHandler &mg_dof, - std::vector > &result, - std::vector target_block = std::vector()); + count_dofs_per_block (const DH& dof_handler, - std::vector >& dofs_per_block, - std::vector target_block = std::vector()); ++ std::vector >& dofs_per_block, ++ std::vector target_block = std::vector()); /** * Count the dofs component-wise diff --cc deal.II/include/deal.II/numerics/vector_tools.h index b770a80853,5e3cd82257..6f489e0511 --- a/deal.II/include/deal.II/numerics/vector_tools.h +++ b/deal.II/include/deal.II/numerics/vector_tools.h @@@ -747,8 -747,8 +747,8 @@@ namespace VectorTool interpolate_boundary_values (const Mapping &mapping, const DH &dof, const typename FunctionMap::type &function_map, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask = std::vector()); + const ComponentMask &component_mask = ComponentMask()); /** * @deprecated This function exists mainly @@@ -769,8 -769,8 +769,8 @@@ const DH &dof, const types::boundary_id boundary_component, const Function &boundary_function, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask = std::vector()); + const ComponentMask &component_mask = ComponentMask()); /** * Calls the other @@@ -783,8 -783,8 +783,8 @@@ interpolate_boundary_values (const DH &dof, const types::boundary_id boundary_component, const Function &boundary_function, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask = std::vector()); + const ComponentMask &component_mask = ComponentMask()); /** @@@ -797,8 -797,8 +797,8 @@@ void interpolate_boundary_values (const DH &dof, const typename FunctionMap::type &function_map, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask = std::vector()); + const ComponentMask &component_mask = ComponentMask()); /** diff --cc deal.II/include/deal.II/numerics/vector_tools.templates.h index ed9ffbac25,aa7adf8d71..c7c30bc364 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@@ -1477,8 -1477,8 +1477,8 @@@ namespace VectorTool void interpolate_boundary_values (const Mapping &, const DH &dof, const typename FunctionMap::type &function_map, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask_, + const ComponentMask &component_mask, const dealii::internal::int2type<1>) { const unsigned int dim = DH::dimension; @@@ -1568,8 -1561,8 +1561,8 @@@ interpolate_boundary_values (const Mapping &mapping, const DH &dof, const typename FunctionMap::type &function_map, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask_, + const ComponentMask &component_mask, const dealii::internal::int2type) { const unsigned int dim = DH::dimension; @@@ -1599,17 -1591,8 +1591,8 @@@ Assert (n_components == i->second->n_components, ExcDimensionMismatch(n_components, i->second->n_components)); - // set the component mask to either - // the original value or a vector - // of trues - const std::vector component_mask ((component_mask_.size() == 0) ? - std::vector (n_components, true) : - component_mask_); - Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, - ExcNoComponentSelected()); - // field to store the indices - std::vector face_dofs; + std::vector face_dofs; face_dofs.reserve (DoFTools::max_dofs_per_face(dof)); std::vector > dof_locations; @@@ -1927,8 -1906,8 +1906,8 @@@ interpolate_boundary_values (const Mapping &mapping, const DH &dof, const typename FunctionMap::type &function_map, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask_) + const ComponentMask &component_mask_) { internal:: interpolate_boundary_values (mapping, dof, function_map, boundary_values, @@@ -1944,8 -1923,8 +1923,8 @@@ const DH &dof, const types::boundary_id boundary_component, const Function &boundary_function, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask) + const ComponentMask &component_mask) { typename FunctionMap::type function_map; function_map[boundary_component] = &boundary_function; @@@ -1960,8 -1939,8 +1939,8 @@@ interpolate_boundary_values (const DH &dof, const types::boundary_id boundary_component, const Function &boundary_function, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask) + const ComponentMask &component_mask) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); interpolate_boundary_values(StaticMappingQ1::mapping, @@@ -1975,8 -1954,8 +1954,8 @@@ void interpolate_boundary_values (const DH &dof, const typename FunctionMap::type &function_map, - std::map &boundary_values, + std::map &boundary_values, - const std::vector &component_mask) + const ComponentMask &component_mask) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); interpolate_boundary_values(StaticMappingQ1::mapping, @@@ -1998,12 -1977,12 +1977,12 @@@ const DH &dof, const typename FunctionMap::type &function_map, ConstraintMatrix &constraints, - const std::vector &component_mask_) + const ComponentMask &component_mask_) { - std::map boundary_values; + std::map boundary_values; interpolate_boundary_values (mapping, dof, function_map, boundary_values, component_mask_); - std::map::const_iterator boundary_value = + std::map::const_iterator boundary_value = boundary_values.begin(); for ( ; boundary_value !=boundary_values.end(); ++boundary_value) { @@@ -3549,8 -3565,9 +3565,9 @@@ update_JxW_values | update_quadrature_points | update_values); + std::vector dofs_processed; std::vector dof_values; - std::vector face_dof_indices; + std::vector face_dof_indices; typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (); switch (dim) diff --cc deal.II/source/dofs/dof_handler.cc index b808d04f5e,1d4cacf434..0a68adaa5c --- a/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/source/dofs/dof_handler.cc @@@ -1137,15 -714,15 +1044,15 @@@ DoFHandler::last_active_f switch (dim) { case 1: - Assert (false, ExcImpossibleInDim(1)); - return raw_face_iterator(); + Assert (false, ExcImpossibleInDim(1)); + return raw_face_iterator(); case 2: - return last_active_line (); + return last_active_line (); case 3: - return last_active_quad (); + return last_active_quad (); default: - Assert (false, ExcNotImplemented()); - return raw_face_iterator (); + Assert (false, ExcNotImplemented()); + return face_iterator (); } } @@@ -1155,30 -741,34 +1071,32 @@@ DoFHandler::begin_raw_lin template - typename DoFHandler::raw_line_iterator - DoFHandler::begin_raw_line (const unsigned int level) const + typename DoFHandler::line_iterator + DoFHandler::begin_line (const unsigned int level) const { - Assert(tria != 0, ExcNotInitialized()); - switch (dim) - { - case 1: - Assert (leveln_levels(), ExcInvalidLevel(level)); + typename Triangulation::line_iterator ti = tria->begin_line(level); + typename DoFHandler::line_iterator ri(*ti, this); + return ri; + } - if (tria->n_raw_lines(level) == 0) - return end_line (); - - + template + typename DoFHandler::active_line_iterator + DoFHandler::begin_active_line (const unsigned int level) const + { + typename Triangulation::active_line_iterator ti = tria->begin_active_line(level); + typename DoFHandler::active_line_iterator ri(*ti, this); + return ri; + } - return raw_line_iterator (tria, - level, - 0, - this); - default: - Assert (level == 0, ExcFacesHaveNoLevel()); - return raw_line_iterator (tria, - 0, - 0, - this); - } + + template + typename DoFHandler::raw_line_iterator + DoFHandler::end_raw_line (const unsigned int level) const + { + typename Triangulation::raw_line_iterator ti = tria->end_raw_line(level); + typename DoFHandler::raw_line_iterator ri(ti.access_any(), this); + return ri; } diff --cc deal.II/source/dofs/dof_tools.cc index fdf9356d4e,a955d6d26b..ab9851811a --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@@ -80,13 -80,13 +80,13 @@@ namespace DoFTool || (subdomain_id == dof.get_tria().locally_owned_subdomain()), ExcMessage ("For parallel::distributed::Triangulation objects and " - "associated DoF handler objects, asking for any subdomain other " - "than the locally owned one does not make sense.")); + "associated DoF handler objects, asking for any subdomain other " + "than the locally owned one does not make sense.")); - std::vector dofs_on_this_cell; + std::vector dofs_on_this_cell; dofs_on_this_cell.reserve (max_dofs_per_cell(dof)); typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); + endc = dof.end(); // In case we work with a distributed // sparsity pattern of Trilinos type, we @@@ -174,46 -174,38 +174,38 @@@ if (need_dof_mask == true) for (unsigned int f=0; f dofs_on_this_cell(fe_collection.max_dofs_per_cell()); + std::vector dofs_on_this_cell(fe_collection.max_dofs_per_cell()); typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); + endc = dof.end(); // In case we work with a distributed // sparsity pattern of Trilinos type, we @@@ -264,7 -256,7 +256,7 @@@ const std::list > ++ typename DH::cell_iterator> > cell_list = GridTools::get_finest_common_cells (dof_row, dof_col); @@@ -276,74 -268,74 +268,74 @@@ for (; cell_iter!=cell_list.end(); ++cell_iter) { - const typename DH::cell_iterator cell_row = cell_iter->first; - const typename DH::cell_iterator cell_col = cell_iter->second; - - if (!cell_row->has_children() && !cell_col->has_children()) + const typename DH::cell_iterator cell_row = cell_iter->first; + const typename DH::cell_iterator cell_col = cell_iter->second; + - if (!cell_row->has_children() && !cell_col->has_children()) - { - const unsigned int dofs_per_cell_row = - cell_row->get_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row->get_dof_indices (local_dof_indices_row); - cell_col->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; ihas_children()) - { - const std::vector - child_cells = GridTools::get_active_child_cells (cell_row); - for (unsigned int i=0; iget_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row_child->get_dof_indices (local_dof_indices_row); - cell_col->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; i - child_cells = GridTools::get_active_child_cells (cell_col); - for (unsigned int i=0; iget_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col_child->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row->get_dof_indices (local_dof_indices_row); - cell_col_child->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; ihas_children() && !cell_col->has_children()) ++ { ++ const unsigned int dofs_per_cell_row = ++ cell_row->get_fe().dofs_per_cell; ++ const unsigned int dofs_per_cell_col = ++ cell_col->get_fe().dofs_per_cell; ++ std::vector ++ local_dof_indices_row(dofs_per_cell_row); ++ std::vector ++ local_dof_indices_col(dofs_per_cell_col); ++ cell_row->get_dof_indices (local_dof_indices_row); ++ cell_col->get_dof_indices (local_dof_indices_col); ++ for (unsigned int i=0; ihas_children()) ++ { ++ const std::vector ++ child_cells = GridTools::get_active_child_cells (cell_row); ++ for (unsigned int i=0; iget_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row->get_dof_indices (local_dof_indices_row); - cell_col->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; iget_fe().dofs_per_cell; ++ const unsigned int dofs_per_cell_col = ++ cell_col->get_fe().dofs_per_cell; ++ std::vector ++ local_dof_indices_row(dofs_per_cell_row); ++ std::vector ++ local_dof_indices_col(dofs_per_cell_col); ++ cell_row_child->get_dof_indices (local_dof_indices_row); ++ cell_col->get_dof_indices (local_dof_indices_col); ++ for (unsigned int i=0; ihas_children()) ++ } ++ else ++ { ++ std::vector ++ child_cells = GridTools::get_active_child_cells (cell_col); ++ for (unsigned int i=0; i - child_cells = GridTools::get_active_child_cells (cell_row); - for (unsigned int i=0; iget_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row_child->get_dof_indices (local_dof_indices_row); - cell_col->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; i - child_cells = GridTools::get_active_child_cells (cell_col); - for (unsigned int i=0; iget_fe().dofs_per_cell; - const unsigned int dofs_per_cell_col = - cell_col_child->get_fe().dofs_per_cell; - std::vector - local_dof_indices_row(dofs_per_cell_row); - std::vector - local_dof_indices_col(dofs_per_cell_col); - cell_row->get_dof_indices (local_dof_indices_row); - cell_col_child->get_dof_indices (local_dof_indices_col); - for (unsigned int i=0; iget_fe().dofs_per_cell; ++ const unsigned int dofs_per_cell_col = ++ cell_col_child->get_fe().dofs_per_cell; ++ std::vector ++ local_dof_indices_row(dofs_per_cell_row); ++ std::vector ++ local_dof_indices_col(dofs_per_cell_col); ++ cell_row->get_dof_indices (local_dof_indices_row); ++ cell_col_child->get_dof_indices (local_dof_indices_col); ++ for (unsigned int i=0; i::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != DH::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - AssertDimension (max_element, sparsity.n_rows()-1); - unsigned int max_element = 0; - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != DH::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - AssertDimension (max_element, sparsity.n_rows()-1); ++ types::global_dof_index max_element = 0; ++ for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); ++ i!=dof_to_boundary_mapping.end(); ++i) ++ if ((*i != DH::invalid_dof_index) && ++ (*i > max_element)) ++ max_element = *i; ++ AssertDimension (max_element, sparsity.n_rows()-1); }; #endif @@@ -446,14 -438,14 +438,14 @@@ // find active cell at that // boundary: first go to // left/right, then to children - typename DH::cell_iterator cell = dof.begin(0); - while (!cell->at_boundary(direction)) - cell = cell->neighbor(direction); - while (!cell->active()) - cell = cell->child(direction); + typename DH::cell_iterator cell = dof.begin(0); + while (!cell->at_boundary(direction)) + cell = cell->neighbor(direction); + while (!cell->active()) + cell = cell->child(direction); - const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex; - std::vector boundary_dof_boundary_indices (dofs_per_vertex); - const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex; - std::vector boundary_dof_boundary_indices (dofs_per_vertex); ++ const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex; ++ std::vector boundary_dof_boundary_indices (dofs_per_vertex); // next get boundary mapped dof // indices of boundary dofs @@@ -481,36 -473,36 +473,36 @@@ #ifdef DEBUG if (sparsity.n_rows() != 0) { - unsigned int max_element = 0; - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != DH::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - AssertDimension (max_element, sparsity.n_rows()-1); - unsigned int max_element = 0; - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != DH::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - AssertDimension (max_element, sparsity.n_rows()-1); ++ types::global_dof_index max_element = 0; ++ for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); ++ i!=dof_to_boundary_mapping.end(); ++i) ++ if ((*i != DH::invalid_dof_index) && ++ (*i > max_element)) ++ max_element = *i; ++ AssertDimension (max_element, sparsity.n_rows()-1); }; #endif - std::vector dofs_on_this_face; + std::vector dofs_on_this_face; dofs_on_this_face.reserve (max_dofs_per_face(dof)); typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); + endc = dof.end(); for (; cell!=endc; ++cell) for (unsigned int f=0; f::faces_per_cell; ++f) - if (boundary_indicators.find(cell->face(f)->boundary_indicator()) != - boundary_indicators.end()) - { - const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; - dofs_on_this_face.resize (dofs_per_face); - cell->face(f)->get_dof_indices (dofs_on_this_face, - cell->active_fe_index()); + if (boundary_indicators.find(cell->face(f)->boundary_indicator()) != + boundary_indicators.end()) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_this_face.resize (dofs_per_face); + cell->face(f)->get_dof_indices (dofs_on_this_face, + cell->active_fe_index()); // make sparsity pattern for this cell - for (unsigned int i=0; i dofs_on_this_cell; - std::vector dofs_on_other_cell; + std::vector dofs_on_this_cell; + std::vector dofs_on_other_cell; dofs_on_this_cell.reserve (max_dofs_per_cell(dof)); dofs_on_other_cell.reserve (max_dofs_per_cell(dof)); typename DH::active_cell_iterator cell = dof.begin_active(), @@@ -742,230 -726,230 +726,230 @@@ template void make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const Table<2,Coupling> &int_mask, - const Table<2,Coupling> &flux_mask) + SparsityPattern &sparsity, + const Table<2,Coupling> &int_mask, + const Table<2,Coupling> &flux_mask) { - const FiniteElement &fe = dof.get_fe(); + const FiniteElement &fe = dof.get_fe(); - std::vector dofs_on_this_cell(fe.dofs_per_cell); - std::vector dofs_on_other_cell(fe.dofs_per_cell); - std::vector dofs_on_this_cell(fe.dofs_per_cell); - std::vector dofs_on_other_cell(fe.dofs_per_cell); ++ std::vector dofs_on_this_cell(fe.dofs_per_cell); ++ std::vector dofs_on_other_cell(fe.dofs_per_cell); - const Table<2,Coupling> - int_dof_mask = dof_couplings_from_component_couplings(fe, int_mask), - flux_dof_mask = dof_couplings_from_component_couplings(fe, flux_mask); + const Table<2,Coupling> + int_dof_mask = dof_couplings_from_component_couplings(fe, int_mask), + flux_dof_mask = dof_couplings_from_component_couplings(fe, flux_mask); - Table<2,bool> support_on_face(fe.dofs_per_cell, - GeometryInfo::faces_per_cell); - for (unsigned int i=0; i::faces_per_cell;++f) - support_on_face(i,f) = fe.has_support_on_face(i,f); + Table<2,bool> support_on_face(fe.dofs_per_cell, + GeometryInfo::faces_per_cell); + for (unsigned int i=0; i::faces_per_cell;++f) + support_on_face(i,f) = fe.has_support_on_face(i,f); - typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) if (cell->is_locally_owned()) - { - cell->get_dof_indices (dofs_on_this_cell); + { + cell->get_dof_indices (dofs_on_this_cell); // make sparsity pattern for this cell - for (unsigned int i=0; i::faces_per_cell; - ++face) - { - const typename DH::face_iterator - cell_face = cell->face(face); - if (cell_face->user_flag_set ()) - continue; - - if (cell->at_boundary (face) ) - { - for (unsigned int i=0; ineighbor(face); + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; + ++face) + { + const typename DH::face_iterator + cell_face = cell->face(face); + if (cell_face->user_flag_set ()) + continue; + + if (cell->at_boundary (face) ) + { + for (unsigned int i=0; ineighbor(face); // Refinement edges are taken care of // by coarser cells - if (cell->neighbor_is_coarser(face)) - continue; - - typename DH::face_iterator cell_face = cell->face(face); - const unsigned int - neighbor_face = cell->neighbor_of_neighbor(face); - - if (cell_face->has_children()) - { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->n_children(); - ++sub_nr) - { - const typename DH::cell_iterator - sub_neighbor - = cell->neighbor_child_on_subface (face, sub_nr); - - sub_neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); - } - } - else - { - neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); - } - } - } + if (cell->neighbor_is_coarser(face)) + continue; + + typename DH::face_iterator cell_face = cell->face(face); + const unsigned int + neighbor_face = cell->neighbor_of_neighbor(face); + + if (cell_face->has_children()) + { + for (unsigned int sub_nr = 0; + sub_nr != cell_face->n_children(); + ++sub_nr) + { + const typename DH::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); + + sub_neighbor->get_dof_indices (dofs_on_other_cell); + for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); + } + } + else + { + neighbor->get_dof_indices (dofs_on_other_cell); + for (unsigned int i=0; iface(neighbor_face)->set_user_flag (); + } + } + } } } @@@ -992,174 -976,174 +976,174 @@@ // probably less efficient but at least // readable... - const dealii::hp::FECollection &fe = dof.get_fe(); + const dealii::hp::FECollection &fe = dof.get_fe(); - std::vector dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof)); - std::vector dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof)); - std::vector dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof)); - std::vector dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof)); ++ std::vector dofs_on_this_cell(DoFTools::max_dofs_per_cell(dof)); ++ std::vector dofs_on_other_cell(DoFTools::max_dofs_per_cell(dof)); - const std::vector > - int_dof_mask - = dof_couplings_from_component_couplings(fe, int_mask); + const std::vector > + int_dof_mask + = dof_couplings_from_component_couplings(fe, int_mask); - typename dealii::hp::DoFHandler::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - { - dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell); - cell->get_dof_indices (dofs_on_this_cell); + typename dealii::hp::DoFHandler::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + { + dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell); + cell->get_dof_indices (dofs_on_this_cell); // make sparsity pattern for this cell - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) - if (int_dof_mask[cell->active_fe_index()](i,j) != none) - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) + if (int_dof_mask[cell->active_fe_index()](i,j) != none) + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); // Loop over all interior neighbors - for (unsigned int face = 0; - face < GeometryInfo::faces_per_cell; - ++face) - { - const typename dealii::hp::DoFHandler::face_iterator - cell_face = cell->face(face); - if (cell_face->user_flag_set ()) - continue; - - if (cell->at_boundary (face) ) - { - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) - if ((flux_mask(cell->get_fe().system_to_component_index(i).first, - cell->get_fe().system_to_component_index(j).first) - == always) - || - (flux_mask(cell->get_fe().system_to_component_index(i).first, - cell->get_fe().system_to_component_index(j).first) - == nonzero)) - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - } - else - { - typename dealii::hp::DoFHandler::cell_iterator - neighbor = cell->neighbor(face); + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; + ++face) + { + const typename dealii::hp::DoFHandler::face_iterator + cell_face = cell->face(face); + if (cell_face->user_flag_set ()) + continue; + + if (cell->at_boundary (face) ) + { + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) + if ((flux_mask(cell->get_fe().system_to_component_index(i).first, + cell->get_fe().system_to_component_index(j).first) + == always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + cell->get_fe().system_to_component_index(j).first) + == nonzero)) + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); + } + else + { + typename dealii::hp::DoFHandler::cell_iterator + neighbor = cell->neighbor(face); // Refinement edges are taken care of // by coarser cells - if (cell->neighbor_is_coarser(face)) - continue; - - typename dealii::hp::DoFHandler::face_iterator - cell_face = cell->face(face); - const unsigned int - neighbor_face = cell->neighbor_of_neighbor(face); - - if (cell_face->has_children()) - { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->n_children(); - ++sub_nr) - { - const typename dealii::hp::DoFHandler::cell_iterator - sub_neighbor - = cell->neighbor_child_on_subface (face, sub_nr); - - dofs_on_other_cell.resize (sub_neighbor->get_fe().dofs_per_cell); - sub_neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - { - for (unsigned int j=0; jget_fe().dofs_per_cell; - ++j) - { - if ((flux_mask(cell->get_fe().system_to_component_index(i).first, - sub_neighbor->get_fe().system_to_component_index(j).first) - == always) - || - (flux_mask(cell->get_fe().system_to_component_index(i).first, - sub_neighbor->get_fe().system_to_component_index(j).first) - == nonzero)) - { - sparsity.add (dofs_on_this_cell[i], - dofs_on_other_cell[j]); - sparsity.add (dofs_on_other_cell[i], - dofs_on_this_cell[j]); - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - sparsity.add (dofs_on_other_cell[i], - dofs_on_other_cell[j]); - } - - if ((flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, - cell->get_fe().system_to_component_index(i).first) - == always) - || - (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, - cell->get_fe().system_to_component_index(i).first) - == nonzero)) - { - sparsity.add (dofs_on_this_cell[j], - dofs_on_other_cell[i]); - sparsity.add (dofs_on_other_cell[j], - dofs_on_this_cell[i]); - sparsity.add (dofs_on_this_cell[j], - dofs_on_this_cell[i]); - sparsity.add (dofs_on_other_cell[j], - dofs_on_other_cell[i]); - } - } - } - sub_neighbor->face(neighbor_face)->set_user_flag (); - } - } - else - { - dofs_on_other_cell.resize (neighbor->get_fe().dofs_per_cell); - neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) - { - for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) - { - if ((flux_mask(cell->get_fe().system_to_component_index(i).first, - neighbor->get_fe().system_to_component_index(j).first) - == always) - || - (flux_mask(cell->get_fe().system_to_component_index(i).first, - neighbor->get_fe().system_to_component_index(j).first) - == nonzero)) - { - sparsity.add (dofs_on_this_cell[i], - dofs_on_other_cell[j]); - sparsity.add (dofs_on_other_cell[i], - dofs_on_this_cell[j]); - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - sparsity.add (dofs_on_other_cell[i], - dofs_on_other_cell[j]); - } - - if ((flux_mask(neighbor->get_fe().system_to_component_index(j).first, - cell->get_fe().system_to_component_index(i).first) - == always) - || - (flux_mask(neighbor->get_fe().system_to_component_index(j).first, - cell->get_fe().system_to_component_index(i).first) - == nonzero)) - { - sparsity.add (dofs_on_this_cell[j], - dofs_on_other_cell[i]); - sparsity.add (dofs_on_other_cell[j], - dofs_on_this_cell[i]); - sparsity.add (dofs_on_this_cell[j], - dofs_on_this_cell[i]); - sparsity.add (dofs_on_other_cell[j], - dofs_on_other_cell[i]); - } - } - } - neighbor->face(neighbor_face)->set_user_flag (); - } - } - } - } + if (cell->neighbor_is_coarser(face)) + continue; + + typename dealii::hp::DoFHandler::face_iterator + cell_face = cell->face(face); + const unsigned int + neighbor_face = cell->neighbor_of_neighbor(face); + + if (cell_face->has_children()) + { + for (unsigned int sub_nr = 0; + sub_nr != cell_face->n_children(); + ++sub_nr) + { + const typename dealii::hp::DoFHandler::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); + + dofs_on_other_cell.resize (sub_neighbor->get_fe().dofs_per_cell); + sub_neighbor->get_dof_indices (dofs_on_other_cell); + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + { + for (unsigned int j=0; jget_fe().dofs_per_cell; + ++j) + { + if ((flux_mask(cell->get_fe().system_to_component_index(i).first, + sub_neighbor->get_fe().system_to_component_index(j).first) + == always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + sub_neighbor->get_fe().system_to_component_index(j).first) + == nonzero)) + { + sparsity.add (dofs_on_this_cell[i], + dofs_on_other_cell[j]); + sparsity.add (dofs_on_other_cell[i], + dofs_on_this_cell[j]); + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); + sparsity.add (dofs_on_other_cell[i], + dofs_on_other_cell[j]); + } + + if ((flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == always) + || + (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == nonzero)) + { + sparsity.add (dofs_on_this_cell[j], + dofs_on_other_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_this_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_other_cell[i]); + } + } + } + sub_neighbor->face(neighbor_face)->set_user_flag (); + } + } + else + { + dofs_on_other_cell.resize (neighbor->get_fe().dofs_per_cell); + neighbor->get_dof_indices (dofs_on_other_cell); + for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) + { + for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) + { + if ((flux_mask(cell->get_fe().system_to_component_index(i).first, + neighbor->get_fe().system_to_component_index(j).first) + == always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + neighbor->get_fe().system_to_component_index(j).first) + == nonzero)) + { + sparsity.add (dofs_on_this_cell[i], + dofs_on_other_cell[j]); + sparsity.add (dofs_on_other_cell[i], + dofs_on_this_cell[j]); + sparsity.add (dofs_on_this_cell[i], + dofs_on_this_cell[j]); + sparsity.add (dofs_on_other_cell[i], + dofs_on_other_cell[j]); + } + + if ((flux_mask(neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == always) + || + (flux_mask(neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == nonzero)) + { + sparsity.add (dofs_on_this_cell[j], + dofs_on_other_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_this_cell[j], + dofs_on_this_cell[i]); + sparsity.add (dofs_on_other_cell[j], + dofs_on_other_cell[i]); + } + } + } + neighbor->face(neighbor_face)->set_user_flag (); + } + } + } + } } } @@@ -1228,14 -1212,14 +1212,14 @@@ inline bool check_master_dof_list (const FullMatrix &face_interpolation_matrix, - const std::vector &master_dof_list) - const std::vector &master_dof_list) ++ const std::vector &master_dof_list) { - const unsigned int N = master_dof_list.size(); + const unsigned int N = master_dof_list.size(); - FullMatrix tmp (N,N); - for (unsigned int i=0; i tmp (N,N); + for (unsigned int i=0; i master_dof_list; - unsigned int index = 0; - for (int v=0; - v(GeometryInfo::vertices_per_face); - ++v) - std::vector master_dof_list; - unsigned int index = 0; - for (int v=0; - v(GeometryInfo::vertices_per_face); - ++v) - { - unsigned int dofs_added = 0; - unsigned int i = 0; - while (dofs_added < fe2.dofs_per_vertex) - { ++ std::vector master_dof_list; ++ unsigned int index = 0; ++ for (int v=0; ++ v(GeometryInfo::vertices_per_face); ++ ++v) ++ { ++ unsigned int dofs_added = 0; ++ unsigned int i = 0; ++ while (dofs_added < fe2.dofs_per_vertex) + { - unsigned int dofs_added = 0; - unsigned int i = 0; - while (dofs_added < fe2.dofs_per_vertex) - { // make sure that we // were able to find // a set of master @@@ -1533,10 -1517,10 +1517,10 @@@ // finally copy the list into the // mask - std::fill (master_dof_mask.begin(), master_dof_mask.end(), false); - for (std::vector::const_iterator i=master_dof_list.begin(); - i!=master_dof_list.end(); ++i) - master_dof_mask[*i] = true; - std::fill (master_dof_mask.begin(), master_dof_mask.end(), false); - for (std::vector::const_iterator i=master_dof_list.begin(); - i!=master_dof_list.end(); ++i) - master_dof_mask[*i] = true; ++ std::fill (master_dof_mask.begin(), master_dof_mask.end(), false); ++ for (std::vector::const_iterator i=master_dof_list.begin(); ++ i!=master_dof_list.end(); ++i) ++ master_dof_mask[*i] = true; } @@@ -1826,20 -1810,20 +1810,20 @@@ * fuller than necessary. */ void - filter_constraints (const std::vector &master_dofs, - const std::vector &slave_dofs, - const FullMatrix &face_constraints, - ConstraintMatrix &constraints) + filter_constraints (const std::vector &master_dofs, - const std::vector &slave_dofs, - const FullMatrix &face_constraints, - ConstraintMatrix &constraints) ++ const std::vector &slave_dofs, ++ const FullMatrix &face_constraints, ++ ConstraintMatrix &constraints) { - Assert (face_constraints.n () == master_dofs.size (), - ExcDimensionMismatch(master_dofs.size (), - face_constraints.n())); - Assert (face_constraints.m () == slave_dofs.size (), - ExcDimensionMismatch(slave_dofs.size (), - face_constraints.m())); + Assert (face_constraints.n () == master_dofs.size (), + ExcDimensionMismatch(master_dofs.size (), + face_constraints.n())); + Assert (face_constraints.m () == slave_dofs.size (), + ExcDimensionMismatch(slave_dofs.size (), + face_constraints.m())); - const unsigned int n_master_dofs = master_dofs.size (); - const unsigned int n_slave_dofs = slave_dofs.size (); + const unsigned int n_master_dofs = master_dofs.size (); + const unsigned int n_slave_dofs = slave_dofs.size (); // check for a couple // conditions that happened @@@ -3328,22 -3312,22 +3312,22 @@@ } return; } -- // .. otherwise we should be in the case -- // were both faces are active and have -- // no children .. ++ // .. otherwise we should be in the case ++ // were both faces are active and have ++ // no children .. Assert (!face_1->has_children() && !face_2->has_children(), ExcNotImplemented()); Assert (face_1->n_active_fe_indices() == 1 && face_2->n_active_fe_indices() == 1, ExcInternalError()); -- // .. then we match the -- // corresponding DoFs of both faces .. ++ // .. then we match the ++ // corresponding DoFs of both faces .. const unsigned int face_1_index = face_1->nth_active_fe_index(0); const unsigned int face_2_index = face_2->nth_active_fe_index(0); Assert ( face_1->get_fe(face_1_index) -- == face_2->get_fe(face_1_index), -- ExcMessage ("Matching periodic cells need to use the same finite element")); ++ == face_2->get_fe(face_1_index), ++ ExcMessage ("Matching periodic cells need to use the same finite element")); const dealii::FiniteElement &fe = face_1->get_fe(face_1_index); @@@ -3361,16 -3344,19 +3344,19 @@@ face_1->get_dof_indices(dofs_1, face_1_index); face_2->get_dof_indices(dofs_2, face_2_index); -- // .. and constrain them (respecting -- // component_mask): ++ // .. and constrain them (respecting ++ // component_mask): for (unsigned int i = 0; i < dofs_per_face; ++i) { - if (component_mask.size() == 0 || - component_mask[fe.face_system_to_component_index(i).first]) { - if(!constraint_matrix.is_constrained(dofs_1[i])) { - constraint_matrix.add_line(dofs_1[i]); - constraint_matrix.add_entry(dofs_1[i], dofs_2[i], 1.0); + if ((component_mask.n_selected_components(fe.n_components()) == fe.n_components()) + || + (component_mask[fe.face_system_to_component_index(i).first] == true)) + { + if(!constraint_matrix.is_constrained(dofs_1[i])) + { + constraint_matrix.add_line(dofs_1[i]); + constraint_matrix.add_entry(dofs_1[i], dofs_2[i], 1.0); + } } - } } } @@@ -3399,9 -3385,9 +3385,9 @@@ const types::boundary_id boundary_component, const int direction, dealii::Tensor<1,DH::space_dimension> -- &offset, ++ &offset, dealii::ConstraintMatrix &constraint_matrix, - const std::vector &component_mask) + const ComponentMask &component_mask) { static const int space_dim = DH::space_dimension; Assert (0<=direction && direction matched_cells = -- GridTools::collect_periodic_cell_pairs(dof_handler.begin(0), -- dof_handler.end(0), -- boundary_component, -- direction, -- offset); -- -- // And apply the low level -- // make_periodicity_constraints function -- // to every matching pair: ++ GridTools::collect_periodic_cell_pairs(dof_handler.begin(0), ++ dof_handler.end(0), ++ boundary_component, ++ direction, ++ offset); ++ ++ // And apply the low level ++ // make_periodicity_constraints function ++ // to every matching pair: for (typename std::map::iterator it = matched_cells.begin(); - it != matched_cells.end(); ++it) + it != matched_cells.end(); ++it) { - typedef typename DH::face_iterator FaceIterator; - FaceIterator face_1 = it->first->face(2*direction); - FaceIterator face_2 = it->second->face(2*direction+1); + typedef typename DH::face_iterator FaceIterator; + FaceIterator face_1 = it->first->face(2*direction); + FaceIterator face_2 = it->second->face(2*direction+1); - Assert(face_1->at_boundary() && face_2->at_boundary(), - ExcInternalError()); + Assert(face_1->at_boundary() && face_2->at_boundary(), + ExcInternalError()); - Assert (face_1->boundary_indicator() == boundary_component && - face_2->boundary_indicator() == boundary_component, - ExcInternalError()); + Assert (face_1->boundary_indicator() == boundary_component && + face_2->boundary_indicator() == boundary_component, + ExcInternalError()); - make_periodicity_constraints(face_1, - face_2, - constraint_matrix, - component_mask); + make_periodicity_constraints(face_1, + face_2, + constraint_matrix, + component_mask); } } @@@ -3501,81 -3557,100 +3557,100 @@@ for (unsigned int f=0; f &fe = - fe_collection[f]; - local_component_association[f].resize(fe.dofs_per_cell); - if (sort_by_blocks == true) - // compute the block each - // local dof belongs to - { - for (unsigned int i=0; i indices; + for (typename DH::active_cell_iterator c=dof.begin_active(); + c!=dof.end(); ++ c) + if (c->is_locally_owned()) + { + const unsigned int fe_index = c->active_fe_index(); + const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell; + indices.resize(dofs_per_cell); + c->get_dof_indices(indices); + for (unsigned int i=0; i + inline + void + get_block_association (const DH &dof, + std::vector &dofs_by_block) + { + const dealii::hp::FECollection + fe_collection (dof.get_fe()); + Assert (fe_collection.n_components() < 256, ExcNotImplemented()); + Assert (dofs_by_block.size() == dof.n_locally_owned_dofs(), + ExcDimensionMismatch(dofs_by_block.size(), + dof.n_locally_owned_dofs())); + + // next set up a table for the degrees + // of freedom on each of the cells + // (regardless of the fact whether it + // is listed in the component_select + // argument or not) + // + // for each element 'f' of the + // FECollection, + // local_block_association[f][d] + // then returns the vector block + // that degree of freedom 'd' belongs + // to + std::vector > local_block_association + (fe_collection.size()); + for (unsigned int f=0; f &fe = + fe_collection[f]; + local_block_association[f].resize(fe.dofs_per_cell, + (unsigned char)(-1)); + for (unsigned int i=0; i indices; + std::vector indices; for (typename DH::active_cell_iterator c=dof.begin_active(); - c!=dof.end(); ++ c) - if (c->is_locally_owned()) - { - const unsigned int fe_index = c->active_fe_index(); - const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell; - indices.resize(dofs_per_cell); - c->get_dof_indices(indices); - for (unsigned int i=0; iis_locally_owned()) + { + const unsigned int fe_index = c->active_fe_index(); + const unsigned int dofs_per_cell = c->get_fe().dofs_per_cell; + indices.resize(dofs_per_cell); + c->get_dof_indices(indices); + for (unsigned int i=0; i touch_count (dof_handler.n_dofs(), 0); typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + endc = dof_handler.end(); - std::vector dof_indices; + std::vector dof_indices; dof_indices.reserve (max_dofs_per_cell(dof_handler)); for (unsigned int present_cell = 0; cell!=endc; ++cell, ++present_cell) @@@ -3799,30 -3868,20 +3868,20 @@@ - template + template void - extract_level_dofs( - const unsigned int level, - const DH &dof, - const std::vector &component_select, - std::vector &selected_dofs, - const bool count_by_blocks) + extract_level_dofs(const unsigned int level, - const MGDoFHandler &dof, ++ const DH &dof, + const ComponentMask &component_mask, + std::vector &selected_dofs) { - const FiniteElement& fe = dof.get_fe(); + const FiniteElement& fe = dof.get_fe(); - if (count_by_blocks == true) - { - Assert(component_select.size() == fe.n_blocks(), - ExcDimensionMismatch(component_select.size(), fe.n_blocks())); - } - else - { - Assert(component_select.size() == fe.n_components(), - ExcDimensionMismatch(component_select.size(), fe.n_components())); - } - + Assert(component_mask.represents_n_components(n_components(dof)), + ExcMessage ("The given component mask is not sized correctly to represent the " + "components of the given finite element.")); Assert(selected_dofs.size() == dof.n_dofs(level), - ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level))); + ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level))); // two special cases: no component // is selected, and all components @@@ -3901,18 -3913,32 +3913,32 @@@ // then loop over all cells and do // work - std::vector indices(fe.dofs_per_cell); - typename MGDoFHandler::cell_iterator c; + std::vector indices(fe.dofs_per_cell); + typename DH::cell_iterator c; for (c = dof.begin(level) ; c != dof.end(level) ; ++ c) { - c->get_mg_dof_indices(indices); - for (unsigned int i=0; iget_mg_dof_indices(indices); + for (unsigned int i=0; i ++ template + void + extract_level_dofs(const unsigned int level, - const MGDoFHandler &dof, ++ const DH &dof, + const BlockMask &block_mask, + std::vector &selected_dofs) + { + // simply defer to the other extract_level_dofs() function + extract_level_dofs (level, dof, dof.get_fe().component_mask(block_mask), + selected_dofs); + } + + + template void extract_boundary_dofs (const DH &dof_handler, @@@ -3928,14 -3954,14 +3954,14 @@@ "See the documentation for more information.")); IndexSet indices; - extract_boundary_dofs (dof_handler, component_select, + extract_boundary_dofs (dof_handler, component_mask, indices, boundary_indicators); -- // clear and reset array by default values ++ // clear and reset array by default values selected_dofs.clear (); selected_dofs.resize (dof_handler.n_dofs(), false); -- // then convert the values computed above to the binary vector ++ // then convert the values computed above to the binary vector indices.fill_binary_vector(selected_dofs); } @@@ -3952,7 -3979,7 +3979,7 @@@ ExcInvalidBoundaryIndicator()); const unsigned int dim=DH::dimension; -- // first reset output argument ++ // first reset output argument selected_dofs.clear (); selected_dofs.set_size(dof_handler.n_dofs()); @@@ -3966,10 -3993,12 +3993,12 @@@ // check whether a certain vector // component is selected, or all const bool check_vector_component - = (component_select != std::vector(component_select.size(), - true)); + = ((component_mask.represents_the_all_selected_mask() == false) + || + (component_mask.n_selected_components(n_components(dof_handler)) != + n_components(dof_handler))); - std::vector face_dof_indices; + std::vector face_dof_indices; face_dof_indices.reserve (max_dofs_per_face(dof_handler)); // now loop over all cells and @@@ -3987,15 -4016,15 +4016,15 @@@ // sooner or later for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) -- // only work on cells that are either locally owned or at -- // least ghost cells ++ // only work on cells that are either locally owned or at ++ // least ghost cells if (cell->is_artificial() == false) for (unsigned int face=0; -- face::faces_per_cell; ++face) ++ face::faces_per_cell; ++face) if (cell->at_boundary(face)) if (! check_boundary_indicator || (boundary_indicators.find (cell->face(face)->boundary_indicator()) -- != boundary_indicators.end())) ++ != boundary_indicators.end())) { const FiniteElement &fe = cell->get_fe(); @@@ -4008,44 -4037,44 +4037,44 @@@ if (!check_vector_component) selected_dofs.add_index (face_dof_indices[i]); else -- // check for -- // component is -- // required. somewhat -- // tricky as usual -- // for the case that -- // the shape function -- // is non-primitive, -- // but use usual -- // convention (see -- // docs) ++ // check for ++ // component is ++ // required. somewhat ++ // tricky as usual ++ // for the case that ++ // the shape function ++ // is non-primitive, ++ // but use usual ++ // convention (see ++ // docs) { -- // first get at the -- // cell-global -- // number of a face -- // dof, to ask the -- // fe certain -- // questions ++ // first get at the ++ // cell-global ++ // number of a face ++ // dof, to ask the ++ // fe certain ++ // questions const unsigned int cell_index -- = (dim == 1 ? -- i -- : -- (dim == 2 ? -- (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex) ++ = (dim == 1 ? ++ i ++ : ++ (dim == 2 ? ++ (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex) ++ : ++ (dim == 3 ? ++ (i<4*fe.dofs_per_vertex ? ++ i : -- (dim == 3 ? -- (i<4*fe.dofs_per_vertex ? -- i -- : -- (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ? -- i+4*fe.dofs_per_vertex -- : -- i+4*fe.dofs_per_vertex+8*fe.dofs_per_line)) -- : -- numbers::invalid_unsigned_int))); ++ (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ? ++ i+4*fe.dofs_per_vertex ++ : ++ i+4*fe.dofs_per_vertex+8*fe.dofs_per_line)) ++ : ++ numbers::invalid_unsigned_int))); if (fe.is_primitive (cell_index)) { - if (component_select[fe.face_system_to_component_index(i).first] - == true) + if (component_mask[fe.face_system_to_component_index(i).first] - == true) ++ == true) selected_dofs.add_index (face_dof_indices[i]); } else // not primitive @@@ -4363,23 -4384,23 +4384,23 @@@ // because we can't be sure that the dof // range of locally_owned_dofs is really // contiguous. - std::vector dof_indices; - std::set global_dof_indices; + std::vector dof_indices; + std::set global_dof_indices; typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + endc = dof_handler.end(); for (; cell!=endc; ++cell) if (cell->is_locally_owned()) - { - dof_indices.resize(cell->get_fe().dofs_per_cell); - cell->get_dof_indices(dof_indices); + { + dof_indices.resize(cell->get_fe().dofs_per_cell); + cell->get_dof_indices(dof_indices); - for (std::vector::iterator it=dof_indices.begin(); - it!=dof_indices.end(); - ++it) - if (!dof_set.is_element(*it)) - global_dof_indices.insert(*it); - } - for (std::vector::iterator it=dof_indices.begin(); - it!=dof_indices.end(); - ++it) - if (!dof_set.is_element(*it)) - global_dof_indices.insert(*it); - } ++ for (std::vector::iterator it=dof_indices.begin(); ++ it!=dof_indices.end(); ++ ++it) ++ if (!dof_set.is_element(*it)) ++ global_dof_indices.insert(*it); ++ } dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end()); @@@ -4402,23 -4423,23 +4423,23 @@@ // because we can't be sure that the dof // range of locally_owned_dofs is really // contiguous. - std::vector dof_indices; - std::set global_dof_indices; + std::vector dof_indices; + std::set global_dof_indices; typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + endc = dof_handler.end(); for (; cell!=endc; ++cell) if (cell->is_ghost()) - { - dof_indices.resize(cell->get_fe().dofs_per_cell); - cell->get_dof_indices(dof_indices); + { + dof_indices.resize(cell->get_fe().dofs_per_cell); + cell->get_dof_indices(dof_indices); - for (std::vector::iterator it=dof_indices.begin(); - it!=dof_indices.end(); - ++it) - if (!dof_set.is_element(*it)) - global_dof_indices.insert(*it); - } - for (std::vector::iterator it=dof_indices.begin(); - it!=dof_indices.end(); - ++it) - if (!dof_set.is_element(*it)) - global_dof_indices.insert(*it); - } ++ for (std::vector::iterator it=dof_indices.begin(); ++ it!=dof_indices.end(); ++ ++it) ++ if (!dof_set.is_element(*it)) ++ global_dof_indices.insert(*it); ++ } dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end()); @@@ -4509,9 -4530,9 +4530,9 @@@ // preset all values by an invalid value std::fill_n (subdomain_association.begin(), dof_handler.n_dofs(), - types::invalid_subdomain_id); + types::invalid_subdomain_id); - std::vector local_dof_indices; + std::vector