const unsigned int component = 0);
/**
- * Extract the indices of the degrees
- * of freedom belonging to certain
- * components. The bit vector @p{select}
- * defines, which components of an
- * @ref{FESystem} are to be extracted
- * from the @ref{DoFHandler} @p{dof}. The
- * respective entries in @p{selected_dofs}
- * are then flagged @p{true}, while all
- * others are set to @p{false}.
+ * Extract the indices of the
+ * degrees of freedom belonging
+ * to certain vector components
+ * of a vector-valued finite
+ * element. The bit vector
+ * @p{select} defines, which
+ * components of an
+ * @ref{FESystem} are to be
+ * extracted from the
+ * @ref{DoFHandler} @p{dof}. The
+ * entries in @p{selected_dofs}
+ * corresponding to degrees of
+ * freedom belonging to these
+ * components are then flagged
+ * @p{true}, while all others are
+ * set to @p{false}.
*
- * The size of @p{component_select}
- * shall equal the number of
- * components in the finite
- * element used by @p{dof}. The
- * size of @p{selected_dofs} shall
- * equal
+ * The size of
+ * @p{component_select} shall
+ * equal the number of components
+ * in the finite element used by
+ * @p{dof}. The size of
+ * @p{selected_dofs} shall equal
* @p{dof_handler.n_dofs()}. Previous
* contents of this array or
* overwritten.
+ *
+ * If the finite element under
+ * consideration is not
+ * primitive, that is some or all
+ * of its shape functions are
+ * non-zero in more than one
+ * vector component (which holds,
+ * for example, for Nedelec or
+ * Raviart-Thomas elements), then
+ * shape functions cannot be
+ * associated with a single
+ * vector component. In this
+ * case, if @em{one} shape vector
+ * component of this element is
+ * flagged in
+ * @p{component_select}, then
+ * this is equivalent to
+ * selecting @em{all} vector
+ * components corresponding to
+ * this non-primitive base
+ * element.
*/
template <int dim>
static void
* number belong to each
* components. If the number of
* components the finite element
- * has (i.e. you only have one
- * scalar variable), then the
+ * has is one (i.e. you only have
+ * one scalar variable), then the
* number in this component
* obviously equals the total
* number of degrees of
* needs to equal the total
* number.
*
+ * However, the last statement
+ * does not hold true if the
+ * finite element is not
+ * primitive, i.e. some or all of
+ * its shape functions are
+ * non-zero in more than one
+ * vector component. This
+ * applies, for example, to the
+ * Nedelec or Raviart-Thomas
+ * elements. In this case, a
+ * degree of freedom is counted
+ * in each component in which it
+ * is non-zero, so that the sum
+ * mentioned above is greater
+ * than the total number of
+ * degrees of freedom.
+ *
* The result is returned in the
* last argument.
*/
Assert(selected_dofs.size() == dof.n_dofs(),
ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs()));
+ // two special cases: no component
+ // is selected, and all components
+ // are selected, both rather
+ // stupid, but easy to catch
+ if (std::count (component_select.begin(), component_select.end(), true)
+ == 0)
+ {
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(), false);
+ return;
+ };
+ if (std::count (component_select.begin(), component_select.end(), true)
+ == static_cast<signed int>(component_select.size()))
+ {
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(), true);
+ return;
+ };
+
+
// preset all values by false
- fill_n (selected_dofs.begin(), dof.n_dofs(), false);
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(), false);
+
+ // next set up a table for the
+ // degrees of freedom on each of
+ // the cells whether it is
+ // something interesting or not
+ std::vector<bool> local_selected_dofs (fe.dofs_per_cell, false);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ if (fe.is_primitive(i))
+ local_selected_dofs[i]
+ = component_select[fe.system_to_component_index(i).first];
+ else
+ // if this shape function is
+ // not primitive, then we have
+ // to work harder. we have to
+ // find out whether _any_ of
+ // the vector components of
+ // this element is selected or
+ // not
+ //
+ // to do so, get the first and
+ // last vector components of
+ // the base element to which
+ // the local dof with index i
+ // belongs
+ {
+ unsigned int first_comp = 0;
+ const unsigned int this_base = fe.system_to_base_index(i).first.first;
+ const unsigned int this_multiplicity
+ = fe.system_to_base_index(i).first.second;
+
+ for (unsigned int b=0; b<this_base; ++b)
+ first_comp += fe.base_element(b).n_components() *
+ fe.element_multiplicity(b);
+ for (unsigned int m=0; m<this_multiplicity; ++m)
+ first_comp += fe.base_element(this_base).n_components();
+ const unsigned int end_comp = first_comp +
+ fe.base_element(this_base).n_components();
+
+ Assert (first_comp < fe.n_components(), ExcInternalError());
+ Assert (end_comp <= fe.n_components(), ExcInternalError());
+
+ // now check whether any of
+ // the components in between
+ // is set
+ for (unsigned int c=first_comp; c<end_comp; ++c)
+ if (component_select[c] == true)
+ {
+ local_selected_dofs[i] = true;
+ break;
+ };
+ };
+ // then loop over all cells and do
+ // the work
std::vector<unsigned int> indices(fe.dofs_per_cell);
-
typename DoFHandler<dim>::active_cell_iterator c;
- for (c = dof.begin_active() ; c != dof.end() ; ++ c)
+ for (c=dof.begin_active(); c!=dof.end(); ++ c)
{
c->get_dof_indices(indices);
- for (unsigned int i=0;i<fe.dofs_per_cell;++i)
- {
- const unsigned int component = fe.system_to_component_index(i).first;
-
- if (component_select[component] == true)
- selected_dofs[indices[i]] = true;
- }
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ selected_dofs[indices[i]] = local_selected_dofs[i];
}
}
Assert(selected_dofs.size() == dof.n_dofs(level),
ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level)));
- // preset all values by false
- fill_n (selected_dofs.begin(), dof.n_dofs(level), false);
+ // two special cases: no component
+ // is selected, and all components
+ // are selected, both rather
+ // stupid, but easy to catch
+ if (std::count (component_select.begin(), component_select.end(), true)
+ == 0)
+ {
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false);
+ return;
+ };
+ if (std::count (component_select.begin(), component_select.end(), true)
+ == static_cast<signed int>(component_select.size()))
+ {
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(level), true);
+ return;
+ };
- std::vector<unsigned int> indices(fe.dofs_per_cell);
-
+ // preset all values by false
+ std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false);
+
+ // next set up a table for the
+ // degrees of freedom on each of
+ // the cells whether it is
+ // something interesting or not
+ std::vector<bool> local_selected_dofs (fe.dofs_per_cell, false);
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ if (fe.is_primitive(i))
+ local_selected_dofs[i]
+ = component_select[fe.system_to_component_index(i).first];
+ else
+ // if this shape function is
+ // not primitive, then we have
+ // to work harder. we have to
+ // find out whether _any_ of
+ // the vector components of
+ // this element is selected or
+ // not
+ //
+ // to do so, get the first and
+ // last vector components of
+ // the base element to which
+ // the local dof with index i
+ // belongs
+ {
+ unsigned int first_comp = 0;
+ const unsigned int this_base = fe.system_to_base_index(i).first.first;
+ const unsigned int this_multiplicity
+ = fe.system_to_base_index(i).first.second;
+
+ for (unsigned int b=0; b<this_base; ++b)
+ first_comp += fe.base_element(b).n_components() *
+ fe.element_multiplicity(b);
+ for (unsigned int m=0; m<this_multiplicity; ++m)
+ first_comp += fe.base_element(this_base).n_components();
+ const unsigned int end_comp = first_comp +
+ fe.base_element(this_base).n_components();
+
+ Assert (first_comp < fe.n_components(), ExcInternalError());
+ Assert (end_comp <= fe.n_components(), ExcInternalError());
+
+ // now check whether any of
+ // the components in between
+ // is set
+ for (unsigned int c=first_comp; c<end_comp; ++c)
+ if (component_select[c] == true)
+ {
+ local_selected_dofs[i] = true;
+ break;
+ };
+ };
+
+ // then loop over all cells and do
+ // work
+ std::vector<unsigned int> indices(fe.dofs_per_cell);
typename MGDoFHandler<dim>::cell_iterator c;
for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)
{
c->get_mg_dof_indices(indices);
- for (unsigned int i=0;i<fe.dofs_per_cell;++i)
- {
- const unsigned int component = fe.system_to_component_index(i).first;
- if (component_select[component] == true)
- selected_dofs[indices[i]] = true;
- }
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ selected_dofs[indices[i]] = local_selected_dofs[i];
}
}
#endif
template <int dim>
void
-DoFTools::count_dofs_per_component (const DoFHandler<dim> &dof_handler,
- std::vector<unsigned int> &dofs_per_component)
+DoFTools::
+count_dofs_per_component (const DoFHandler<dim> &dof_handler,
+ std::vector<unsigned int> &dofs_per_component)
{
const unsigned int n_components = dof_handler.get_fe().n_components();
dofs_per_component.resize (n_components);
// otherwise determine the number
// of dofs in each component
// separately. do so in parallel
- std::vector<std::vector<bool> > dofs_in_component (n_components,
- std::vector<bool>(dof_handler.n_dofs(),
- false));
- std::vector<std::vector<bool> > component_select (n_components,
- std::vector<bool>(n_components, false));
+ std::vector<std::vector<bool> >
+ dofs_in_component (n_components,
+ std::vector<bool>(dof_handler.n_dofs(), false));
+ std::vector<std::vector<bool> >
+ component_select (n_components,
+ std::vector<bool>(n_components, false));
Threads::ThreadManager thread_manager;
for (unsigned int i=0; i<n_components; ++i)
{
void (*fun_ptr) (const DoFHandler<dim> &,
const std::vector<bool> &,
- std::vector<bool> &) = &DoFTools::template extract_dofs<dim>;
+ std::vector<bool> &)
+ = &DoFTools::template extract_dofs<dim>;
component_select[i][i] = true;
Threads::spawn (thread_manager,
Threads::encapsulate (fun_ptr)
dofs_in_component[i].end(),
true);
- // finally sanity check
- Assert (std::accumulate (dofs_per_component.begin(), dofs_per_component.end(), 0U)
- ==
- dof_handler.n_dofs(),
+ // finally sanity check. this is
+ // 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());
};
--- /dev/null
+//---------------------------- anna_3.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002 by the deal.II authors and Anna Schneebeli
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- anna_3.cc ---------------------------
+
+
+// check some things about Nedelec elements, here that
+// DoFTools::component_select and DoFTools::count_dofs_per_component
+// works
+//
+// this program is a modified version of one by Anna Schneebeli,
+// University of Basel
+
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_system.h>
+#include <fe/fe_q.h>
+#include <fe/fe_nedelec.h>
+#include <fe/fe_base.h>
+#include <dofs/dof_tools.h>
+#include <numerics/dof_renumbering.h>
+#include <iostream>
+#include <fstream>
+
+
+template <int dim>
+class SystemTest
+{
+ public:
+ SystemTest ();
+ void run ();
+
+ private:
+ void make_grid_and_dofs ();
+ void check ();
+
+
+ Triangulation<dim> triangulation;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+
+};
+
+template <int dim>
+SystemTest<dim>::SystemTest () :
+ fe (FE_Nedelec<dim>(1), 2,
+ FE_Q<dim>(1), 1),
+ dof_handler (triangulation)
+{};
+
+
+template <int dim>
+void SystemTest<dim>::make_grid_and_dofs ()
+{
+
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (0);
+ deallog << "Number of active cells: " << triangulation.n_active_cells()
+ << std::endl;
+ deallog << "Total number of cells: " << triangulation.n_cells()
+ << std::endl;
+
+ dof_handler.distribute_dofs (fe);
+ deallog << "Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
+
+};
+
+
+template <int dim>
+void SystemTest<dim>::check ()
+{
+ for (unsigned int c=0; c<fe.n_components(); ++c)
+ {
+ deallog << "Checking for component " << c << std::endl;
+ std::vector<bool> x(fe.n_components(), false);
+ x[c] = true;
+ std::vector<bool> sel(dof_handler.n_dofs());
+ DoFTools::extract_dofs (dof_handler, x, sel);
+
+ for (unsigned int i=0; i<sel.size(); ++i)
+ if (sel[i])
+ deallog << " DoF " << i << std::endl;
+ };
+
+ std::vector<unsigned int> dofs_per_component (fe.n_components(), 0U);
+ DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
+ deallog << "DoFs per component: ";
+ for (unsigned int i=0; i<fe.n_components(); ++i)
+ deallog << dofs_per_component[i] << ' ';
+ deallog << std::endl;
+};
+
+
+template <int dim>
+void SystemTest<dim>::run ()
+{
+ deallog << "************* " << dim << "D *************" << std::endl;
+ make_grid_and_dofs ();
+ check();
+
+ // renumber degrees of freedom and
+ // try again
+ deallog << std::endl << "*** Renumbering ***" << std::endl;
+ DoFRenumbering::component_wise (dof_handler);
+ check ();
+};
+
+
+
+int main ()
+{
+ std::ofstream logfile("anna_3.output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ SystemTest<2>().run();
+ SystemTest<3>().run();
+ return 0;
+};
--- /dev/null
+
+DEAL::************* 2D *************
+DEAL::Number of active cells: 1
+DEAL::Total number of cells: 1
+DEAL::Number of degrees of freedom: 12
+DEAL::Checking for component 0
+DEAL:: DoF 4
+DEAL:: DoF 6
+DEAL:: DoF 8
+DEAL:: DoF 10
+DEAL::Checking for component 1
+DEAL:: DoF 4
+DEAL:: DoF 6
+DEAL:: DoF 8
+DEAL:: DoF 10
+DEAL::Checking for component 2
+DEAL:: DoF 5
+DEAL:: DoF 7
+DEAL:: DoF 9
+DEAL:: DoF 11
+DEAL::Checking for component 3
+DEAL:: DoF 5
+DEAL:: DoF 7
+DEAL:: DoF 9
+DEAL:: DoF 11
+DEAL::Checking for component 4
+DEAL:: DoF 0
+DEAL:: DoF 1
+DEAL:: DoF 2
+DEAL:: DoF 3
+DEAL::DoFs per component: 4 4 4 4 4
+DEAL::
+DEAL::*** Renumbering ***
+DEAL::Checking for component 0
+DEAL:: DoF 0
+DEAL:: DoF 1
+DEAL:: DoF 2
+DEAL:: DoF 3
+DEAL::Checking for component 1
+DEAL:: DoF 0
+DEAL:: DoF 1
+DEAL:: DoF 2
+DEAL:: DoF 3
+DEAL::Checking for component 2
+DEAL:: DoF 4
+DEAL:: DoF 5
+DEAL:: DoF 6
+DEAL:: DoF 7
+DEAL::Checking for component 3
+DEAL:: DoF 4
+DEAL:: DoF 5
+DEAL:: DoF 6
+DEAL:: DoF 7
+DEAL::Checking for component 4
+DEAL:: DoF 8
+DEAL:: DoF 9
+DEAL:: DoF 10
+DEAL:: DoF 11
+DEAL::DoFs per component: 4 4 4 4 4
+DEAL::************* 3D *************
+DEAL::Number of active cells: 1
+DEAL::Total number of cells: 1
+DEAL::Number of degrees of freedom: 32
+DEAL::Checking for component 0
+DEAL:: DoF 8
+DEAL:: DoF 10
+DEAL:: DoF 12
+DEAL:: DoF 14
+DEAL:: DoF 16
+DEAL:: DoF 18
+DEAL:: DoF 20
+DEAL:: DoF 22
+DEAL:: DoF 24
+DEAL:: DoF 26
+DEAL:: DoF 28
+DEAL:: DoF 30
+DEAL::Checking for component 1
+DEAL:: DoF 8
+DEAL:: DoF 10
+DEAL:: DoF 12
+DEAL:: DoF 14
+DEAL:: DoF 16
+DEAL:: DoF 18
+DEAL:: DoF 20
+DEAL:: DoF 22
+DEAL:: DoF 24
+DEAL:: DoF 26
+DEAL:: DoF 28
+DEAL:: DoF 30
+DEAL::Checking for component 2
+DEAL:: DoF 8
+DEAL:: DoF 10
+DEAL:: DoF 12
+DEAL:: DoF 14
+DEAL:: DoF 16
+DEAL:: DoF 18
+DEAL:: DoF 20
+DEAL:: DoF 22
+DEAL:: DoF 24
+DEAL:: DoF 26
+DEAL:: DoF 28
+DEAL:: DoF 30
+DEAL::Checking for component 3
+DEAL:: DoF 9
+DEAL:: DoF 11
+DEAL:: DoF 13
+DEAL:: DoF 15
+DEAL:: DoF 17
+DEAL:: DoF 19
+DEAL:: DoF 21
+DEAL:: DoF 23
+DEAL:: DoF 25
+DEAL:: DoF 27
+DEAL:: DoF 29
+DEAL:: DoF 31
+DEAL::Checking for component 4
+DEAL:: DoF 9
+DEAL:: DoF 11
+DEAL:: DoF 13
+DEAL:: DoF 15
+DEAL:: DoF 17
+DEAL:: DoF 19
+DEAL:: DoF 21
+DEAL:: DoF 23
+DEAL:: DoF 25
+DEAL:: DoF 27
+DEAL:: DoF 29
+DEAL:: DoF 31
+DEAL::Checking for component 5
+DEAL:: DoF 9
+DEAL:: DoF 11
+DEAL:: DoF 13
+DEAL:: DoF 15
+DEAL:: DoF 17
+DEAL:: DoF 19
+DEAL:: DoF 21
+DEAL:: DoF 23
+DEAL:: DoF 25
+DEAL:: DoF 27
+DEAL:: DoF 29
+DEAL:: DoF 31
+DEAL::Checking for component 6
+DEAL:: DoF 0
+DEAL:: DoF 1
+DEAL:: DoF 2
+DEAL:: DoF 3
+DEAL:: DoF 4
+DEAL:: DoF 5
+DEAL:: DoF 6
+DEAL:: DoF 7
+DEAL::DoFs per component: 12 12 12 12 12 12 8
+DEAL::
+DEAL::*** Renumbering ***
+DEAL::Checking for component 0
+DEAL:: DoF 0
+DEAL:: DoF 1
+DEAL:: DoF 2
+DEAL:: DoF 3
+DEAL:: DoF 4
+DEAL:: DoF 5
+DEAL:: DoF 6
+DEAL:: DoF 7
+DEAL:: DoF 8
+DEAL:: DoF 9
+DEAL:: DoF 10
+DEAL:: DoF 11
+DEAL::Checking for component 1
+DEAL:: DoF 0
+DEAL:: DoF 1
+DEAL:: DoF 2
+DEAL:: DoF 3
+DEAL:: DoF 4
+DEAL:: DoF 5
+DEAL:: DoF 6
+DEAL:: DoF 7
+DEAL:: DoF 8
+DEAL:: DoF 9
+DEAL:: DoF 10
+DEAL:: DoF 11
+DEAL::Checking for component 2
+DEAL:: DoF 0
+DEAL:: DoF 1
+DEAL:: DoF 2
+DEAL:: DoF 3
+DEAL:: DoF 4
+DEAL:: DoF 5
+DEAL:: DoF 6
+DEAL:: DoF 7
+DEAL:: DoF 8
+DEAL:: DoF 9
+DEAL:: DoF 10
+DEAL:: DoF 11
+DEAL::Checking for component 3
+DEAL:: DoF 12
+DEAL:: DoF 13
+DEAL:: DoF 14
+DEAL:: DoF 15
+DEAL:: DoF 16
+DEAL:: DoF 17
+DEAL:: DoF 18
+DEAL:: DoF 19
+DEAL:: DoF 20
+DEAL:: DoF 21
+DEAL:: DoF 22
+DEAL:: DoF 23
+DEAL::Checking for component 4
+DEAL:: DoF 12
+DEAL:: DoF 13
+DEAL:: DoF 14
+DEAL:: DoF 15
+DEAL:: DoF 16
+DEAL:: DoF 17
+DEAL:: DoF 18
+DEAL:: DoF 19
+DEAL:: DoF 20
+DEAL:: DoF 21
+DEAL:: DoF 22
+DEAL:: DoF 23
+DEAL::Checking for component 5
+DEAL:: DoF 12
+DEAL:: DoF 13
+DEAL:: DoF 14
+DEAL:: DoF 15
+DEAL:: DoF 16
+DEAL:: DoF 17
+DEAL:: DoF 18
+DEAL:: DoF 19
+DEAL:: DoF 20
+DEAL:: DoF 21
+DEAL:: DoF 22
+DEAL:: DoF 23
+DEAL::Checking for component 6
+DEAL:: DoF 24
+DEAL:: DoF 25
+DEAL:: DoF 26
+DEAL:: DoF 27
+DEAL:: DoF 28
+DEAL:: DoF 29
+DEAL:: DoF 30
+DEAL:: DoF 31
+DEAL::DoFs per component: 12 12 12 12 12 12 8