From b2fe09235b558ebae1d56b59f0402c753e6dcc54 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 27 Sep 2012 17:02:48 +0000 Subject: [PATCH] Reindent a bunch of places to make merging from mainline simpler. git-svn-id: https://svn.dealii.org/branches/branch_merge_mg_into_dof_handler@26799 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/dofs/dof_tools.h | 202 +- deal.II/include/deal.II/multigrid/mg_tools.h | 6 +- deal.II/source/dofs/dof_tools.cc | 5036 +++++++++--------- deal.II/source/dofs/dof_tools.inst.in | 28 +- 4 files changed, 2636 insertions(+), 2636 deletions(-) diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index 09479d5afb..1c568de5b5 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -201,11 +201,11 @@ namespace DoFTools * Two components do not * couple. */ - none, + none, /** * Two components do couple. */ - always, + always, /** * Two components couple only * if their shape functions are @@ -214,7 +214,7 @@ namespace DoFTools * when computing integrals over * faces of cells. */ - nonzero + nonzero }; /** @@ -493,10 +493,10 @@ namespace DoFTools template void make_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints = ConstraintMatrix(), - const bool keep_constrained_dofs = true, - const types::subdomain_id subdomain_id = types::invalid_subdomain_id); + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints = ConstraintMatrix(), + const bool keep_constrained_dofs = true, + const types::subdomain_id subdomain_id = types::invalid_subdomain_id); /** * Locate non-zero entries for @@ -654,11 +654,11 @@ namespace DoFTools template void make_sparsity_pattern (const DH &dof, - const Table<2, Coupling> &coupling, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints = ConstraintMatrix(), - const bool keep_constrained_dofs = true, - const types::subdomain_id subdomain_id = types::invalid_subdomain_id); + const Table<2, Coupling> &coupling, + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints = ConstraintMatrix(), + const bool keep_constrained_dofs = true, + const types::subdomain_id subdomain_id = types::invalid_subdomain_id); /** * @deprecated This is the old @@ -675,8 +675,8 @@ namespace DoFTools template void make_sparsity_pattern (const DH &dof, - const std::vector > &mask, - SparsityPattern &sparsity_pattern); + const std::vector > &mask, + SparsityPattern &sparsity_pattern); /** * Construct a sparsity pattern that @@ -714,8 +714,8 @@ namespace DoFTools template void make_sparsity_pattern (const DH &dof_row, - const DH &dof_col, - SparsityPattern &sparsity); + const DH &dof_col, + SparsityPattern &sparsity); /** * Create the sparsity pattern for @@ -737,8 +737,8 @@ namespace DoFTools 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); /** * Write the sparsity structure of the @@ -773,9 +773,9 @@ namespace DoFTools 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); /** * Generate sparsity pattern for @@ -796,7 +796,7 @@ namespace DoFTools template void make_flux_sparsity_pattern (const DH &dof_handler, - SparsityPattern &sparsity_pattern); + SparsityPattern &sparsity_pattern); /** * This function does the same as @@ -812,10 +812,10 @@ namespace DoFTools template void make_flux_sparsity_pattern (const DH &dof_handler, - SparsityPattern &sparsity_pattern, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs = true, - const types::subdomain_id subdomain_id = numbers::invalid_unsigned_int); + SparsityPattern &sparsity_pattern, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs = true, + const types::subdomain_id subdomain_id = numbers::invalid_unsigned_int); /** * This function does the same as @@ -840,9 +840,9 @@ namespace DoFTools 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); //@} /** @@ -925,7 +925,7 @@ namespace DoFTools template void make_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints); + ConstraintMatrix &constraints); //@} @@ -1170,9 +1170,9 @@ namespace DoFTools template void distribute_cell_to_dof_vector (const DH &dof_handler, - const Vector &cell_data, - Vector &dof_data, - const unsigned int component = 0); + const Vector &cell_data, + Vector &dof_data, + const unsigned int component = 0); /** * Extract the indices of the @@ -1223,9 +1223,9 @@ namespace DoFTools template void extract_dofs (const DoFHandler &dof_handler, - const std::vector &select, - std::vector &selected_dofs, - const bool blocks = false); + const std::vector &select, + std::vector &selected_dofs, + const bool blocks = false); /** * The same function as above, @@ -1234,9 +1234,9 @@ namespace DoFTools template void extract_dofs (const hp::DoFHandler &dof_handler, - const std::vector &select, - std::vector &selected_dofs, - const bool blocks = false); + const std::vector &select, + std::vector &selected_dofs, + const bool blocks = false); /** * Do the same thing as @@ -1246,10 +1246,10 @@ namespace DoFTools 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 DH &dof, + const std::vector &select, + std::vector &selected_dofs, + const bool blocks = false); /** * Extract all degrees of freedom @@ -1442,9 +1442,9 @@ namespace DoFTools template void extract_dofs_with_support_on_boundary (const DH &dof_handler, - const std::vector &component_mask, - std::vector &selected_dofs, - const std::set &boundary_indicators = std::set()); + const std::vector &component_mask, + std::vector &selected_dofs, + const std::set &boundary_indicators = std::set()); /** * @name Hanging Nodes @@ -1466,7 +1466,7 @@ namespace DoFTools template void extract_hanging_node_dofs (const DoFHandler &dof_handler, - std::vector &selected_dofs); + std::vector &selected_dofs); //@} /** @@ -1494,8 +1494,8 @@ namespace DoFTools template void extract_subdomain_dofs (const DH &dof_handler, - const types::subdomain_id subdomain_id, - std::vector &selected_dofs); + const types::subdomain_id subdomain_id, + std::vector &selected_dofs); /** @@ -1512,7 +1512,7 @@ namespace DoFTools template void extract_locally_owned_dofs (const DH & dof_handler, - IndexSet & dof_set); + IndexSet & dof_set); /** @@ -1546,7 +1546,7 @@ namespace DoFTools template void extract_locally_active_dofs (const DH & dof_handler, - IndexSet & dof_set); + IndexSet & dof_set); /** * Extract the set of global DoF @@ -1567,7 +1567,7 @@ namespace DoFTools template void extract_locally_relevant_dofs (const DH & dof_handler, - IndexSet & dof_set); + IndexSet & dof_set); /** * For each DoF, return in the output @@ -1625,7 +1625,7 @@ namespace DoFTools template void get_subdomain_association (const DH &dof_handler, - std::vector &subdomain); + std::vector &subdomain); /** * Count how many degrees of freedom are @@ -1669,7 +1669,7 @@ namespace DoFTools template unsigned int count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id subdomain); + const types::subdomain_id subdomain); /** * Count how many degrees of freedom are @@ -1706,8 +1706,8 @@ namespace DoFTools template void count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id subdomain, - std::vector &n_dofs_on_subdomain); + const types::subdomain_id subdomain, + std::vector &n_dofs_on_subdomain); /** * Return a set of indices that denotes @@ -1748,7 +1748,7 @@ namespace DoFTools 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 @@ -1992,8 +1992,8 @@ namespace DoFTools template void extract_constant_modes (const DH &dof_handler, - const std::vector &component_mask, - std::vector > &constant_modes); + const std::vector &component_mask, + std::vector > &constant_modes); /** * For each active cell of a DoFHandler @@ -2013,7 +2013,7 @@ namespace DoFTools template void get_active_fe_indices (const DH &dof_handler, - std::vector &active_fe_indices); + std::vector &active_fe_indices); /** * Count how many degrees of @@ -2095,10 +2095,10 @@ namespace DoFTools 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()); /** * Count the degrees of freedom @@ -2149,8 +2149,8 @@ namespace DoFTools 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); /** * This function can be used when @@ -2369,11 +2369,11 @@ namespace DoFTools template void compute_intergrid_constraints (const DoFHandler &coarse_grid, - const unsigned int coarse_component, - const DoFHandler &fine_grid, - const unsigned int fine_component, - const InterGridMap > &coarse_to_fine_grid_map, - ConstraintMatrix &constraints); + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap > &coarse_to_fine_grid_map, + ConstraintMatrix &constraints); /** @@ -2402,11 +2402,11 @@ namespace DoFTools template void compute_intergrid_transfer_representation (const DoFHandler &coarse_grid, - const unsigned int coarse_component, - const DoFHandler &fine_grid, - const unsigned int fine_component, - const InterGridMap > &coarse_to_fine_grid_map, - std::vector > &transfer_representation); + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap > &coarse_to_fine_grid_map, + std::vector > &transfer_representation); /** * Create a mapping from degree @@ -2438,7 +2438,7 @@ namespace DoFTools template void map_dof_to_boundary_indices (const DH &dof_handler, - std::vector &mapping); + std::vector &mapping); /** * Same as the previous function, @@ -2454,8 +2454,8 @@ namespace DoFTools 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); /** * Return a list of support @@ -2479,8 +2479,8 @@ namespace DoFTools 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 above for the hp case. @@ -2524,8 +2524,8 @@ namespace DoFTools template void map_support_points_to_dofs (const Mapping &mapping, - const DH &dof_handler, - std::map, unsigned int, Comp> &point_to_index_map); + const DH &dof_handler, + std::map, unsigned int, Comp> &point_to_index_map); /** * Map a coupling table from the @@ -2545,8 +2545,8 @@ namespace DoFTools template void convert_couplings_to_blocks (const hp::DoFHandler& dof_handler, - const Table<2, Coupling>& table_by_component, - std::vector >& tables_by_block); + const Table<2, Coupling>& table_by_component, + std::vector >& tables_by_block); /** * Make a constraint matrix for the @@ -2581,8 +2581,8 @@ namespace DoFTools template class DH> void make_zero_boundary_constraints (const DH &dof, - ConstraintMatrix &zero_boundary_constraints, - const std::vector &component_mask_=std::vector()); + ConstraintMatrix &zero_boundary_constraints, + const std::vector &component_mask_=std::vector()); /** * Map a coupling table from the @@ -2602,8 +2602,8 @@ namespace DoFTools template void convert_couplings_to_blocks (const DoFHandler& dof_handler, - const Table<2, Coupling>& table_by_component, - std::vector >& tables_by_block); + const Table<2, Coupling>& table_by_component, + std::vector >& tables_by_block); /** * Given a finite element and a table how @@ -2616,7 +2616,7 @@ namespace DoFTools template Table<2,Coupling> dof_couplings_from_component_couplings (const FiniteElement &fe, - const Table<2,Coupling> &component_couplings); + const Table<2,Coupling> &component_couplings); /** * Same function as above for a @@ -2630,7 +2630,7 @@ namespace DoFTools template std::vector > dof_couplings_from_component_couplings (const hp::FECollection &fe, - const Table<2,Coupling> &component_couplings); + const Table<2,Coupling> &component_couplings); /** * Exception * @ingroup Exceptions @@ -2673,7 +2673,7 @@ namespace DoFTools */ inline Coupling operator |= (Coupling& c1, - const Coupling c2) + const Coupling c2) { if (c2 == always) c1 = always; @@ -2690,7 +2690,7 @@ namespace DoFTools */ inline Coupling operator | (const Coupling c1, - const Coupling c2) + const Coupling c2) { if (c1 == always || c2 == always) return always; @@ -2797,24 +2797,24 @@ namespace DoFTools inline void make_sparsity_pattern (const DH &dof, - const std::vector > &mask, - SparsityPattern &sparsity_pattern) + const std::vector > &mask, + SparsityPattern &sparsity_pattern) { const unsigned int ncomp = dof.get_fe().n_components(); Assert (mask.size() == ncomp, - ExcDimensionMismatch(mask.size(), ncomp)); + ExcDimensionMismatch(mask.size(), ncomp)); for (unsigned int i=0; i couplings(ncomp, ncomp); for (unsigned int i=0;i void 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 @@ -257,7 +257,7 @@ namespace MGTools apply_boundary_values (const std::set &boundary_dofs, SparseMatrix& matrix, const bool preserve_symmetry, - const bool ignore_zeros = false); + const bool ignore_zeros = false); template void diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index a72d614b2e..a10b59bcf2 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -57,17 +57,17 @@ namespace DoFTools template void make_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, - const types::subdomain_id subdomain_id) + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, + const types::subdomain_id subdomain_id) { const unsigned int n_dofs = dof.n_dofs(); Assert (sparsity.n_rows() == n_dofs, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); Assert (sparsity.n_cols() == n_dofs, - ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); + ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); // If we have a distributed::Triangulation only // allow locally_owned subdomain. Not setting a @@ -80,13 +80,13 @@ namespace DoFTools || (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; 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 @@ -95,23 +95,23 @@ namespace DoFTools // processor. Otherwise, just continue. for (; cell!=endc; ++cell) if (((subdomain_id == types::invalid_subdomain_id) - || - (subdomain_id == cell->subdomain_id())) - && + || + (subdomain_id == cell->subdomain_id())) + && cell->is_locally_owned()) - { - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - dofs_on_this_cell.resize (dofs_per_cell); - cell->get_dof_indices (dofs_on_this_cell); + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + dofs_on_this_cell.resize (dofs_per_cell); + cell->get_dof_indices (dofs_on_this_cell); // make sparsity pattern for this // cell. if no constraints pattern was // given, then the following call acts // as if simply no constraints existed - constraints.add_entries_local_to_global (dofs_on_this_cell, - sparsity, - keep_constrained_dofs); - } + constraints.add_entries_local_to_global (dofs_on_this_cell, + sparsity, + keep_constrained_dofs); + } } @@ -119,22 +119,22 @@ namespace DoFTools template void make_sparsity_pattern (const DH &dof, - const Table<2,Coupling> &couplings, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, - const types::subdomain_id subdomain_id) + const Table<2,Coupling> &couplings, + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, + const types::subdomain_id subdomain_id) { const unsigned int n_dofs = dof.n_dofs(); Assert (sparsity.n_rows() == n_dofs, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); Assert (sparsity.n_cols() == n_dofs, - ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); + ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); Assert (couplings.n_rows() == dof.get_fe().n_components(), - ExcDimensionMismatch(couplings.n_rows(), dof.get_fe().n_components())); + ExcDimensionMismatch(couplings.n_rows(), dof.get_fe().n_components())); Assert (couplings.n_cols() == dof.get_fe().n_components(), - ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components())); + ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components())); // If we have a distributed::Triangulation only // allow locally_owned subdomain. Not setting a @@ -147,8 +147,8 @@ namespace DoFTools || (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.")); const hp::FECollection fe_collection (dof.get_fe()); @@ -169,51 +169,51 @@ namespace DoFTools bool need_dof_mask = false; for (unsigned int i=0; i 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 @@ -222,27 +222,27 @@ namespace DoFTools // processor. Otherwise, just continue. for (; cell!=endc; ++cell) if (((subdomain_id == types::invalid_subdomain_id) - || - (subdomain_id == cell->subdomain_id())) - && + || + (subdomain_id == cell->subdomain_id())) + && cell->is_locally_owned()) - { - const unsigned int fe_index = cell->active_fe_index(); - const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell; + { + const unsigned int fe_index = cell->active_fe_index(); + const unsigned int dofs_per_cell =fe_collection[fe_index].dofs_per_cell; - dofs_on_this_cell.resize (dofs_per_cell); - cell->get_dof_indices (dofs_on_this_cell); + dofs_on_this_cell.resize (dofs_per_cell); + cell->get_dof_indices (dofs_on_this_cell); // make sparsity pattern for this // cell. if no constraints pattern was // given, then the following call acts // as if simply no constraints existed - constraints.add_entries_local_to_global (dofs_on_this_cell, - sparsity, - keep_constrained_dofs, - dof_mask[fe_index]); - } + constraints.add_entries_local_to_global (dofs_on_this_cell, + sparsity, + keep_constrained_dofs, + dof_mask[fe_index]); + } } @@ -258,13 +258,13 @@ namespace DoFTools const unsigned int n_dofs_col = dof_col.n_dofs(); Assert (sparsity.n_rows() == n_dofs_row, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs_row)); + ExcDimensionMismatch (sparsity.n_rows(), n_dofs_row)); Assert (sparsity.n_cols() == n_dofs_col, - ExcDimensionMismatch (sparsity.n_cols(), n_dofs_col)); + ExcDimensionMismatch (sparsity.n_cols(), n_dofs_col)); const std::list > + typename DH::cell_iterator> > cell_list = GridTools::get_finest_common_cells (dof_row, dof_col); @@ -276,74 +276,74 @@ namespace DoFTools 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 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; ifirst; + 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; i (dof, - boundary_indicators, - dof_to_boundary_mapping, - sparsity); - return; + typename DH::FunctionMap boundary_indicators; + boundary_indicators[0] = 0; + boundary_indicators[1] = 0; + make_boundary_sparsity_pattern (dof, + boundary_indicators, + dof_to_boundary_mapping, + sparsity); + return; } const unsigned int n_dofs = dof.n_dofs(); @@ -380,13 +380,13 @@ namespace DoFTools #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); }; #endif @@ -404,22 +404,22 @@ namespace DoFTools // and so every boundary line is // also part of a boundary face. 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 (cell->at_boundary(f)) - { - 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 (cell->at_boundary(f)) + { + 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; iat_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); // next get boundary mapped dof // indices of boundary dofs - for (unsigned int i=0; ivertex_dof_index(direction,i)]; - - for (unsigned int i=0; ivertex_dof_index(direction,i)]; + + 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); }; #endif 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 void make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, const types::subdomain_id subdomain_id) { const unsigned int n_dofs = dof.n_dofs(); @@ -539,15 +539,15 @@ namespace DoFTools || (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_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(), - endc = dof.end(); + endc = dof.end(); // TODO: in an old implementation, we used // user flags before to tag faces that were @@ -563,57 +563,57 @@ namespace DoFTools // processor. Otherwise, just continue. for (; cell!=endc; ++cell) if (((subdomain_id == types::invalid_subdomain_id) - || - (subdomain_id == cell->subdomain_id())) + || + (subdomain_id == cell->subdomain_id())) && cell->is_locally_owned()) - { - const unsigned int n_dofs_on_this_cell = cell->get_fe().dofs_per_cell; - dofs_on_this_cell.resize (n_dofs_on_this_cell); - cell->get_dof_indices (dofs_on_this_cell); + { + const unsigned int n_dofs_on_this_cell = cell->get_fe().dofs_per_cell; + dofs_on_this_cell.resize (n_dofs_on_this_cell); + cell->get_dof_indices (dofs_on_this_cell); // make sparsity pattern for this // cell. if no constraints pattern was // given, then the following call acts // as if simply no constraints existed - constraints.add_entries_local_to_global (dofs_on_this_cell, - sparsity, - keep_constrained_dofs); - - for (unsigned int face = 0; - face < GeometryInfo::faces_per_cell; - ++face) - { - typename DH::face_iterator cell_face = cell->face(face); - if (! cell->at_boundary(face) ) - { - typename DH::cell_iterator neighbor = cell->neighbor(face); - - if (cell_face->has_children()) - { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->number_of_children(); - ++sub_nr) - { - const typename DH::cell_iterator - sub_neighbor - = cell->neighbor_child_on_subface (face, sub_nr); - - const unsigned int n_dofs_on_neighbor - = sub_neighbor->get_fe().dofs_per_cell; - dofs_on_other_cell.resize (n_dofs_on_neighbor); - sub_neighbor->get_dof_indices (dofs_on_other_cell); - - constraints.add_entries_local_to_global - (dofs_on_this_cell, dofs_on_other_cell, - sparsity, keep_constrained_dofs); - constraints.add_entries_local_to_global - (dofs_on_other_cell, dofs_on_this_cell, - sparsity, keep_constrained_dofs); - } - } - else - { + constraints.add_entries_local_to_global (dofs_on_this_cell, + sparsity, + keep_constrained_dofs); + + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; + ++face) + { + typename DH::face_iterator cell_face = cell->face(face); + if (! cell->at_boundary(face) ) + { + typename DH::cell_iterator neighbor = cell->neighbor(face); + + if (cell_face->has_children()) + { + for (unsigned int sub_nr = 0; + sub_nr != cell_face->number_of_children(); + ++sub_nr) + { + const typename DH::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); + + const unsigned int n_dofs_on_neighbor + = sub_neighbor->get_fe().dofs_per_cell; + dofs_on_other_cell.resize (n_dofs_on_neighbor); + sub_neighbor->get_dof_indices (dofs_on_other_cell); + + constraints.add_entries_local_to_global + (dofs_on_this_cell, dofs_on_other_cell, + sparsity, keep_constrained_dofs); + constraints.add_entries_local_to_global + (dofs_on_other_cell, dofs_on_this_cell, + sparsity, keep_constrained_dofs); + } + } + else + { // Refinement edges are taken care of // by coarser cells @@ -621,32 +621,32 @@ namespace DoFTools // the constraints when the neighbor cell is // coarser, but only the current cell is owned // locally! - if (cell->neighbor_is_coarser(face)) - continue; + if (cell->neighbor_is_coarser(face)) + continue; - const unsigned int n_dofs_on_neighbor - = neighbor->get_fe().dofs_per_cell; - dofs_on_other_cell.resize (n_dofs_on_neighbor); + const unsigned int n_dofs_on_neighbor + = neighbor->get_fe().dofs_per_cell; + dofs_on_other_cell.resize (n_dofs_on_neighbor); - neighbor->get_dof_indices (dofs_on_other_cell); + neighbor->get_dof_indices (dofs_on_other_cell); - constraints.add_entries_local_to_global - (dofs_on_this_cell, dofs_on_other_cell, - sparsity, keep_constrained_dofs); + constraints.add_entries_local_to_global + (dofs_on_this_cell, dofs_on_other_cell, + sparsity, keep_constrained_dofs); // only need to add these in case the neighbor // cell is not locally owned - otherwise, we // touch each face twice and hence put the // indices the other way around - if (cell->neighbor(face)->subdomain_id() != - cell->subdomain_id()) - constraints.add_entries_local_to_global - (dofs_on_other_cell, dofs_on_this_cell, - sparsity, keep_constrained_dofs); - } - } - } - } + if (cell->neighbor(face)->subdomain_id() != + cell->subdomain_id()) + constraints.add_entries_local_to_global + (dofs_on_other_cell, dofs_on_this_cell, + sparsity, keep_constrained_dofs); + } + } + } + } } @@ -669,11 +669,11 @@ namespace DoFTools const Table<2,Coupling> &component_couplings) { Assert(component_couplings.n_rows() == fe.n_components(), - ExcDimensionMismatch(component_couplings.n_rows(), - fe.n_components())); + ExcDimensionMismatch(component_couplings.n_rows(), + fe.n_components())); Assert(component_couplings.n_cols() == fe.n_components(), - ExcDimensionMismatch(component_couplings.n_cols(), - fe.n_components())); + ExcDimensionMismatch(component_couplings.n_cols(), + fe.n_components())); const unsigned int n_dofs = fe.dofs_per_cell; @@ -681,34 +681,34 @@ namespace DoFTools for (unsigned int i=0; i > return_value (fe.size()); for (unsigned int i=0; i 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); - 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 (); + } + } + } } } @@ -976,9 +976,9 @@ namespace DoFTools template void make_flux_sparsity_pattern (const dealii::hp::DoFHandler &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) { // while the implementation above is // quite optimized and caches a lot of @@ -992,174 +992,174 @@ namespace DoFTools // 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)); - 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 (); + } + } + } + } } } @@ -1171,9 +1171,9 @@ namespace DoFTools 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) { // do the error checking and frame code // here, and then pass on to more @@ -1183,17 +1183,17 @@ namespace DoFTools const unsigned int n_comp = dof.get_fe().n_components(); Assert (sparsity.n_rows() == n_dofs, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); Assert (sparsity.n_cols() == n_dofs, - ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); + ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); Assert (int_mask.n_rows() == n_comp, - ExcDimensionMismatch (int_mask.n_rows(), n_comp)); + ExcDimensionMismatch (int_mask.n_rows(), n_comp)); Assert (int_mask.n_cols() == n_comp, - ExcDimensionMismatch (int_mask.n_cols(), n_comp)); + ExcDimensionMismatch (int_mask.n_cols(), n_comp)); Assert (flux_mask.n_rows() == n_comp, - ExcDimensionMismatch (flux_mask.n_rows(), n_comp)); + ExcDimensionMismatch (flux_mask.n_rows(), n_comp)); Assert (flux_mask.n_cols() == n_comp, - ExcDimensionMismatch (flux_mask.n_cols(), n_comp)); + ExcDimensionMismatch (flux_mask.n_cols(), n_comp)); // Clear user flags because we will @@ -1209,7 +1209,7 @@ namespace DoFTools const_cast &>(dof.get_tria()).clear_user_flags (); internal::make_flux_sparsity_pattern (dof, sparsity, - int_mask, flux_mask); + int_mask, flux_mask); // finally restore the user flags const_cast &>(dof.get_tria()).load_user_flags(user_flags); @@ -1228,14 +1228,14 @@ namespace DoFTools inline bool check_master_dof_list (const FullMatrix &face_interpolation_matrix, - 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 p(N); - for (unsigned int i=0; i p(N); + for (unsigned int i=0; i max) - { - max = std::fabs(tmp(i,j)); - r = i; - } - } + double max = std::fabs(tmp(j,j)); + unsigned int r = j; + for (unsigned int i=j+1; i max) + { + max = std::fabs(tmp(i,j)); + r = i; + } + } // check whether the // pivot is too small. if // that is the case, then @@ -1295,44 +1295,44 @@ namespace DoFTools // and we shouldn't use // this set of master // dofs - if (max < 1.e-12*typical_diagonal_element) - return false; + if (max < 1.e-12*typical_diagonal_element) + return false; // row interchange - if (r>j) - { - for (unsigned int k=0; kj) + { + for (unsigned int k=0; k void select_master_dofs_for_face_restriction (const FiniteElement &fe1, - const FiniteElement &fe2, - const FullMatrix &face_interpolation_matrix, - std::vector &master_dof_mask) + const FiniteElement &fe2, + const FullMatrix &face_interpolation_matrix, + std::vector &master_dof_mask) { - Assert (fe1.dofs_per_face >= fe2.dofs_per_face, - ExcInternalError()); - AssertDimension (master_dof_mask.size(), fe1.dofs_per_face); - - Assert (fe2.dofs_per_vertex <= fe1.dofs_per_vertex, - ExcInternalError()); - Assert (fe2.dofs_per_line <= fe1.dofs_per_line, - ExcInternalError()); - Assert ((dim < 3) - || - (fe2.dofs_per_quad <= fe1.dofs_per_quad), - ExcInternalError()); + Assert (fe1.dofs_per_face >= fe2.dofs_per_face, + ExcInternalError()); + AssertDimension (master_dof_mask.size(), fe1.dofs_per_face); + + Assert (fe2.dofs_per_vertex <= fe1.dofs_per_vertex, + ExcInternalError()); + Assert (fe2.dofs_per_line <= fe1.dofs_per_line, + ExcInternalError()); + Assert ((dim < 3) + || + (fe2.dofs_per_quad <= fe1.dofs_per_quad), + ExcInternalError()); // the idea here is to designate as // many DoFs in fe1 per object @@ -1431,16 +1431,16 @@ namespace DoFTools // the other one. this latter case // shows up when running // hp/hp_constraints_q_system_06 - 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) + { // make sure that we // were able to find // a set of master @@ -1448,20 +1448,20 @@ namespace DoFTools // code down below // didn't just reject // all our efforts - Assert (i < fe1.dofs_per_vertex, - ExcInternalError()); + Assert (i < fe1.dofs_per_vertex, + ExcInternalError()); // tentatively push // this vertex dof - master_dof_list.push_back (index+i); + master_dof_list.push_back (index+i); // then see what // happens. if it // succeeds, fine - if (check_master_dof_list (face_interpolation_matrix, - master_dof_list) - == true) - ++dofs_added; - else + if (check_master_dof_list (face_interpolation_matrix, + master_dof_list) + == true) + ++dofs_added; + else // well, it // didn't. simply // pop that dof @@ -1469,74 +1469,74 @@ namespace DoFTools // again and try // with the next // dof - master_dof_list.pop_back (); + master_dof_list.pop_back (); // forward counter by // one - ++i; - } - index += fe1.dofs_per_vertex; - } - - for (int l=0; - l(GeometryInfo::lines_per_face); - ++l) - { + ++i; + } + index += fe1.dofs_per_vertex; + } + + for (int l=0; + l(GeometryInfo::lines_per_face); + ++l) + { // same algorithm as above - unsigned int dofs_added = 0; - unsigned int i = 0; - while (dofs_added < fe2.dofs_per_line) - { - Assert (i < fe1.dofs_per_line, - ExcInternalError()); - - master_dof_list.push_back (index+i); - if (check_master_dof_list (face_interpolation_matrix, - master_dof_list) - == true) - ++dofs_added; - else - master_dof_list.pop_back (); - - ++i; - } - index += fe1.dofs_per_line; - } - - for (int q=0; - q(GeometryInfo::quads_per_face); - ++q) - { + unsigned int dofs_added = 0; + unsigned int i = 0; + while (dofs_added < fe2.dofs_per_line) + { + Assert (i < fe1.dofs_per_line, + ExcInternalError()); + + master_dof_list.push_back (index+i); + if (check_master_dof_list (face_interpolation_matrix, + master_dof_list) + == true) + ++dofs_added; + else + master_dof_list.pop_back (); + + ++i; + } + index += fe1.dofs_per_line; + } + + for (int q=0; + q(GeometryInfo::quads_per_face); + ++q) + { // same algorithm as above - unsigned int dofs_added = 0; - unsigned int i = 0; - while (dofs_added < fe2.dofs_per_quad) - { - Assert (i < fe1.dofs_per_quad, - ExcInternalError()); - - master_dof_list.push_back (index+i); - if (check_master_dof_list (face_interpolation_matrix, - master_dof_list) - == true) - ++dofs_added; - else - master_dof_list.pop_back (); - - ++i; - } - index += fe1.dofs_per_quad; - } - - AssertDimension (index, fe1.dofs_per_face); - AssertDimension (master_dof_list.size(), fe2.dofs_per_face); + unsigned int dofs_added = 0; + unsigned int i = 0; + while (dofs_added < fe2.dofs_per_quad) + { + Assert (i < fe1.dofs_per_quad, + ExcInternalError()); + + master_dof_list.push_back (index+i); + if (check_master_dof_list (face_interpolation_matrix, + master_dof_list) + == true) + ++dofs_added; + else + master_dof_list.pop_back (); + + ++i; + } + index += fe1.dofs_per_quad; + } + + AssertDimension (index, fe1.dofs_per_face); + AssertDimension (master_dof_list.size(), fe2.dofs_per_face); // 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; } @@ -1550,19 +1550,19 @@ namespace DoFTools template void ensure_existence_of_master_dof_mask (const FiniteElement &fe1, - const FiniteElement &fe2, - const FullMatrix &face_interpolation_matrix, - std_cxx1x::shared_ptr > &master_dof_mask) + const FiniteElement &fe2, + const FullMatrix &face_interpolation_matrix, + std_cxx1x::shared_ptr > &master_dof_mask) { - if (master_dof_mask == std_cxx1x::shared_ptr >()) - { - master_dof_mask = std_cxx1x::shared_ptr > - (new std::vector (fe1.dofs_per_face)); - select_master_dofs_for_face_restriction (fe1, - fe2, - face_interpolation_matrix, - *master_dof_mask); - } + if (master_dof_mask == std_cxx1x::shared_ptr >()) + { + master_dof_mask = std_cxx1x::shared_ptr > + (new std::vector (fe1.dofs_per_face)); + select_master_dofs_for_face_restriction (fe1, + fe2, + face_interpolation_matrix, + *master_dof_mask); + } } @@ -1578,17 +1578,17 @@ namespace DoFTools template void ensure_existence_of_face_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - std_cxx1x::shared_ptr > &matrix) + const FiniteElement &fe2, + std_cxx1x::shared_ptr > &matrix) { - if (matrix == std_cxx1x::shared_ptr >()) - { - matrix = std_cxx1x::shared_ptr > - (new FullMatrix (fe2.dofs_per_face, - fe1.dofs_per_face)); - fe1.get_face_interpolation_matrix (fe2, - *matrix); - } + if (matrix == std_cxx1x::shared_ptr >()) + { + matrix = std_cxx1x::shared_ptr > + (new FullMatrix (fe2.dofs_per_face, + fe1.dofs_per_face)); + fe1.get_face_interpolation_matrix (fe2, + *matrix); + } } @@ -1600,19 +1600,19 @@ namespace DoFTools template void ensure_existence_of_subface_matrix (const FiniteElement &fe1, - const FiniteElement &fe2, - const unsigned int subface, - std_cxx1x::shared_ptr > &matrix) + const FiniteElement &fe2, + const unsigned int subface, + std_cxx1x::shared_ptr > &matrix) { - if (matrix == std_cxx1x::shared_ptr >()) - { - matrix = std_cxx1x::shared_ptr > - (new FullMatrix (fe2.dofs_per_face, - fe1.dofs_per_face)); - fe1.get_subface_interpolation_matrix (fe2, - subface, - *matrix); - } + if (matrix == std_cxx1x::shared_ptr >()) + { + matrix = std_cxx1x::shared_ptr > + (new FullMatrix (fe2.dofs_per_face, + fe1.dofs_per_face)); + fe1.get_subface_interpolation_matrix (fe2, + subface, + *matrix); + } } @@ -1626,57 +1626,57 @@ namespace DoFTools */ void ensure_existence_of_split_face_matrix (const FullMatrix &face_interpolation_matrix, - const std::vector &master_dof_mask, - std_cxx1x::shared_ptr,FullMatrix > > &split_matrix) + const std::vector &master_dof_mask, + std_cxx1x::shared_ptr,FullMatrix > > &split_matrix) { - AssertDimension (master_dof_mask.size(), face_interpolation_matrix.m()); - Assert (std::count (master_dof_mask.begin(), master_dof_mask.end(), true) == - static_cast(face_interpolation_matrix.n()), - ExcInternalError()); + AssertDimension (master_dof_mask.size(), face_interpolation_matrix.m()); + Assert (std::count (master_dof_mask.begin(), master_dof_mask.end(), true) == + static_cast(face_interpolation_matrix.n()), + ExcInternalError()); - if (split_matrix == - std_cxx1x::shared_ptr,FullMatrix > >()) - { - split_matrix - = std_cxx1x::shared_ptr,FullMatrix > > - (new std::pair,FullMatrix >()); + if (split_matrix == + std_cxx1x::shared_ptr,FullMatrix > >()) + { + split_matrix + = std_cxx1x::shared_ptr,FullMatrix > > + (new std::pair,FullMatrix >()); - const unsigned int n_master_dofs = face_interpolation_matrix.n(); - const unsigned int n_dofs = face_interpolation_matrix.m(); + const unsigned int n_master_dofs = face_interpolation_matrix.n(); + const unsigned int n_dofs = face_interpolation_matrix.m(); - Assert (n_master_dofs <= n_dofs, ExcInternalError()); + Assert (n_master_dofs <= n_dofs, ExcInternalError()); // copy and invert the master // component, copy the slave // component - split_matrix->first.reinit (n_master_dofs, n_master_dofs); - split_matrix->second.reinit (n_dofs-n_master_dofs, n_master_dofs); - - unsigned int nth_master_dof = 0, - nth_slave_dof = 0; - - for (unsigned int i=0; ifirst(nth_master_dof,j) - = face_interpolation_matrix(i,j); - ++nth_master_dof; - } - else - { - for (unsigned int j=0; jsecond(nth_slave_dof,j) - = face_interpolation_matrix(i,j); - ++nth_slave_dof; - } - - AssertDimension (nth_master_dof, n_master_dofs); - AssertDimension (nth_slave_dof, n_dofs-n_master_dofs); + split_matrix->first.reinit (n_master_dofs, n_master_dofs); + split_matrix->second.reinit (n_dofs-n_master_dofs, n_master_dofs); + + unsigned int nth_master_dof = 0, + nth_slave_dof = 0; + + for (unsigned int i=0; ifirst(nth_master_dof,j) + = face_interpolation_matrix(i,j); + ++nth_master_dof; + } + else + { + for (unsigned int j=0; jsecond(nth_slave_dof,j) + = face_interpolation_matrix(i,j); + ++nth_slave_dof; + } + + AssertDimension (nth_master_dof, n_master_dofs); + AssertDimension (nth_slave_dof, n_dofs-n_master_dofs); //TODO[WB]: We should make sure very small entries are removed after inversion - split_matrix->first.gauss_jordan (); - } + split_matrix->first.gauss_jordan (); + } } @@ -1688,21 +1688,21 @@ namespace DoFTools template struct DoFHandlerSupportsDifferentFEs { - static const bool value = true; + static const bool value = true; }; template struct DoFHandlerSupportsDifferentFEs< dealii::DoFHandler > { - static const bool value = false; + static const bool value = false; }; template struct DoFHandlerSupportsDifferentFEs< dealii::MGDoFHandler > { - static const bool value = false; + static const bool value = false; }; @@ -1719,7 +1719,7 @@ namespace DoFTools unsigned int n_finite_elements (const dealii::hp::DoFHandler &dof_handler) { - return dof_handler.get_fe().size(); + return dof_handler.get_fe().size(); } @@ -1727,7 +1727,7 @@ namespace DoFTools unsigned int n_finite_elements (const DH &) { - return 1; + return 1; } @@ -1744,44 +1744,44 @@ namespace DoFTools unsigned int get_most_dominating_subface_fe_index (const face_iterator &face) { - const unsigned int dim - = face_iterator::AccessorType::dimension; - const unsigned int spacedim - = face_iterator::AccessorType::space_dimension; - - unsigned int dominating_subface_no = 0; - for (; dominating_subface_non_children(); - ++dominating_subface_no) - { + const unsigned int dim + = face_iterator::AccessorType::dimension; + const unsigned int spacedim + = face_iterator::AccessorType::space_dimension; + + unsigned int dominating_subface_no = 0; + for (; dominating_subface_non_children(); + ++dominating_subface_no) + { // each of the subfaces // can have only a single // fe_index associated // with them, since there // is no cell on the // other side - Assert (face->child(dominating_subface_no) - ->n_active_fe_indices() - == 1, - ExcInternalError()); - - const FiniteElement & - this_subface_fe = (face->child(dominating_subface_no) - ->get_fe (face->child(dominating_subface_no) - ->nth_active_fe_index(0))); - - FiniteElementDomination::Domination - domination = FiniteElementDomination::either_element_can_dominate; - for (unsigned int sf=0; sfn_children(); ++sf) - if (sf != dominating_subface_no) - { - const FiniteElement & - that_subface_fe = (face->child(sf) - ->get_fe (face->child(sf) - ->nth_active_fe_index(0))); - - domination = domination & - this_subface_fe.compare_for_face_domination(that_subface_fe); - } + Assert (face->child(dominating_subface_no) + ->n_active_fe_indices() + == 1, + ExcInternalError()); + + const FiniteElement & + this_subface_fe = (face->child(dominating_subface_no) + ->get_fe (face->child(dominating_subface_no) + ->nth_active_fe_index(0))); + + FiniteElementDomination::Domination + domination = FiniteElementDomination::either_element_can_dominate; + for (unsigned int sf=0; sfn_children(); ++sf) + if (sf != dominating_subface_no) + { + const FiniteElement & + that_subface_fe = (face->child(sf) + ->get_fe (face->child(sf) + ->nth_active_fe_index(0))); + + domination = domination & + this_subface_fe.compare_for_face_domination(that_subface_fe); + } // see if the element // on this subface is @@ -1789,22 +1789,22 @@ namespace DoFTools // the ones on all // other subfaces, // and if so take it - if ((domination == FiniteElementDomination::this_element_dominates) - || - (domination == FiniteElementDomination::either_element_can_dominate)) - break; - } + if ((domination == FiniteElementDomination::this_element_dominates) + || + (domination == FiniteElementDomination::either_element_can_dominate)) + break; + } // check that we have // found one such subface - Assert (dominating_subface_no < face->n_children(), - ExcNotImplemented()); + Assert (dominating_subface_no < face->n_children(), + ExcNotImplemented()); // return the finite element // index used on it. note // that only a single fe can // be active on such subfaces - return face->child (dominating_subface_no)->nth_active_fe_index(0); + return face->child (dominating_subface_no)->nth_active_fe_index(0); } @@ -1827,60 +1827,60 @@ namespace DoFTools */ void 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 // in parallel distributed // mode - for (unsigned int row=0; row!=n_slave_dofs; ++row) - Assert (slave_dofs[row] != numbers::invalid_unsigned_int, - ExcInternalError()); - for (unsigned int col=0; col!=n_master_dofs; ++col) - Assert (master_dofs[col] != numbers::invalid_unsigned_int, - ExcInternalError()); + for (unsigned int row=0; row!=n_slave_dofs; ++row) + Assert (slave_dofs[row] != numbers::invalid_unsigned_int, + ExcInternalError()); + for (unsigned int col=0; col!=n_master_dofs; ++col) + Assert (master_dofs[col] != numbers::invalid_unsigned_int, + ExcInternalError()); - for (unsigned int row=0; row!=n_slave_dofs; ++row) - if (constraints.is_constrained (slave_dofs[row]) == false) - { - bool constraint_already_satisfied = false; + for (unsigned int row=0; row!=n_slave_dofs; ++row) + if (constraints.is_constrained (slave_dofs[row]) == false) + { + bool constraint_already_satisfied = false; // Check if we have an identity // constraint, which is already // satisfied by unification of // the corresponding global dof // indices - for (unsigned int i=0; i= 1e-14*abs_sum)) - constraints.add_entry (slave_dofs[row], - master_dofs[i], - face_constraints (row,i)); - constraints.set_inhomogeneity (slave_dofs[row], 0.); - } - } + constraints.add_line (slave_dofs[row]); + for (unsigned int i=0; i= 1e-14*abs_sum)) + constraints.add_entry (slave_dofs[row], + master_dofs[i], + face_constraints (row,i)); + constraints.set_inhomogeneity (slave_dofs[row], 0.); + } + } } } @@ -1916,7 +1916,7 @@ namespace DoFTools void make_hp_hanging_node_constraints (const dealii::DoFHandler<1> &, - ConstraintMatrix &) + ConstraintMatrix &) { // nothing to do for regular // dof handlers in 1d @@ -1926,8 +1926,8 @@ namespace DoFTools void make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1> &, - ConstraintMatrix &, - dealii::internal::int2type<1>) + ConstraintMatrix &, + dealii::internal::int2type<1>) { // nothing to do for regular // dof handlers in 1d @@ -1936,7 +1936,7 @@ namespace DoFTools void make_hp_hanging_node_constraints (const dealii::MGDoFHandler<1> &, - ConstraintMatrix &) + ConstraintMatrix &) { // nothing to do for regular // dof handlers in 1d @@ -1946,8 +1946,8 @@ namespace DoFTools void make_oldstyle_hanging_node_constraints (const dealii::MGDoFHandler<1> &, - ConstraintMatrix &, - dealii::internal::int2type<1>) + ConstraintMatrix &, + dealii::internal::int2type<1>) { // nothing to do for regular // dof handlers in 1d @@ -1956,7 +1956,7 @@ namespace DoFTools void make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, - ConstraintMatrix &/*constraints*/) + ConstraintMatrix &/*constraints*/) { // we may have to compute // constraints for @@ -1969,8 +1969,8 @@ namespace DoFTools void make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, - ConstraintMatrix &/*constraints*/, - dealii::internal::int2type<1>) + ConstraintMatrix &/*constraints*/, + dealii::internal::int2type<1>) { // we may have to compute // constraints for @@ -1982,7 +1982,7 @@ namespace DoFTools void make_hp_hanging_node_constraints (const dealii::DoFHandler<1,2> &, - ConstraintMatrix &) + ConstraintMatrix &) { // nothing to do for regular // dof handlers in 1d @@ -1992,8 +1992,8 @@ namespace DoFTools void make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,2> &, - ConstraintMatrix &, - dealii::internal::int2type<1>) + ConstraintMatrix &, + dealii::internal::int2type<1>) { // nothing to do for regular // dof handlers in 1d @@ -2058,8 +2058,8 @@ namespace DoFTools template void make_oldstyle_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints, - dealii::internal::int2type<2>) + ConstraintMatrix &constraints, + dealii::internal::int2type<2>) { const unsigned int dim = 2; @@ -2085,16 +2085,16 @@ namespace DoFTools // refined, we can only visit each // face with hanging nodes once typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + endc = dof_handler.end(); for (; cell!=endc; ++cell) // artificial cells can at // best neighbor ghost cells, // but we're not interested // in these interfaces - if (!cell->is_artificial ()) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->has_children()) - { + if (!cell->is_artificial ()) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->has_children()) + { // in any case, faces // can have at most two // active fe indices, @@ -2107,84 +2107,84 @@ namespace DoFTools // children can have // only one as // well. check this - Assert (cell->face(face)->n_active_fe_indices() == 1, - ExcInternalError()); - Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) - == true, - ExcInternalError()); - for (unsigned int c=0; cface(face)->n_children(); ++c) - if (!cell->neighbor_child_on_subface(face,c)->is_artificial()) - Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, - ExcInternalError()); + Assert (cell->face(face)->n_active_fe_indices() == 1, + ExcInternalError()); + Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) + == true, + ExcInternalError()); + for (unsigned int c=0; cface(face)->n_children(); ++c) + if (!cell->neighbor_child_on_subface(face,c)->is_artificial()) + Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, + ExcInternalError()); // right now, all that // is implemented is // the case that both // sides use the same // fe - for (unsigned int c=0; cface(face)->n_children(); ++c) - if (!cell->neighbor_child_on_subface(face,c)->is_artificial()) - Assert (cell->face(face)->child(c) - ->fe_index_is_active(cell->active_fe_index()) == true, - ExcNotImplemented()); + for (unsigned int c=0; cface(face)->n_children(); ++c) + if (!cell->neighbor_child_on_subface(face,c)->is_artificial()) + Assert (cell->face(face)->child(c) + ->fe_index_is_active(cell->active_fe_index()) == true, + ExcNotImplemented()); // ok, start up the work - const FiniteElement &fe = cell->get_fe(); - const unsigned int fe_index = cell->active_fe_index(); + const FiniteElement &fe = cell->get_fe(); + const unsigned int fe_index = cell->active_fe_index(); - const unsigned int - n_dofs_on_mother = 2*fe.dofs_per_vertex + fe.dofs_per_line, - n_dofs_on_children = fe.dofs_per_vertex + 2*fe.dofs_per_line; + const unsigned int + n_dofs_on_mother = 2*fe.dofs_per_vertex + fe.dofs_per_line, + n_dofs_on_children = fe.dofs_per_vertex + 2*fe.dofs_per_line; - dofs_on_mother.resize (n_dofs_on_mother); - dofs_on_children.resize (n_dofs_on_children); + dofs_on_mother.resize (n_dofs_on_mother); + dofs_on_children.resize (n_dofs_on_children); - Assert(n_dofs_on_mother == fe.constraints().n(), - ExcDimensionMismatch(n_dofs_on_mother, - fe.constraints().n())); - Assert(n_dofs_on_children == fe.constraints().m(), - ExcDimensionMismatch(n_dofs_on_children, - fe.constraints().m())); + Assert(n_dofs_on_mother == fe.constraints().n(), + ExcDimensionMismatch(n_dofs_on_mother, + fe.constraints().n())); + Assert(n_dofs_on_children == fe.constraints().m(), + ExcDimensionMismatch(n_dofs_on_children, + fe.constraints().m())); - const typename DH::line_iterator this_face = cell->face(face); + const typename DH::line_iterator this_face = cell->face(face); // fill the dofs indices. Use same // enumeration scheme as in // @p{FiniteElement::constraints()} - unsigned int next_index = 0; - for (unsigned int vertex=0; vertex<2; ++vertex) - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof, - fe_index); - for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) - dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index); - AssertDimension (next_index, dofs_on_mother.size()); - - next_index = 0; - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - dofs_on_children[next_index++] - = this_face->child(0)->vertex_dof_index(1,dof,fe_index); - for (unsigned int child=0; child<2; ++child) - for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) - dofs_on_children[next_index++] - = this_face->child(child)->dof_index(dof, fe_index); - AssertDimension (next_index, dofs_on_children.size()); + unsigned int next_index = 0; + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof, + fe_index); + for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) + dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index); + AssertDimension (next_index, dofs_on_mother.size()); + + next_index = 0; + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + dofs_on_children[next_index++] + = this_face->child(0)->vertex_dof_index(1,dof,fe_index); + for (unsigned int child=0; child<2; ++child) + for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) + dofs_on_children[next_index++] + = this_face->child(child)->dof_index(dof, fe_index); + AssertDimension (next_index, dofs_on_children.size()); // for each row in the constraint // matrix for this line: - for (unsigned int row=0; row!=dofs_on_children.size(); ++row) - { - constraints.add_line (dofs_on_children[row]); - for (unsigned int i=0; i!=dofs_on_mother.size(); ++i) - constraints.add_entry (dofs_on_children[row], - dofs_on_mother[i], - fe.constraints()(row,i)); - - constraints.set_inhomogeneity (dofs_on_children[row], 0.); - } - } - else - { + for (unsigned int row=0; row!=dofs_on_children.size(); ++row) + { + constraints.add_line (dofs_on_children[row]); + for (unsigned int i=0; i!=dofs_on_mother.size(); ++i) + constraints.add_entry (dofs_on_children[row], + dofs_on_mother[i], + fe.constraints()(row,i)); + + constraints.set_inhomogeneity (dofs_on_children[row], 0.); + } + } + else + { // this face has no // children, but it // could still be @@ -2197,16 +2197,16 @@ namespace DoFTools // case that the // neighbor is an // artificial cell - if (!cell->at_boundary(face) && - !cell->neighbor(face)->is_artificial()) - { - Assert (cell->face(face)->n_active_fe_indices() == 1, - ExcNotImplemented()); - Assert (cell->face(face) - ->fe_index_is_active(cell->active_fe_index()) == true, - ExcInternalError()); - } - } + if (!cell->at_boundary(face) && + !cell->neighbor(face)->is_artificial()) + { + Assert (cell->face(face)->n_active_fe_indices() == 1, + ExcNotImplemented()); + Assert (cell->face(face) + ->fe_index_is_active(cell->active_fe_index()) == true, + ExcInternalError()); + } + } } @@ -2214,8 +2214,8 @@ namespace DoFTools template void make_oldstyle_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints, - dealii::internal::int2type<3>) + ConstraintMatrix &constraints, + dealii::internal::int2type<3>) { const unsigned int dim = 3; @@ -2239,26 +2239,26 @@ namespace DoFTools // refined, we can only visit each // face with hanging nodes once typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + endc = dof_handler.end(); for (; cell!=endc; ++cell) // artificial cells can at // best neighbor ghost cells, // but we're not interested // in these interfaces - if (!cell->is_artificial ()) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->has_children()) - { + if (!cell->is_artificial ()) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->has_children()) + { // first of all, make sure that // we treat a case which is // possible, i.e. either no dofs // on the face at all or no // anisotropic refinement - if (cell->get_fe().dofs_per_face == 0) - continue; + if (cell->get_fe().dofs_per_face == 0) + continue; - Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, - ExcNotImplemented()); + Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, + ExcNotImplemented()); // in any case, faces // can have at most two @@ -2272,12 +2272,12 @@ namespace DoFTools // children can have // only one as // well. check this - AssertDimension (cell->face(face)->n_active_fe_indices(), 1); - Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) - == true, - ExcInternalError()); - for (unsigned int c=0; cface(face)->n_children(); ++c) - AssertDimension (cell->face(face)->child(c)->n_active_fe_indices(), 1); + AssertDimension (cell->face(face)->n_active_fe_indices(), 1); + Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) + == true, + ExcInternalError()); + for (unsigned int c=0; cface(face)->n_children(); ++c) + AssertDimension (cell->face(face)->child(c)->n_active_fe_indices(), 1); // right now, all that // is implemented is @@ -2289,141 +2289,141 @@ namespace DoFTools // this face and the // children have the // same fe - for (unsigned int c=0; cface(face)->n_children(); ++c) - if (!cell->neighbor_child_on_subface(face,c)->is_artificial()) - { - Assert (cell->face(face)->child(c) - ->fe_index_is_active(cell->active_fe_index()) == true, - ExcNotImplemented()); - for (unsigned int e=0; e<4; ++e) - { - Assert (cell->face(face)->child(c)->line(e) - ->n_active_fe_indices() == 1, - ExcNotImplemented()); - Assert (cell->face(face)->child(c)->line(e) - ->fe_index_is_active(cell->active_fe_index()) == true, - ExcNotImplemented()); - } - } - for (unsigned int e=0; e<4; ++e) - { - Assert (cell->face(face)->line(e) - ->n_active_fe_indices() == 1, - ExcNotImplemented()); - Assert (cell->face(face)->line(e) - ->fe_index_is_active(cell->active_fe_index()) == true, - ExcNotImplemented()); - } + for (unsigned int c=0; cface(face)->n_children(); ++c) + if (!cell->neighbor_child_on_subface(face,c)->is_artificial()) + { + Assert (cell->face(face)->child(c) + ->fe_index_is_active(cell->active_fe_index()) == true, + ExcNotImplemented()); + for (unsigned int e=0; e<4; ++e) + { + Assert (cell->face(face)->child(c)->line(e) + ->n_active_fe_indices() == 1, + ExcNotImplemented()); + Assert (cell->face(face)->child(c)->line(e) + ->fe_index_is_active(cell->active_fe_index()) == true, + ExcNotImplemented()); + } + } + for (unsigned int e=0; e<4; ++e) + { + Assert (cell->face(face)->line(e) + ->n_active_fe_indices() == 1, + ExcNotImplemented()); + Assert (cell->face(face)->line(e) + ->fe_index_is_active(cell->active_fe_index()) == true, + ExcNotImplemented()); + } // ok, start up the work - const FiniteElement &fe = cell->get_fe(); - const unsigned int fe_index = cell->active_fe_index(); + const FiniteElement &fe = cell->get_fe(); + const unsigned int fe_index = cell->active_fe_index(); - const unsigned int n_dofs_on_mother = fe.dofs_per_face; - const unsigned int n_dofs_on_children = (5*fe.dofs_per_vertex+ - 12*fe.dofs_per_line+ - 4*fe.dofs_per_quad); + const unsigned int n_dofs_on_mother = fe.dofs_per_face; + const unsigned int n_dofs_on_children = (5*fe.dofs_per_vertex+ + 12*fe.dofs_per_line+ + 4*fe.dofs_per_quad); //TODO[TL]: think about this and the following in case of anisotropic refinement - dofs_on_mother.resize (n_dofs_on_mother); - dofs_on_children.resize (n_dofs_on_children); + dofs_on_mother.resize (n_dofs_on_mother); + dofs_on_children.resize (n_dofs_on_children); - Assert(n_dofs_on_mother == fe.constraints().n(), - ExcDimensionMismatch(n_dofs_on_mother, - fe.constraints().n())); - Assert(n_dofs_on_children == fe.constraints().m(), - ExcDimensionMismatch(n_dofs_on_children, - fe.constraints().m())); + Assert(n_dofs_on_mother == fe.constraints().n(), + ExcDimensionMismatch(n_dofs_on_mother, + fe.constraints().n())); + Assert(n_dofs_on_children == fe.constraints().m(), + ExcDimensionMismatch(n_dofs_on_children, + fe.constraints().m())); - const typename DH::face_iterator this_face = cell->face(face); + const typename DH::face_iterator this_face = cell->face(face); // fill the dofs indices. Use same // enumeration scheme as in // @p{FiniteElement::constraints()} - unsigned int next_index = 0; - for (unsigned int vertex=0; vertex<4; ++vertex) - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof, - fe_index); - for (unsigned int line=0; line<4; ++line) - for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) - dofs_on_mother[next_index++] - = this_face->line(line)->dof_index(dof, fe_index); - for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) - dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index); - AssertDimension (next_index, dofs_on_mother.size()); - - next_index = 0; + unsigned int next_index = 0; + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof, + fe_index); + for (unsigned int line=0; line<4; ++line) + for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) + dofs_on_mother[next_index++] + = this_face->line(line)->dof_index(dof, fe_index); + for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) + dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index); + AssertDimension (next_index, dofs_on_mother.size()); + + next_index = 0; // assert some consistency // assumptions //TODO[TL]: think about this in case of anisotropic refinement - Assert (dof_handler.get_tria().get_anisotropic_refinement_flag() || - ((this_face->child(0)->vertex_index(3) == - this_face->child(1)->vertex_index(2)) && - (this_face->child(0)->vertex_index(3) == - this_face->child(2)->vertex_index(1)) && - (this_face->child(0)->vertex_index(3) == - this_face->child(3)->vertex_index(0))), - ExcInternalError()); - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - dofs_on_children[next_index++] - = this_face->child(0)->vertex_dof_index(3,dof); + Assert (dof_handler.get_tria().get_anisotropic_refinement_flag() || + ((this_face->child(0)->vertex_index(3) == + this_face->child(1)->vertex_index(2)) && + (this_face->child(0)->vertex_index(3) == + this_face->child(2)->vertex_index(1)) && + (this_face->child(0)->vertex_index(3) == + this_face->child(3)->vertex_index(0))), + ExcInternalError()); + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + dofs_on_children[next_index++] + = this_face->child(0)->vertex_dof_index(3,dof); // dof numbers on the centers of // the lines bounding this face - for (unsigned int line=0; line<4; ++line) - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - dofs_on_children[next_index++] - = this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index); + for (unsigned int line=0; line<4; ++line) + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + dofs_on_children[next_index++] + = this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index); // next the dofs on the lines interior // to the face; the order of these // lines is laid down in the // FiniteElement class documentation - for (unsigned int dof=0; dofchild(0)->line(1)->dof_index(dof, fe_index); - for (unsigned int dof=0; dofchild(2)->line(1)->dof_index(dof, fe_index); - for (unsigned int dof=0; dofchild(0)->line(3)->dof_index(dof, fe_index); - for (unsigned int dof=0; dofchild(1)->line(3)->dof_index(dof, fe_index); + for (unsigned int dof=0; dofchild(0)->line(1)->dof_index(dof, fe_index); + for (unsigned int dof=0; dofchild(2)->line(1)->dof_index(dof, fe_index); + for (unsigned int dof=0; dofchild(0)->line(3)->dof_index(dof, fe_index); + for (unsigned int dof=0; dofchild(1)->line(3)->dof_index(dof, fe_index); // dofs on the bordering lines - for (unsigned int line=0; line<4; ++line) - for (unsigned int child=0; child<2; ++child) - for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) - dofs_on_children[next_index++] - = this_face->line(line)->child(child)->dof_index(dof, fe_index); + for (unsigned int line=0; line<4; ++line) + for (unsigned int child=0; child<2; ++child) + for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) + dofs_on_children[next_index++] + = this_face->line(line)->child(child)->dof_index(dof, fe_index); // finally, for the dofs interior // to the four child faces - for (unsigned int child=0; child<4; ++child) - for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) - dofs_on_children[next_index++] - = this_face->child(child)->dof_index(dof, fe_index); - AssertDimension (next_index, dofs_on_children.size()); + for (unsigned int child=0; child<4; ++child) + for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) + dofs_on_children[next_index++] + = this_face->child(child)->dof_index(dof, fe_index); + AssertDimension (next_index, dofs_on_children.size()); // for each row in the constraint // matrix for this line: - for (unsigned int row=0; row!=dofs_on_children.size(); ++row) - { - constraints.add_line (dofs_on_children[row]); - for (unsigned int i=0; i!=dofs_on_mother.size(); ++i) - constraints.add_entry (dofs_on_children[row], - dofs_on_mother[i], - fe.constraints()(row,i)); - - constraints.set_inhomogeneity(dofs_on_children[row], 0.); - } - } - else - { + for (unsigned int row=0; row!=dofs_on_children.size(); ++row) + { + constraints.add_line (dofs_on_children[row]); + for (unsigned int i=0; i!=dofs_on_mother.size(); ++i) + constraints.add_entry (dofs_on_children[row], + dofs_on_mother[i], + fe.constraints()(row,i)); + + constraints.set_inhomogeneity(dofs_on_children[row], 0.); + } + } + else + { // this face has no // children, but it // could still be @@ -2436,23 +2436,23 @@ namespace DoFTools // case that the // neighbor is an // artificial cell - if (!cell->at_boundary(face) && - !cell->neighbor(face)->is_artificial()) - { - Assert (cell->face(face)->n_active_fe_indices() == 1, - ExcNotImplemented()); - Assert (cell->face(face) - ->fe_index_is_active(cell->active_fe_index()) == true, - ExcInternalError()); - } - } + if (!cell->at_boundary(face) && + !cell->neighbor(face)->is_artificial()) + { + Assert (cell->face(face)->n_active_fe_indices() == 1, + ExcNotImplemented()); + Assert (cell->face(face) + ->fe_index_is_active(cell->active_fe_index()) == true, + ExcInternalError()); + } + } } template void make_hp_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints) + ConstraintMatrix &constraints) { // note: this function is going // to be hard to understand if @@ -2492,12 +2492,12 @@ namespace DoFTools // time they are needed, and // then just reuse them Table<2,std_cxx1x::shared_ptr > > - face_interpolation_matrices (n_finite_elements (dof_handler), - n_finite_elements (dof_handler)); + face_interpolation_matrices (n_finite_elements (dof_handler), + n_finite_elements (dof_handler)); Table<3,std_cxx1x::shared_ptr > > - subface_interpolation_matrices (n_finite_elements (dof_handler), - n_finite_elements (dof_handler), - GeometryInfo::max_children_per_face); + subface_interpolation_matrices (n_finite_elements (dof_handler), + n_finite_elements (dof_handler), + GeometryInfo::max_children_per_face); // similarly have a cache for // the matrices that are split @@ -2509,8 +2509,8 @@ namespace DoFTools // interpolation matrix as // described in the @ref hp_paper "hp paper" Table<2,std_cxx1x::shared_ptr,FullMatrix > > > - split_face_interpolation_matrices (n_finite_elements (dof_handler), - n_finite_elements (dof_handler)); + split_face_interpolation_matrices (n_finite_elements (dof_handler), + n_finite_elements (dof_handler)); // finally, for each pair of finite // elements, have a mask that states @@ -2518,8 +2518,8 @@ namespace DoFTools // the coarse side of a refined face // will act as master dofs. Table<2,std_cxx1x::shared_ptr > > - master_dof_masks (n_finite_elements (dof_handler), - n_finite_elements (dof_handler)); + master_dof_masks (n_finite_elements (dof_handler), + n_finite_elements (dof_handler)); // loop over all faces // @@ -2529,26 +2529,26 @@ namespace DoFTools // refined, we can only visit each // face with hanging nodes once typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + endc = dof_handler.end(); for (; cell!=endc; ++cell) // artificial cells can at // best neighbor ghost cells, // but we're not interested // in these interfaces - if (!cell->is_artificial ()) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->has_children()) - { + if (!cell->is_artificial ()) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->has_children()) + { // first of all, make sure that // we treat a case which is // possible, i.e. either no dofs // on the face at all or no // anisotropic refinement - if (cell->get_fe().dofs_per_face == 0) - continue; + if (cell->get_fe().dofs_per_face == 0) + continue; - Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, - ExcNotImplemented()); + Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, + ExcNotImplemented()); // so now we've found a // face of an active @@ -2569,14 +2569,14 @@ namespace DoFTools // each of the children // can have only one as // well. check this - Assert (cell->face(face)->n_active_fe_indices() == 1, - ExcInternalError()); - Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) - == true, - ExcInternalError()); - for (unsigned int c=0; cface(face)->n_children(); ++c) - Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, - ExcInternalError()); + Assert (cell->face(face)->n_active_fe_indices() == 1, + ExcInternalError()); + Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) + == true, + ExcInternalError()); + for (unsigned int c=0; cface(face)->n_children(); ++c) + Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, + ExcInternalError()); // first find out // whether we can @@ -2595,21 +2595,21 @@ namespace DoFTools // ignore all // interfaces with // artificial cells - FiniteElementDomination::Domination - mother_face_dominates = FiniteElementDomination::either_element_can_dominate; - - if (DoFHandlerSupportsDifferentFEs::value == true) - for (unsigned int c=0; cface(face)->number_of_children(); ++c) - if (!cell->neighbor_child_on_subface (face, c)->is_artificial()) - mother_face_dominates = mother_face_dominates & - (cell->get_fe().compare_for_face_domination - (cell->neighbor_child_on_subface (face, c)->get_fe())); - - switch (mother_face_dominates) - { - case FiniteElementDomination::this_element_dominates: - case FiniteElementDomination::either_element_can_dominate: - { + FiniteElementDomination::Domination + mother_face_dominates = FiniteElementDomination::either_element_can_dominate; + + if (DoFHandlerSupportsDifferentFEs::value == true) + for (unsigned int c=0; cface(face)->number_of_children(); ++c) + if (!cell->neighbor_child_on_subface (face, c)->is_artificial()) + mother_face_dominates = mother_face_dominates & + (cell->get_fe().compare_for_face_domination + (cell->neighbor_child_on_subface (face, c)->get_fe())); + + switch (mother_face_dominates) + { + case FiniteElementDomination::this_element_dominates: + case FiniteElementDomination::either_element_can_dominate: + { // Case 1 (the // simple case // and the only @@ -2634,10 +2634,10 @@ namespace DoFTools // against the // DoFs on the // face itself - master_dofs.resize (cell->get_fe().dofs_per_face); + master_dofs.resize (cell->get_fe().dofs_per_face); - cell->face(face)->get_dof_indices (master_dofs, - cell->active_fe_index ()); + cell->face(face)->get_dof_indices (master_dofs, + cell->active_fe_index ()); // Now create // constraint matrix @@ -2675,21 +2675,21 @@ namespace DoFTools // node is only // detected if we also // look between ghosts - for (unsigned int c=0; cface(face)->n_children(); ++c) - { - if (cell->neighbor_child_on_subface (face, c)->is_artificial() - || - (dim == 2 && cell->neighbor_child_on_subface (face, c)->is_ghost())) - continue; + for (unsigned int c=0; cface(face)->n_children(); ++c) + { + if (cell->neighbor_child_on_subface (face, c)->is_artificial() + || + (dim == 2 && cell->neighbor_child_on_subface (face, c)->is_ghost())) + continue; - const typename DH::active_face_iterator - subface = cell->face(face)->child(c); + const typename DH::active_face_iterator + subface = cell->face(face)->child(c); - Assert (subface->n_active_fe_indices() == 1, - ExcInternalError()); + Assert (subface->n_active_fe_indices() == 1, + ExcInternalError()); - const unsigned int - subface_fe_index = subface->nth_active_fe_index(0); + const unsigned int + subface_fe_index = subface->nth_active_fe_index(0); // we sometime run // into the @@ -2719,22 +2719,22 @@ namespace DoFTools // care of. in that // case, just // continue - if (cell->get_fe().compare_for_face_domination - (subface->get_fe(subface_fe_index)) - == - FiniteElementDomination::no_requirements) - continue; + if (cell->get_fe().compare_for_face_domination + (subface->get_fe(subface_fe_index)) + == + FiniteElementDomination::no_requirements) + continue; // Same procedure as for the // mother cell. Extract the face // DoFs from the cell DoFs. - slave_dofs.resize (subface->get_fe(subface_fe_index) - .dofs_per_face); - subface->get_dof_indices (slave_dofs, subface_fe_index); + slave_dofs.resize (subface->get_fe(subface_fe_index) + .dofs_per_face); + subface->get_dof_indices (slave_dofs, subface_fe_index); - for (unsigned int i=0; iget_fe(), - subface->get_fe(subface_fe_index), - c, - subface_interpolation_matrices - [cell->active_fe_index()][subface_fe_index][c]); + ensure_existence_of_subface_matrix + (cell->get_fe(), + subface->get_fe(subface_fe_index), + c, + subface_interpolation_matrices + [cell->active_fe_index()][subface_fe_index][c]); // Add constraints to global constraint // matrix. - filter_constraints (master_dofs, - slave_dofs, - *(subface_interpolation_matrices - [cell->active_fe_index()][subface_fe_index][c]), - constraints); - } - - break; - } - - case FiniteElementDomination::other_element_dominates: - case FiniteElementDomination::neither_element_dominates: - { + filter_constraints (master_dofs, + slave_dofs, + *(subface_interpolation_matrices + [cell->active_fe_index()][subface_fe_index][c]), + constraints); + } + + break; + } + + case FiniteElementDomination::other_element_dominates: + case FiniteElementDomination::neither_element_dominates: + { // Case 2 (the "complex" // case): at least one // (the neither_... case) @@ -2830,8 +2830,8 @@ namespace DoFTools // dof handlers, // add a check // here... - Assert (DoFHandlerSupportsDifferentFEs::value == true, - ExcInternalError()); + Assert (DoFHandlerSupportsDifferentFEs::value == true, + ExcInternalError()); // we first have // to find the @@ -2842,11 +2842,11 @@ namespace DoFTools // the other ones // can be // constrained to - const unsigned int dominating_fe_index - = get_most_dominating_subface_fe_index (cell->face(face)); + const unsigned int dominating_fe_index + = get_most_dominating_subface_fe_index (cell->face(face)); - const FiniteElement &dominating_fe - = dof_handler.get_fe()[dominating_fe_index]; + const FiniteElement &dominating_fe + = dof_handler.get_fe()[dominating_fe_index]; // check also // that it is @@ -2860,59 +2860,59 @@ namespace DoFTools // into the // branch for the // 'complex' case - Assert ((dominating_fe.compare_for_face_domination - (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0))) - == FiniteElementDomination::this_element_dominates) - || - (dominating_fe.compare_for_face_domination - (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0))) - == FiniteElementDomination::either_element_can_dominate), - ExcInternalError()); + Assert ((dominating_fe.compare_for_face_domination + (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0))) + == FiniteElementDomination::this_element_dominates) + || + (dominating_fe.compare_for_face_domination + (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0))) + == FiniteElementDomination::either_element_can_dominate), + ExcInternalError()); // first get the // interpolation matrix // from the mother to the // virtual dofs - Assert (dominating_fe.dofs_per_face <= - cell->get_fe().dofs_per_face, - ExcInternalError()); + Assert (dominating_fe.dofs_per_face <= + cell->get_fe().dofs_per_face, + ExcInternalError()); - ensure_existence_of_face_matrix - (dominating_fe, - cell->get_fe(), - face_interpolation_matrices - [dominating_fe_index][cell->active_fe_index()]); + ensure_existence_of_face_matrix + (dominating_fe, + cell->get_fe(), + face_interpolation_matrices + [dominating_fe_index][cell->active_fe_index()]); // split this matrix into // master and slave // components. invert the // master component - ensure_existence_of_master_dof_mask - (cell->get_fe(), - dominating_fe, - (*face_interpolation_matrices - [dominating_fe_index] - [cell->active_fe_index()]), - master_dof_masks - [dominating_fe_index] - [cell->active_fe_index()]); - - ensure_existence_of_split_face_matrix - (*face_interpolation_matrices - [dominating_fe_index][cell->active_fe_index()], - (*master_dof_masks - [dominating_fe_index][cell->active_fe_index()]), - split_face_interpolation_matrices - [dominating_fe_index][cell->active_fe_index()]); - - const FullMatrix &restrict_mother_to_virtual_master_inv - = (split_face_interpolation_matrices - [dominating_fe_index][cell->active_fe_index()]->first); - - const FullMatrix &restrict_mother_to_virtual_slave - = (split_face_interpolation_matrices - [dominating_fe_index][cell->active_fe_index()]->second); + ensure_existence_of_master_dof_mask + (cell->get_fe(), + dominating_fe, + (*face_interpolation_matrices + [dominating_fe_index] + [cell->active_fe_index()]), + master_dof_masks + [dominating_fe_index] + [cell->active_fe_index()]); + + ensure_existence_of_split_face_matrix + (*face_interpolation_matrices + [dominating_fe_index][cell->active_fe_index()], + (*master_dof_masks + [dominating_fe_index][cell->active_fe_index()]), + split_face_interpolation_matrices + [dominating_fe_index][cell->active_fe_index()]); + + const FullMatrix &restrict_mother_to_virtual_master_inv + = (split_face_interpolation_matrices + [dominating_fe_index][cell->active_fe_index()]->first); + + const FullMatrix &restrict_mother_to_virtual_slave + = (split_face_interpolation_matrices + [dominating_fe_index][cell->active_fe_index()]->second); // now compute // the constraint @@ -2922,12 +2922,12 @@ namespace DoFTools // inverse matrix // and the slave // part - constraint_matrix.reinit (cell->get_fe().dofs_per_face - - dominating_fe.dofs_per_face, - dominating_fe.dofs_per_face); - restrict_mother_to_virtual_slave - .mmult (constraint_matrix, - restrict_mother_to_virtual_master_inv); + constraint_matrix.reinit (cell->get_fe().dofs_per_face - + dominating_fe.dofs_per_face, + dominating_fe.dofs_per_face); + restrict_mother_to_virtual_slave + .mmult (constraint_matrix, + restrict_mother_to_virtual_master_inv); // then figure // out the global @@ -2936,29 +2936,29 @@ namespace DoFTools // slave dofs and // apply // constraints - scratch_dofs.resize (cell->get_fe().dofs_per_face); - cell->face(face)->get_dof_indices (scratch_dofs, - cell->active_fe_index ()); + scratch_dofs.resize (cell->get_fe().dofs_per_face); + cell->face(face)->get_dof_indices (scratch_dofs, + cell->active_fe_index ()); // split dofs into master // and slave components - master_dofs.clear (); - slave_dofs.clear (); - for (unsigned int i=0; iget_fe().dofs_per_face; ++i) - if ((*master_dof_masks - [dominating_fe_index][cell->active_fe_index()])[i] == true) - master_dofs.push_back (scratch_dofs[i]); - else - slave_dofs.push_back (scratch_dofs[i]); + master_dofs.clear (); + slave_dofs.clear (); + for (unsigned int i=0; iget_fe().dofs_per_face; ++i) + if ((*master_dof_masks + [dominating_fe_index][cell->active_fe_index()])[i] == true) + master_dofs.push_back (scratch_dofs[i]); + else + slave_dofs.push_back (scratch_dofs[i]); - AssertDimension (master_dofs.size(), dominating_fe.dofs_per_face); - AssertDimension (slave_dofs.size(), - cell->get_fe().dofs_per_face - dominating_fe.dofs_per_face); + AssertDimension (master_dofs.size(), dominating_fe.dofs_per_face); + AssertDimension (slave_dofs.size(), + cell->get_fe().dofs_per_face - dominating_fe.dofs_per_face); - filter_constraints (master_dofs, - slave_dofs, - constraint_matrix, - constraints); + filter_constraints (master_dofs, + slave_dofs, + constraint_matrix, + constraints); @@ -2969,9 +2969,9 @@ namespace DoFTools // as discussed // in the hp // paper - for (unsigned int sf=0; - sfface(face)->n_children(); ++sf) - { + for (unsigned int sf=0; + sfface(face)->n_children(); ++sf) + { // ignore // interfaces // with @@ -2982,62 +2982,62 @@ namespace DoFTools // between // ghost // cells in 2d - if (cell->neighbor_child_on_subface (face, sf)->is_artificial() - || - (dim==2 && cell->is_ghost() - && - cell->neighbor_child_on_subface (face, sf)->is_ghost())) - continue; - - Assert (cell->face(face)->child(sf) - ->n_active_fe_indices() == 1, - ExcInternalError()); - - const unsigned int subface_fe_index - = cell->face(face)->child(sf)->nth_active_fe_index(0); - const FiniteElement &subface_fe - = dof_handler.get_fe()[subface_fe_index]; + if (cell->neighbor_child_on_subface (face, sf)->is_artificial() + || + (dim==2 && cell->is_ghost() + && + cell->neighbor_child_on_subface (face, sf)->is_ghost())) + continue; + + Assert (cell->face(face)->child(sf) + ->n_active_fe_indices() == 1, + ExcInternalError()); + + const unsigned int subface_fe_index + = cell->face(face)->child(sf)->nth_active_fe_index(0); + const FiniteElement &subface_fe + = dof_handler.get_fe()[subface_fe_index]; // first get the // interpolation // matrix from the // subface to the // virtual dofs - Assert (dominating_fe.dofs_per_face <= - subface_fe.dofs_per_face, - ExcInternalError()); - ensure_existence_of_subface_matrix - (dominating_fe, - subface_fe, - sf, - subface_interpolation_matrices - [dominating_fe_index][subface_fe_index][sf]); - - const FullMatrix &restrict_subface_to_virtual - = *(subface_interpolation_matrices - [dominating_fe_index][subface_fe_index][sf]); - - constraint_matrix.reinit (subface_fe.dofs_per_face, - dominating_fe.dofs_per_face); - - restrict_subface_to_virtual - .mmult (constraint_matrix, - restrict_mother_to_virtual_master_inv); - - slave_dofs.resize (subface_fe.dofs_per_face); - cell->face(face)->child(sf)->get_dof_indices (slave_dofs, - subface_fe_index); - - filter_constraints (master_dofs, - slave_dofs, - constraint_matrix, - constraints); - } - - break; - } - - case FiniteElementDomination::no_requirements: + Assert (dominating_fe.dofs_per_face <= + subface_fe.dofs_per_face, + ExcInternalError()); + ensure_existence_of_subface_matrix + (dominating_fe, + subface_fe, + sf, + subface_interpolation_matrices + [dominating_fe_index][subface_fe_index][sf]); + + const FullMatrix &restrict_subface_to_virtual + = *(subface_interpolation_matrices + [dominating_fe_index][subface_fe_index][sf]); + + constraint_matrix.reinit (subface_fe.dofs_per_face, + dominating_fe.dofs_per_face); + + restrict_subface_to_virtual + .mmult (constraint_matrix, + restrict_mother_to_virtual_master_inv); + + slave_dofs.resize (subface_fe.dofs_per_face); + cell->face(face)->child(sf)->get_dof_indices (slave_dofs, + subface_fe_index); + + filter_constraints (master_dofs, + slave_dofs, + constraint_matrix, + constraints); + } + + break; + } + + case FiniteElementDomination::no_requirements: // there // are no // continuity @@ -3047,24 +3047,24 @@ namespace DoFTools // elements. record // no // constraints - break; + break; - default: + default: // we shouldn't get here - Assert (false, ExcInternalError()); - } - } - else - { + Assert (false, ExcInternalError()); + } + } + else + { // this face has no // children, but it // could still be that // it is shared by two // cells that use a // different fe index - Assert (cell->face(face) - ->fe_index_is_active(cell->active_fe_index()) == true, - ExcInternalError()); + Assert (cell->face(face) + ->fe_index_is_active(cell->active_fe_index()) == true, + ExcInternalError()); // see if there is a // neighbor that is @@ -3078,10 +3078,10 @@ namespace DoFTools // cells may not have // an active_fe_index // set, etc - if (!cell->at_boundary(face) - && - cell->neighbor(face)->is_artificial()) - continue; + if (!cell->at_boundary(face) + && + cell->neighbor(face)->is_artificial()) + continue; // Only if there is // a neighbor with @@ -3090,42 +3090,42 @@ namespace DoFTools // and the same h-level, // some action has // to be taken. - if ((DoFHandlerSupportsDifferentFEs::value == true) - && - !cell->face(face)->at_boundary () - && - (cell->neighbor(face)->active_fe_index () != - cell->active_fe_index ()) - && - (!cell->face(face)->has_children() && - !cell->neighbor_is_coarser(face) )) - { - const typename DH::cell_iterator neighbor = cell->neighbor (face); + if ((DoFHandlerSupportsDifferentFEs::value == true) + && + !cell->face(face)->at_boundary () + && + (cell->neighbor(face)->active_fe_index () != + cell->active_fe_index ()) + && + (!cell->face(face)->has_children() && + !cell->neighbor_is_coarser(face) )) + { + const typename DH::cell_iterator neighbor = cell->neighbor (face); // see which side of the // face we have to // constrain - switch (cell->get_fe().compare_for_face_domination (neighbor->get_fe ())) - { - case FiniteElementDomination::this_element_dominates: - { + switch (cell->get_fe().compare_for_face_domination (neighbor->get_fe ())) + { + case FiniteElementDomination::this_element_dominates: + { // Get DoFs on // dominating and // dominated side of // the face - master_dofs.resize (cell->get_fe().dofs_per_face); - cell->face(face)->get_dof_indices (master_dofs, - cell->active_fe_index ()); + master_dofs.resize (cell->get_fe().dofs_per_face); + cell->face(face)->get_dof_indices (master_dofs, + cell->active_fe_index ()); - slave_dofs.resize (neighbor->get_fe().dofs_per_face); - cell->face(face)->get_dof_indices (slave_dofs, - neighbor->active_fe_index ()); + slave_dofs.resize (neighbor->get_fe().dofs_per_face); + cell->face(face)->get_dof_indices (slave_dofs, + neighbor->active_fe_index ()); // break if the n_master_dofs == 0, // because we are attempting to // constrain to an element that has // has no face dofs - if(master_dofs.size() == 0) break; + if(master_dofs.size() == 0) break; // make sure // the element @@ -3133,26 +3133,26 @@ namespace DoFTools // for this // face are // available - ensure_existence_of_face_matrix - (cell->get_fe(), - neighbor->get_fe(), - face_interpolation_matrices - [cell->active_fe_index()][neighbor->active_fe_index()]); + ensure_existence_of_face_matrix + (cell->get_fe(), + neighbor->get_fe(), + face_interpolation_matrices + [cell->active_fe_index()][neighbor->active_fe_index()]); // Add constraints to global constraint // matrix. - filter_constraints (master_dofs, - slave_dofs, - *(face_interpolation_matrices - [cell->active_fe_index()] - [neighbor->active_fe_index()]), - constraints); - - break; - } - - case FiniteElementDomination::other_element_dominates: - { + filter_constraints (master_dofs, + slave_dofs, + *(face_interpolation_matrices + [cell->active_fe_index()] + [neighbor->active_fe_index()]), + constraints); + + break; + } + + case FiniteElementDomination::other_element_dominates: + { // we don't do anything // here since we will // come back to this @@ -3161,11 +3161,11 @@ namespace DoFTools // we will fall into // the first case // clause above - break; - } + break; + } - case FiniteElementDomination::either_element_can_dominate: - { + case FiniteElementDomination::either_element_can_dominate: + { // it appears as if // neither element has // any constraints on @@ -3208,38 +3208,38 @@ namespace DoFTools // finite elements // actually do have // dofs on the face - if ((cell->get_fe().dofs_per_face != 0) - || - (cell->neighbor(face)->get_fe().dofs_per_face != 0)) - { - Assert (cell->get_fe().dofs_per_face - == - cell->neighbor(face)->get_fe().dofs_per_face, - ExcNotImplemented()); + if ((cell->get_fe().dofs_per_face != 0) + || + (cell->neighbor(face)->get_fe().dofs_per_face != 0)) + { + Assert (cell->get_fe().dofs_per_face + == + cell->neighbor(face)->get_fe().dofs_per_face, + ExcNotImplemented()); // (ab)use the master // and slave dofs // arrays for a // moment here - master_dofs.resize (cell->get_fe().dofs_per_face); - cell->face(face) - ->get_dof_indices (master_dofs, - cell->active_fe_index ()); + master_dofs.resize (cell->get_fe().dofs_per_face); + cell->face(face) + ->get_dof_indices (master_dofs, + cell->active_fe_index ()); - slave_dofs.resize (cell->neighbor(face)->get_fe().dofs_per_face); - cell->face(face) - ->get_dof_indices (slave_dofs, - cell->neighbor(face)->active_fe_index ()); + slave_dofs.resize (cell->neighbor(face)->get_fe().dofs_per_face); + cell->face(face) + ->get_dof_indices (slave_dofs, + cell->neighbor(face)->active_fe_index ()); - for (unsigned int i=0; iget_fe().dofs_per_face; ++i) - AssertDimension (master_dofs[i], slave_dofs[i]); - } + for (unsigned int i=0; iget_fe().dofs_per_face; ++i) + AssertDimension (master_dofs[i], slave_dofs[i]); + } - break; - } + break; + } - case FiniteElementDomination::neither_element_dominates: - { + case FiniteElementDomination::neither_element_dominates: + { // we don't presently // know what exactly to // do here. it isn't quite @@ -3250,23 +3250,23 @@ namespace DoFTools // following statement // and see what exactly // is going on - Assert (false, ExcNotImplemented()); - break; - } + Assert (false, ExcNotImplemented()); + break; + } - case FiniteElementDomination::no_requirements: - { + case FiniteElementDomination::no_requirements: + { // nothing to do here - break; - } + break; + } - default: + default: // we shouldn't get // here - Assert (false, ExcInternalError()); - } - } - } + Assert (false, ExcInternalError()); + } + } + } } } @@ -3276,7 +3276,7 @@ namespace DoFTools template void make_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints) + ConstraintMatrix &constraints) { // Decide whether to use the // new or old make_hanging_node_constraints @@ -3288,13 +3288,13 @@ namespace DoFTools // for the moment. if (dof_handler.get_fe().hp_constraints_are_implemented ()) internal:: - make_hp_hanging_node_constraints (dof_handler, - constraints); + make_hp_hanging_node_constraints (dof_handler, + constraints); else internal:: - make_oldstyle_hanging_node_constraints (dof_handler, - constraints, - dealii::internal::int2type()); + make_oldstyle_hanging_node_constraints (dof_handler, + constraints, + dealii::internal::int2type()); } @@ -3311,10 +3311,10 @@ namespace DoFTools ExcMessage ("Faces for periodicity constraints must be on the boundary")); - // In the case that both faces have - // children, we loop over all children - // and applu make_periodicty_constrains - // recursively: + // In the case that both faces have + // children, we loop over all children + // and applu make_periodicty_constrains + // recursively: if (face_1->has_children() && face_2->has_children()) { Assert(face_1->n_children() == GeometryInfo::max_children_per_face && face_2->n_children() == GeometryInfo::max_children_per_face, @@ -3328,22 +3328,22 @@ namespace DoFTools } 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,8 +3361,8 @@ namespace DoFTools 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]) { @@ -3399,7 +3399,7 @@ namespace DoFTools 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) { @@ -3409,37 +3409,37 @@ namespace DoFTools typedef typename DH::cell_iterator CellIterator; - // We collect matching periodic cells on - // the coarsest level: + // We collect matching periodic cells on + // the coarsest level: std::map 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); } } @@ -3459,61 +3459,61 @@ namespace DoFTools // component can be assigned. Then, we sort // them to the first selected component of the // vector system. - // - // the output array dofs_by_component - // lists for each dof the corresponding - // vector component. if the DoFHandler is - // based on a parallel distributed - // triangulation then the output array is - // index by - // dof.locally_owned_dofs().index_within_set(indices[i]) + // + // the output array dofs_by_component + // lists for each dof the corresponding + // vector component. if the DoFHandler is + // based on a parallel distributed + // triangulation then the output array is + // index by + // dof.locally_owned_dofs().index_within_set(indices[i]) template inline void extract_dofs_by_component (const DH &dof, - const std::vector &component_select, - const bool sort_by_blocks, - std::vector &dofs_by_component) + const std::vector &component_select, + const bool sort_by_blocks, + std::vector &dofs_by_component) { const dealii::hp::FECollection - fe_collection (dof.get_fe()); + fe_collection (dof.get_fe()); Assert (fe_collection.n_components() < 256, ExcNotImplemented()); Assert (dofs_by_component.size() == dof.n_locally_owned_dofs(), - ExcDimensionMismatch(dofs_by_component.size(), - dof.n_locally_owned_dofs())); + ExcDimensionMismatch(dofs_by_component.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_component_association[f][d] - // then returns the vector component - // that degree of freedom 'd' belongs - // to (or, in the case of - // sort_by_blocks, the block that it - // corresponds to) + // + // for each element 'f' of the + // FECollection, + // local_component_association[f][d] + // then returns the vector component + // that degree of freedom 'd' belongs + // to (or, in the case of + // sort_by_blocks, the block that it + // corresponds to) std::vector > local_component_association - (fe_collection.size()); + (fe_collection.size()); 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 + // 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; 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::ExcFENotPrimitive()); + typename FiniteElement::ExcFENotPrimitive()); // store a flag whether we should care // about different components. this is @@ -3611,16 +3611,16 @@ namespace DoFTools dof_data = 0; else { - std::vector component_dofs (dof_handler.n_locally_owned_dofs()); - std::vector component_mask (dof_handler.get_fe().n_components(), - false); - component_mask[component] = true; - internal::extract_dofs_by_component (dof_handler, component_mask, - false, component_dofs); - - for (unsigned int i=0; i(component)) - dof_data(i) = 0; + std::vector component_dofs (dof_handler.n_locally_owned_dofs()); + std::vector component_mask (dof_handler.get_fe().n_components(), + false); + component_mask[component] = true; + internal::extract_dofs_by_component (dof_handler, component_mask, + false, component_dofs); + + for (unsigned int i=0; i(component)) + dof_data(i) = 0; } // count how often we have added a value @@ -3628,31 +3628,31 @@ namespace DoFTools std::vector 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; dof_indices.reserve (max_dofs_per_cell(dof_handler)); for (unsigned int present_cell = 0; cell!=endc; ++cell, ++present_cell) { - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (dof_indices); + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (dof_indices); - for (unsigned int i=0; iget_fe().system_to_component_index(i).first == component)) - { + if (!consider_components || + (cell->get_fe().system_to_component_index(i).first == component)) + { // sum up contribution of the // present_cell to this dof - dof_data(dof_indices[i]) += cell_data(present_cell); + dof_data(dof_indices[i]) += cell_data(present_cell); // note that we added another // summand - ++touch_count[dof_indices[i]]; - } + ++touch_count[dof_indices[i]]; + } } // compute the mean value on all the @@ -3664,10 +3664,10 @@ namespace DoFTools // at least once. this needs not be // the case if the vector has more than // one component - Assert (consider_components || (touch_count[i]!=0), - ExcInternalError()); - if (touch_count[i] != 0) - dof_data(i) /= touch_count[i]; + Assert (consider_components || (touch_count[i]!=0), + ExcInternalError()); + if (touch_count[i] != 0) + dof_data(i) /= touch_count[i]; } } @@ -3685,33 +3685,33 @@ namespace DoFTools if (count_by_blocks == true) { - Assert(component_select.size() == fe.n_blocks(), - ExcDimensionMismatch(component_select.size(), fe.n_blocks())); + Assert(component_select.size() == fe.n_blocks(), + ExcDimensionMismatch(component_select.size(), fe.n_blocks())); } else { - Assert(component_select.size() == n_components(dof), - ExcDimensionMismatch(component_select.size(), n_components(dof))); + Assert(component_select.size() == n_components(dof), + ExcDimensionMismatch(component_select.size(), n_components(dof))); } Assert(selected_dofs.size() == dof.n_locally_owned_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof.n_locally_owned_dofs())); + ExcDimensionMismatch(selected_dofs.size(), dof.n_locally_owned_dofs())); // two special cases: no component // is selected, and all components // are selected; both rather // stupid, but easy to catch if (std::count (component_select.begin(), component_select.end(), true) - == 0) + == 0) { - std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false); - return; + std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), false); + return; } else if (std::count (component_select.begin(), component_select.end(), true) - == static_cast(component_select.size())) + == static_cast(component_select.size())) { - std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), true); - return; + std::fill_n (selected_dofs.begin(), dof.n_locally_owned_dofs(), true); + return; } @@ -3727,11 +3727,11 @@ namespace DoFTools // job. std::vector dofs_by_component (dof.n_locally_owned_dofs()); internal::extract_dofs_by_component (dof, component_select, count_by_blocks, - dofs_by_component); + dofs_by_component); for (unsigned int i=0; i(component_select.size())) + == static_cast(component_select.size())) { - std::fill_n (selected_dofs.begin(), dof.n_dofs(), true); - return; + std::fill_n (selected_dofs.begin(), dof.n_dofs(), true); + return; }; @@ -3790,11 +3790,11 @@ namespace DoFTools // job. std::vector dofs_by_component (dof.n_dofs()); internal::extract_dofs_by_component (dof, component_select, count_by_blocks, - dofs_by_component); + dofs_by_component); for (unsigned int i=0; i(component_select.size())) + == static_cast(component_select.size())) { - std::fill_n (selected_dofs.begin(), dof.n_dofs(level), true); - return; + std::fill_n (selected_dofs.begin(), dof.n_dofs(level), true); + return; }; // preset all values by false @@ -3851,13 +3851,13 @@ namespace DoFTools std::vector local_selected_dofs (fe.dofs_per_cell, false); for (unsigned int i=0; iget_mg_dof_indices(indices); - for (unsigned int i=0; iget_mg_dof_indices(indices); + for (unsigned int i=0; iis_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,54 +4008,54 @@ namespace DoFTools 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) + == true) selected_dofs.add_index (face_dof_indices[i]); } else // not primitive { const unsigned int first_nonzero_comp - = (std::find (fe.get_nonzero_components(cell_index).begin(), - fe.get_nonzero_components(cell_index).end(), - true) - - - fe.get_nonzero_components(cell_index).begin()); + = (std::find (fe.get_nonzero_components(cell_index).begin(), + fe.get_nonzero_components(cell_index).end(), + true) + - + fe.get_nonzero_components(cell_index).begin()); Assert (first_nonzero_comp < fe.n_components(), ExcInternalError()); @@ -4071,13 +4071,13 @@ namespace DoFTools template void extract_dofs_with_support_on_boundary (const DH &dof_handler, - const std::vector &component_select, - std::vector &selected_dofs, - const std::set &boundary_indicators) + const std::vector &component_select, + std::vector &selected_dofs, + const std::set &boundary_indicators) { AssertDimension (component_select.size(), n_components(dof_handler)); Assert (boundary_indicators.find (numbers::internal_face_boundary_id) == boundary_indicators.end(), - ExcInvalidBoundaryIndicator()); + ExcInvalidBoundaryIndicator()); // let's see whether we have to // check for certain boundary @@ -4090,7 +4090,7 @@ namespace DoFTools // component is selected, or all const bool check_vector_component = (component_select != std::vector(component_select.size(), - true)); + true)); // clear and reset array by default // values @@ -4113,26 +4113,26 @@ namespace DoFTools // face which we will be visiting // sooner or later for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); - cell!=dof_handler.end(); ++cell) + cell!=dof_handler.end(); ++cell) for (unsigned int face=0; - 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())) - { - const FiniteElement &fe = cell->get_fe(); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - cell_dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (cell_dof_indices); - - for (unsigned int i=0; i::faces_per_cell; ++face) + if (cell->at_boundary(face)) + if (! check_boundary_indicator || + (boundary_indicators.find (cell->face(face)->boundary_indicator()) + != boundary_indicators.end())) + { + const FiniteElement &fe = cell->get_fe(); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + cell_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (cell_dof_indices); + + for (unsigned int i=0; i void extract_hanging_node_dofs (const dealii::DoFHandler<1,spacedim> &dof_handler, - std::vector &selected_dofs) + std::vector &selected_dofs) { - Assert(selected_dofs.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + Assert(selected_dofs.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); // preset all values by false - std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); + std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); // there are no hanging nodes in 1d } @@ -4190,101 +4190,101 @@ namespace DoFTools template void extract_hanging_node_dofs (const dealii::DoFHandler<2,spacedim> &dof_handler, - std::vector &selected_dofs) + std::vector &selected_dofs) { - const unsigned int dim = 2; + const unsigned int dim = 2; - Assert(selected_dofs.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + Assert(selected_dofs.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); // preset all values by false - std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); + std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); - const FiniteElement &fe = dof_handler.get_fe(); + const FiniteElement &fe = dof_handler.get_fe(); // this function is similar to the // make_sparsity_pattern function, // see there for more information - typename dealii::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->has_children()) - { - const typename dealii::DoFHandler::line_iterator - line = cell->face(face); - - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - selected_dofs[line->child(0)->vertex_dof_index(1,dof)] = true; - - for (unsigned int child=0; child<2; ++child) - for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) - selected_dofs[line->child(child)->dof_index(dof)] = true; - } + typename dealii::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->has_children()) + { + const typename dealii::DoFHandler::line_iterator + line = cell->face(face); + + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + selected_dofs[line->child(0)->vertex_dof_index(1,dof)] = true; + + for (unsigned int child=0; child<2; ++child) + for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) + selected_dofs[line->child(child)->dof_index(dof)] = true; + } } template void extract_hanging_node_dofs (const dealii::DoFHandler<3,spacedim> &dof_handler, - std::vector &selected_dofs) + std::vector &selected_dofs) { - const unsigned int dim = 3; + const unsigned int dim = 3; - Assert(selected_dofs.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + Assert(selected_dofs.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); // preset all values by false - std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); + std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); - const FiniteElement &fe = dof_handler.get_fe(); + const FiniteElement &fe = dof_handler.get_fe(); // this function is similar to the // make_sparsity_pattern function, // see there for more information - typename dealii::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->has_children()) - { - const typename dealii::DoFHandler::face_iterator - face = cell->face(f); + typename dealii::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->has_children()) + { + const typename dealii::DoFHandler::face_iterator + face = cell->face(f); - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true; + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true; // dof numbers on the centers of // the lines bounding this face - for (unsigned int line=0; line<4; ++line) - for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) - selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true; + for (unsigned int line=0; line<4; ++line) + for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) + selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true; // next the dofs on the lines interior // to the face; the order of these // lines is laid down in the // FiniteElement class documentation - for (unsigned int dof=0; dofchild(0)->line(1)->dof_index(dof)] = true; - for (unsigned int dof=0; dofchild(1)->line(2)->dof_index(dof)] = true; - for (unsigned int dof=0; dofchild(2)->line(3)->dof_index(dof)] = true; - for (unsigned int dof=0; dofchild(3)->line(0)->dof_index(dof)] = true; + for (unsigned int dof=0; dofchild(0)->line(1)->dof_index(dof)] = true; + for (unsigned int dof=0; dofchild(1)->line(2)->dof_index(dof)] = true; + for (unsigned int dof=0; dofchild(2)->line(3)->dof_index(dof)] = true; + for (unsigned int dof=0; dofchild(3)->line(0)->dof_index(dof)] = true; // dofs on the bordering lines - for (unsigned int line=0; line<4; ++line) - for (unsigned int child=0; child<2; ++child) - for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) - selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true; + for (unsigned int line=0; line<4; ++line) + for (unsigned int child=0; child<2; ++child) + for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof) + selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true; // finally, for the dofs interior // to the four child faces - for (unsigned int child=0; child<4; ++child) - for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) - selected_dofs[face->child(child)->dof_index(dof)] = true; - } + for (unsigned int child=0; child<4; ++child) + for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) + selected_dofs[face->child(child)->dof_index(dof)] = true; + } } } } @@ -4295,10 +4295,10 @@ namespace DoFTools void extract_hanging_node_dofs (const DoFHandler &dof_handler, - std::vector &selected_dofs) + std::vector &selected_dofs) { internal::extract_hanging_node_dofs (dof_handler, - selected_dofs); + selected_dofs); } @@ -4306,11 +4306,11 @@ namespace DoFTools template void extract_subdomain_dofs (const DH &dof_handler, - const types::subdomain_id subdomain_id, - std::vector &selected_dofs) + const types::subdomain_id subdomain_id, + std::vector &selected_dofs) { Assert(selected_dofs.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); + ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs())); // preset all values by false std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); @@ -4326,13 +4326,13 @@ namespace DoFTools endc = dof_handler.end(); for (; cell!=endc; ++cell) if (cell->subdomain_id() == subdomain_id) - { - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - local_dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; iget_fe().dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i void extract_locally_owned_dofs (const DH & dof_handler, - IndexSet & dof_set) + IndexSet & dof_set) { // collect all the locally owned dofs dof_set = dof_handler.locally_owned_dofs(); @@ -4352,7 +4352,7 @@ namespace DoFTools template void extract_locally_active_dofs (const DH & dof_handler, - IndexSet & dof_set) + IndexSet & dof_set) { // collect all the locally owned dofs dof_set = dof_handler.locally_owned_dofs(); @@ -4367,19 +4367,19 @@ namespace DoFTools 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); - - 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_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); + } dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end()); @@ -4391,7 +4391,7 @@ namespace DoFTools template void extract_locally_relevant_dofs (const DH & dof_handler, - IndexSet & dof_set) + IndexSet & dof_set) { // collect all the locally owned dofs dof_set = dof_handler.locally_owned_dofs(); @@ -4406,19 +4406,19 @@ namespace DoFTools 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); - - 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_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); + } dof_set.add_indices(global_dof_indices.begin(), global_dof_indices.end()); @@ -4430,33 +4430,33 @@ namespace DoFTools template void extract_constant_modes (const DH &dof_handler, - const std::vector &component_select, - std::vector > &constant_modes) + const std::vector &component_select, + std::vector > &constant_modes) { const unsigned int n_components = dof_handler.get_fe().n_components(); Assert (n_components == component_select.size(), - ExcDimensionMismatch(n_components, - component_select.size())); + ExcDimensionMismatch(n_components, + component_select.size())); std::vector localized_component (n_components, - numbers::invalid_unsigned_int); + numbers::invalid_unsigned_int); unsigned int n_components_selected = 0; for (unsigned int i=0; i dofs_by_component (dof_handler.n_locally_owned_dofs()); internal::extract_dofs_by_component (dof_handler, component_select, false, - dofs_by_component); + dofs_by_component); unsigned int n_selected_dofs = 0; for (unsigned int i=0; i(n_selected_dofs, - false)); + false)); std::vector component_list (n_components, 0); for (unsigned int d=0; d void get_active_fe_indices (const DH &dof_handler, - std::vector &active_fe_indices) + std::vector &active_fe_indices) { AssertDimension (active_fe_indices.size(), dof_handler.get_tria().n_active_cells()); @@ -4491,25 +4491,25 @@ namespace DoFTools template void get_subdomain_association (const DH &dof_handler, - std::vector &subdomain_association) + std::vector &subdomain_association) { // if the Triangulation is distributed, the // only thing we can usefully ask is for // its locally owned subdomain Assert ((dynamic_cast*> - (&dof_handler.get_tria()) == 0), - 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.")); + Triangulation*> + (&dof_handler.get_tria()) == 0), + 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.")); Assert(subdomain_association.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(subdomain_association.size(), - dof_handler.n_dofs())); + ExcDimensionMismatch(subdomain_association.size(), + dof_handler.n_dofs())); // 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; local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); @@ -4529,14 +4529,14 @@ namespace DoFTools endc = dof_handler.end(); for (; cell!=endc; ++cell) { - Assert (cell->is_artificial() == false, - ExcMessage ("You can't call this function for meshes that " - "have artificial cells.")); + Assert (cell->is_artificial() == false, + ExcMessage ("You can't call this function for meshes that " + "have artificial cells.")); - const types::subdomain_id subdomain_id = cell->subdomain_id(); - const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; - local_dof_indices.resize (dofs_per_cell); - cell->get_dof_indices (local_dof_indices); + const types::subdomain_id subdomain_id = cell->subdomain_id(); + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); + cell->get_dof_indices (local_dof_indices); // set subdomain ids. if dofs // already have their values @@ -4548,23 +4548,23 @@ namespace DoFTools // one, where we take "random" // to be "once this way once // that way" - for (unsigned int i=0; i unsigned int count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id subdomain) + const types::subdomain_id subdomain) { std::vector subdomain_association (dof_handler.n_dofs()); get_subdomain_association (dof_handler, subdomain_association); return std::count (subdomain_association.begin(), - subdomain_association.end(), - subdomain); + subdomain_association.end(), + subdomain); } @@ -4587,7 +4587,7 @@ namespace DoFTools template IndexSet dof_indices_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id subdomain) + const types::subdomain_id subdomain) { // If we have a distributed::Triangulation only @@ -4597,22 +4597,22 @@ namespace DoFTools || (subdomain == dof_handler.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.")); IndexSet index_set (dof_handler.n_dofs()); std::vector local_dof_indices; local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - // first generate an unsorted list of all - // indices which we fill from the back. could - // also insert them directly into the - // IndexSet, but that inserts indices in the - // middle, which is an O(n^2) algorithm and - // hence too expensive. Could also use - // std::set, but that is in general more - // expensive than a vector + // first generate an unsorted list of all + // indices which we fill from the back. could + // also insert them directly into the + // IndexSet, but that inserts indices in the + // middle, which is an O(n^2) algorithm and + // hence too expensive. Could also use + // std::set, but that is in general more + // expensive than a vector std::vector subdomain_indices; typename DH::active_cell_iterator @@ -4630,13 +4630,13 @@ namespace DoFTools local_dof_indices.begin(), local_dof_indices.end()); } - // sort indices and remove duplicates + // sort indices and remove duplicates std::sort (subdomain_indices.begin(), subdomain_indices.end()); subdomain_indices.erase (std::unique(subdomain_indices.begin(), subdomain_indices.end()), subdomain_indices.end()); - // insert into IndexSet + // insert into IndexSet index_set.add_indices (subdomain_indices.begin(), subdomain_indices.end()); index_set.compress (); @@ -4648,12 +4648,12 @@ namespace DoFTools template void count_dofs_with_subdomain_association (const DH &dof_handler, - const types::subdomain_id subdomain, - std::vector &n_dofs_on_subdomain) + const types::subdomain_id subdomain, + std::vector &n_dofs_on_subdomain) { Assert (n_dofs_on_subdomain.size() == dof_handler.get_fe().n_components(), - ExcDimensionMismatch (n_dofs_on_subdomain.size(), - dof_handler.get_fe().n_components())); + ExcDimensionMismatch (n_dofs_on_subdomain.size(), + dof_handler.get_fe().n_components())); std::fill (n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0); // in debug mode, make sure that there are @@ -4663,15 +4663,15 @@ namespace DoFTools { bool found = false; for (typename Triangulation::active_cell_iterator - cell=dof_handler.get_tria().begin_active(); - cell!=dof_handler.get_tria().end(); ++cell) - if (cell->subdomain_id() == subdomain) - { - found = true; - break; - } + cell=dof_handler.get_tria().begin_active(); + cell!=dof_handler.get_tria().end(); ++cell) + if (cell->subdomain_id() == subdomain) + { + found = true; + break; + } Assert (found == true, - ExcMessage ("There are no cells for the given subdomain!")); + ExcMessage ("There are no cells for the given subdomain!")); } #endif @@ -4680,14 +4680,14 @@ namespace DoFTools std::vector component_association (dof_handler.n_dofs()); internal::extract_dofs_by_component (dof_handler, std::vector(), false, - component_association); + component_association); for (unsigned int c=0; c(c))) - ++n_dofs_on_subdomain[c]; + for (unsigned int i=0; i(c))) + ++n_dofs_on_subdomain[c]; } } @@ -4702,41 +4702,41 @@ namespace DoFTools 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) + 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); + { + const FiniteElement& 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 1) + resolve_components(base, dofs_by_component, target_component, + only_once, dofs_per_component, component); + else + { + for (unsigned int dd=0;dd @@ -4748,23 +4748,23 @@ namespace DoFTools std::vector &dofs_per_component, unsigned int &component) { - // assert that all elements in the collection have the same - // structure (base elements and multiplicity, components per base - // element) and then simply call the function above + // assert that all elements in the collection have the same + // structure (base elements and multiplicity, components per base + // element) and then simply call the function above for (unsigned int fe=1; fe bool all_elements_are_primitive (const FiniteElement &fe) { @@ -4786,9 +4786,9 @@ namespace DoFTools } - /** - * Return true if each element of the given element collection is primitive. - */ + /** + * Return true if each element of the given element collection is primitive. + */ template bool all_elements_are_primitive (const dealii::hp::FECollection &fe_collection) { @@ -4818,19 +4818,19 @@ namespace DoFTools // vector as identity. if (target_component.size()==0) { - target_component.resize(n_components); - for (unsigned int i=0; i dofs_by_component (dof_handler.n_locally_owned_dofs()); internal::extract_dofs_by_component (dof_handler, std::vector(), false, - dofs_by_component); + dofs_by_component); // next count what we got unsigned int component = 0; internal::resolve_components(dof_handler.get_fe(), dofs_by_component, target_component, - only_once, dofs_per_component, component); + only_once, dofs_per_component, component); Assert (n_components == component, ExcInternalError()); // finally sanity check. this is @@ -4865,11 +4865,11 @@ namespace DoFTools // is actually primitive, so // exclude other elements from this Assert ((internal::all_elements_are_primitive(dof_handler.get_fe()) == false) - || - (std::accumulate (dofs_per_component.begin(), - dofs_per_component.end(), 0U) - == dof_handler.n_locally_owned_dofs()), - ExcInternalError()); + || + (std::accumulate (dofs_per_component.begin(), + dofs_per_component.end(), 0U) + == dof_handler.n_locally_owned_dofs()), + ExcInternalError()); // reduce information from all CPUs #if defined(DEAL_II_USE_P4EST) && defined(DEAL_II_COMPILER_SUPPORTS_MPI) @@ -4877,13 +4877,13 @@ namespace DoFTools const unsigned int spacedim = DH::space_dimension; if (const parallel::distributed::Triangulation * tria - = (dynamic_cast*> - (&dof_handler.get_tria()))) + = (dynamic_cast*> + (&dof_handler.get_tria()))) { - std::vector local_dof_count = dofs_per_component; + std::vector local_dof_count = dofs_per_component; - MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components, - MPI_UNSIGNED, MPI_SUM, tria->get_communicator()); + MPI_Allreduce ( &local_dof_count[0], &dofs_per_component[0], n_target_components, + MPI_UNSIGNED, MPI_SUM, tria->get_communicator()); } #endif } @@ -4977,11 +4977,11 @@ namespace DoFTools 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) { count_dofs_per_component (dof_handler, dofs_per_component, - false, target_component); + false, target_component); } @@ -5004,18 +5004,18 @@ namespace DoFTools template void compute_intergrid_weights_3 ( - const dealii::DoFHandler &coarse_grid, - const unsigned int coarse_component, - const InterGridMap > &coarse_to_fine_grid_map, - const std::vector > ¶meter_dofs, - const std::vector &weight_mapping, - std::vector > &weights, - const typename dealii::DoFHandler::active_cell_iterator &begin, - const typename dealii::DoFHandler::active_cell_iterator &end) + const dealii::DoFHandler &coarse_grid, + const unsigned int coarse_component, + const InterGridMap > &coarse_to_fine_grid_map, + const std::vector > ¶meter_dofs, + const std::vector &weight_mapping, + std::vector > &weights, + const typename dealii::DoFHandler::active_cell_iterator &begin, + const typename dealii::DoFHandler::active_cell_iterator &end) { // aliases to the finite elements // used by the dof handlers: - const FiniteElement &coarse_fe = coarse_grid.get_fe(); + const FiniteElement &coarse_fe = coarse_grid.get_fe(); // for each cell on the parameter grid: // find out which degrees of freedom on the @@ -5086,44 +5086,44 @@ namespace DoFTools // a single degree of freedom on the // coarse grid (for the selected fe) // on the fine grid - const unsigned int n_fine_dofs = weight_mapping.size(); - dealii::Vector global_parameter_representation (n_fine_dofs); + const unsigned int n_fine_dofs = weight_mapping.size(); + dealii::Vector global_parameter_representation (n_fine_dofs); - typename dealii::DoFHandler::active_cell_iterator cell; - std::vector parameter_dof_indices (coarse_fe.dofs_per_cell); + typename dealii::DoFHandler::active_cell_iterator cell; + std::vector parameter_dof_indices (coarse_fe.dofs_per_cell); - for (cell=begin; cell!=end; ++cell) - { + for (cell=begin; cell!=end; ++cell) + { // get the global indices of the // parameter dofs on this parameter // grid cell - cell->get_dof_indices (parameter_dof_indices); + cell->get_dof_indices (parameter_dof_indices); // loop over all dofs on this // cell and check whether they // are interesting for us - for (unsigned int local_dof=0; - local_dof - set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof], - global_parameter_representation); + coarse_to_fine_grid_map[cell]-> + set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof], + global_parameter_representation); // now that we've got the global // representation of each parameter // dof, we've only got to clobber the @@ -5160,28 +5160,28 @@ namespace DoFTools // different intergrid // weights, have only one // mutex for all of them - static Threads::ThreadMutex mutex; - Threads::ThreadMutex::ScopedLock lock (mutex); - for (unsigned int i=0; i void compute_intergrid_weights_2 ( - const dealii::DoFHandler &coarse_grid, - const unsigned int coarse_component, - const InterGridMap > &coarse_to_fine_grid_map, - const std::vector > ¶meter_dofs, - const std::vector &weight_mapping, - std::vector > &weights) + const dealii::DoFHandler &coarse_grid, + const unsigned int coarse_component, + const InterGridMap > &coarse_to_fine_grid_map, + const std::vector > ¶meter_dofs, + const std::vector &weight_mapping, + std::vector > &weights) { // simply distribute the range of // cells to different threads - typedef typename dealii::DoFHandler::active_cell_iterator active_cell_iterator; - std::vector > - cell_intervals = Threads::split_range (coarse_grid.begin_active(), - coarse_grid.end(), - multithread_info.n_default_threads); + typedef typename dealii::DoFHandler::active_cell_iterator active_cell_iterator; + std::vector > + cell_intervals = Threads::split_range (coarse_grid.begin_active(), + coarse_grid.end(), + multithread_info.n_default_threads); //TODO: use WorkStream here - Threads::TaskGroup<> tasks; - void (*fun_ptr) (const dealii::DoFHandler &, - const unsigned int , - const InterGridMap > &, - const std::vector > &, - const std::vector &, - std::vector > &, - const typename dealii::DoFHandler::active_cell_iterator &, - const typename dealii::DoFHandler::active_cell_iterator &) - = &compute_intergrid_weights_3; - for (unsigned int i=0; i tasks; + void (*fun_ptr) (const dealii::DoFHandler &, + const unsigned int , + const InterGridMap > &, + const std::vector > &, + const std::vector &, + std::vector > &, + const typename dealii::DoFHandler::active_cell_iterator &, + const typename dealii::DoFHandler::active_cell_iterator &) + = &compute_intergrid_weights_3; + for (unsigned int i=0; i unsigned int compute_intergrid_weights_1 ( - const dealii::DoFHandler &coarse_grid, - const unsigned int coarse_component, - const dealii::DoFHandler &fine_grid, - const unsigned int fine_component, - const InterGridMap > &coarse_to_fine_grid_map, - std::vector > &weights, - std::vector &weight_mapping) + const dealii::DoFHandler &coarse_grid, + const unsigned int coarse_component, + const dealii::DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap > &coarse_to_fine_grid_map, + std::vector > &weights, + std::vector &weight_mapping) { // aliases to the finite elements // used by the dof handlers: - const FiniteElement &coarse_fe = coarse_grid.get_fe(), - &fine_fe = fine_grid.get_fe(); + const FiniteElement &coarse_fe = coarse_grid.get_fe(), + &fine_fe = fine_grid.get_fe(); // global numbers of dofs - const unsigned int n_coarse_dofs = coarse_grid.n_dofs(), - n_fine_dofs = fine_grid.n_dofs(); + const unsigned int n_coarse_dofs = coarse_grid.n_dofs(), + n_fine_dofs = fine_grid.n_dofs(); // local numbers of dofs - const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell; + const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell; // alias the number of dofs per // cell belonging to the // coarse_component which is to be // the restriction of the fine // grid: - const unsigned int coarse_dofs_per_cell_component - = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell; + const unsigned int coarse_dofs_per_cell_component + = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell; // Try to find out whether the // grids stem from the same coarse // grid. This is a rather crude // test, but better than nothing - Assert (coarse_grid.get_tria().n_cells(0) == fine_grid.get_tria().n_cells(0), - ExcGridsDontMatch()); + Assert (coarse_grid.get_tria().n_cells(0) == fine_grid.get_tria().n_cells(0), + ExcGridsDontMatch()); // check whether the map correlates // the right objects - Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid, - ExcGridsDontMatch ()); - Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid, - ExcGridsDontMatch ()); + Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid, + ExcGridsDontMatch ()); + Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid, + ExcGridsDontMatch ()); // check whether component numbers @@ -5295,21 +5295,21 @@ namespace DoFTools AssertIndexRange (fine_component, fine_fe.n_components()); // check whether respective finite // elements are equal - Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first) - == - fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first), - ExcFiniteElementsDontMatch()); + Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first) + == + fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first), + ExcFiniteElementsDontMatch()); #ifdef DEBUG // if in debug mode, check whether // the coarse grid is indeed // coarser everywhere than the fine // grid - for (typename dealii::DoFHandler::active_cell_iterator - cell=coarse_grid.begin_active(); - cell != coarse_grid.end(); ++cell) - Assert (cell->level() <= coarse_to_fine_grid_map[cell]->level(), - ExcGridNotCoarser()); + for (typename dealii::DoFHandler::active_cell_iterator + cell=coarse_grid.begin_active(); + cell != coarse_grid.end(); ++cell) + Assert (cell->level() <= coarse_to_fine_grid_map[cell]->level(), + ExcGridNotCoarser()); #endif @@ -5338,89 +5338,89 @@ namespace DoFTools // degree of freedom of the // coarse-grid variable in the // fine-grid element - std::vector > - parameter_dofs (coarse_dofs_per_cell_component, - dealii::Vector(fine_dofs_per_cell)); + std::vector > + parameter_dofs (coarse_dofs_per_cell_component, + dealii::Vector(fine_dofs_per_cell)); // for each coarse dof: find its // position within the fine element // and set this value to one in the // respective vector (all other values // are zero by construction) - for (unsigned int local_coarse_dof=0; - local_coarse_dof dof_is_interesting (fine_grid.n_dofs(), false); - std::vector local_dof_indices (fine_fe.dofs_per_cell); + std::vector dof_is_interesting (fine_grid.n_dofs(), false); + std::vector local_dof_indices (fine_fe.dofs_per_cell); - for (typename dealii::DoFHandler::active_cell_iterator - cell=fine_grid.begin_active(); - cell!=fine_grid.end(); ++cell) - { - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i::active_cell_iterator + cell=fine_grid.begin_active(); + cell!=fine_grid.end(); ++cell) + { + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i local_dof_indices(fine_fe.dofs_per_cell); - unsigned int next_free_index=0; - for (typename dealii::DoFHandler::active_cell_iterator - cell=fine_grid.begin_active(); - cell != fine_grid.end(); ++cell) - { - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i local_dof_indices(fine_fe.dofs_per_cell); + unsigned int next_free_index=0; + for (typename dealii::DoFHandler::active_cell_iterator + cell=fine_grid.begin_active(); + cell != fine_grid.end(); ++cell) + { + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i1) && (sum==0)), ExcInternalError()); - }; + for (unsigned int col=0; col1) && (sum==0)), ExcInternalError()); + }; #endif - return n_parameters_on_fine_grid; + return n_parameters_on_fine_grid; } @@ -5543,13 +5543,13 @@ namespace DoFTools const unsigned int n_parameters_on_fine_grid = internal::compute_intergrid_weights_1 (coarse_grid, coarse_component, - fine_grid, fine_component, - coarse_to_fine_grid_map, - weights, weight_mapping); + fine_grid, fine_component, + coarse_to_fine_grid_map, + weights, weight_mapping); // global numbers of dofs const unsigned int n_coarse_dofs = coarse_grid.n_dofs(), - n_fine_dofs = fine_grid.n_dofs(); + n_fine_dofs = fine_grid.n_dofs(); // get an array in which we store @@ -5558,10 +5558,10 @@ namespace DoFTools std::vector coarse_dof_is_parameter (coarse_grid.n_dofs()); if (true) { - std::vector mask (coarse_grid.get_fe().n_components(), - false); - mask[coarse_component] = true; - extract_dofs (coarse_grid, mask, coarse_dof_is_parameter); + std::vector mask (coarse_grid.get_fe().n_components(), + false); + mask[coarse_component] = true; + extract_dofs (coarse_grid, mask, coarse_dof_is_parameter); }; // now we know that the weights in @@ -5585,50 +5585,50 @@ namespace DoFTools // possibly others) std::vector representants(n_coarse_dofs, -1); for (unsigned int parameter_dof=0; parameter_dof 0, ExcInternalError()); + Assert (weights[parameter_dof].size() > 0, ExcInternalError()); // find the column where the // representant is mentioned - std::map::const_iterator i = weights[parameter_dof].begin(); - for (; i!=weights[parameter_dof].end(); ++i) - if (i->second == 1) - break; - Assert (i!=weights[parameter_dof].end(), ExcInternalError()); - const unsigned int column = i->first; + std::map::const_iterator i = weights[parameter_dof].begin(); + for (; i!=weights[parameter_dof].end(); ++i) + if (i->second == 1) + break; + Assert (i!=weights[parameter_dof].end(), ExcInternalError()); + const unsigned int column = i->first; // now we know in which column of // weights the representant is, but // we don't know its global index. get // it using the inverse operation of // the weight_mapping - unsigned int global_dof=0; - for (; global_dof(column)) - break; - Assert (global_dof < weight_mapping.size(), ExcInternalError()); + unsigned int global_dof=0; + for (; global_dof(column)) + break; + Assert (global_dof < weight_mapping.size(), ExcInternalError()); // now enter the representants global // index into our list - representants[parameter_dof] = global_dof; - } + representants[parameter_dof] = global_dof; + } else - { + { // consistency check: if this // is no parameter dof on the // coarse grid, then the // respective row must be // empty! - Assert (weights[parameter_dof].size() == 0, ExcInternalError()); - }; + Assert (weights[parameter_dof].size() == 0, ExcInternalError()); + }; @@ -5658,52 +5658,52 @@ namespace DoFTools // dof, then we consider this dof // to be unconstrained. otherwise, // all other dofs are constrained - { - const unsigned int col = weight_mapping[global_dof]; - Assert (col < n_parameters_on_fine_grid, ExcInternalError()); - - unsigned int first_used_row=0; - - { - Assert (weights.size() > 0, ExcInternalError()); - std::map::const_iterator - col_entry = weights[0].end(); - for (; first_used_rowsecond == 1) && - (representants[first_used_row] == static_cast(global_dof))) + { + const unsigned int col = weight_mapping[global_dof]; + Assert (col < n_parameters_on_fine_grid, ExcInternalError()); + + unsigned int first_used_row=0; + + { + Assert (weights.size() > 0, ExcInternalError()); + std::map::const_iterator + col_entry = weights[0].end(); + for (; first_used_rowsecond == 1) && + (representants[first_used_row] == static_cast(global_dof))) // dof unconstrained or // constrained to itself // (in case this cell is // mapped to itself, rather // than to children of // itself) - continue; - } + continue; + } // otherwise enter all constraints - constraints.add_line (global_dof); - - constraint_line.clear (); - for (unsigned int row=first_used_row; row::const_iterator - j = weights[row].find(col); - if ((j != weights[row].end()) && (j->second != 0)) - constraint_line.push_back (std::pair(representants[row], - j->second)); - }; - - constraints.add_entries (global_dof, constraint_line); - }; + constraints.add_line (global_dof); + + constraint_line.clear (); + for (unsigned int row=first_used_row; row::const_iterator + j = weights[row].find(col); + if ((j != weights[row].end()) && (j->second != 0)) + constraint_line.push_back (std::pair(representants[row], + j->second)); + }; + + constraints.add_entries (global_dof, constraint_line); + }; } @@ -5767,33 +5767,33 @@ namespace DoFTools std::vector weight_mapping; internal::compute_intergrid_weights_1 (coarse_grid, coarse_component, - fine_grid, fine_component, - coarse_to_fine_grid_map, - weights, weight_mapping); + fine_grid, fine_component, + coarse_to_fine_grid_map, + weights, weight_mapping); // now compute the requested // representation const unsigned int n_global_parm_dofs = std::count_if (weight_mapping.begin(), weight_mapping.end(), - std::bind2nd (std::not_equal_to (), -1)); + std::bind2nd (std::not_equal_to (), -1)); // first construct the inverse // mapping of weight_mapping std::vector inverse_weight_mapping (n_global_parm_dofs, - DoFHandler::invalid_dof_index); + DoFHandler::invalid_dof_index); for (unsigned int i=0; i::invalid_dof_index), - ExcInternalError()); - - inverse_weight_mapping[parameter_dof] = i; - }; + if (parameter_dof != numbers::invalid_unsigned_int) + { + Assert (parameter_dof < n_global_parm_dofs, ExcInternalError()); + Assert ((inverse_weight_mapping[parameter_dof] == DoFHandler::invalid_dof_index), + ExcInternalError()); + + inverse_weight_mapping[parameter_dof] = i; + }; }; // next copy over weights array @@ -5807,14 +5807,14 @@ namespace DoFTools const unsigned int n_coarse_dofs = coarse_grid.n_dofs(); for (unsigned int i=0; i::const_iterator j = weights[i].begin(); - for (; j!=weights[i].end(); ++j) - { - const unsigned int p = inverse_weight_mapping[j->first]; - Assert (psecond; - }; + std::map::const_iterator j = weights[i].begin(); + for (; j!=weights[i].end(); ++j) + { + const unsigned int p = inverse_weight_mapping[j->first]; + Assert (psecond; + }; }; } @@ -5823,13 +5823,13 @@ namespace DoFTools template void map_dof_to_boundary_indices (const DH &dof_handler, - std::vector &mapping) + std::vector &mapping) { Assert (&dof_handler.get_fe() != 0, ExcNoFESelected()); mapping.clear (); mapping.insert (mapping.end(), dof_handler.n_dofs(), - DH::invalid_dof_index); + DH::invalid_dof_index); std::vector dofs_on_face; dofs_on_face.reserve (max_dofs_per_face(dof_handler)); @@ -5849,19 +5849,19 @@ namespace DoFTools // face which we will be visiting // sooner or later typename DH::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + endc = dof_handler.end(); for (; cell!=endc; ++cell) for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->at_boundary(f)) - { - const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; - dofs_on_face.resize (dofs_per_face); - cell->face(f)->get_dof_indices (dofs_on_face, - cell->active_fe_index()); - for (unsigned int i=0; iat_boundary(f)) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_face.resize (dofs_per_face); + cell->face(f)->get_dof_indices (dofs_on_face, + cell->active_fe_index()); + for (unsigned int i=0; i::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_face.resize (dofs_per_face); - cell->face(f)->get_dof_indices (dofs_on_face, cell->active_fe_index()); - for (unsigned int i=0; iface(f)->boundary_indicator()) != + boundary_indicators.end()) + { + const unsigned int dofs_per_face = cell->get_fe().dofs_per_face; + dofs_on_face.resize (dofs_per_face); + cell->face(f)->get_dof_indices (dofs_on_face, cell->active_fe_index()); + for (unsigned int i=0; i > &support_points) { - const unsigned int dim = DH::dimension; - const unsigned int spacedim = DH::space_dimension; + const unsigned int dim = DH::dimension; + const unsigned int spacedim = DH::space_dimension; - hp::FECollection fe_collection(dof_handler.get_fe()); - hp::QCollection q_coll_dummy; + hp::FECollection fe_collection(dof_handler.get_fe()); + hp::QCollection q_coll_dummy; - for (unsigned int fe_index = 0; fe_index < fe_collection.size(); ++fe_index) - { - // check whether every fe in the collection - // has support points - Assert(fe_collection[fe_index].has_support_points(), - typename FiniteElement::ExcFEHasNoSupportPoints()); - q_coll_dummy.push_back( - Quadrature ( - fe_collection[fe_index].get_unit_support_points())); - } + for (unsigned int fe_index = 0; fe_index < fe_collection.size(); ++fe_index) + { + // check whether every fe in the collection + // has support points + Assert(fe_collection[fe_index].has_support_points(), + typename FiniteElement::ExcFEHasNoSupportPoints()); + q_coll_dummy.push_back( + Quadrature ( + fe_collection[fe_index].get_unit_support_points())); + } - // now loop over all cells and - // enquire the support points on - // each of these. we use dummy - // quadrature formulas where the - // quadrature points are located at - // the unit support points to - // enquire the location of the - // support points in real space - // - // the weights of the quadrature - // rule have been set to invalid values - // by the used constructor. - hp::FEValues hp_fe_values(mapping, fe_collection, - q_coll_dummy, update_quadrature_points); - typename DH::active_cell_iterator cell = - dof_handler.begin_active(), endc = dof_handler.end(); - - std::vector local_dof_indices; - for (; cell != endc; ++cell) - // only work on locally relevant cells - if (cell->is_artificial() == false) - { - hp_fe_values.reinit(cell); - const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + // now loop over all cells and + // enquire the support points on + // each of these. we use dummy + // quadrature formulas where the + // quadrature points are located at + // the unit support points to + // enquire the location of the + // support points in real space + // + // the weights of the quadrature + // rule have been set to invalid values + // by the used constructor. + hp::FEValues hp_fe_values(mapping, fe_collection, + q_coll_dummy, update_quadrature_points); + typename DH::active_cell_iterator cell = + dof_handler.begin_active(), endc = dof_handler.end(); + + std::vector local_dof_indices; + for (; cell != endc; ++cell) + // only work on locally relevant cells + if (cell->is_artificial() == false) + { + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); - local_dof_indices.resize(cell->get_fe().dofs_per_cell); - cell->get_dof_indices(local_dof_indices); + local_dof_indices.resize(cell->get_fe().dofs_per_cell); + cell->get_dof_indices(local_dof_indices); - const std::vector > & points = - fe_values.get_quadrature_points(); - for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i) - // insert the values into the map - support_points[local_dof_indices[i]] = points[i]; - } + const std::vector > & points = + fe_values.get_quadrature_points(); + for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i) + // insert the values into the map + support_points[local_dof_indices[i]] = points[i]; + } } @@ -5978,28 +5978,28 @@ namespace DoFTools map_dofs_to_support_points(const hp::MappingCollection & mapping, const DH &dof_handler, std::vector > &support_points) - { - // get the data in the form of the map as above - std::map > x_support_points; - map_dofs_to_support_points(mapping, dof_handler, x_support_points); + { + // get the data in the form of the map as above + std::map > x_support_points; + map_dofs_to_support_points(mapping, dof_handler, x_support_points); - // now convert from the map to the linear vector. make sure every - // entry really appeared in the map - for (unsigned int i=0; i 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) { AssertDimension(support_points.size(), dof_handler.n_dofs()); Assert ((dynamic_cast*> @@ -6009,8 +6009,8 @@ namespace DoFTools ExcMessage ("This function can not be used with distributed triangulations." "See the documentation for more information.")); - //Let the internal function do all the work, - //just make sure that it gets a MappingCollection + //Let the internal function do all the work, + //just make sure that it gets a MappingCollection const hp::MappingCollection mapping_collection(mapping); internal::map_dofs_to_support_points (mapping_collection, @@ -6033,8 +6033,8 @@ namespace DoFTools ExcMessage ("This function can not be used with distributed triangulations." "See the documentation for more information.")); - //Let the internal function do all the work, - //just make sure that it gets a MappingCollection + //Let the internal function do all the work, + //just make sure that it gets a MappingCollection internal::map_dofs_to_support_points (mapping, dof_handler, support_points); @@ -6049,8 +6049,8 @@ namespace DoFTools { support_points.clear(); - //Let the internal function do all the work, - //just make sure that it gets a MappingCollection + //Let the internal function do all the work, + //just make sure that it gets a MappingCollection const hp::MappingCollection mapping_collection(mapping); internal::map_dofs_to_support_points (mapping_collection, @@ -6067,8 +6067,8 @@ namespace DoFTools { support_points.clear(); - //Let the internal function do all the work, - //just make sure that it gets a MappingCollection + //Let the internal function do all the work, + //just make sure that it gets a MappingCollection internal::map_dofs_to_support_points (mapping, dof_handler, support_points); @@ -6091,12 +6091,12 @@ namespace DoFTools for (unsigned int i=0;i& fe = fe_collection[f]; - - const unsigned int nb = fe.n_blocks(); - tables_by_block[f].reinit(nb, nb); - tables_by_block[f].fill(none); - for (unsigned int i=0;i& fe = fe_collection[f]; + + const unsigned int nb = fe.n_blocks(); + tables_by_block[f].reinit(nb, nb); + tables_by_block[f].fill(none); + for (unsigned int i=0;i class DH> void make_zero_boundary_constraints (const DH &dof, - ConstraintMatrix &zero_boundary_constraints, - const std::vector &component_mask_) + ConstraintMatrix &zero_boundary_constraints, + const std::vector &component_mask_) { Assert ((component_mask_.size() == 0) || - (component_mask_.size() == dof.get_fe().n_components()), - ExcMessage ("The number of components in the mask has to be either " - "zero or equal to the number of components in the finite " - "element.")); + (component_mask_.size() == dof.get_fe().n_components()), + ExcMessage ("The number of components in the mask has to be either " + "zero or equal to the number of components in the finite " + "element.")); const unsigned int n_components = DoFTools::n_components (dof); @@ -6150,10 +6150,10 @@ namespace DoFTools // the original value or a vector // of trues const std::vector component_mask ((component_mask_.size() == 0) ? - std::vector (n_components, true) : - component_mask_); + std::vector (n_components, true) : + component_mask_); Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0, - VectorTools::ExcNoComponentSelected()); + VectorTools::ExcNoComponentSelected()); // a field to store the indices std::vector face_dofs; @@ -6164,54 +6164,54 @@ namespace DoFTools endc = dof.end(); for (; cell!=endc; ++cell) for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) - { - const FiniteElement &fe = cell->get_fe(); + ++face_no) + { + const FiniteElement &fe = cell->get_fe(); - typename DH::face_iterator face = cell->face(face_no); + typename DH::face_iterator face = cell->face(face_no); - // if face is on the boundary - if (face->at_boundary ()) - { + // if face is on the boundary + if (face->at_boundary ()) + { // get indices and physical // location on this face - face_dofs.resize (fe.dofs_per_face); - face->get_dof_indices (face_dofs, cell->active_fe_index()); + face_dofs.resize (fe.dofs_per_face); + face->get_dof_indices (face_dofs, cell->active_fe_index()); // enter those dofs into the list // that match the component // signature. - for (unsigned int i=0; i &nonzero_component_array - = cell->get_fe().get_nonzero_components (i); - bool nonzero = false; - for (unsigned int c=0; c &nonzero_component_array + = cell->get_fe().get_nonzero_components (i); + bool nonzero = false; + for (unsigned int c=0; c void make_cell_patches( - Sparsity& block_list, - const DH& dof_handler, - const unsigned int level, - const std::vector& selected_dofs, - unsigned int offset) + Sparsity& block_list, + const DH& dof_handler, + const unsigned int level, + const std::vector& selected_dofs, + unsigned int offset) { typename DH::cell_iterator cell; typename DH::cell_iterator endc = dof_handler.end(level); @@ -6227,15 +6227,15 @@ namespace DoFTools AssertDimension(indices.size(), selected_dofs.size()); for (unsigned int j=0;j& fe = dof_handler.get_fe(); - block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level)); - typename DH::cell_iterator cell; - typename DH::cell_iterator endc = dof_handler.end(level); + { + const FiniteElement& fe = dof_handler.get_fe(); + block_list.reinit(1, dof_handler.n_dofs(level), dof_handler.n_dofs(level)); + typename DH::cell_iterator cell; + typename DH::cell_iterator endc = dof_handler.end(level); - std::vector indices; - std::vector exclude; + std::vector indices; + std::vector exclude; - for (cell=dof_handler.begin(level); cell != endc;++cell) + for (cell=dof_handler.begin(level); cell != endc;++cell) { indices.resize(cell->get_fe().dofs_per_cell); cell->get_mg_dof_indices(indices); @@ -6283,7 +6283,7 @@ namespace DoFTools block_list.add(0, indices[j]); } } - } + } template diff --git a/deal.II/source/dofs/dof_tools.inst.in b/deal.II/source/dofs/dof_tools.inst.in index 8a2c627d71..206641e784 100644 --- a/deal.II/source/dofs/dof_tools.inst.in +++ b/deal.II/source/dofs/dof_tools.inst.in @@ -312,7 +312,7 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS) template void DoFTools::make_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints); + ConstraintMatrix &constraints); } @@ -364,7 +364,7 @@ for (deal_II_dimension : DIMENSIONS) template void DoFTools::make_single_patch (SparsityPattern&, const MGDoFHandler&, unsigned int, bool); - + template void DoFTools::make_child_patches(SparsityPattern&, const MGDoFHandler&, unsigned int, bool, bool); @@ -374,7 +374,7 @@ template void DoFTools:: make_hanging_node_constraints (const DoFHandler &dof_handler, - ConstraintMatrix &constraints); + ConstraintMatrix &constraints); #endif #if deal_II_dimension == 3 @@ -382,7 +382,7 @@ template void DoFTools:: make_hanging_node_constraints (const DoFHandler<1,3> &dof_handler, - ConstraintMatrix &constraints); + ConstraintMatrix &constraints); #endif @@ -605,13 +605,13 @@ unsigned int DoFTools::count_dofs_with_subdomain_association > (const hp::DoFHandler &, const types::subdomain_id); - + template IndexSet DoFTools::dof_indices_with_subdomain_association > (const hp::DoFHandler &, const types::subdomain_id); - + template void DoFTools::get_subdomain_association > @@ -731,15 +731,15 @@ void DoFTools::count_dofs_per_component > ( const hp::DoFHandler&, std::vector&, bool, std::vector); - + template void DoFTools::extract_level_dofs > (const unsigned int, const DoFHandler&, const std::vector&, std::vector&, bool); - -template -void + +template +void DoFTools::extract_level_dofs > (const unsigned int, const MGDoFHandler&, const std::vector&, std::vector&, bool); @@ -765,9 +765,9 @@ void DoFTools::extract_level_dofs > (const unsigned int, const DoFHandler<1, 3>&, const std::vector&, std::vector&, bool); - -template -void + +template +void DoFTools::extract_level_dofs > (const unsigned int, const MGDoFHandler<1, 3>&, const std::vector&, std::vector&, bool); @@ -1027,7 +1027,7 @@ DoFTools::extract_hanging_node_dofs DoFTools::map_dof_to_boundary_indices > (const hp::DoFHandler &, std::vector &); - + #endif #if deal_II_dimension == 3 -- 2.39.5