From f19881c824803c1084fe715478f4a3dfb85d5643 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 27 Apr 2011 15:27:52 +0000 Subject: [PATCH] Implement DoFTools::make_flux_sparsity_pattern for hp DoFHandlers. git-svn-id: https://svn.dealii.org/trunk@23654 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 5 + deal.II/include/deal.II/dofs/dof_tools.h | 48 +- deal.II/source/dofs/dof_tools.cc | 577 ++++++++++++++++------- deal.II/source/dofs/dof_tools.inst.in | 8 +- 4 files changed, 459 insertions(+), 179 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 79ff92045d..8a19ba25bb 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -85,6 +85,11 @@ should be fixed now.

Specific improvements

    +
  1. New: The version of DoFTools::make_flux_sparsity_pattern that takes +the coupling masks is now also available for hp::DoFHandler objects. +
    +(Wolfgang Bangerth, 2011/04/27) +
  2. Fixed: If Triangulation::create_triangulation is called after an hp::DoFHandler object is attached to the triangulation object, setting active FE indices leads to a crash. The problem did not happen if the mesh was diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index 940033f53e..9c2fa7407b 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -178,27 +178,40 @@ class DoFTools public: /** * The flags used in tables by certain - * make_*_pattern - * functions to determine whether - * two components of the solution - * couple in the interior of mesh - * cells or at the boundary. + * make_*_pattern functions to + * describe whether two components of the + * solution couple in the bilinear forms + * corresponding to cell or face + * terms. An example of using these flags + * is shown in the introduction of + * step-46. + * + * In the descriptions of the individual + * elements below, remember that these + * flags are used as elements of tables + * of size FiniteElement::n_components + * times FiniteElement::n_components + * where each element indicates whether + * two components do or do not couple. */ enum Coupling { /** - * The two components do not + * Two components do not * couple. */ none, /** - * The two components do couple. + * Two components do couple. */ always, /** - * The two components couple only - * if their shape functions can be - * nonzero on this face. + * Two components couple only + * if their shape functions are + * both nonzero on a given + * face. This flag is only used + * when computing integrals over + * faces of cells. */ nonzero }; @@ -2084,6 +2097,21 @@ class DoFTools dof_couplings_from_component_couplings (const FiniteElement &fe, const Table<2,Coupling> &component_couplings); + /** + * Same function as above for a + * collection of finite elements, + * returning a collection of tables. + * + * The function currently treats + * DoFTools::Couplings::nonzero the same + * as DoFTools::Couplings::always . + */ + template + static + std::vector > + dof_couplings_from_component_couplings (const hp::FECollection &fe, + const Table<2,Coupling> &component_couplings); + /** * Exception */ diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index c2c6df733a..001dac30c4 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -669,135 +669,198 @@ DoFTools::dof_couplings_from_component_couplings } - -// TODO: look whether one can employ collective add operations in sparsity -// pattern in this function. -template -void -DoFTools::make_flux_sparsity_pattern ( - const DH &dof, - SparsityPattern &sparsity, - const Table<2,Coupling> &int_mask, - const Table<2,Coupling> &flux_mask) +template +std::vector > +DoFTools::dof_couplings_from_component_couplings +(const hp::FECollection &fe, + const Table<2,Coupling> &component_couplings) { - const unsigned int n_dofs = dof.n_dofs(); - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_comp = fe.n_components(); + std::vector > return_value (fe.size()); + for (unsigned int i=0; i dofs_on_this_cell(total_dofs); - std::vector dofs_on_other_cell(total_dofs); - Table<2,bool> support_on_face( - total_dofs, GeometryInfo::faces_per_cell); + - typename DH::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - 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); +namespace internal +{ + namespace DoFTools + { + // implementation of the same function in + // namespace DoFTools for non-hp + // DoFHandlers + template + void + make_flux_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity, + const Table<2,dealii::DoFTools::Coupling> &int_mask, + const Table<2,dealii::DoFTools::Coupling> &flux_mask) + { + const FiniteElement &fe = dof.get_fe(); - for (unsigned int i=0; i::faces_per_cell;++f) - support_on_face(i,f) = fe.has_support_on_face(i,f); + std::vector dofs_on_this_cell(fe.dofs_per_cell); + std::vector dofs_on_other_cell(fe.dofs_per_cell); - // Clear user flags because we will - // need them. But first we save - // them and make sure that we - // restore them later such that at - // the end of this function the - // Triangulation will be in the - // same state as it was at the - // beginning of this function. - std::vector user_flags; - dof.get_tria().save_user_flags(user_flags); - const_cast &>(dof.get_tria()).clear_user_flags (); + const Table<2,dealii::DoFTools::Coupling> + int_dof_mask = dealii::DoFTools::dof_couplings_from_component_couplings(fe, int_mask), + flux_dof_mask = dealii::DoFTools::dof_couplings_from_component_couplings(fe, flux_mask); - for (; cell!=endc; ++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; + 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); - if (cell->at_boundary (face) ) + typename DH::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++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) { - for (unsigned int i=0; iface(face); + if (cell_face->user_flag_set ()) + continue; + + if (cell->at_boundary (face) ) { - const bool i_non_zero_i = support_on_face (i, face); - for (unsigned int j=0; jneighbor(face); - // Refinement edges are taken care of - // by coarser cells - if (cell->neighbor_is_coarser(face)) - continue; + else + { + typename DH::cell_iterator + neighbor = cell->neighbor(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); + 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) + if (cell_face->has_children()) { - const typename DH::cell_iterator - sub_neighbor - = cell->neighbor_child_on_subface (face, sub_nr); + 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; iget_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 (); + neighbor->face(neighbor_face)->set_user_flag (); } } - else - { - neighbor->get_dof_indices (dofs_on_other_cell); - for (unsigned int i=0; i + void + make_flux_sparsity_pattern (const dealii::hp::DoFHandler &dof, + SparsityPattern &sparsity, + const Table<2,dealii::DoFTools::Coupling> &int_mask, + const Table<2,dealii::DoFTools::Coupling> &flux_mask) + { + // while the implementation above is + // quite optimized and caches a lot of + // data (see e.g. the int/flux_dof_mask + // tables), this is no longer practical + // for the hp version since we would + // have to have it for all combinations + // of elements in the + // hp::FECollection. consequently, the + // implementation here is simpler and + // probably less efficient but at least + // readable... + + const dealii::hp::FECollection &fe = dof.get_fe(); + + std::vector dofs_on_this_cell(dealii::DoFTools::max_dofs_per_cell(dof)); + std::vector dofs_on_other_cell(dealii::DoFTools::max_dofs_per_cell(dof)); + + const std::vector > + int_dof_mask + = dealii::DoFTools::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); + // 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) != dealii::DoFTools::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) + == dealii::DoFTools::always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + cell->get_fe().system_to_component_index(j).first) + == dealii::DoFTools::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()) { - const bool i_non_zero_i = support_on_face (i, face); - const bool i_non_zero_e = support_on_face (i, neighbor_face); - for (unsigned int j=0; jn_children(); + ++sub_nr) { - const bool j_non_zero_i = support_on_face (j, face); - const bool j_non_zero_e = support_on_face (j, neighbor_face); - if (flux_dof_mask(i,j) == always) - { - 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_dof_mask(i,j) == nonzero) - { - if (i_non_zero_i && j_non_zero_e) - sparsity.add (dofs_on_this_cell[i], - dofs_on_other_cell[j]); - if (i_non_zero_e && j_non_zero_i) - sparsity.add (dofs_on_other_cell[i], - dofs_on_this_cell[j]); - if (i_non_zero_i && j_non_zero_i) - sparsity.add (dofs_on_this_cell[i], - dofs_on_this_cell[j]); - if (i_non_zero_e && j_non_zero_e) - sparsity.add (dofs_on_other_cell[i], - dofs_on_other_cell[j]); - } + const typename dealii::hp::DoFHandler::cell_iterator + sub_neighbor + = cell->neighbor_child_on_subface (face, sub_nr); - if (flux_dof_mask(j,i) == always) + 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) { - 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]); + 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) + == dealii::DoFTools::always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + sub_neighbor->get_fe().system_to_component_index(j).first) + == dealii::DoFTools::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) + == dealii::DoFTools::always) + || + (flux_mask(sub_neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == dealii::DoFTools::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]); + } + } } - if (flux_dof_mask(j,i) == nonzero) + 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 (j_non_zero_i && i_non_zero_e) - sparsity.add (dofs_on_this_cell[j], - dofs_on_other_cell[i]); - if (j_non_zero_e && i_non_zero_i) - sparsity.add (dofs_on_other_cell[j], - dofs_on_this_cell[i]); - if (j_non_zero_i && i_non_zero_i) - sparsity.add (dofs_on_this_cell[j], - dofs_on_this_cell[i]); - if (j_non_zero_e && i_non_zero_e) - sparsity.add (dofs_on_other_cell[j], - dofs_on_other_cell[i]); + if ((flux_mask(cell->get_fe().system_to_component_index(i).first, + neighbor->get_fe().system_to_component_index(j).first) + == dealii::DoFTools::always) + || + (flux_mask(cell->get_fe().system_to_component_index(i).first, + neighbor->get_fe().system_to_component_index(j).first) + == dealii::DoFTools::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) + == dealii::DoFTools::always) + || + (flux_mask(neighbor->get_fe().system_to_component_index(j).first, + cell->get_fe().system_to_component_index(i).first) + == dealii::DoFTools::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 (); } - neighbor->face(neighbor_face)->set_user_flag (); } } } } + + } +} + + + +template +void +DoFTools:: +make_flux_sparsity_pattern (const DH &dof, + 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 + // specialized functions in the internal + // namespace + const unsigned int n_dofs = dof.n_dofs(); + const unsigned int n_comp = dof.get_fe().n_components(); + + Assert (sparsity.n_rows() == n_dofs, + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + Assert (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)); + Assert (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)); + Assert (flux_mask.n_cols() == n_comp, + ExcDimensionMismatch (flux_mask.n_cols(), n_comp)); - // finally restore the user flags + + // Clear user flags because we will + // need them. But first we save + // them and make sure that we + // restore them later such that at + // the end of this function the + // Triangulation will be in the + // same state as it was at the + // beginning of this function. + std::vector user_flags; + dof.get_tria().save_user_flags(user_flags); + const_cast &>(dof.get_tria()).clear_user_flags (); + + internal::DoFTools::make_flux_sparsity_pattern (dof, sparsity, + int_mask, flux_mask); + + // finally restore the user flags const_cast &>(dof.get_tria()).load_user_flags(user_flags); } + + + + namespace internal { namespace DoFTools diff --git a/deal.II/source/dofs/dof_tools.inst.in b/deal.II/source/dofs/dof_tools.inst.in index bc4707ccd2..adbc04d862 100644 --- a/deal.II/source/dofs/dof_tools.inst.in +++ b/deal.II/source/dofs/dof_tools.inst.in @@ -137,6 +137,12 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS) const ConstraintMatrix &constraints, const bool, const unsigned int); + template void + DoFTools::make_flux_sparsity_pattern,SP> + (const hp::DoFHandler &dof, + SP &, + const Table<2,Coupling>&, + const Table<2,Coupling>&); #endif #if deal_II_dimension < 3 @@ -637,4 +643,4 @@ DoFTools::extract_hanging_node_dofs #endif -} \ No newline at end of file +} -- 2.39.5