From: kronbichler Date: Thu, 22 Aug 2013 04:19:24 +0000 (+0000) Subject: Split two more files into three files each. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a3c89b72fc59fb01880d0e5854e5a578b534ec6b;p=dealii-svn.git Split two more files into three files each. git-svn-id: https://svn.dealii.org/trunk@30407 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/fe/component_mask.h b/deal.II/include/deal.II/fe/component_mask.h index 0661aa2c63..bb60046fc4 100644 --- a/deal.II/include/deal.II/fe/component_mask.h +++ b/deal.II/include/deal.II/fe/component_mask.h @@ -209,6 +209,11 @@ public: std::size_t memory_consumption () const; + /** + * Exception + */ + DeclException0 (ExcNoComponentSelected); + private: /** * The actual component mask. diff --git a/deal.II/include/deal.II/numerics/vector_tools.h b/deal.II/include/deal.II/numerics/vector_tools.h index d7cf41a6ed..1121e1ff66 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.h +++ b/deal.II/include/deal.II/numerics/vector_tools.h @@ -2418,11 +2418,6 @@ namespace VectorTools * Exception */ DeclException0 (ExcNonInterpolatingFE); - - /** - * Exception - */ - DeclException0 (ExcNoComponentSelected); } diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 0c0799fc5d..32bd467905 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -1524,7 +1524,7 @@ namespace VectorTools boundary_function.n_components)); Assert (component_mask.n_selected_components(fe.n_components()) > 0, - ExcNoComponentSelected()); + ComponentMask::ExcNoComponentSelected()); // now set the value of // the vertex degree of @@ -5482,7 +5482,7 @@ namespace VectorTools } // Error out if we have not constrained anything. Note that in this // case the vector v is always nonempty. - Assert (n == 0 || counter > 0, ExcNoComponentSelected()); + Assert (n == 0 || counter > 0, ComponentMask::ExcNoComponentSelected()); s /= counter; diff --git a/deal.II/source/dofs/CMakeLists.txt b/deal.II/source/dofs/CMakeLists.txt index f791bd398e..e3bf9d5631 100644 --- a/deal.II/source/dofs/CMakeLists.txt +++ b/deal.II/source/dofs/CMakeLists.txt @@ -23,23 +23,31 @@ FILE(GLOB _header SET(_src block_info.cc dof_accessor.cc + dof_accessor_get.cc + dof_accessor_set.cc dof_faces.cc dof_handler.cc dof_handler_policy.cc dof_objects.cc dof_renumbering.cc dof_tools.cc + dof_tools_constraints.cc + dof_tools_sparsity.cc number_cache.cc ) SET(_inst block_info.inst.in dof_accessor.inst.in + dof_accessor_get.inst.in + dof_accessor_set.inst.in dof_handler.inst.in dof_handler_policy.inst.in dof_objects.inst.in dof_renumbering.inst.in dof_tools.inst.in + dof_tools_constraints.inst.in + dof_tools_sparsity.inst.in ) DEAL_II_ADD_LIBRARY(obj_dofs OBJECT ${_src} ${_header} ${_inst}) diff --git a/deal.II/source/dofs/dof_accessor.cc b/deal.II/source/dofs/dof_accessor.cc index 53b67ac4d1..34dbb34906 100644 --- a/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/source/dofs/dof_accessor.cc @@ -14,16 +14,6 @@ // // --------------------------------------------------------------------- -#include -#include -#include -#include -#include -#include -#include -#include -#include - #include #include #include @@ -97,134 +87,6 @@ DoFCellAccessor::neighbor_child_on_subface (const unsigned int face, -template -template -void -DoFCellAccessor:: -get_interpolated_dof_values (const InputVector &values, - Vector &interpolated_values) const -{ - const FiniteElement &fe = this->get_fe(); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - - Assert (this->dof_handler != 0, - typename BaseClass::ExcInvalidObject()); - Assert (&fe != 0, - typename BaseClass::ExcInvalidObject()); - Assert (interpolated_values.size() == dofs_per_cell, - typename BaseClass::ExcVectorDoesNotMatch()); - Assert (values.size() == this->dof_handler->n_dofs(), - typename BaseClass::ExcVectorDoesNotMatch()); - - if (!this->has_children()) - // if this cell has no children: simply return the exact values on this - // cell - this->get_dof_values (values, interpolated_values); - else - // otherwise clobber them from the children - { - Vector tmp1(dofs_per_cell); - Vector tmp2(dofs_per_cell); - - interpolated_values = 0; - - // later on we will have to push the values interpolated from the child - // to the mother cell into the output vector. unfortunately, there are - // two types of elements: ones where you add up the contributions from - // the different child cells, and ones where you overwrite. - // - // an example for the first is piecewise constant (and discontinuous) - // elements, where we build the value on the coarse cell by averaging - // the values from the cell (i.e. by adding up a fraction of the values - // of their values) - // - // an example for the latter are the usual continuous elements. the - // value on a vertex of a coarse cell must there be the same, - // irrespective of the adjacent cell we are presently on. so we always - // overwrite. in fact, we must, since we cannot know in advance how many - // neighbors there will be, so there is no way to compute the average - // with fixed factors - // - // so we have to find out to which type this element belongs. the - // difficulty is: the finite element may be a composed one, so we can - // only hope to do this for each shape function individually. in fact, - // there are even weird finite elements (for example the Raviart-Thomas - // element) which have shape functions that are additive (interior ones) - // and others that are overwriting (face degrees of freedom that need to - // be continuous across the face). to avoid checking this over and over - // again, we do this once now and cache the results - std::vector restriction_is_additive (dofs_per_cell); - for (unsigned int i=0; in_children(); ++child) - { - // get the values from the present child, if necessary by - // interpolation itself - this->child(child)->get_interpolated_dof_values (values, - tmp1); - // interpolate these to the mother cell - fe.get_restriction_matrix(child, this->refinement_case()).vmult (tmp2, tmp1); - - // and add up or set them in the output vector - for (unsigned int i=0; i -template -void -DoFCellAccessor:: -set_dof_values_by_interpolation (const Vector &local_values, - OutputVector &values) const -{ - const unsigned int dofs_per_cell = this->get_fe().dofs_per_cell; - - Assert (this->dof_handler != 0, - typename BaseClass::ExcInvalidObject()); - Assert (&this->get_fe() != 0, - typename BaseClass::ExcInvalidObject()); - Assert (local_values.size() == dofs_per_cell, - typename BaseClass::ExcVectorDoesNotMatch()); - Assert (values.size() == this->dof_handler->n_dofs(), - typename BaseClass::ExcVectorDoesNotMatch()); - - if (!this->has_children()) - // if this cell has no children: simply set the values on this cell - this->set_dof_values (local_values, values); - else - // otherwise distribute them to the children - { - Vector tmp(dofs_per_cell); - - for (unsigned int child=0; childn_children(); ++child) - { - Assert (this->child(child)->get_fe().dofs_per_cell == dofs_per_cell, - ExcNotImplemented()); - - // prolong the given data to the present cell. FullMatrix only wants - // us to call vmult if the matrix size is actually non-zero, so - // check that case - if (tmp.size() > 0) - { - this->get_fe().get_prolongation_matrix(child, this->refinement_case()) - .vmult (tmp, local_values); - - this->child(child)->set_dof_values_by_interpolation (tmp, values); - } - } - } -} - - - // -------------------------------------------------------------------------- // explicit instantiations #include "dof_accessor.inst" diff --git a/deal.II/source/dofs/dof_accessor.inst.in b/deal.II/source/dofs/dof_accessor.inst.in index 103290a711..297e98c829 100644 --- a/deal.II/source/dofs/dof_accessor.inst.in +++ b/deal.II/source/dofs/dof_accessor.inst.in @@ -165,82 +165,3 @@ for (deal_II_dimension : DIMENSIONS; lda : BOOL) #endif } - -for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS; lda : BOOL) - { - template - void - DoFCellAccessor, lda>::get_interpolated_dof_values - (const VEC&, Vector&) const; - template - void - DoFCellAccessor, lda>::set_dof_values_by_interpolation - (const Vector&, VEC&) const; - -#if deal_II_dimension != 3 - - template - void - DoFCellAccessor, lda>::get_interpolated_dof_values - (const VEC&, Vector&) const; - template - void - DoFCellAccessor, lda>::set_dof_values_by_interpolation - (const Vector&, VEC&) const; - -#endif - -#if deal_II_dimension == 3 - - template - void - DoFCellAccessor, lda>::get_interpolated_dof_values - (const VEC&, Vector&) const; - template - void - DoFCellAccessor, lda>::set_dof_values_by_interpolation - (const Vector&, VEC&) const; - -#endif - - } - - -for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS; lda : BOOL) - { - template - void - DoFCellAccessor, lda>::get_interpolated_dof_values - (const VEC&, Vector&) const; - template - void - DoFCellAccessor, lda>::set_dof_values_by_interpolation - (const Vector&, VEC&) const; - -#if deal_II_dimension != 3 - - template - void - DoFCellAccessor, lda>::get_interpolated_dof_values - (const VEC&, Vector&) const; - template - void - DoFCellAccessor, lda>::set_dof_values_by_interpolation - (const Vector&, VEC&) const; - -#endif - -#if deal_II_dimension == 3 - - template - void - DoFCellAccessor, lda>::get_interpolated_dof_values - (const VEC&, Vector&) const; - template - void - DoFCellAccessor, lda>::set_dof_values_by_interpolation - (const Vector&, VEC&) const; - -#endif - } - diff --git a/deal.II/source/dofs/dof_accessor_get.cc b/deal.II/source/dofs/dof_accessor_get.cc new file mode 100644 index 0000000000..9b9ca8b7b8 --- /dev/null +++ b/deal.II/source/dofs/dof_accessor_get.cc @@ -0,0 +1,125 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + +template +template +void +DoFCellAccessor:: +get_interpolated_dof_values (const InputVector &values, + Vector &interpolated_values) const +{ + const FiniteElement &fe = this->get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + Assert (this->dof_handler != 0, + typename BaseClass::ExcInvalidObject()); + Assert (&fe != 0, + typename BaseClass::ExcInvalidObject()); + Assert (interpolated_values.size() == dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (values.size() == this->dof_handler->n_dofs(), + typename BaseClass::ExcVectorDoesNotMatch()); + + if (!this->has_children()) + // if this cell has no children: simply return the exact values on this + // cell + this->get_dof_values (values, interpolated_values); + else + // otherwise clobber them from the children + { + Vector tmp1(dofs_per_cell); + Vector tmp2(dofs_per_cell); + + interpolated_values = 0; + + // later on we will have to push the values interpolated from the child + // to the mother cell into the output vector. unfortunately, there are + // two types of elements: ones where you add up the contributions from + // the different child cells, and ones where you overwrite. + // + // an example for the first is piecewise constant (and discontinuous) + // elements, where we build the value on the coarse cell by averaging + // the values from the cell (i.e. by adding up a fraction of the values + // of their values) + // + // an example for the latter are the usual continuous elements. the + // value on a vertex of a coarse cell must there be the same, + // irrespective of the adjacent cell we are presently on. so we always + // overwrite. in fact, we must, since we cannot know in advance how many + // neighbors there will be, so there is no way to compute the average + // with fixed factors + // + // so we have to find out to which type this element belongs. the + // difficulty is: the finite element may be a composed one, so we can + // only hope to do this for each shape function individually. in fact, + // there are even weird finite elements (for example the Raviart-Thomas + // element) which have shape functions that are additive (interior ones) + // and others that are overwriting (face degrees of freedom that need to + // be continuous across the face). to avoid checking this over and over + // again, we do this once now and cache the results + std::vector restriction_is_additive (dofs_per_cell); + for (unsigned int i=0; in_children(); ++child) + { + // get the values from the present child, if necessary by + // interpolation itself + this->child(child)->get_interpolated_dof_values (values, + tmp1); + // interpolate these to the mother cell + fe.get_restriction_matrix(child, this->refinement_case()).vmult (tmp2, tmp1); + + // and add up or set them in the output vector + for (unsigned int i=0; i, lda>::get_interpolated_dof_values + (const VEC&, Vector&) const; + +#if deal_II_dimension != 3 + + template + void + DoFCellAccessor, lda>::get_interpolated_dof_values + (const VEC&, Vector&) const; + +#endif + +#if deal_II_dimension == 3 + + template + void + DoFCellAccessor, lda>::get_interpolated_dof_values + (const VEC&, Vector&) const; + +#endif + + } + + +for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS; lda : BOOL) + { + template + void + DoFCellAccessor, lda>::get_interpolated_dof_values + (const VEC&, Vector&) const; + +#if deal_II_dimension != 3 + + template + void + DoFCellAccessor, lda>::get_interpolated_dof_values + (const VEC&, Vector&) const; + +#endif + +#if deal_II_dimension == 3 + + template + void + DoFCellAccessor, lda>::get_interpolated_dof_values + (const VEC&, Vector&) const; + +#endif + } + diff --git a/deal.II/source/dofs/dof_accessor_set.cc b/deal.II/source/dofs/dof_accessor_set.cc new file mode 100644 index 0000000000..f11fc845a0 --- /dev/null +++ b/deal.II/source/dofs/dof_accessor_set.cc @@ -0,0 +1,93 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + + + +template +template +void +DoFCellAccessor:: +set_dof_values_by_interpolation (const Vector &local_values, + OutputVector &values) const +{ + const unsigned int dofs_per_cell = this->get_fe().dofs_per_cell; + + Assert (this->dof_handler != 0, + typename BaseClass::ExcInvalidObject()); + Assert (&this->get_fe() != 0, + typename BaseClass::ExcInvalidObject()); + Assert (local_values.size() == dofs_per_cell, + typename BaseClass::ExcVectorDoesNotMatch()); + Assert (values.size() == this->dof_handler->n_dofs(), + typename BaseClass::ExcVectorDoesNotMatch()); + + if (!this->has_children()) + // if this cell has no children: simply set the values on this cell + this->set_dof_values (local_values, values); + else + // otherwise distribute them to the children + { + Vector tmp(dofs_per_cell); + + for (unsigned int child=0; childn_children(); ++child) + { + Assert (this->child(child)->get_fe().dofs_per_cell == dofs_per_cell, + ExcNotImplemented()); + + // prolong the given data to the present cell. FullMatrix only wants + // us to call vmult if the matrix size is actually non-zero, so + // check that case + if (tmp.size() > 0) + { + this->get_fe().get_prolongation_matrix(child, this->refinement_case()) + .vmult (tmp, local_values); + + this->child(child)->set_dof_values_by_interpolation (tmp, values); + } + } + } +} + + + +// -------------------------------------------------------------------------- +// explicit instantiations +#include "dof_accessor_set.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/dofs/dof_accessor_set.inst.in b/deal.II/source/dofs/dof_accessor_set.inst.in new file mode 100644 index 0000000000..74948318e8 --- /dev/null +++ b/deal.II/source/dofs/dof_accessor_set.inst.in @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS; lda : BOOL) + { + template + void + DoFCellAccessor, lda>::set_dof_values_by_interpolation + (const Vector&, VEC&) const; + +#if deal_II_dimension != 3 + + template + void + DoFCellAccessor, lda>::set_dof_values_by_interpolation + (const Vector&, VEC&) const; + +#endif + +#if deal_II_dimension == 3 + + template + void + DoFCellAccessor, lda>::set_dof_values_by_interpolation + (const Vector&, VEC&) const; + +#endif + + } + + +for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS; lda : BOOL) + { + template + void + DoFCellAccessor, lda>::set_dof_values_by_interpolation + (const Vector&, VEC&) const; + +#if deal_II_dimension != 3 + + template + void + DoFCellAccessor, lda>::set_dof_values_by_interpolation + (const Vector&, VEC&) const; + +#endif + +#if deal_II_dimension == 3 + + template + void + DoFCellAccessor, lda>::set_dof_values_by_interpolation + (const Vector&, VEC&) const; + +#endif + } + diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index ad41b28892..51b5f2ec55 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -41,7 +41,6 @@ #include #include #include -#include #include @@ -54,3212 +53,6 @@ DEAL_II_NAMESPACE_OPEN 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) - { - const types::global_dof_index n_dofs = dof.n_dofs(); - - Assert (sparsity.n_rows() == n_dofs, - ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); - Assert (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 subdomain is also okay, because we skip - // ghost cells in the loop below. - Assert ( - (dof.get_tria().locally_owned_subdomain() == numbers::invalid_subdomain_id) - || - (subdomain_id == numbers::invalid_subdomain_id) - || - (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.")); - - 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(); - - // In case we work with a distributed sparsity pattern of Trilinos - // type, we only have to do the work if the current cell is owned by - // the calling processor. Otherwise, just continue. - for (; cell!=endc; ++cell) - if (((subdomain_id == numbers::invalid_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); - - // 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); - } - } - - - - 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 types::global_dof_index n_dofs = dof.n_dofs(); - - 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 (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())); - - // If we have a distributed::Triangulation only allow locally_owned - // subdomain. Not setting a subdomain is also okay, because we skip - // ghost cells in the loop below. - Assert ( - (dof.get_tria().locally_owned_subdomain() == numbers::invalid_subdomain_id) - || - (subdomain_id == numbers::invalid_subdomain_id) - || - (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.")); - - const hp::FECollection fe_collection (dof.get_fe()); - - // first, for each finite element, build a mask for each dof, not like - // the one given which represents components. make sure we do the right - // thing also with respect to non-primitive shape functions, which - // takes some additional thought - std::vector > dof_mask(fe_collection.size()); - - // check whether the table of couplings contains only true arguments, - // i.e., we do not exclude any index. that is the easy case, since we - // don't have to set up the tables - 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(); - - // In case we work with a distributed sparsity pattern of Trilinos - // type, we only have to do the work if the current cell is owned by - // the calling processor. Otherwise, just continue. - for (; cell!=endc; ++cell) - if (((subdomain_id == numbers::invalid_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; - - 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]); - } - } - - - - template - void - make_sparsity_pattern ( - const DH &dof_row, - const DH &dof_col, - SparsityPattern &sparsity) - { - const types::global_dof_index n_dofs_row = dof_row.n_dofs(); - const types::global_dof_index n_dofs_col = dof_col.n_dofs(); - - Assert (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)); - -//TODO: Looks like wasteful memory management here - - const std::list > - cell_list - = GridTools::get_finest_common_cells (dof_row, dof_col); - - - typename std::list >::const_iterator - cell_iter = cell_list.begin(); - - 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; i - void - make_boundary_sparsity_pattern ( - const DH &dof, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &sparsity) - { - if (DH::dimension == 1) - { - // there are only 2 boundary indicators in 1d, so it is no - // performance problem to call the other function - 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 types::global_dof_index n_dofs = dof.n_dofs(); - - AssertDimension (dof_to_boundary_mapping.size(), n_dofs); - AssertDimension (sparsity.n_rows(), dof.n_boundary_dofs()); - AssertDimension (sparsity.n_cols(), dof.n_boundary_dofs()); -#ifdef DEBUG - if (sparsity.n_rows() != 0) - { - types::global_dof_index max_element = 0; - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != DH::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - AssertDimension (max_element, sparsity.n_rows()-1); - }; -#endif - - std::vector dofs_on_this_face; - dofs_on_this_face.reserve (max_dofs_per_face(dof)); - - // loop over all faces to check whether they are at a boundary. note - // that we need not take special care of single lines (using - // @p{cell->has_boundary_lines}), since we do not support boundaries of - // dimension dim-2, and so every boundary line is also part of a - // boundary face. - typename DH::active_cell_iterator cell = dof.begin_active(), - 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()); - - // make sparsity pattern for this cell - for (unsigned int i=0; i - void make_boundary_sparsity_pattern ( - const DH &dof, - const typename FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SparsityPattern &sparsity) - { - if (DH::dimension == 1) - { - // first check left, then right boundary point - for (unsigned int direction=0; direction<2; ++direction) - { - // if this boundary is not requested, then go on with next one - if (boundary_indicators.find(direction) == - boundary_indicators.end()) - continue; - - // find active cell at that boundary: first go to left/right, - // then to children - typename DH::level_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); - - // next get boundary mapped dof indices of boundary dofs - 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); - }; -#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(); - 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()); - - // 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, - const types::subdomain_id subdomain_id) - - // TODO: QA: reduce the indentation level of this method..., Maier 2012 - - { - const types::global_dof_index n_dofs = dof.n_dofs(); - - AssertDimension (sparsity.n_rows(), n_dofs); - AssertDimension (sparsity.n_cols(), n_dofs); - - // If we have a distributed::Triangulation only allow locally_owned - // subdomain. Not setting a subdomain is also okay, because we skip - // ghost cells in the loop below. - Assert ( - (dof.get_tria().locally_owned_subdomain() == numbers::invalid_subdomain_id) - || - (subdomain_id == numbers::invalid_subdomain_id) - || - (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.")); - - 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(); - - // TODO: in an old implementation, we used user flags before to tag - // faces that were already touched. this way, we could reduce the work - // a little bit. now, we instead add only data from one side. this - // should be OK, but we need to actually verify it. - - // In case we work with a distributed sparsity pattern of Trilinos - // type, we only have to do the work if the current cell is owned by - // the calling processor. Otherwise, just continue. - for (; cell!=endc; ++cell) - if (((subdomain_id == numbers::invalid_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); - - // 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::level_cell_iterator neighbor = cell->neighbor(face); - - // in 1d, we do not need to worry whether the neighbor - // might have children and then loop over those children. - // rather, we may as well go straight to to cell behind - // this particular cell's most terminal child - if (DH::dimension==1) - while (neighbor->has_children()) - neighbor = neighbor->child(face==0 ? 1 : 0); - - if (neighbor->has_children()) - { - for (unsigned int sub_nr = 0; - sub_nr != cell_face->number_of_children(); - ++sub_nr) - { - const typename DH::level_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 - - // TODO: in the distributed case, we miss out the - // constraints when the neighbor cell is coarser, but - // only the current cell is owned locally! - 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); - - 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); - - // 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)->active() - || - (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); - } - } - } - } - } - - - - template - void - make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity) - { - ConstraintMatrix constraints; - make_flux_sparsity_pattern (dof, sparsity, constraints); - } - - - - template - Table<2,Coupling> - dof_couplings_from_component_couplings (const FiniteElement &fe, - const Table<2,Coupling> &component_couplings) - { - Assert(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())); - - const unsigned int n_dofs = fe.dofs_per_cell; - - Table<2,Coupling> dof_couplings (n_dofs, n_dofs); - - for (unsigned int i=0; i - std::vector > - dof_couplings_from_component_couplings - (const hp::FECollection &fe, - const Table<2,Coupling> &component_couplings) - { - std::vector > 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) - { - 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); - - 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); - - 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); - // 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); - // 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::level_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 (); - } - } - } - } - } - - - // implementation of the same function in namespace DoFTools for - // non-hp DoFHandlers - 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) - { - // 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(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); - - 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]); - - // 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::level_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::level_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 (); - } - } - } - } - } - } - - } - - - - - template - void - 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 types::global_dof_index 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)); - - // 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::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 - { - inline bool - check_master_dof_list (const FullMatrix &face_interpolation_matrix, - const std::vector &master_dof_list) - { - const unsigned int N = master_dof_list.size(); - - FullMatrix tmp (N,N); - for (unsigned int i=0; i p(N); - for (unsigned int i=0; i max) - { - max = std::fabs(tmp(i,j)); - r = i; - } - } - // check whether the pivot is too small. if that is the case, - // then the matrix is singular and we shouldn't use this set of - // master dofs - if (max < 1.e-12*typical_diagonal_element) - return false; - - // row interchange - if (r>j) - { - 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) - { - 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 - // (vertex, line, quad) as master as there are such dofs in fe2 - // (indices are int, because we want to avoid the 'unsigned int < 0 - // is always false warning for the cases at the bottom in 1d and - // 2d) - // - // as mentioned in the paper, it is not always easy to find a set - // of master dofs that produces an invertible matrix. to this end, - // we check in each step whether the matrix is still invertible and - // simply discard this dof if the matrix is not invertible anymore. - // - // the cases where we did have trouble in the past were with adding - // more quad dofs when Q3 and Q4 elements meet at a refined face in - // 3d (see the hp/crash_12 test that tests that we can do exactly - // this, and failed before we had code to compensate for this - // case). the other case are system elements: if we have say a Q1Q2 - // vs a Q2Q3 element, then we can't just take all master dofs on a - // line from a single base element, since the shape functions of - // that base element are independent of that of 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) - { - // make sure that we were able to find a set of master dofs - // and that the code down below didn't just reject all our - // efforts - Assert (i < fe1.dofs_per_vertex, - ExcInternalError()); - - // tentatively push this vertex dof - 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 - // well, it didn't. simply pop that dof from the list - // again and try with the next dof - 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) - { - // 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) - { - // 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); - - // 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; - } - - - - /** - * Make sure that the mask exists that determines which dofs will be - * the masters on refined faces where an fe1 and a fe2 meet. - */ - 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) - { - 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); - } - } - - - - /** - * Make sure that the given @p face_interpolation_matrix pointer - * points to a valid matrix. If the pointer is zero beforehand, - * create an entry with the correct data. If it is nonzero, don't - * touch it. - */ - template - void - ensure_existence_of_face_matrix (const FiniteElement &fe1, - 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); - } - } - - - - /** - * Same, but for subface interpolation matrices. - */ - template - void - ensure_existence_of_subface_matrix (const FiniteElement &fe1, - 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); - } - } - - - - /** - * Given the face interpolation matrix between two elements, split it - * into its master and slave parts and invert the master part as - * explained in the @ref hp_paper "hp paper". - */ - 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) - { - 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 >()); - - 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()); - - // 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); - - //TODO[WB]: We should make sure very small entries are removed after inversion - split_matrix->first.gauss_jordan (); - } - } - - - // a template that can determine statically whether a given - // DoFHandler class supports different finite element elements - template - struct DoFHandlerSupportsDifferentFEs - { - static const bool value = true; - }; - - - template - struct DoFHandlerSupportsDifferentFEs< dealii::DoFHandler > - { - static const bool value = false; - }; - - - template - struct DoFHandlerSupportsDifferentFEs< dealii::MGDoFHandler > - { - static const bool value = false; - }; - - - /** - * A function that returns how many different finite elements a dof - * handler uses. This is one for non-hp DoFHandlers and - * dof_handler.get_fe().size() for the hp-versions. - */ - template - unsigned int - n_finite_elements (const dealii::hp::DoFHandler &dof_handler) - { - return dof_handler.get_fe().size(); - } - - - template - unsigned int - n_finite_elements (const DH &) - { - return 1; - } - - - /** - * For a given face belonging to an active cell that borders to a - * more refined cell, return the fe_index of the most dominating - * finite element used on any of the face's subfaces. - */ - template - 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) - { - // 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); - } - - // see if the element on this subface is able to dominate the - // ones on all other subfaces, and if so take it - 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()); - - // 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); - } - - - - /** - * Copy constraints into a constraint matrix object. - * - * This function removes zero constraints and those, which constrain - * a DoF which was already eliminated in one of the previous steps of - * the hp hanging node procedure. - * - * It also suppresses very small entries in the constraint matrix to - * avoid making the sparsity pattern fuller than necessary. - */ - void - filter_constraints (const std::vector &master_dofs, - 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())); - - 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_dof_index, - ExcInternalError()); - for (unsigned int col=0; col!=n_master_dofs; ++col) - Assert (master_dofs[col] != numbers::invalid_dof_index, - ExcInternalError()); - - - 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.); - } - } - } - - } - - - - void - make_hp_hanging_node_constraints (const dealii::DoFHandler<1> &, - ConstraintMatrix &) - { - // nothing to do for regular dof handlers in 1d - } - - - - void - make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1> &, - ConstraintMatrix &, - dealii::internal::int2type<1>) - { - // nothing to do for regular dof handlers in 1d - } - - - void - make_hp_hanging_node_constraints (const dealii::MGDoFHandler<1> &, - ConstraintMatrix &) - { - // nothing to do for regular dof handlers in 1d - } - - - - void - make_oldstyle_hanging_node_constraints (const dealii::MGDoFHandler<1> &, - ConstraintMatrix &, - dealii::internal::int2type<1>) - { - // nothing to do for regular dof handlers in 1d - } - - - void - make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, - ConstraintMatrix &/*constraints*/) - { - // we may have to compute constraints for vertices. gotta think about - // that a bit more - - //TODO[WB]: think about what to do here... - } - - - - void - make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, - ConstraintMatrix &/*constraints*/, - dealii::internal::int2type<1>) - { - // we may have to compute constraints for vertices. gotta think about - // that a bit more - - //TODO[WB]: think about what to do here... - } - - - void - make_hp_hanging_node_constraints (const dealii::DoFHandler<1,2> &, - ConstraintMatrix &) - { - // nothing to do for regular dof handlers in 1d - } - - - - void - make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,2> &, - ConstraintMatrix &, - dealii::internal::int2type<1>) - { - // nothing to do for regular dof handlers in 1d - } - - - void - make_hp_hanging_node_constraints (const dealii::DoFHandler<1,3> &, - ConstraintMatrix &) - { - // nothing to do for regular dof handlers in 1d - } - - void - make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,3> &, - ConstraintMatrix &, - dealii::internal::int2type<1>) - { - // nothing to do for regular dof handlers in 1d - } - - -// currently not used but may be in the future: - -// void -// make_hp_hanging_node_constraints (const dealii::MDoFHandler<1,2> &, -// ConstraintMatrix &) -// { -// // nothing to do for regular -// // dof handlers in 1d -// } - - - -// void -// make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,2> &, -// ConstraintMatrix &, -// dealii::internal::int2type<1>) -// { -// // nothing to do for regular -// // dof handlers in 1d -// } - - -// void -// make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1,2> &/*dof_handler*/, -// ConstraintMatrix &/*constraints*/, -// dealii::internal::int2type<1>) -// { -// // we may have to compute -// // constraints for -// // vertices. gotta think about -// // that a bit more -// -// //TODO[WB]: think about what to do here... -// } -//#endif - - - - template - void - make_oldstyle_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints, - dealii::internal::int2type<2>) - { - const unsigned int dim = 2; - - const unsigned int spacedim = DH::space_dimension; - - std::vector dofs_on_mother; - std::vector dofs_on_children; - - // loop over all lines; only on lines there can be constraints. We do - // so by looping over all active cells and checking whether any of - // the faces are refined which can only be from the neighboring cell - // because this one is active. In that case, the face is subject to - // constraints - // - // note that even though we may visit a face twice if the neighboring - // cells are equally 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(); - 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()) - { - // in any case, faces can have at most two active fe - // indices, but here the face can have only one (namely the - // same as that from the cell we're sitting on), and 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) - 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()); - - // ok, start up the work - 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; - - 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())); - - 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()); - - // 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 - { - // this face has no children, but it could still be that it - // is shared by two cells that use a different fe index. - // check a couple of things, but ignore the 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()); - } - } - } - - - - template - void - make_oldstyle_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints, - dealii::internal::int2type<3>) - { - const unsigned int dim = 3; - - std::vector dofs_on_mother; - std::vector dofs_on_children; - - // loop over all quads; only on quads there can be constraints. We do - // so by looping over all active cells and checking whether any of - // the faces are refined which can only be from the neighboring cell - // because this one is active. In that case, the face is subject to - // constraints - // - // note that even though we may visit a face twice if the neighboring - // cells are equally 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(); - 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()) - { - // 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; - - Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, - ExcNotImplemented()); - - // in any case, faces can have at most two active fe - // indices, but here the face can have only one (namely the - // same as that from the cell we're sitting on), and each - // of the 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); - - // right now, all that is implemented is the case that both - // sides use the same fe, and not only that but also that - // all lines bounding 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()); - } - - // ok, start up the work - 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); - - //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); - - 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); - - // 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; - - // 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); - - // 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); - - // 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); - - // 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); - - // 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 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 - { - // this face has no children, but it could still be that it - // is shared by two cells that use a different fe index. - // check a couple of things, but ignore the 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()); - } - } - } - - - template - void - make_hp_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints) - { - // note: this function is going to be hard to understand if you - // haven't read the hp paper. however, we try to follow the notation - // laid out there, so go read the paper before you try to understand - // what is going on here - - const unsigned int dim = DH::dimension; - - const unsigned int spacedim = DH::space_dimension; - - - // a matrix to be used for constraints below. declared here and - // simply resized down below to avoid permanent re-allocation of - // memory - FullMatrix constraint_matrix; - - // similarly have arrays that will hold master and slave dof numbers, - // as well as a scratch array needed for the complicated case below - std::vector master_dofs; - std::vector slave_dofs; - std::vector scratch_dofs; - - // caches for the face and subface interpolation matrices between - // different (or the same) finite elements. we compute them only - // once, namely the first 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)); - Table<3,std_cxx1x::shared_ptr > > - 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 into their - // master and slave parts, and for which the master part is inverted. - // these two matrices are derived from the face 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)); - - // finally, for each pair of finite elements, have a mask that states - // which of the degrees of freedom on 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)); - - // loop over all faces - // - // note that even though we may visit a face twice if the neighboring - // cells are equally 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(); - 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()) - { - // 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; - - Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, - ExcNotImplemented()); - - // so now we've found a face of an active cell that has - // children. that means that there are hanging nodes here. - - // in any case, faces can have at most two sets of active - // fe indices, but here the face can have only one (namely - // the same as that from the cell we're sitting on), and - // 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()); - - // first find out whether we can constrain each of the - // subfaces to the mother face. in the lingo of the hp - // paper, this would be the simple case. note that we can - // short-circuit this decision if the dof_handler doesn't - // support hp at all - // - // 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: - { - // Case 1 (the simple case and the only case that can - // happen for non-hp DoFHandlers): The coarse element - // dominates the elements on the subfaces (or they are - // all the same) - // - // so we are going to constrain the DoFs on the face - // children against the DoFs on the face itself - master_dofs.resize (cell->get_fe().dofs_per_face); - - cell->face(face)->get_dof_indices (master_dofs, - cell->active_fe_index ()); - - // Now create constraint matrix for the subfaces and - // assemble it. ignore all interfaces with artificial - // cells because we can only get to such interfaces if - // the current cell is a ghost cell - for (unsigned int c=0; cface(face)->n_children(); ++c) - { - if (cell->neighbor_child_on_subface (face, c)->is_artificial()) - continue; - - const typename DH::active_face_iterator - subface = cell->face(face)->child(c); - - Assert (subface->n_active_fe_indices() == 1, - ExcInternalError()); - - const unsigned int - subface_fe_index = subface->nth_active_fe_index(0); - - // we sometime run into the situation where for - // example on one big cell we have a FE_Q(1) and on - // the subfaces we have a mixture of FE_Q(1) and - // FE_Nothing. In that case, the face domination is - // either_element_can_dominate for the whole - // collection of subfaces, but on the particular - // subface between FE_Q(1) and FE_Nothing, there - // are no constraints that we need to take 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; - - // 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); - - 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]); - - // 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: - { - // Case 2 (the "complex" case): at least one (the - // neither_... case) of the finer elements or all of - // them (the other_... case) is dominating. See the hp - // paper for a way how to deal with this situation - // - // since this is something that can only happen for hp - // dof handlers, add a check here... - Assert (DoFHandlerSupportsDifferentFEs::value == true, - ExcInternalError()); - - // we first have to find the finite element that is - // able to generate a space that all the other ones can - // be constrained to - 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]; - - // check also that it is able to constrain the mother - // face. it should be, or we wouldn't have gotten 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()); - - - // 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()); - - 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); - - // now compute the constraint matrix as the product - // between the 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); - - // then figure out the global numbers of master and - // 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 ()); - - // 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]); - - 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); - - - - // next we have to deal with the subfaces. do as - // discussed in the hp paper - for (unsigned int sf=0; - sfface(face)->n_children(); ++sf) - { - // ignore interfaces with artificial cells as well - // as interfaces 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]; - - // 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: - // there are no continuity requirements between the two - // elements. record no constraints - break; - - default: - // we shouldn't get here - 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()); - - // see if there is a neighbor that is an artificial cell. - // in that case, we're not interested in this interface. we - // test this case first since artificial cells may not have - // an active_fe_index set, etc - if (!cell->at_boundary(face) - && - cell->neighbor(face)->is_artificial()) - continue; - - // Only if there is a neighbor with a different - // active_fe_index 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::level_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: - { - // 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 ()); - - 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 - // no face dofs - if (master_dofs.size() == 0) break; - - // make sure the element constraints 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()]); - - // 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: - { - // we don't do anything here since we will come - // back to this face from the other cell, at which - // time we will fall into the first case clause - // above - break; - } - - case FiniteElementDomination::either_element_can_dominate: - { - // it appears as if neither element has any - // constraints on its neighbor. this may be because - // neither element has any DoFs on faces at all. or - // that the two elements are actually the same, - // although they happen to run under different - // fe_indices (this is what happens in - // hp/hp_hanging_nodes_01 for example). - // - // another possibility is what happens in crash_13. - // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs. - // FESystem(FE_Q(1),FE_DGQ(1)). neither of them - // dominates the other. - // - // a final possibility is that we have something like - // FESystem(FE_Q(1),FE_Q(1)) vs - // FESystem(FE_Q(1),FE_Nothing()), see - // hp/fe_nothing_18/19. - // - // in any case, the point is that it doesn't - // matter. there is nothing to do here. - break; - } - - case FiniteElementDomination::neither_element_dominates: - { - // we don't presently know what exactly to do here. - // it isn't quite clear what exactly we would have - // to do here. sit tight until someone trips over - // the following statement and see what exactly is - // going on - Assert (false, ExcNotImplemented()); - break; - } - - case FiniteElementDomination::no_requirements: - { - // nothing to do here - break; - } - - default: - // we shouldn't get here - Assert (false, ExcInternalError()); - } - } - } - } - } - - - - - template - void - make_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints) - { - // Decide whether to use the new or old make_hanging_node_constraints - // function. If all the FiniteElement or all elements in a FECollection - // support the new face constraint matrix, the new code will be used. - // Otherwise, the old implementation is used for the moment. - if (dof_handler.get_fe().hp_constraints_are_implemented ()) - internal:: - make_hp_hanging_node_constraints (dof_handler, - constraints); - else - internal:: - make_oldstyle_hanging_node_constraints (dof_handler, - constraints, - dealii::internal::int2type()); - } - - - - namespace - { - // enter constraints for periodicity into the given ConstraintMatrix object. - // this function is called when at least one of the two face iterators corresponds - // to an active object without further children - // - // @param transformation A matrix that maps degrees of freedom from one face - // to another. If the DoFs on the two faces are supposed to match exactly, then - // the matrix so provided will be the identity matrix. if face 2 is once refined - // from face 1, then the matrix needs to be the interpolation matrix from a face - // to this particular child - // - // @precondition: face_1 is supposed to be active - // - // @note As bug #82 ((http://code.google.com/p/dealii/issues/detail?id=82) and the - // corresponding testcase bits/periodicity_05 demonstrate, we can occasionally - // get into trouble if we already have the constraint x1=x2 and want to insert - // x2=x1. we avoid this by skipping an identity constraint if the opposite - // constraint already exists - template - void - set_periodicity_constraints (const FaceIterator &face_1, - const typename identity::type &face_2, - const FullMatrix &transformation, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask, - const bool face_orientation, - const bool face_flip, - const bool face_rotation) - { - static const int dim = FaceIterator::AccessorType::dimension; - static const int spacedim = FaceIterator::AccessorType::space_dimension; - - // we should be in the case where face_1 is active, i.e. has no children: - Assert (!face_1->has_children(), - ExcInternalError()); - - Assert (face_1->n_active_fe_indices() == 1, - ExcInternalError()); - - // if face_2 does have children, then we need to iterate over them - if (face_2->has_children()) - { - Assert (face_2->n_children() == GeometryInfo::max_children_per_face, - ExcNotImplemented()); - const unsigned int dofs_per_face - = face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face; - FullMatrix child_transformation (dofs_per_face, dofs_per_face); - FullMatrix subface_interpolation (dofs_per_face, dofs_per_face); - for (unsigned int c=0; cn_children(); ++c) - { - // get the interpolation matrix recursively from the one that - // interpolated from face_1 to face_2 by multiplying from the - // left with the one that interpolates from face_2 to - // its child - face_1->get_fe(face_1->nth_active_fe_index(0)) - .get_subface_interpolation_matrix (face_1->get_fe(face_1->nth_active_fe_index(0)), - c, - subface_interpolation); - subface_interpolation.mmult (child_transformation, transformation); - set_periodicity_constraints(face_1, face_2->child(c), - child_transformation, - constraint_matrix, component_mask, - face_orientation, face_flip, face_rotation); - } - } - else - // both faces are active. we need to 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")); - - const FiniteElement &fe = face_1->get_fe(face_1_index); - - Assert(component_mask.represents_n_components(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 dofs_per_face = fe.dofs_per_face; - - std::vector dofs_1(dofs_per_face); - std::vector dofs_2(dofs_per_face); - - face_1->get_dof_indices(dofs_1, face_1_index); - face_2->get_dof_indices(dofs_2, face_2_index); - - // Well, this is a hack: - // - // There is no - // face_to_face_index(face_index, - // face_orientation, - // face_flip, - // face_rotation) - // function in FiniteElementData, so we have to use - // face_to_cell_index(face_index, face - // face_orientation, - // face_flip, - // face_rotation) - // But this will give us an index on a cell - something we cannot work - // with directly. But luckily we can match them back :-] - - std::map cell_to_rotated_face_index; - - // Build up a cell to face index for face_2: - for (unsigned int i = 0; i < dofs_per_face; ++i) - { - const unsigned int cell_index = fe.face_to_cell_index(i, 0, /* It doesn't really matter, just assume - * we're on the first face... - */ - true, false, false // default orientation - ); - cell_to_rotated_face_index[cell_index] = i; - } - - // loop over all dofs on face 2 and constrain them again the ones on face 1 - for (unsigned int i=0; i - void - make_periodicity_constraints (const FaceIterator &face_1, - const typename identity::type &face_2, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask, - const bool face_orientation, - const bool face_flip, - const bool face_rotation) - { - static const int dim = FaceIterator::AccessorType::dimension; - - Assert( (dim != 1) || - (face_orientation == true && - face_flip == false && - face_rotation == false), - ExcMessage ("The supplied orientation " - "(face_orientation, face_flip, face_rotation) " - "is invalid for 1D")); - - Assert( (dim != 2) || - (face_orientation == true && - face_rotation == false), - ExcMessage ("The supplied orientation " - "(face_orientation, face_flip, face_rotation) " - "is invalid for 2D")); - - Assert(face_1 != face_2, - ExcMessage ("face_1 and face_2 are equal! Cannot constrain DoFs " - "on the very same face")); - - Assert(face_1->at_boundary() && face_2->at_boundary(), - ExcMessage ("Faces for periodicity constraints must be on the boundary")); - - - // A lookup table on how to go through the child faces depending on the - // orientation: - - static const int lookup_table_2d[2][2] = - { - // flip: - {0, 1}, // false - {1, 0}, // true - }; - - static const int lookup_table_3d[2][2][2][4] = - { - // orientation flip rotation - { { {0, 2, 1, 3}, // false false false - {2, 3, 0, 1}, // false false true - }, - { {3, 1, 2, 0}, // false true false - {1, 0, 3, 2}, // false true true - }, - }, - { { {0, 1, 2, 3}, // true false false - {1, 3, 0, 2}, // true false true - }, - { {3, 2, 1, 0}, // true true false - {2, 0, 3, 1}, // true true true - }, - }, - }; - - // In the case that both faces have children, we loop over all - // children and apply 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, - ExcNotImplemented()); - - for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - { - // Lookup the index for the second face - unsigned int j; - switch (dim) - { - case 2: - j = lookup_table_2d[face_flip][i]; - break; - case 3: - j = lookup_table_3d[face_orientation][face_flip][face_rotation][i]; - break; - default: - AssertThrow(false, ExcNotImplemented()); - } - - make_periodicity_constraints (face_1->child(i), - face_2->child(j), - constraint_matrix, - component_mask, - face_orientation, - face_flip, - face_rotation); - } - } - else - // otherwise at least one of the two faces is active and - // we need to enter the constraints - { - if (face_2->has_children() == false) - set_periodicity_constraints(face_2, face_1, - FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), - constraint_matrix, - component_mask, - face_orientation, face_flip, face_rotation); - else - set_periodicity_constraints(face_1, face_2, - FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), - constraint_matrix, - component_mask, - face_orientation, face_flip, face_rotation); - } - } - - - - template - void - make_periodicity_constraints (const DH &dof_handler, - const types::boundary_id b_id1, - const types::boundary_id b_id2, - const int direction, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask) - { - Tensor<1,DH::space_dimension> dummy; - make_periodicity_constraints (dof_handler, - b_id1, - b_id2, - direction, - dummy, - constraint_matrix, - component_mask); - } - - - - template - void - make_periodicity_constraints (const DH &dof_handler, - const types::boundary_id b_id1, - const types::boundary_id b_id2, - const int direction, - dealii::Tensor<1,DH::space_dimension> &offset, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask) - { - static const int space_dim = DH::space_dimension; - Assert (0<=direction && direction PTRIA; - const PTRIA *ptria_p = dynamic_cast (&dof_handler.get_tria()); - Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1), - ExcMessage ("This function can not be used with distributed triangulations." - "See the documentation for more information.")); - } -#endif - - Assert (b_id1 != b_id2, - ExcMessage ("The boundary indicators b_id1 and b_id2 must be" - "different to denote different boundaries.")); - - typedef typename DH::face_iterator FaceIterator; - typedef std::map > > FaceMap; - - // Collect matching periodic cells on the coarsest level: - FaceMap matched_cells = - GridTools::collect_periodic_face_pairs(dof_handler, - b_id1, b_id2, - direction, offset); - - // And apply the low level make_periodicity_constraints function to - // every matching pair: - for (typename FaceMap::iterator it = matched_cells.begin(); - it != matched_cells.end(); ++it) - { - typedef typename DH::face_iterator FaceIterator; - const FaceIterator &face_1 = it->first; - const FaceIterator &face_2 = it->second.first; - const std::bitset<3> &orientation = it->second.second; - - Assert(face_1->at_boundary() && face_2->at_boundary(), - ExcInternalError()); - - Assert (face_1->boundary_indicator() == b_id1 && - face_2->boundary_indicator() == b_id2, - ExcInternalError()); - - Assert (face_1 != face_2, - ExcInternalError()); - - make_periodicity_constraints(face_1, - face_2, - constraint_matrix, - component_mask, - orientation[0], - orientation[1], - orientation[2]); - } - } - - - - template - void - make_periodicity_constraints (const DH &dof_handler, - const types::boundary_id b_id, - const int direction, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask) - { - Tensor<1,DH::space_dimension> dummy; - make_periodicity_constraints (dof_handler, - b_id, - direction, - dummy, - constraint_matrix, - component_mask); - } - - - - template - void - make_periodicity_constraints (const DH &dof_handler, - const types::boundary_id b_id, - const int direction, - dealii::Tensor<1,DH::space_dimension> &offset, - dealii::ConstraintMatrix &constraint_matrix, - const ComponentMask &component_mask) - { - static const int dim = DH::dimension; - static const int space_dim = DH::space_dimension; - - Assert (0<=direction && direction PTRIA; - const PTRIA *ptria_p = dynamic_cast (&dof_handler.get_tria()); - Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1), - ExcMessage ("This function can not be used with distributed triangulations." - "See the documentation for more information.")); - } -#endif - - typedef typename DH::face_iterator FaceIterator; - typedef std::map FaceMap; - - // Collect matching periodic cells on the coarsest level: - FaceMap matched_cells = - GridTools::collect_periodic_face_pairs(dof_handler, - b_id, - direction, offset); - - // And apply the low level make_periodicity_constraints function to - // every matching pair: - for (typename FaceMap::iterator it = matched_cells.begin(); - it != matched_cells.end(); ++it) - { - typedef typename DH::face_iterator FaceIterator; - const FaceIterator &face_1 = it->first; - const FaceIterator &face_2 = it->second; - - Assert(face_1->at_boundary() && face_2->at_boundary(), - ExcInternalError()); - - Assert (face_1->boundary_indicator() == b_id && - face_2->boundary_indicator() == b_id, - ExcInternalError()); - - Assert (face_1 != face_2, - ExcInternalError()); - - make_periodicity_constraints(face_1, - face_2, - constraint_matrix, - component_mask - /* standard orientation */); - } - } - - - namespace internal { // return an array that for each dof on the reference cell @@ -4719,677 +1512,6 @@ namespace DoFTools - - namespace internal - { - namespace - { - /** - * This is a function that is called by the _2 function and that - * operates on a range of cells only. It is used to split up the - * whole range of cells into chunks which are then worked on in - * parallel, if multithreading is available. - */ - 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) - { - // aliases to the finite elements used by the dof handlers: - const FiniteElement &coarse_fe = coarse_grid.get_fe(); - - // for each cell on the parameter grid: find out which degrees of - // freedom on the fine grid correspond in which way to the degrees - // of freedom on the parameter grid - // - // since for continuous FEs some dofs exist on more than one cell, - // we have to track which ones were already visited. the problem is - // that if we visit a dof first on one cell and compute its weight - // with respect to some global dofs to be non-zero, and later visit - // the dof again on another cell and (since we are on another cell) - // recompute the weights with respect to the same dofs as above to - // be zero now, we have to preserve them. we therefore overwrite - // all weights if they are nonzero and do not enforce zero weights - // since that might be only due to the fact that we are on another - // cell. - // - // example: - // coarse grid - // | | | - // *-----*-----* - // | cell|cell | - // | 1 | 2 | - // | | | - // 0-----1-----* - // - // fine grid - // | | | | | - // *--*--*--*--* - // | | | | | - // *--*--*--*--* - // | | | | | - // *--x--y--*--* - // - // when on cell 1, we compute the weights of dof 'x' to be 1/2 from - // parameter dofs 0 and 1, respectively. however, when later we are - // on cell 2, we again compute the prolongation of shape function 1 - // restricted to cell 2 to the globla grid and find that the weight - // of global dof 'x' now is zero. however, we should not overwrite - // the old value. - // - // we therefore always only set nonzero values. why adding up is - // not useful: dof 'y' would get weight 1 from parameter dof 1 on - // both cells 1 and 2, but the correct weight is nevertheless only - // 1. - - // vector to hold the representation of a single degree of freedom - // on the coarse grid (for the selected fe) on the fine grid - const types::global_dof_index 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); - - 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); - - // 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); - // now that we've got the global representation of each - // parameter dof, we've only got to clobber the non-zero - // entries in that vector and store the result - // - // what we have learned: if entry @p{i} of the global - // vector holds the value @p{v[i]}, then this is the - // weight with which the present dof contributes to - // @p{i}. there may be several such @p{i}s and their - // weights' sum should be one. Then, @p{v[i]} should be - // equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the - // values of the degrees of freedom on the coarse grid. - // we can thus compute constraints which link the degrees - // of freedom @p{v[i]} on the fine grid to those on the - // coarse grid, @p{p[j]}. Now to use these as real - // constraints, rather than as additional equations, we - // have to identify representants among the @p{i} for - // each @p{j}. this will be done by simply taking the - // first @p{i} for which @p{w_{ij}==1}. - // - // guard modification of the weights array by a Mutex. - // since it should happen rather rarely that there are - // several threads operating on different intergrid - // weights, have only one mutex for all of them - static Threads::Mutex mutex; - Threads::Mutex::ScopedLock lock (mutex); - for (types::global_dof_index 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) - { - // 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); - - // 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 - 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) - { - // aliases to the finite elements used by the dof handlers: - const FiniteElement &coarse_fe = coarse_grid.get_fe(), - &fine_fe = fine_grid.get_fe(); - - // global numbers of dofs - const types::global_dof_index 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; - - // 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; - - - // 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()); - - // 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 ()); - - - // check whether component numbers are valid - AssertIndexRange (coarse_component,coarse_fe.n_components()); - 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()); - -#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()); -#endif - - /* - * From here on: the term `parameter' refers to the selected - * component on the coarse grid and its analogon on the fine grid. - * The naming of variables containing this term is due to the fact - * that `selected_component' is longer, but also due to the fact - * that the code of this function was initially written for a - * program where the component which we wanted to match between - * grids was actually the `parameter' variable. - * - * Likewise, the terms `parameter grid' and `state grid' refer to - * the coarse and fine grids, respectively. - * - * Changing the names of variables would in principle be a good - * idea, but would not make things simpler and would be another - * source of errors. If anyone feels like doing so: patches would - * be welcome! - */ - - - - // set up vectors of cell-local data; each vector represents one - // 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)); - // 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); - - 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()); - }; -#endif - - - return n_parameters_on_fine_grid; - } - - - } - } - - - - 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) - { - // store the weights with which a dof on the parameter grid contributes - // to a dof on the fine grid. see the long doc below for more info - // - // allocate as many rows as there are parameter dofs on the coarse grid - // and as many columns as there are parameter dofs on the fine grid. - // - // weight_mapping is used to map the global (fine grid) parameter dof - // indices to the columns - // - // in the original implementation, the weights array was actually of - // FullMatrix type. this wasted huge amounts of memory, but was - // fast. nonetheless, since the memory consumption was quadratic in the - // number of degrees of freedom, this was not very practical, so we now - // use a vector of rows of the matrix, and in each row a vector of - // pairs (colnum,value). this seems like the best tradeoff between - // memory and speed, as it is now linear in memory and still fast - // enough. - // - // to save some memory and since the weights are usually (negative) - // powers of 2, we choose the value type of the matrix to be @p{float} - // rather than @p{double}. - std::vector > weights; - - // this is this mapping. there is one entry for each dof on the fine - // grid; if it is a parameter dof, then its value is the column in - // weights for that parameter dof, if it is any other dof, then its - // value is -1, indicating an error - std::vector weight_mapping; - - 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); - - // global numbers of dofs - const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(), - n_fine_dofs = fine_grid.n_dofs(); - - - // get an array in which we store which dof on the coarse grid is a - // parameter and which is not - 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, ComponentMask(mask), coarse_dof_is_parameter); - } - - // now we know that the weights in each row constitute a constraint. - // enter this into the constraints object - // - // first task: for each parameter dof on the parameter grid, find a - // representant on the fine, global grid. this is possible since we use - // conforming finite element. we take this representant to be the first - // element in this row with weight identical to one. the representant - // will become an unconstrained degree of freedom, while all others - // will be constrained to this dof (and possibly others) - std::vector representants(n_coarse_dofs, numbers::invalid_dof_index); - for (types::global_dof_index parameter_dof=0; parameter_dof 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 types::global_dof_index 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 - types::global_dof_index 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; - } - 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()); - }; - - - - // note for people that want to optimize this function: the largest - // part of the computing time is spent in the following, rather - // innocent block of code. basically, it must be the - // ConstraintMatrix::add_entry call which takes the bulk of the time, - // but it is not known to the author how to make it faster... - std::vector > constraint_line; - for (types::global_dof_index global_dof=0; global_dof 0, ExcInternalError()); - std::map::const_iterator - col_entry = weights[0].end(); - for (; first_used_rowsecond == 1) && - (representants[first_used_row] == global_dof)) - // dof unconstrained or constrained to itself (in case this - // cell is mapped to itself, rather than to children of - // itself) - continue; - } - - - // otherwise enter all constraints - constraints.add_line (global_dof); - - constraint_line.clear (); - for (types::global_dof_index 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); - }; - } - - - - 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) - { - // store the weights with which a dof on the parameter grid contributes - // to a dof on the fine grid. see the long doc below for more info - // - // allocate as many rows as there are parameter dofs on the coarse grid - // and as many columns as there are parameter dofs on the fine grid. - // - // weight_mapping is used to map the global (fine grid) parameter dof - // indices to the columns - // - // in the original implementation, the weights array was actually of - // FullMatrix type. this wasted huge amounts of memory, but was - // fast. nonetheless, since the memory consumption was quadratic in the - // number of degrees of freedom, this was not very practical, so we now - // use a vector of rows of the matrix, and in each row a vector of - // pairs (colnum,value). this seems like the best tradeoff between - // memory and speed, as it is now linear in memory and still fast - // enough. - // - // to save some memory and since the weights are usually (negative) - // powers of 2, we choose the value type of the matrix to be @p{float} - // rather than @p{double}. - std::vector > weights; - - // this is this mapping. there is one entry for each dof on the fine - // grid; if it is a parameter dof, then its value is the column in - // weights for that parameter dof, if it is any other dof, then its - // value is -1, indicating an error - 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); - - // now compute the requested representation - const types::global_dof_index n_global_parm_dofs - = std::count_if (weight_mapping.begin(), weight_mapping.end(), - std::bind2nd (std::not_equal_to (), numbers::invalid_dof_index)); - - // first construct the inverse mapping of weight_mapping - std::vector inverse_weight_mapping (n_global_parm_dofs, - DoFHandler::invalid_dof_index); - for (types::global_dof_index i=0; i::invalid_dof_index), - ExcInternalError()); - - inverse_weight_mapping[parameter_dof] = i; - }; - }; - - // next copy over weights array and replace respective numbers - const types::global_dof_index n_rows = weight_mapping.size(); - - transfer_representation.clear (); - transfer_representation.resize (n_rows); - - const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(); - for (types::global_dof_index i=0; i::const_iterator j = weights[i].begin(); - for (; j!=weights[i].end(); ++j) - { - const types::global_dof_index p = inverse_weight_mapping[j->first]; - Assert (psecond; - }; - }; - } - - - template void map_dof_to_boundary_indices (const DH &dof_handler, @@ -5689,87 +1811,6 @@ namespace DoFTools - template class DH> - void - make_zero_boundary_constraints (const DH &dof, - const types::boundary_id boundary_indicator, - ConstraintMatrix &zero_boundary_constraints, - const ComponentMask &component_mask) - { - Assert (component_mask.represents_n_components(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); - - Assert (component_mask.n_selected_components(n_components) > 0, - VectorTools::ExcNoComponentSelected()); - - // a field to store the indices - std::vector face_dofs; - face_dofs.reserve (max_dofs_per_face(dof)); - - typename DH::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - if (!cell->is_artificial()) - for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; - ++face_no) - { - const FiniteElement &fe = cell->get_fe(); - - const typename DH::face_iterator face = cell->face(face_no); - - // if face is on the boundary and satisfies the correct - // boundary id property - if (face->at_boundary () - && - ((boundary_indicator == numbers::invalid_boundary_id) - || - (face->boundary_indicator() == boundary_indicator))) - { - // 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()); - - // enter those dofs into the list that match the component - // signature. - for (unsigned int i=0; iget_fe().get_nonzero_components (i); - bool nonzero = false; - for (unsigned int c=0; c class DH> - void - make_zero_boundary_constraints (const DH &dof, - ConstraintMatrix &zero_boundary_constraints, - const ComponentMask &component_mask) - { - make_zero_boundary_constraints(dof, numbers::invalid_boundary_id, - zero_boundary_constraints, component_mask); - } - - template void make_cell_patches( Sparsity &block_list, @@ -6029,7 +2070,9 @@ namespace DoFTools } } } -} + + +} // end of namespace DoFTools diff --git a/deal.II/source/dofs/dof_tools.inst.in b/deal.II/source/dofs/dof_tools.inst.in index 17d1beaa6c..f4121a9d63 100644 --- a/deal.II/source/dofs/dof_tools.inst.in +++ b/deal.II/source/dofs/dof_tools.inst.in @@ -15,394 +15,21 @@ // --------------------------------------------------------------------- - for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS) { - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler &dof, - SP &sparsity, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler &dof, - SP &sparsity, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const MGDoFHandler &dof, - SP &sparsity, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler&, - const Table<2,Coupling>&, - SP &, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler&, - const Table<2,Coupling>&, - SP &, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const MGDoFHandler&, - const Table<2,Coupling>&, - SP &, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler &dof_row, - const DoFHandler &dof_col, - SP &sparsity); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler &dof_row, - const hp::DoFHandler &dof_col, - SP &sparsity); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const DoFHandler& dof, - const std::vector &, - SP &); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const hp::DoFHandler& dof, - const std::vector &, - SP &); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const DoFHandler& dof, - const FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const hp::DoFHandler& dof, - const FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - -#if deal_II_dimension < 3 - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const hp::DoFHandler& dof, - const FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - #endif - - template void - DoFTools::make_flux_sparsity_pattern,SP> - (const DoFHandler &dof, - SP &sparsity); - - template void - DoFTools::make_flux_sparsity_pattern,SP> - (const hp::DoFHandler &dof, - SP &sparsity); - - template void - DoFTools::make_cell_patches,SP> - (SP&, const DoFHandler&, const unsigned int, const std::vector&, types::global_dof_index); - - template void - DoFTools::make_cell_patches,SP> - (SP&, const MGDoFHandler&, const unsigned int, const std::vector&, types::global_dof_index); - -#if deal_II_dimension > 1 - - template void - DoFTools::make_flux_sparsity_pattern,SP> - (const DoFHandler &dof, - SP &, - const Table<2,Coupling>&, - const Table<2,Coupling>&); - - template void - DoFTools::make_flux_sparsity_pattern,SP> - (const DoFHandler &dof, - SP &sparsity, - const ConstraintMatrix &constraints, - const bool, const unsigned int); - - template void - DoFTools::make_flux_sparsity_pattern,SP> - (const hp::DoFHandler &dof, - SP &sparsity, - 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 - - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler &dof, - SP &sparsity, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler &dof, - SP &sparsity, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler&, - const Table<2,Coupling>&, - SP &, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler&, - const Table<2,Coupling>&, - SP &, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler &dof_row, - const DoFHandler &dof_col, - SP &sparsity); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler &dof_row, - const hp::DoFHandler &dof_col, - SP &sparsity); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const DoFHandler& dof, - const std::vector &, - SP &); - - //template void - //DoFTools::make_boundary_sparsity_pattern,SP> - //(const hp::DoFHandler& dof, - // const std::vector &, - // SP &); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const DoFHandler& dof, - const FunctionMap::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - - //template void - //DoFTools::make_boundary_sparsity_pattern,SP> - //(const hp::DoFHandler& dof, - // const FunctionMap::type &boundary_indicators, - // const std::vector &dof_to_boundary_mapping, - // SP &sparsity); - -#endif - - -#if deal_II_dimension == 3 - - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler<1,3> &dof, - SP &sparsity, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler<1,3> &dof, - SP &sparsity, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler<1,3>&, - const Table<2,Coupling>&, - SP &, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler<1,3>&, - const Table<2,Coupling>&, - SP &, - const ConstraintMatrix &, - const bool, - const unsigned int); - - template void - DoFTools::make_sparsity_pattern, SP> - (const DoFHandler<1,3> &dof_row, - const DoFHandler<1,3> &dof_col, - SP &sparsity); - - template void - DoFTools::make_sparsity_pattern, SP> - (const hp::DoFHandler<1,3> &dof_row, - const hp::DoFHandler<1,3> &dof_col, - SP &sparsity); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const DoFHandler<1,3>& dof, - const std::vector &, - SP &); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const hp::DoFHandler<1,3>& dof, - const std::vector &, - SP &); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const DoFHandler<1,3>& dof, - const FunctionMap<3>::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - - template void - DoFTools::make_boundary_sparsity_pattern,SP> - (const hp::DoFHandler<1,3>& dof, - const FunctionMap<3>::type &boundary_indicators, - const std::vector &dof_to_boundary_mapping, - SP &sparsity); - -#endif - - } - + template void + DoFTools::make_cell_patches,SP> + (SP&, const DoFHandler&, const unsigned int, const std::vector&, types::global_dof_index); - - -for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS) -{ - template - void - DoFTools::make_hanging_node_constraints (const DH &dof_handler, - ConstraintMatrix &constraints); -} - - - -for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS) -{ - -#if deal_II_dimension != 1 - template - void - DoFTools::make_periodicity_constraints (const DH::face_iterator &, - const DH::face_iterator &, - dealii::ConstraintMatrix &, - const ComponentMask &, - bool, bool, bool); - - template - void - DoFTools::make_periodicity_constraints(const DH &, - const types::boundary_id, - const types::boundary_id, - const int, - dealii::ConstraintMatrix &, - const ComponentMask &); - - template - void - DoFTools::make_periodicity_constraints(const DH &, - const types::boundary_id, - const types::boundary_id, - const int, - dealii::Tensor<1,DH::space_dimension> &, - dealii::ConstraintMatrix &, - const ComponentMask &); - - template - void - DoFTools::make_periodicity_constraints(const DH &, - const types::boundary_id, - const int, - dealii::ConstraintMatrix &, - const ComponentMask &); - - template - void - DoFTools::make_periodicity_constraints(const DH &, - const types::boundary_id, - const int, - dealii::Tensor<1,DH::space_dimension> &, - dealii::ConstraintMatrix &, - const ComponentMask &); -#endif + template void + DoFTools::make_cell_patches,SP> + (SP&, const MGDoFHandler&, const unsigned int, const std::vector&, types::global_dof_index); } -for (deal_II_dimension : DIMENSIONS) -{ -#if deal_II_dimension != 1 - template - void - DoFTools::make_hanging_node_constraints (const MGDoFHandler &dof_handler, - ConstraintMatrix &constraints); - -#endif - } for (deal_II_dimension : DIMENSIONS) { - template - Table<2,DoFTools::Coupling> - DoFTools::dof_couplings_from_component_couplings - (const FiniteElement &fe, - const Table<2,DoFTools::Coupling> &component_couplings); - template void DoFTools::make_vertex_patches (SparsityPattern&, const DoFHandler&, unsigned int, bool, bool, bool, bool); @@ -427,22 +54,6 @@ for (deal_II_dimension : DIMENSIONS) void DoFTools::make_child_patches(SparsityPattern&, const MGDoFHandler&, unsigned int, bool, bool); -#if deal_II_dimension < 3 -template -void -DoFTools:: -make_hanging_node_constraints (const DoFHandler &dof_handler, - ConstraintMatrix &constraints); -#endif - -#if deal_II_dimension == 3 -template -void -DoFTools:: -make_hanging_node_constraints (const DoFHandler<1,3> &dof_handler, - ConstraintMatrix &constraints); -#endif - // TODO: can cleanup a bit more to fit into the scheme used above @@ -912,23 +523,6 @@ DoFTools::count_dofs_per_component ( const DoFHandler&, std::vector&, std::vector); -template -void -DoFTools::compute_intergrid_constraints ( - const DoFHandler &, const unsigned int, - const DoFHandler &, const unsigned int, - const InterGridMap > &, - ConstraintMatrix&); - -template -void -DoFTools::compute_intergrid_transfer_representation -(const DoFHandler &, const unsigned int, - const DoFHandler &, const unsigned int, - const InterGridMap > &, - std::vector > &); - - template void DoFTools::map_dof_to_boundary_indices > @@ -1089,36 +683,6 @@ DoFTools::convert_couplings_to_blocks ( const hp::DoFHandler&, const Table<2, Coupling>&, std::vector >&); -template -void -DoFTools::make_zero_boundary_constraints -(const DoFHandler &, - ConstraintMatrix &, - const ComponentMask &); - -template -void -DoFTools::make_zero_boundary_constraints -(const DoFHandler &, - const types::boundary_id , - ConstraintMatrix &, - const ComponentMask &); - -template -void -DoFTools::make_zero_boundary_constraints -(const hp::DoFHandler &, - ConstraintMatrix &, - const ComponentMask &); - -template -void -DoFTools::make_zero_boundary_constraints -(const hp::DoFHandler &, - const types::boundary_id , - ConstraintMatrix &, - const ComponentMask &); - #if deal_II_dimension < 3 diff --git a/deal.II/source/dofs/dof_tools_constraints.cc b/deal.II/source/dofs/dof_tools_constraints.cc new file mode 100644 index 0000000000..6e21770c30 --- /dev/null +++ b/deal.II/source/dofs/dof_tools_constraints.cc @@ -0,0 +1,2891 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1999 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + + +namespace DoFTools +{ + namespace internal + { + namespace + { + inline bool + check_master_dof_list (const FullMatrix &face_interpolation_matrix, + const std::vector &master_dof_list) + { + const unsigned int N = master_dof_list.size(); + + FullMatrix tmp (N,N); + for (unsigned int i=0; i p(N); + for (unsigned int i=0; i max) + { + max = std::fabs(tmp(i,j)); + r = i; + } + } + // check whether the pivot is too small. if that is the case, + // then the matrix is singular and we shouldn't use this set of + // master dofs + if (max < 1.e-12*typical_diagonal_element) + return false; + + // row interchange + if (r>j) + { + 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) + { + 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 + // (vertex, line, quad) as master as there are such dofs in fe2 + // (indices are int, because we want to avoid the 'unsigned int < 0 + // is always false warning for the cases at the bottom in 1d and + // 2d) + // + // as mentioned in the paper, it is not always easy to find a set + // of master dofs that produces an invertible matrix. to this end, + // we check in each step whether the matrix is still invertible and + // simply discard this dof if the matrix is not invertible anymore. + // + // the cases where we did have trouble in the past were with adding + // more quad dofs when Q3 and Q4 elements meet at a refined face in + // 3d (see the hp/crash_12 test that tests that we can do exactly + // this, and failed before we had code to compensate for this + // case). the other case are system elements: if we have say a Q1Q2 + // vs a Q2Q3 element, then we can't just take all master dofs on a + // line from a single base element, since the shape functions of + // that base element are independent of that of 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) + { + // make sure that we were able to find a set of master dofs + // and that the code down below didn't just reject all our + // efforts + Assert (i < fe1.dofs_per_vertex, + ExcInternalError()); + + // tentatively push this vertex dof + 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 + // well, it didn't. simply pop that dof from the list + // again and try with the next dof + 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) + { + // 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) + { + // 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); + + // 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; + } + + + + /** + * Make sure that the mask exists that determines which dofs will be + * the masters on refined faces where an fe1 and a fe2 meet. + */ + 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) + { + 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); + } + } + + + + /** + * Make sure that the given @p face_interpolation_matrix pointer + * points to a valid matrix. If the pointer is zero beforehand, + * create an entry with the correct data. If it is nonzero, don't + * touch it. + */ + template + void + ensure_existence_of_face_matrix (const FiniteElement &fe1, + 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); + } + } + + + + /** + * Same, but for subface interpolation matrices. + */ + template + void + ensure_existence_of_subface_matrix (const FiniteElement &fe1, + 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); + } + } + + + + /** + * Given the face interpolation matrix between two elements, split it + * into its master and slave parts and invert the master part as + * explained in the @ref hp_paper "hp paper". + */ + 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) + { + 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 >()); + + 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()); + + // 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); + + //TODO[WB]: We should make sure very small entries are removed after inversion + split_matrix->first.gauss_jordan (); + } + } + + + // a template that can determine statically whether a given + // DoFHandler class supports different finite element elements + template + struct DoFHandlerSupportsDifferentFEs + { + static const bool value = true; + }; + + + template + struct DoFHandlerSupportsDifferentFEs< dealii::DoFHandler > + { + static const bool value = false; + }; + + + /** + * A function that returns how many different finite elements a dof + * handler uses. This is one for non-hp DoFHandlers and + * dof_handler.get_fe().size() for the hp-versions. + */ + template + unsigned int + n_finite_elements (const dealii::hp::DoFHandler &dof_handler) + { + return dof_handler.get_fe().size(); + } + + + template + unsigned int + n_finite_elements (const DH &) + { + return 1; + } + + + /** + * For a given face belonging to an active cell that borders to a + * more refined cell, return the fe_index of the most dominating + * finite element used on any of the face's subfaces. + */ + template + 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) + { + // 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); + } + + // see if the element on this subface is able to dominate the + // ones on all other subfaces, and if so take it + 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()); + + // 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); + } + + + + /** + * Copy constraints into a constraint matrix object. + * + * This function removes zero constraints and those, which constrain + * a DoF which was already eliminated in one of the previous steps of + * the hp hanging node procedure. + * + * It also suppresses very small entries in the constraint matrix to + * avoid making the sparsity pattern fuller than necessary. + */ + void + filter_constraints (const std::vector &master_dofs, + 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())); + + 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_dof_index, + ExcInternalError()); + for (unsigned int col=0; col!=n_master_dofs; ++col) + Assert (master_dofs[col] != numbers::invalid_dof_index, + ExcInternalError()); + + + 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.); + } + } + } + + } + + + + void + make_hp_hanging_node_constraints (const dealii::DoFHandler<1> &, + ConstraintMatrix &) + { + // nothing to do for regular dof handlers in 1d + } + + + + void + make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1> &, + ConstraintMatrix &, + dealii::internal::int2type<1>) + { + // nothing to do for regular dof handlers in 1d + } + + + void + make_hp_hanging_node_constraints (const dealii::MGDoFHandler<1> &, + ConstraintMatrix &) + { + // nothing to do for regular dof handlers in 1d + } + + + + void + make_oldstyle_hanging_node_constraints (const dealii::MGDoFHandler<1> &, + ConstraintMatrix &, + dealii::internal::int2type<1>) + { + // nothing to do for regular dof handlers in 1d + } + + + void + make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, + ConstraintMatrix &/*constraints*/) + { + // we may have to compute constraints for vertices. gotta think about + // that a bit more + + //TODO[WB]: think about what to do here... + } + + + + void + make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, + ConstraintMatrix &/*constraints*/, + dealii::internal::int2type<1>) + { + // we may have to compute constraints for vertices. gotta think about + // that a bit more + + //TODO[WB]: think about what to do here... + } + + + void + make_hp_hanging_node_constraints (const dealii::DoFHandler<1,2> &, + ConstraintMatrix &) + { + // nothing to do for regular dof handlers in 1d + } + + + + void + make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,2> &, + ConstraintMatrix &, + dealii::internal::int2type<1>) + { + // nothing to do for regular dof handlers in 1d + } + + + void + make_hp_hanging_node_constraints (const dealii::DoFHandler<1,3> &, + ConstraintMatrix &) + { + // nothing to do for regular dof handlers in 1d + } + + void + make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,3> &, + ConstraintMatrix &, + dealii::internal::int2type<1>) + { + // nothing to do for regular dof handlers in 1d + } + + +// currently not used but may be in the future: + +// void +// make_hp_hanging_node_constraints (const dealii::MDoFHandler<1,2> &, +// ConstraintMatrix &) +// { +// // nothing to do for regular +// // dof handlers in 1d +// } + + + +// void +// make_oldstyle_hanging_node_constraints (const dealii::DoFHandler<1,2> &, +// ConstraintMatrix &, +// dealii::internal::int2type<1>) +// { +// // nothing to do for regular +// // dof handlers in 1d +// } + + +// void +// make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1,2> &/*dof_handler*/, +// ConstraintMatrix &/*constraints*/, +// dealii::internal::int2type<1>) +// { +// // we may have to compute +// // constraints for +// // vertices. gotta think about +// // that a bit more +// +// //TODO[WB]: think about what to do here... +// } +//#endif + + + + template + void + make_oldstyle_hanging_node_constraints (const DH &dof_handler, + ConstraintMatrix &constraints, + dealii::internal::int2type<2>) + { + const unsigned int dim = 2; + + const unsigned int spacedim = DH::space_dimension; + + std::vector dofs_on_mother; + std::vector dofs_on_children; + + // loop over all lines; only on lines there can be constraints. We do + // so by looping over all active cells and checking whether any of + // the faces are refined which can only be from the neighboring cell + // because this one is active. In that case, the face is subject to + // constraints + // + // note that even though we may visit a face twice if the neighboring + // cells are equally 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(); + 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()) + { + // in any case, faces can have at most two active fe + // indices, but here the face can have only one (namely the + // same as that from the cell we're sitting on), and 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) + 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()); + + // ok, start up the work + 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; + + 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())); + + 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()); + + // 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 + { + // this face has no children, but it could still be that it + // is shared by two cells that use a different fe index. + // check a couple of things, but ignore the 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()); + } + } + } + + + + template + void + make_oldstyle_hanging_node_constraints (const DH &dof_handler, + ConstraintMatrix &constraints, + dealii::internal::int2type<3>) + { + const unsigned int dim = 3; + + std::vector dofs_on_mother; + std::vector dofs_on_children; + + // loop over all quads; only on quads there can be constraints. We do + // so by looping over all active cells and checking whether any of + // the faces are refined which can only be from the neighboring cell + // because this one is active. In that case, the face is subject to + // constraints + // + // note that even though we may visit a face twice if the neighboring + // cells are equally 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(); + 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()) + { + // 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; + + Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, + ExcNotImplemented()); + + // in any case, faces can have at most two active fe + // indices, but here the face can have only one (namely the + // same as that from the cell we're sitting on), and each + // of the 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); + + // right now, all that is implemented is the case that both + // sides use the same fe, and not only that but also that + // all lines bounding 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()); + } + + // ok, start up the work + 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); + + //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); + + 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); + + // 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; + + // 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); + + // 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); + + // 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); + + // 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); + + // 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 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 + { + // this face has no children, but it could still be that it + // is shared by two cells that use a different fe index. + // check a couple of things, but ignore the 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()); + } + } + } + + + template + void + make_hp_hanging_node_constraints (const DH &dof_handler, + ConstraintMatrix &constraints) + { + // note: this function is going to be hard to understand if you + // haven't read the hp paper. however, we try to follow the notation + // laid out there, so go read the paper before you try to understand + // what is going on here + + const unsigned int dim = DH::dimension; + + const unsigned int spacedim = DH::space_dimension; + + + // a matrix to be used for constraints below. declared here and + // simply resized down below to avoid permanent re-allocation of + // memory + FullMatrix constraint_matrix; + + // similarly have arrays that will hold master and slave dof numbers, + // as well as a scratch array needed for the complicated case below + std::vector master_dofs; + std::vector slave_dofs; + std::vector scratch_dofs; + + // caches for the face and subface interpolation matrices between + // different (or the same) finite elements. we compute them only + // once, namely the first 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)); + Table<3,std_cxx1x::shared_ptr > > + 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 into their + // master and slave parts, and for which the master part is inverted. + // these two matrices are derived from the face 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)); + + // finally, for each pair of finite elements, have a mask that states + // which of the degrees of freedom on 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)); + + // loop over all faces + // + // note that even though we may visit a face twice if the neighboring + // cells are equally 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(); + 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()) + { + // 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; + + Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, + ExcNotImplemented()); + + // so now we've found a face of an active cell that has + // children. that means that there are hanging nodes here. + + // in any case, faces can have at most two sets of active + // fe indices, but here the face can have only one (namely + // the same as that from the cell we're sitting on), and + // 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()); + + // first find out whether we can constrain each of the + // subfaces to the mother face. in the lingo of the hp + // paper, this would be the simple case. note that we can + // short-circuit this decision if the dof_handler doesn't + // support hp at all + // + // 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: + { + // Case 1 (the simple case and the only case that can + // happen for non-hp DoFHandlers): The coarse element + // dominates the elements on the subfaces (or they are + // all the same) + // + // so we are going to constrain the DoFs on the face + // children against the DoFs on the face itself + master_dofs.resize (cell->get_fe().dofs_per_face); + + cell->face(face)->get_dof_indices (master_dofs, + cell->active_fe_index ()); + + // Now create constraint matrix for the subfaces and + // assemble it. ignore all interfaces with artificial + // cells because we can only get to such interfaces if + // the current cell is a ghost cell + for (unsigned int c=0; cface(face)->n_children(); ++c) + { + if (cell->neighbor_child_on_subface (face, c)->is_artificial()) + continue; + + const typename DH::active_face_iterator + subface = cell->face(face)->child(c); + + Assert (subface->n_active_fe_indices() == 1, + ExcInternalError()); + + const unsigned int + subface_fe_index = subface->nth_active_fe_index(0); + + // we sometime run into the situation where for + // example on one big cell we have a FE_Q(1) and on + // the subfaces we have a mixture of FE_Q(1) and + // FE_Nothing. In that case, the face domination is + // either_element_can_dominate for the whole + // collection of subfaces, but on the particular + // subface between FE_Q(1) and FE_Nothing, there + // are no constraints that we need to take 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; + + // 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); + + 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]); + + // 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: + { + // Case 2 (the "complex" case): at least one (the + // neither_... case) of the finer elements or all of + // them (the other_... case) is dominating. See the hp + // paper for a way how to deal with this situation + // + // since this is something that can only happen for hp + // dof handlers, add a check here... + Assert (DoFHandlerSupportsDifferentFEs::value == true, + ExcInternalError()); + + // we first have to find the finite element that is + // able to generate a space that all the other ones can + // be constrained to + 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]; + + // check also that it is able to constrain the mother + // face. it should be, or we wouldn't have gotten 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()); + + + // 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()); + + 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); + + // now compute the constraint matrix as the product + // between the 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); + + // then figure out the global numbers of master and + // 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 ()); + + // 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]); + + 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); + + + + // next we have to deal with the subfaces. do as + // discussed in the hp paper + for (unsigned int sf=0; + sfface(face)->n_children(); ++sf) + { + // ignore interfaces with artificial cells as well + // as interfaces 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]; + + // 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: + // there are no continuity requirements between the two + // elements. record no constraints + break; + + default: + // we shouldn't get here + 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()); + + // see if there is a neighbor that is an artificial cell. + // in that case, we're not interested in this interface. we + // test this case first since artificial cells may not have + // an active_fe_index set, etc + if (!cell->at_boundary(face) + && + cell->neighbor(face)->is_artificial()) + continue; + + // Only if there is a neighbor with a different + // active_fe_index 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::level_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: + { + // 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 ()); + + 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 + // no face dofs + if (master_dofs.size() == 0) break; + + // make sure the element constraints 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()]); + + // 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: + { + // we don't do anything here since we will come + // back to this face from the other cell, at which + // time we will fall into the first case clause + // above + break; + } + + case FiniteElementDomination::either_element_can_dominate: + { + // it appears as if neither element has any + // constraints on its neighbor. this may be because + // neither element has any DoFs on faces at all. or + // that the two elements are actually the same, + // although they happen to run under different + // fe_indices (this is what happens in + // hp/hp_hanging_nodes_01 for example). + // + // another possibility is what happens in crash_13. + // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs. + // FESystem(FE_Q(1),FE_DGQ(1)). neither of them + // dominates the other. + // + // a final possibility is that we have something like + // FESystem(FE_Q(1),FE_Q(1)) vs + // FESystem(FE_Q(1),FE_Nothing()), see + // hp/fe_nothing_18/19. + // + // in any case, the point is that it doesn't + // matter. there is nothing to do here. + break; + } + + case FiniteElementDomination::neither_element_dominates: + { + // we don't presently know what exactly to do here. + // it isn't quite clear what exactly we would have + // to do here. sit tight until someone trips over + // the following statement and see what exactly is + // going on + Assert (false, ExcNotImplemented()); + break; + } + + case FiniteElementDomination::no_requirements: + { + // nothing to do here + break; + } + + default: + // we shouldn't get here + Assert (false, ExcInternalError()); + } + } + } + } + } + + + + + template + void + make_hanging_node_constraints (const DH &dof_handler, + ConstraintMatrix &constraints) + { + // Decide whether to use the new or old make_hanging_node_constraints + // function. If all the FiniteElement or all elements in a FECollection + // support the new face constraint matrix, the new code will be used. + // Otherwise, the old implementation is used for the moment. + if (dof_handler.get_fe().hp_constraints_are_implemented ()) + internal:: + make_hp_hanging_node_constraints (dof_handler, + constraints); + else + internal:: + make_oldstyle_hanging_node_constraints (dof_handler, + constraints, + dealii::internal::int2type()); + } + + + + namespace + { + // enter constraints for periodicity into the given ConstraintMatrix object. + // this function is called when at least one of the two face iterators corresponds + // to an active object without further children + // + // @param transformation A matrix that maps degrees of freedom from one face + // to another. If the DoFs on the two faces are supposed to match exactly, then + // the matrix so provided will be the identity matrix. if face 2 is once refined + // from face 1, then the matrix needs to be the interpolation matrix from a face + // to this particular child + // + // @precondition: face_1 is supposed to be active + // + // @note As bug #82 ((http://code.google.com/p/dealii/issues/detail?id=82) and the + // corresponding testcase bits/periodicity_05 demonstrate, we can occasionally + // get into trouble if we already have the constraint x1=x2 and want to insert + // x2=x1. we avoid this by skipping an identity constraint if the opposite + // constraint already exists + template + void + set_periodicity_constraints (const FaceIterator &face_1, + const typename identity::type &face_2, + const FullMatrix &transformation, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask, + const bool face_orientation, + const bool face_flip, + const bool face_rotation) + { + static const int dim = FaceIterator::AccessorType::dimension; + static const int spacedim = FaceIterator::AccessorType::space_dimension; + + // we should be in the case where face_1 is active, i.e. has no children: + Assert (!face_1->has_children(), + ExcInternalError()); + + Assert (face_1->n_active_fe_indices() == 1, + ExcInternalError()); + + // if face_2 does have children, then we need to iterate over them + if (face_2->has_children()) + { + Assert (face_2->n_children() == GeometryInfo::max_children_per_face, + ExcNotImplemented()); + const unsigned int dofs_per_face + = face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face; + FullMatrix child_transformation (dofs_per_face, dofs_per_face); + FullMatrix subface_interpolation (dofs_per_face, dofs_per_face); + for (unsigned int c=0; cn_children(); ++c) + { + // get the interpolation matrix recursively from the one that + // interpolated from face_1 to face_2 by multiplying from the + // left with the one that interpolates from face_2 to + // its child + face_1->get_fe(face_1->nth_active_fe_index(0)) + .get_subface_interpolation_matrix (face_1->get_fe(face_1->nth_active_fe_index(0)), + c, + subface_interpolation); + subface_interpolation.mmult (child_transformation, transformation); + set_periodicity_constraints(face_1, face_2->child(c), + child_transformation, + constraint_matrix, component_mask, + face_orientation, face_flip, face_rotation); + } + } + else + // both faces are active. we need to 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")); + + const FiniteElement &fe = face_1->get_fe(face_1_index); + + Assert(component_mask.represents_n_components(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 dofs_per_face = fe.dofs_per_face; + + std::vector dofs_1(dofs_per_face); + std::vector dofs_2(dofs_per_face); + + face_1->get_dof_indices(dofs_1, face_1_index); + face_2->get_dof_indices(dofs_2, face_2_index); + + // Well, this is a hack: + // + // There is no + // face_to_face_index(face_index, + // face_orientation, + // face_flip, + // face_rotation) + // function in FiniteElementData, so we have to use + // face_to_cell_index(face_index, face + // face_orientation, + // face_flip, + // face_rotation) + // But this will give us an index on a cell - something we cannot work + // with directly. But luckily we can match them back :-] + + std::map cell_to_rotated_face_index; + + // Build up a cell to face index for face_2: + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int cell_index = fe.face_to_cell_index(i, 0, /* It doesn't really matter, just assume + * we're on the first face... + */ + true, false, false // default orientation + ); + cell_to_rotated_face_index[cell_index] = i; + } + + // loop over all dofs on face 2 and constrain them again the ones on face 1 + for (unsigned int i=0; i + void + make_periodicity_constraints (const FaceIterator &face_1, + const typename identity::type &face_2, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask, + const bool face_orientation, + const bool face_flip, + const bool face_rotation) + { + static const int dim = FaceIterator::AccessorType::dimension; + + Assert( (dim != 1) || + (face_orientation == true && + face_flip == false && + face_rotation == false), + ExcMessage ("The supplied orientation " + "(face_orientation, face_flip, face_rotation) " + "is invalid for 1D")); + + Assert( (dim != 2) || + (face_orientation == true && + face_rotation == false), + ExcMessage ("The supplied orientation " + "(face_orientation, face_flip, face_rotation) " + "is invalid for 2D")); + + Assert(face_1 != face_2, + ExcMessage ("face_1 and face_2 are equal! Cannot constrain DoFs " + "on the very same face")); + + Assert(face_1->at_boundary() && face_2->at_boundary(), + ExcMessage ("Faces for periodicity constraints must be on the boundary")); + + + // A lookup table on how to go through the child faces depending on the + // orientation: + + static const int lookup_table_2d[2][2] = + { + // flip: + {0, 1}, // false + {1, 0}, // true + }; + + static const int lookup_table_3d[2][2][2][4] = + { + // orientation flip rotation + { { {0, 2, 1, 3}, // false false false + {2, 3, 0, 1}, // false false true + }, + { {3, 1, 2, 0}, // false true false + {1, 0, 3, 2}, // false true true + }, + }, + { { {0, 1, 2, 3}, // true false false + {1, 3, 0, 2}, // true false true + }, + { {3, 2, 1, 0}, // true true false + {2, 0, 3, 1}, // true true true + }, + }, + }; + + // In the case that both faces have children, we loop over all + // children and apply 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, + ExcNotImplemented()); + + for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) + { + // Lookup the index for the second face + unsigned int j; + switch (dim) + { + case 2: + j = lookup_table_2d[face_flip][i]; + break; + case 3: + j = lookup_table_3d[face_orientation][face_flip][face_rotation][i]; + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + + make_periodicity_constraints (face_1->child(i), + face_2->child(j), + constraint_matrix, + component_mask, + face_orientation, + face_flip, + face_rotation); + } + } + else + // otherwise at least one of the two faces is active and + // we need to enter the constraints + { + if (face_2->has_children() == false) + set_periodicity_constraints(face_2, face_1, + FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), + constraint_matrix, + component_mask, + face_orientation, face_flip, face_rotation); + else + set_periodicity_constraints(face_1, face_2, + FullMatrix(IdentityMatrix(face_1->get_fe(face_1->nth_active_fe_index(0)).dofs_per_face)), + constraint_matrix, + component_mask, + face_orientation, face_flip, face_rotation); + } + } + + + + template + void + make_periodicity_constraints (const DH &dof_handler, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + const int direction, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask) + { + Tensor<1,DH::space_dimension> dummy; + make_periodicity_constraints (dof_handler, + b_id1, + b_id2, + direction, + dummy, + constraint_matrix, + component_mask); + } + + + + template + void + make_periodicity_constraints (const DH &dof_handler, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + const int direction, + dealii::Tensor<1,DH::space_dimension> &offset, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask) + { + static const int space_dim = DH::space_dimension; + Assert (0<=direction && direction PTRIA; + const PTRIA *ptria_p = dynamic_cast (&dof_handler.get_tria()); + Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1), + ExcMessage ("This function can not be used with distributed triangulations." + "See the documentation for more information.")); + } +#endif + + Assert (b_id1 != b_id2, + ExcMessage ("The boundary indicators b_id1 and b_id2 must be" + "different to denote different boundaries.")); + + typedef typename DH::face_iterator FaceIterator; + typedef std::map > > FaceMap; + + // Collect matching periodic cells on the coarsest level: + FaceMap matched_cells = + GridTools::collect_periodic_face_pairs(dof_handler, + b_id1, b_id2, + direction, offset); + + // And apply the low level make_periodicity_constraints function to + // every matching pair: + for (typename FaceMap::iterator it = matched_cells.begin(); + it != matched_cells.end(); ++it) + { + typedef typename DH::face_iterator FaceIterator; + const FaceIterator &face_1 = it->first; + const FaceIterator &face_2 = it->second.first; + const std::bitset<3> &orientation = it->second.second; + + Assert(face_1->at_boundary() && face_2->at_boundary(), + ExcInternalError()); + + Assert (face_1->boundary_indicator() == b_id1 && + face_2->boundary_indicator() == b_id2, + ExcInternalError()); + + Assert (face_1 != face_2, + ExcInternalError()); + + make_periodicity_constraints(face_1, + face_2, + constraint_matrix, + component_mask, + orientation[0], + orientation[1], + orientation[2]); + } + } + + + + template + void + make_periodicity_constraints (const DH &dof_handler, + const types::boundary_id b_id, + const int direction, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask) + { + Tensor<1,DH::space_dimension> dummy; + make_periodicity_constraints (dof_handler, + b_id, + direction, + dummy, + constraint_matrix, + component_mask); + } + + + + template + void + make_periodicity_constraints (const DH &dof_handler, + const types::boundary_id b_id, + const int direction, + dealii::Tensor<1,DH::space_dimension> &offset, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask) + { + static const int dim = DH::dimension; + static const int space_dim = DH::space_dimension; + + Assert (0<=direction && direction PTRIA; + const PTRIA *ptria_p = dynamic_cast (&dof_handler.get_tria()); + Assert ((ptria_p == 0 || Utilities::MPI::n_mpi_processes(ptria_p->get_communicator()) == 1), + ExcMessage ("This function can not be used with distributed triangulations." + "See the documentation for more information.")); + } +#endif + + typedef typename DH::face_iterator FaceIterator; + typedef std::map FaceMap; + + // Collect matching periodic cells on the coarsest level: + FaceMap matched_cells = + GridTools::collect_periodic_face_pairs(dof_handler, + b_id, + direction, offset); + + // And apply the low level make_periodicity_constraints function to + // every matching pair: + for (typename FaceMap::iterator it = matched_cells.begin(); + it != matched_cells.end(); ++it) + { + typedef typename DH::face_iterator FaceIterator; + const FaceIterator &face_1 = it->first; + const FaceIterator &face_2 = it->second; + + Assert(face_1->at_boundary() && face_2->at_boundary(), + ExcInternalError()); + + Assert (face_1->boundary_indicator() == b_id && + face_2->boundary_indicator() == b_id, + ExcInternalError()); + + Assert (face_1 != face_2, + ExcInternalError()); + + make_periodicity_constraints(face_1, + face_2, + constraint_matrix, + component_mask + /* standard orientation */); + } + } + + + + namespace internal + { + namespace + { + /** + * This is a function that is called by the _2 function and that + * operates on a range of cells only. It is used to split up the + * whole range of cells into chunks which are then worked on in + * parallel, if multithreading is available. + */ + 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) + { + // aliases to the finite elements used by the dof handlers: + const FiniteElement &coarse_fe = coarse_grid.get_fe(); + + // for each cell on the parameter grid: find out which degrees of + // freedom on the fine grid correspond in which way to the degrees + // of freedom on the parameter grid + // + // since for continuous FEs some dofs exist on more than one cell, + // we have to track which ones were already visited. the problem is + // that if we visit a dof first on one cell and compute its weight + // with respect to some global dofs to be non-zero, and later visit + // the dof again on another cell and (since we are on another cell) + // recompute the weights with respect to the same dofs as above to + // be zero now, we have to preserve them. we therefore overwrite + // all weights if they are nonzero and do not enforce zero weights + // since that might be only due to the fact that we are on another + // cell. + // + // example: + // coarse grid + // | | | + // *-----*-----* + // | cell|cell | + // | 1 | 2 | + // | | | + // 0-----1-----* + // + // fine grid + // | | | | | + // *--*--*--*--* + // | | | | | + // *--*--*--*--* + // | | | | | + // *--x--y--*--* + // + // when on cell 1, we compute the weights of dof 'x' to be 1/2 from + // parameter dofs 0 and 1, respectively. however, when later we are + // on cell 2, we again compute the prolongation of shape function 1 + // restricted to cell 2 to the globla grid and find that the weight + // of global dof 'x' now is zero. however, we should not overwrite + // the old value. + // + // we therefore always only set nonzero values. why adding up is + // not useful: dof 'y' would get weight 1 from parameter dof 1 on + // both cells 1 and 2, but the correct weight is nevertheless only + // 1. + + // vector to hold the representation of a single degree of freedom + // on the coarse grid (for the selected fe) on the fine grid + const types::global_dof_index 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); + + 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); + + // 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); + // now that we've got the global representation of each + // parameter dof, we've only got to clobber the non-zero + // entries in that vector and store the result + // + // what we have learned: if entry @p{i} of the global + // vector holds the value @p{v[i]}, then this is the + // weight with which the present dof contributes to + // @p{i}. there may be several such @p{i}s and their + // weights' sum should be one. Then, @p{v[i]} should be + // equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the + // values of the degrees of freedom on the coarse grid. + // we can thus compute constraints which link the degrees + // of freedom @p{v[i]} on the fine grid to those on the + // coarse grid, @p{p[j]}. Now to use these as real + // constraints, rather than as additional equations, we + // have to identify representants among the @p{i} for + // each @p{j}. this will be done by simply taking the + // first @p{i} for which @p{w_{ij}==1}. + // + // guard modification of the weights array by a Mutex. + // since it should happen rather rarely that there are + // several threads operating on different intergrid + // weights, have only one mutex for all of them + static Threads::Mutex mutex; + Threads::Mutex::ScopedLock lock (mutex); + for (types::global_dof_index 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) + { + // 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); + + // 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 + 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) + { + // aliases to the finite elements used by the dof handlers: + const FiniteElement &coarse_fe = coarse_grid.get_fe(), + &fine_fe = fine_grid.get_fe(); + + // global numbers of dofs + const types::global_dof_index 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; + + // 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; + + + // 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()); + + // 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 ()); + + + // check whether component numbers are valid + AssertIndexRange (coarse_component,coarse_fe.n_components()); + 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()); + +#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()); +#endif + + /* + * From here on: the term `parameter' refers to the selected + * component on the coarse grid and its analogon on the fine grid. + * The naming of variables containing this term is due to the fact + * that `selected_component' is longer, but also due to the fact + * that the code of this function was initially written for a + * program where the component which we wanted to match between + * grids was actually the `parameter' variable. + * + * Likewise, the terms `parameter grid' and `state grid' refer to + * the coarse and fine grids, respectively. + * + * Changing the names of variables would in principle be a good + * idea, but would not make things simpler and would be another + * source of errors. If anyone feels like doing so: patches would + * be welcome! + */ + + + + // set up vectors of cell-local data; each vector represents one + // 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)); + // 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); + + 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()); + }; +#endif + + + return n_parameters_on_fine_grid; + } + + + } + } + + + + 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) + { + // store the weights with which a dof on the parameter grid contributes + // to a dof on the fine grid. see the long doc below for more info + // + // allocate as many rows as there are parameter dofs on the coarse grid + // and as many columns as there are parameter dofs on the fine grid. + // + // weight_mapping is used to map the global (fine grid) parameter dof + // indices to the columns + // + // in the original implementation, the weights array was actually of + // FullMatrix type. this wasted huge amounts of memory, but was + // fast. nonetheless, since the memory consumption was quadratic in the + // number of degrees of freedom, this was not very practical, so we now + // use a vector of rows of the matrix, and in each row a vector of + // pairs (colnum,value). this seems like the best tradeoff between + // memory and speed, as it is now linear in memory and still fast + // enough. + // + // to save some memory and since the weights are usually (negative) + // powers of 2, we choose the value type of the matrix to be @p{float} + // rather than @p{double}. + std::vector > weights; + + // this is this mapping. there is one entry for each dof on the fine + // grid; if it is a parameter dof, then its value is the column in + // weights for that parameter dof, if it is any other dof, then its + // value is -1, indicating an error + std::vector weight_mapping; + + 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); + + // global numbers of dofs + const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(), + n_fine_dofs = fine_grid.n_dofs(); + + + // get an array in which we store which dof on the coarse grid is a + // parameter and which is not + 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, ComponentMask(mask), coarse_dof_is_parameter); + } + + // now we know that the weights in each row constitute a constraint. + // enter this into the constraints object + // + // first task: for each parameter dof on the parameter grid, find a + // representant on the fine, global grid. this is possible since we use + // conforming finite element. we take this representant to be the first + // element in this row with weight identical to one. the representant + // will become an unconstrained degree of freedom, while all others + // will be constrained to this dof (and possibly others) + std::vector representants(n_coarse_dofs, numbers::invalid_dof_index); + for (types::global_dof_index parameter_dof=0; parameter_dof 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 types::global_dof_index 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 + types::global_dof_index 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; + } + 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()); + }; + + + + // note for people that want to optimize this function: the largest + // part of the computing time is spent in the following, rather + // innocent block of code. basically, it must be the + // ConstraintMatrix::add_entry call which takes the bulk of the time, + // but it is not known to the author how to make it faster... + std::vector > constraint_line; + for (types::global_dof_index global_dof=0; global_dof 0, ExcInternalError()); + std::map::const_iterator + col_entry = weights[0].end(); + for (; first_used_rowsecond == 1) && + (representants[first_used_row] == global_dof)) + // dof unconstrained or constrained to itself (in case this + // cell is mapped to itself, rather than to children of + // itself) + continue; + } + + + // otherwise enter all constraints + constraints.add_line (global_dof); + + constraint_line.clear (); + for (types::global_dof_index 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); + }; + } + + + + 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) + { + // store the weights with which a dof on the parameter grid contributes + // to a dof on the fine grid. see the long doc below for more info + // + // allocate as many rows as there are parameter dofs on the coarse grid + // and as many columns as there are parameter dofs on the fine grid. + // + // weight_mapping is used to map the global (fine grid) parameter dof + // indices to the columns + // + // in the original implementation, the weights array was actually of + // FullMatrix type. this wasted huge amounts of memory, but was + // fast. nonetheless, since the memory consumption was quadratic in the + // number of degrees of freedom, this was not very practical, so we now + // use a vector of rows of the matrix, and in each row a vector of + // pairs (colnum,value). this seems like the best tradeoff between + // memory and speed, as it is now linear in memory and still fast + // enough. + // + // to save some memory and since the weights are usually (negative) + // powers of 2, we choose the value type of the matrix to be @p{float} + // rather than @p{double}. + std::vector > weights; + + // this is this mapping. there is one entry for each dof on the fine + // grid; if it is a parameter dof, then its value is the column in + // weights for that parameter dof, if it is any other dof, then its + // value is -1, indicating an error + 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); + + // now compute the requested representation + const types::global_dof_index n_global_parm_dofs + = std::count_if (weight_mapping.begin(), weight_mapping.end(), + std::bind2nd (std::not_equal_to (), numbers::invalid_dof_index)); + + // first construct the inverse mapping of weight_mapping + std::vector inverse_weight_mapping (n_global_parm_dofs, + DoFHandler::invalid_dof_index); + for (types::global_dof_index i=0; i::invalid_dof_index), + ExcInternalError()); + + inverse_weight_mapping[parameter_dof] = i; + }; + }; + + // next copy over weights array and replace respective numbers + const types::global_dof_index n_rows = weight_mapping.size(); + + transfer_representation.clear (); + transfer_representation.resize (n_rows); + + const types::global_dof_index n_coarse_dofs = coarse_grid.n_dofs(); + for (types::global_dof_index i=0; i::const_iterator j = weights[i].begin(); + for (; j!=weights[i].end(); ++j) + { + const types::global_dof_index p = inverse_weight_mapping[j->first]; + Assert (psecond; + }; + }; + } + + + + template class DH> + void + make_zero_boundary_constraints (const DH &dof, + const types::boundary_id boundary_indicator, + ConstraintMatrix &zero_boundary_constraints, + const ComponentMask &component_mask) + { + Assert (component_mask.represents_n_components(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); + + Assert (component_mask.n_selected_components(n_components) > 0, + ComponentMask::ExcNoComponentSelected()); + + // a field to store the indices + std::vector face_dofs; + face_dofs.reserve (max_dofs_per_face(dof)); + + typename DH::active_cell_iterator + cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + if (!cell->is_artificial()) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + const FiniteElement &fe = cell->get_fe(); + + const typename DH::face_iterator face = cell->face(face_no); + + // if face is on the boundary and satisfies the correct + // boundary id property + if (face->at_boundary () + && + ((boundary_indicator == numbers::invalid_boundary_id) + || + (face->boundary_indicator() == boundary_indicator))) + { + // 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()); + + // enter those dofs into the list that match the component + // signature. + for (unsigned int i=0; iget_fe().get_nonzero_components (i); + bool nonzero = false; + for (unsigned int c=0; c class DH> + void + make_zero_boundary_constraints (const DH &dof, + ConstraintMatrix &zero_boundary_constraints, + const ComponentMask &component_mask) + { + make_zero_boundary_constraints(dof, numbers::invalid_boundary_id, + zero_boundary_constraints, component_mask); + } + + +} // end of namespace DoFTools + + + +// explicit instantiations + +#include "dof_tools_constraints.inst" + + + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/dofs/dof_tools_constraints.inst.in b/deal.II/source/dofs/dof_tools_constraints.inst.in new file mode 100644 index 0000000000..28b9fc86ce --- /dev/null +++ b/deal.II/source/dofs/dof_tools_constraints.inst.in @@ -0,0 +1,134 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2009 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + + +for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS) +{ + template + void + DoFTools::make_hanging_node_constraints (const DH &dof_handler, + ConstraintMatrix &constraints); + +#if deal_II_dimension != 1 + template + void + DoFTools::make_periodicity_constraints (const DH::face_iterator &, + const DH::face_iterator &, + dealii::ConstraintMatrix &, + const ComponentMask &, + bool, bool, bool); + + template + void + DoFTools::make_periodicity_constraints(const DH &, + const types::boundary_id, + const types::boundary_id, + const int, + dealii::ConstraintMatrix &, + const ComponentMask &); + + template + void + DoFTools::make_periodicity_constraints(const DH &, + const types::boundary_id, + const types::boundary_id, + const int, + dealii::Tensor<1,DH::space_dimension> &, + dealii::ConstraintMatrix &, + const ComponentMask &); + + template + void + DoFTools::make_periodicity_constraints(const DH &, + const types::boundary_id, + const int, + dealii::ConstraintMatrix &, + const ComponentMask &); + + template + void + DoFTools::make_periodicity_constraints(const DH &, + const types::boundary_id, + const int, + dealii::Tensor<1,DH::space_dimension> &, + dealii::ConstraintMatrix &, + const ComponentMask &); +#endif +} + +for (deal_II_dimension : DIMENSIONS) +{ +#if deal_II_dimension != 1 + template + void + DoFTools::make_hanging_node_constraints (const MGDoFHandler &dof_handler, + ConstraintMatrix &constraints); + +#endif + } + + +for (deal_II_dimension : DIMENSIONS) +{ +#if deal_II_dimension < 3 +template +void +DoFTools:: +make_hanging_node_constraints (const DoFHandler &dof_handler, + ConstraintMatrix &constraints); +#endif + +#if deal_II_dimension == 3 +template +void +DoFTools:: +make_hanging_node_constraints (const DoFHandler<1,3> &dof_handler, + ConstraintMatrix &constraints); +#endif + +template +void +DoFTools::make_zero_boundary_constraints +(const DoFHandler &, + ConstraintMatrix &, + const ComponentMask &); + +template +void +DoFTools::make_zero_boundary_constraints +(const DoFHandler &, + const types::boundary_id , + ConstraintMatrix &, + const ComponentMask &); + +template +void +DoFTools::make_zero_boundary_constraints +(const hp::DoFHandler &, + ConstraintMatrix &, + const ComponentMask &); + +template +void +DoFTools::make_zero_boundary_constraints +(const hp::DoFHandler &, + const types::boundary_id , + ConstraintMatrix &, + const ComponentMask &); + +} diff --git a/deal.II/source/dofs/dof_tools_sparsity.cc b/deal.II/source/dofs/dof_tools_sparsity.cc new file mode 100644 index 0000000000..c48bf2f677 --- /dev/null +++ b/deal.II/source/dofs/dof_tools_sparsity.cc @@ -0,0 +1,1179 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1999 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + + +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) + { + const types::global_dof_index n_dofs = dof.n_dofs(); + + Assert (sparsity.n_rows() == n_dofs, + ExcDimensionMismatch (sparsity.n_rows(), n_dofs)); + Assert (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 subdomain is also okay, because we skip + // ghost cells in the loop below. + Assert ( + (dof.get_tria().locally_owned_subdomain() == numbers::invalid_subdomain_id) + || + (subdomain_id == numbers::invalid_subdomain_id) + || + (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.")); + + 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(); + + // In case we work with a distributed sparsity pattern of Trilinos + // type, we only have to do the work if the current cell is owned by + // the calling processor. Otherwise, just continue. + for (; cell!=endc; ++cell) + if (((subdomain_id == numbers::invalid_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); + + // 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); + } + } + + + + 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 types::global_dof_index n_dofs = dof.n_dofs(); + + 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 (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())); + + // If we have a distributed::Triangulation only allow locally_owned + // subdomain. Not setting a subdomain is also okay, because we skip + // ghost cells in the loop below. + Assert ( + (dof.get_tria().locally_owned_subdomain() == numbers::invalid_subdomain_id) + || + (subdomain_id == numbers::invalid_subdomain_id) + || + (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.")); + + const hp::FECollection fe_collection (dof.get_fe()); + + // first, for each finite element, build a mask for each dof, not like + // the one given which represents components. make sure we do the right + // thing also with respect to non-primitive shape functions, which + // takes some additional thought + std::vector > dof_mask(fe_collection.size()); + + // check whether the table of couplings contains only true arguments, + // i.e., we do not exclude any index. that is the easy case, since we + // don't have to set up the tables + 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(); + + // In case we work with a distributed sparsity pattern of Trilinos + // type, we only have to do the work if the current cell is owned by + // the calling processor. Otherwise, just continue. + for (; cell!=endc; ++cell) + if (((subdomain_id == numbers::invalid_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; + + 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]); + } + } + + + + template + void + make_sparsity_pattern ( + const DH &dof_row, + const DH &dof_col, + SparsityPattern &sparsity) + { + const types::global_dof_index n_dofs_row = dof_row.n_dofs(); + const types::global_dof_index n_dofs_col = dof_col.n_dofs(); + + Assert (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)); + +//TODO: Looks like wasteful memory management here + + const std::list > + cell_list + = GridTools::get_finest_common_cells (dof_row, dof_col); + + + typename std::list >::const_iterator + cell_iter = cell_list.begin(); + + 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; i + void + make_boundary_sparsity_pattern ( + const DH &dof, + const std::vector &dof_to_boundary_mapping, + SparsityPattern &sparsity) + { + if (DH::dimension == 1) + { + // there are only 2 boundary indicators in 1d, so it is no + // performance problem to call the other function + 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 types::global_dof_index n_dofs = dof.n_dofs(); + + AssertDimension (dof_to_boundary_mapping.size(), n_dofs); + AssertDimension (sparsity.n_rows(), dof.n_boundary_dofs()); + AssertDimension (sparsity.n_cols(), dof.n_boundary_dofs()); +#ifdef DEBUG + if (sparsity.n_rows() != 0) + { + types::global_dof_index max_element = 0; + for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != DH::invalid_dof_index) && + (*i > max_element)) + max_element = *i; + AssertDimension (max_element, sparsity.n_rows()-1); + }; +#endif + + std::vector dofs_on_this_face; + dofs_on_this_face.reserve (max_dofs_per_face(dof)); + + // loop over all faces to check whether they are at a boundary. note + // that we need not take special care of single lines (using + // @p{cell->has_boundary_lines}), since we do not support boundaries of + // dimension dim-2, and so every boundary line is also part of a + // boundary face. + typename DH::active_cell_iterator cell = dof.begin_active(), + 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()); + + // make sparsity pattern for this cell + for (unsigned int i=0; i + void make_boundary_sparsity_pattern ( + const DH &dof, + const typename FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SparsityPattern &sparsity) + { + if (DH::dimension == 1) + { + // first check left, then right boundary point + for (unsigned int direction=0; direction<2; ++direction) + { + // if this boundary is not requested, then go on with next one + if (boundary_indicators.find(direction) == + boundary_indicators.end()) + continue; + + // find active cell at that boundary: first go to left/right, + // then to children + typename DH::level_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); + + // next get boundary mapped dof indices of boundary dofs + 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); + }; +#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(); + 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()); + + // 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, + const types::subdomain_id subdomain_id) + + // TODO: QA: reduce the indentation level of this method..., Maier 2012 + + { + const types::global_dof_index n_dofs = dof.n_dofs(); + + AssertDimension (sparsity.n_rows(), n_dofs); + AssertDimension (sparsity.n_cols(), n_dofs); + + // If we have a distributed::Triangulation only allow locally_owned + // subdomain. Not setting a subdomain is also okay, because we skip + // ghost cells in the loop below. + Assert ( + (dof.get_tria().locally_owned_subdomain() == numbers::invalid_subdomain_id) + || + (subdomain_id == numbers::invalid_subdomain_id) + || + (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.")); + + 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(); + + // TODO: in an old implementation, we used user flags before to tag + // faces that were already touched. this way, we could reduce the work + // a little bit. now, we instead add only data from one side. this + // should be OK, but we need to actually verify it. + + // In case we work with a distributed sparsity pattern of Trilinos + // type, we only have to do the work if the current cell is owned by + // the calling processor. Otherwise, just continue. + for (; cell!=endc; ++cell) + if (((subdomain_id == numbers::invalid_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); + + // 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::level_cell_iterator neighbor = cell->neighbor(face); + + // in 1d, we do not need to worry whether the neighbor + // might have children and then loop over those children. + // rather, we may as well go straight to to cell behind + // this particular cell's most terminal child + if (DH::dimension==1) + while (neighbor->has_children()) + neighbor = neighbor->child(face==0 ? 1 : 0); + + if (neighbor->has_children()) + { + for (unsigned int sub_nr = 0; + sub_nr != cell_face->number_of_children(); + ++sub_nr) + { + const typename DH::level_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 + + // TODO: in the distributed case, we miss out the + // constraints when the neighbor cell is coarser, but + // only the current cell is owned locally! + 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); + + 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); + + // 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)->active() + || + (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); + } + } + } + } + } + + + + template + void + make_flux_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity) + { + ConstraintMatrix constraints; + make_flux_sparsity_pattern (dof, sparsity, constraints); + } + + template + Table<2,Coupling> + dof_couplings_from_component_couplings (const FiniteElement &fe, + const Table<2,Coupling> &component_couplings) + { + Assert(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())); + + const unsigned int n_dofs = fe.dofs_per_cell; + + Table<2,Coupling> dof_couplings (n_dofs, n_dofs); + + for (unsigned int i=0; i + std::vector > + dof_couplings_from_component_couplings + (const hp::FECollection &fe, + const Table<2,Coupling> &component_couplings) + { + std::vector > 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) + { + 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); + + 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); + + 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); + // 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); + // 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::level_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 (); + } + } + } + } + } + + + // implementation of the same function in namespace DoFTools for + // non-hp DoFHandlers + 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) + { + // 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(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); + + 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]); + + // 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::level_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::level_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 (); + } + } + } + } + } + } + + } + + + + + template + void + 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 types::global_dof_index 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)); + + // 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::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); + } + + +} // end of namespace DoFTools + + +// --------------------------------------------------- explicit instantiations + +#include "dof_tools_sparsity.inst" + + + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/dofs/dof_tools_sparsity.inst.in b/deal.II/source/dofs/dof_tools_sparsity.inst.in new file mode 100644 index 0000000000..6781dbeab2 --- /dev/null +++ b/deal.II/source/dofs/dof_tools_sparsity.inst.in @@ -0,0 +1,323 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2009 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS) + { + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler &dof, + SP &sparsity, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler &dof, + SP &sparsity, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const MGDoFHandler &dof, + SP &sparsity, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler&, + const Table<2,Coupling>&, + SP &, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler&, + const Table<2,Coupling>&, + SP &, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const MGDoFHandler&, + const Table<2,Coupling>&, + SP &, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler &dof_row, + const DoFHandler &dof_col, + SP &sparsity); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler &dof_row, + const hp::DoFHandler &dof_col, + SP &sparsity); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const DoFHandler& dof, + const std::vector &, + SP &); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const hp::DoFHandler& dof, + const std::vector &, + SP &); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const hp::DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + +#if deal_II_dimension < 3 + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const hp::DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + #endif + + template void + DoFTools::make_flux_sparsity_pattern,SP> + (const DoFHandler &dof, + SP &sparsity); + + template void + DoFTools::make_flux_sparsity_pattern,SP> + (const hp::DoFHandler &dof, + SP &sparsity); + +#if deal_II_dimension > 1 + + template void + DoFTools::make_flux_sparsity_pattern,SP> + (const DoFHandler &dof, + SP &, + const Table<2,Coupling>&, + const Table<2,Coupling>&); + + template void + DoFTools::make_flux_sparsity_pattern,SP> + (const DoFHandler &dof, + SP &sparsity, + const ConstraintMatrix &constraints, + const bool, const unsigned int); + + template void + DoFTools::make_flux_sparsity_pattern,SP> + (const hp::DoFHandler &dof, + SP &sparsity, + 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 + + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler &dof, + SP &sparsity, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler &dof, + SP &sparsity, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler&, + const Table<2,Coupling>&, + SP &, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler&, + const Table<2,Coupling>&, + SP &, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler &dof_row, + const DoFHandler &dof_col, + SP &sparsity); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler &dof_row, + const hp::DoFHandler &dof_col, + SP &sparsity); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const DoFHandler& dof, + const std::vector &, + SP &); + + //template void + //DoFTools::make_boundary_sparsity_pattern,SP> + //(const hp::DoFHandler& dof, + // const std::vector &, + // SP &); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const DoFHandler& dof, + const FunctionMap::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + + //template void + //DoFTools::make_boundary_sparsity_pattern,SP> + //(const hp::DoFHandler& dof, + // const FunctionMap::type &boundary_indicators, + // const std::vector &dof_to_boundary_mapping, + // SP &sparsity); + +#endif + + +#if deal_II_dimension == 3 + + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler<1,3> &dof, + SP &sparsity, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler<1,3> &dof, + SP &sparsity, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler<1,3>&, + const Table<2,Coupling>&, + SP &, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler<1,3>&, + const Table<2,Coupling>&, + SP &, + const ConstraintMatrix &, + const bool, + const unsigned int); + + template void + DoFTools::make_sparsity_pattern, SP> + (const DoFHandler<1,3> &dof_row, + const DoFHandler<1,3> &dof_col, + SP &sparsity); + + template void + DoFTools::make_sparsity_pattern, SP> + (const hp::DoFHandler<1,3> &dof_row, + const hp::DoFHandler<1,3> &dof_col, + SP &sparsity); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const DoFHandler<1,3>& dof, + const std::vector &, + SP &); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const hp::DoFHandler<1,3>& dof, + const std::vector &, + SP &); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const DoFHandler<1,3>& dof, + const FunctionMap<3>::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + + template void + DoFTools::make_boundary_sparsity_pattern,SP> + (const hp::DoFHandler<1,3>& dof, + const FunctionMap<3>::type &boundary_indicators, + const std::vector &dof_to_boundary_mapping, + SP &sparsity); + +#endif + + } + + +for (deal_II_dimension : DIMENSIONS) +{ + template + Table<2,DoFTools::Coupling> + DoFTools::dof_couplings_from_component_couplings + (const FiniteElement &fe, + const Table<2,DoFTools::Coupling> &component_couplings); +}