make_sparsity_pattern (const DH &dof,
const std::vector<std::vector<bool> > &mask,
SparsityPattern &sparsity_pattern);
-
- /**
- * Construct a sparsity pattern that
- * allows coupling degrees of freedom on
- * two different but related meshes.
- *
- * The idea is that if the two given
- * DoFHandler objects correspond to two
- * different meshes (and potentially to
- * different finite elements used on
- * these cells), but that if the two
- * triangulations they are based on are
- * derived from the same coarse mesh
- * through hierarchical refinement, then
- * one may set up a problem where one
- * would like to test shape functions
- * from one mesh against the shape
- * functions from another mesh. In
- * particular, this means that shape
- * functions from a cell on the first
- * mesh are tested against those on the
- * second cell that are located on the
- * corresponding cell; this
- * correspondence is something that the
- * IntergridMap class can determine.
- *
- * This function then constructs a
- * sparsity pattern for which the degrees
- * of freedom that represent the rows
- * come from the first given DoFHandler,
- * whereas the ones that correspond to
- * columns come from the second
- * DoFHandler.
- */
- template <class DH, class SparsityPattern>
- static
- void
- make_sparsity_pattern (const DH &dof_row,
- const DH &dof_col,
- SparsityPattern &sparsity);
/**
* Create the sparsity pattern for
* velocities into one block, and
* the pressure into another).
*
- * The result is returned in
- * @p dofs_per_component.
+ * The result is returned in @p
+ * dofs_per_component. Note that
+ * the size of @p
+ * dofs_per_component needs to be
+ * enough to hold all the indices
+ * specified in @p
+ * target_component. If this is
+ * not the case, an assertion is
+ * thrown. The indices not
+ * targetted by target_components
+ * are left untouched.
*/
template <int dim>
static void
* counting is done by
* blocks. See @ref GlossBlock
* "blocks" in the glossary for
- * details.
+ * details. Again the vectors are
+ * assumed to have the correct
+ * size before calling this
+ * function. If this is not the
+ * case, an assertion is thrown.
*/
template <int dim>
static void
* This function can be used when
* different variables shall be
* discretized on different
- * grids, where one grid is
+ * grids, where one grid is
* coarser than the other. This
* idea might seem nonsensical at
* first, but has reasonable
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <grid/intergrid_map.h>
-#include <grid/grid_tools.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_constraints.h>
-template <class DH, class SparsityPattern>
-void
-DoFTools::make_sparsity_pattern (
- const DH &dof_row,
- const DH &dof_col,
- SparsityPattern &sparsity)
-{
- const unsigned int n_dofs_row = dof_row.n_dofs();
- const unsigned int 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));
-
-
- const std::list<std::pair<typename DH::cell_iterator,
- typename DH::cell_iterator> >
- cell_list
- = GridTools::get_finest_common_cells (dof_row, dof_col);
-
-
- typename std::list<std::pair<typename DH::cell_iterator,
- typename DH::cell_iterator> >
- ::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<unsigned int>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<unsigned int>
- 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; i<dofs_per_cell_row; ++i)
- for (unsigned int j=0; j<dofs_per_cell_col; ++j)
- sparsity.add (local_dof_indices_row[i],
- local_dof_indices_col[j]);
- }
- else if (cell_row->has_children())
- {
- const std::vector<typename DH::active_cell_iterator >
- child_cells = GridTools::get_active_child_cells<DH> (cell_row);
- for (unsigned int i=0; i<child_cells.size(); i++)
- {
- const typename DH::active_cell_iterator
- cell_row_child = child_cells[i];
- const unsigned int dofs_per_cell_row =
- cell_row_child->get_fe().dofs_per_cell;
- const unsigned int dofs_per_cell_col =
- cell_col->get_fe().dofs_per_cell;
- std::vector<unsigned int>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<unsigned int>
- 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<dofs_per_cell_row; ++i)
- for (unsigned int j=0; j<dofs_per_cell_col; ++j)
- sparsity.add (local_dof_indices_row[i],
- local_dof_indices_col[j]);
- }
- }
- else
- {
- std::vector<typename DH::active_cell_iterator>
- child_cells = GridTools::get_active_child_cells<DH> (cell_col);
- for (unsigned int i=0; i<child_cells.size(); i++)
- {
- const typename DH::active_cell_iterator
- cell_col_child = child_cells[i];
- const unsigned int dofs_per_cell_row =
- cell_row->get_fe().dofs_per_cell;
- const unsigned int dofs_per_cell_col =
- cell_col_child->get_fe().dofs_per_cell;
- std::vector<unsigned int>
- local_dof_indices_row(dofs_per_cell_row);
- std::vector<unsigned int>
- 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<dofs_per_cell_row; ++i)
- for (unsigned int j=0; j<dofs_per_cell_col; ++j)
- sparsity.add (local_dof_indices_row[i],
- local_dof_indices_col[j]);
- }
- }
- }
-}
-
-
-
#if deal_II_dimension == 1
template <class DH, class SparsityPattern>
ExcInternalError());
Assert (fe2.dofs_per_line <= fe1.dofs_per_line,
ExcInternalError());
- Assert ((dim < 3)
- ||
- (fe2.dofs_per_quad <= fe1.dofs_per_quad),
+ Assert (fe2.dofs_per_quad <= fe1.dofs_per_quad,
ExcInternalError());
// the idea here is to designate as
/**
- * Copy constraints into a constraint
- * matrix object.
- *
- * This function removes zero
+ * This method 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.
*/
#ifdef DEAL_II_ANON_NAMESPACE_BUG
static
}
if (constraint_already_satisfied == false)
- {
- // add up the absolute
- // values of all
- // constraints in this line
- // to get a measure of
- // their absolute size
- double abs_sum = 0;
- for (unsigned int i=0; i<n_master_dofs; ++i)
- abs_sum += std::abs (face_constraints(row,i));
-
- // then enter those
- // constraints that are
- // larger than
- // 1e-14*abs_sum. everything
- // else probably originated
- // from inexact inversion
- // of matrices and similar
- // effects. having those
- // constraints in here will
- // only lead to problems
- // because it makes
- // sparsity patterns fuller
- // than necessary without
- // producing any
- // significant effect
+ {
constraints.add_line (slave_dofs[row]);
for (unsigned int i=0; i<n_master_dofs; ++i)
- if ((face_constraints(row,i) != 0)
- &&
- (std::fabs(face_constraints(row,i)) >= 1e-14*abs_sum))
+ if (face_constraints(row,i) != 0)
{
#ifdef WOLFGANG
std::cout << " " << slave_dofs[row]
// (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. the point is
- // that it doesn't
- // matter, since
- // hopefully for this
- // case, both sides'
- // dofs should have
- // been unified.
- //
- // make sure this is
- // actually true. this
- // actually only
- // matters, of course,
- // if either of the two
- // finite elements
- // actually do have
- // dofs on the face
- if ((cell->get_fe().dofs_per_face != 0)
- ||
- (cell->neighbor(face)->get_fe().dofs_per_face != 0))
+ // for example); check
+ // the latter somewhat
+ // crudely by comparing
+ // fe names
+ if (cell->get_fe().get_name() !=
+ cell->neighbor(face)->get_fe().get_name())
{
- Assert (cell->get_fe().dofs_per_face
- ==
- cell->neighbor(face)->get_fe().dofs_per_face,
+ Assert (cell->get_fe().dofs_per_face == 0,
+ ExcNotImplemented());
+ Assert (cell->neighbor(face)->get_fe().dofs_per_face == 0,
ExcNotImplemented());
-
- // (ab)use the master
- // and slave dofs
- // arrays for a
- // moment here
- master_dofs.resize (cell->get_fe().dofs_per_face);
- cell->face(face)
- ->get_dof_indices (master_dofs,
- cell->active_fe_index ());
-
- slave_dofs.resize (cell->neighbor(face)->get_fe().dofs_per_face);
- cell->face(face)
- ->get_dof_indices (slave_dofs,
- cell->neighbor(face)->active_fe_index ());
-
- for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
- Assert (master_dofs[i] == slave_dofs[i],
- ExcInternalError());
}
break;
Assert(selected_dofs.size() == dof_handler.n_dofs(),
ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
// preset all values by false
- std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
+ fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
const FiniteElement<dim> &fe = dof_handler.get_fe();
{
const FiniteElement<dim>& fe = dof_handler.get_fe();
const unsigned int n_components = fe.n_components();
- dofs_per_component.resize (n_components);
- std::fill (dofs_per_component.begin(), dofs_per_component.end(), 0U);
-
// If the empty vector was given as
// default argument, set up this
// vector as identity.
Assert(target_component.size()==n_components,
ExcDimensionMismatch(target_component.size(),n_components));
+
+
+ // Check that target components
+ // contains sensible information
+ // and reset the counters for
+ // dofs_per_component
+ unsigned int size_dst = dofs_per_component.size();
+ for(unsigned int i=0; i<n_components; ++i) {
+ Assert(target_component[i] < size_dst,
+ ExcIndexRange(target_component[i],
+ 0,
+ size_dst));
+ dofs_per_component[target_component[i]] = 0;
+ }
+
// special case for only one
// component. treat this first
// computations
if (n_components == 1)
{
- dofs_per_component[0] = dof_handler.n_dofs();
+ dofs_per_component[target_component[0]] = dof_handler.n_dofs();
return;
}
-
+
// otherwise determine the number
// of dofs in each component
// separately. do so in parallel
dofs_in_component[i], false);
};
threads.join_all ();
-
+
+ unsigned int total_count = 0;
// next count what we got
unsigned int component = 0;
for (unsigned int b=0;b<fe.n_base_elements();++b)
{
for (unsigned int dd=0;dd<d;++dd)
{
- if (base.is_primitive() || (!only_once || dd==0))
- dofs_per_component[target_component[component]]
- += std::count(dofs_in_component[component].begin(),
- dofs_in_component[component].end(),
- true);
+ if (base.is_primitive() || (!only_once || dd==0)) {
+ unsigned int this_count =
+ std::count(dofs_in_component[component].begin(),
+ dofs_in_component[component].end(),
+ true);
+ dofs_per_component[target_component[component]]
+ += this_count;
+ total_count += this_count;
+ }
++component;
}
}
// only valid if the finite element
// is actually primitive, so
// exclude other elements from this
+
Assert (!dof_handler.get_fe().is_primitive()
||
- (std::accumulate (dofs_per_component.begin(),
- dofs_per_component.end(), 0U)
- == dof_handler.n_dofs()),
- ExcInternalError());
+ (total_count == dof_handler.n_dofs()),
+ ExcDimensionMismatch(total_count, dof_handler.n_dofs()));
}
{
const FiniteElement<dim>& fe = dof_handler.get_fe();
const unsigned int n_blocks = fe.n_blocks();
- dofs_per_block.resize (n_blocks);
- std::fill (dofs_per_block.begin(), dofs_per_block.end(), 0U);
+ unsigned int size_dst = dofs_per_block.size();
// If the empty vector was given as
// default argument, set up this
Assert(target_block.size()==n_blocks,
ExcDimensionMismatch(target_block.size(),n_blocks));
-
+
+ // Check constistency and reset
+ // counters
+ for(unsigned int i=0; i<n_blocks; ++i) {
+ Assert(target_block[i] < size_dst,
+ ExcIndexRange(target_block[i], 0, size_dst));
+ dofs_per_block[target_block[i]] = 0;
+ }
+
// special case for only one
// block. treat this first
// since it does not require any
// computations
if (n_blocks == 1)
{
- dofs_per_block[0] = dof_handler.n_dofs();
+ dofs_per_block[target_block[0]] = dof_handler.n_dofs();
return;
}
// otherwise determine the number
const Table<2,Coupling>&,
CompressedBlockSparsityPattern&);
-
-template void
-DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
- SparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof_row,
- const DoFHandler<deal_II_dimension> &dof_col,
- SparsityPattern &sparsity);
-template void
-DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
- CompressedSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof_row,
- const DoFHandler<deal_II_dimension> &dof_col,
- CompressedSparsityPattern &sparsity);
-template void
-DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
- BlockSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof_row,
- const DoFHandler<deal_II_dimension> &dof_col,
- BlockSparsityPattern &sparsity);
-template void
-DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
- CompressedBlockSparsityPattern>
-(const DoFHandler<deal_II_dimension> &dof_row,
- const DoFHandler<deal_II_dimension> &dof_col,
- CompressedBlockSparsityPattern &sparsity);
-
-template void
-DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
- SparsityPattern>
-(const hp::DoFHandler<deal_II_dimension> &dof_row,
- const hp::DoFHandler<deal_II_dimension> &dof_col,
- SparsityPattern &sparsity);
-template void
-DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
- CompressedSparsityPattern>
-(const hp::DoFHandler<deal_II_dimension> &dof_row,
- const hp::DoFHandler<deal_II_dimension> &dof_col,
- CompressedSparsityPattern &sparsity);
-template void
-DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
- BlockSparsityPattern>
-(const hp::DoFHandler<deal_II_dimension> &dof_row,
- const hp::DoFHandler<deal_II_dimension> &dof_col,
- BlockSparsityPattern &sparsity);
-template void
-DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
- CompressedBlockSparsityPattern>
-(const hp::DoFHandler<deal_II_dimension> &dof_row,
- const hp::DoFHandler<deal_II_dimension> &dof_col,
- CompressedBlockSparsityPattern &sparsity);
-
-
// #if deal_II_dimension > 1
template void
DoFTools::make_boundary_sparsity_pattern<DoFHandler<deal_II_dimension>,SparsityPattern>