From 9e0809c354f9c78360713e6cb888990977980cec Mon Sep 17 00:00:00 2001
From: bangerth
Date: Tue, 25 Sep 2012 03:04:45 +0000
Subject: [PATCH] Merge branch_component_mask into mainline.
git-svn-id: https://svn.dealii.org/trunk@26700 0785d39b-7218-0410-832d-ea1e28bc413d
---
deal.II/doc/doxygen/headers/glossary.h | 146 +-
deal.II/doc/news/changes.h | 30 +-
deal.II/include/deal.II/base/event.h | 1 +
deal.II/include/deal.II/dofs/dof_tools.h | 135 +-
deal.II/include/deal.II/fe/block_mask.h | 417 +
deal.II/include/deal.II/fe/component_mask.h | 410 +
deal.II/include/deal.II/fe/fe.h | 166 +-
.../deal.II/fe/fe_dg_vector.templates.h | 4 +-
deal.II/include/deal.II/fe/fe_poly.h | 2 +-
.../include/deal.II/fe/fe_poly.templates.h | 2 +-
.../deal.II/fe/fe_poly_face.templates.h | 2 +-
deal.II/include/deal.II/fe/fe_poly_tensor.h | 2 +-
deal.II/include/deal.II/fe/fe_system.h | 12 +-
deal.II/include/deal.II/fe/fe_values.h | 482 +-
.../include/deal.II/fe/fe_values_extractors.h | 175 +
deal.II/include/deal.II/hp/dof_handler.h | 5 +-
deal.II/include/deal.II/hp/fe_collection.h | 193 +
.../deal.II/lac/constraint_matrix.templates.h | 2 +-
.../deal.II/multigrid/mg_constrained_dofs.h | 4 +-
deal.II/include/deal.II/multigrid/mg_tools.h | 6 +-
.../deal.II/multigrid/mg_transfer_block.h | 2 +
.../deal.II/multigrid/mg_transfer_component.h | 4 +-
.../deal.II/numerics/error_estimator.h | 32 +-
.../deal.II/numerics/point_value_history.h | 22 +-
.../include/deal.II/numerics/vector_tools.h | 16 +-
.../deal.II/numerics/vector_tools.templates.h | 101 +-
deal.II/source/dofs/dof_renumbering.cc | 6 +-
deal.II/source/dofs/dof_tools.cc | 751 +-
deal.II/source/dofs/dof_tools.inst.in | 82 +-
deal.II/source/fe/block_mask.cc | 50 +
deal.II/source/fe/component_mask.cc | 50 +
deal.II/source/fe/fe.cc | 184 +-
deal.II/source/fe/fe_abf.cc | 4 +-
deal.II/source/fe/fe_bdm.cc | 4 +-
deal.II/source/fe/fe_dgp.cc | 2 +-
deal.II/source/fe/fe_dgp_monomial.cc | 2 +-
deal.II/source/fe/fe_dgp_nonparametric.cc | 2 +-
deal.II/source/fe/fe_dgq.cc | 4 +-
deal.II/source/fe/fe_nedelec.cc | 2 +-
deal.II/source/fe/fe_nothing.cc | 2 +-
deal.II/source/fe/fe_poly_tensor.cc | 11 +-
deal.II/source/fe/fe_q.cc | 4 +-
deal.II/source/fe/fe_q_hierarchical.cc | 2 +-
deal.II/source/fe/fe_raviart_thomas.cc | 4 +-
deal.II/source/fe/fe_raviart_thomas_nodal.cc | 4 +-
deal.II/source/fe/fe_system.cc | 21 +-
deal.II/source/fe/fe_values.cc | 236 +-
deal.II/source/hp/fe_collection.cc | 175 +
deal.II/source/multigrid/mg_tools.cc | 60 +-
deal.II/source/multigrid/mg_tools.inst.in | 10 +-
.../source/multigrid/mg_transfer_component.cc | 164 +-
deal.II/source/multigrid/multigrid.cc | 4 +-
deal.II/source/numerics/error_estimator.cc | 70 +-
.../source/numerics/error_estimator.inst.in | 16 +-
.../source/numerics/point_value_history.cc | 125 +-
deal.II/source/numerics/vector_tools.inst.in | 16 +-
tests/Makefile.paths.in | 2 +-
tests/bits/anna_3.cc | 34 +-
tests/bits/dof_tools_12.cc | 8 +-
tests/bits/dof_tools_13a.cc | 22 +-
tests/deal.II/extract_dofs_by_component_01.cc | 2 +-
.../extract_dofs_by_component_01_hp.cc | 2 +-
.../extract_dofs_by_component_01_mg.cc | 2 +-
.../deal.II/extract_dofs_by_component_01x.cc | 84 +
.../extract_dofs_by_component_01x/cmp/generic | 2 +
.../deal.II/extract_dofs_by_component_01y.cc | 88 +
.../extract_dofs_by_component_01y/cmp/generic | 5 +
tests/deal.II/extract_dofs_by_component_02.cc | 2 +-
.../extract_dofs_by_component_02_hp.cc | 2 +-
.../extract_dofs_by_component_02_mg.cc | 4 +-
.../deal.II/memory_consumption_01/cmp/generic | 24 +-
tests/fe/block_mask_01.cc | 50 +
tests/fe/block_mask_01/cmp/generic | 2 +
tests/fe/block_mask_02.cc | 59 +
tests/fe/block_mask_02/cmp/generic | 3 +
tests/fe/block_mask_03.cc | 53 +
tests/fe/block_mask_03/cmp/generic | 2 +
tests/fe/block_mask_04.cc | 48 +
tests/fe/block_mask_04/cmp/generic | 2 +
tests/fe/block_mask_05.cc | 56 +
tests/fe/block_mask_05/cmp/generic | 2 +
tests/fe/block_mask_06.cc | 64 +
tests/fe/block_mask_06/cmp/generic | 4 +
tests/fe/block_mask_07.cc | 68 +
tests/fe/block_mask_07/cmp/generic | 5 +
tests/fe/block_mask_08.cc | 56 +
tests/fe/block_mask_08/cmp/generic | 2 +
tests/fe/block_mask_09.cc | 60 +
tests/fe/block_mask_09/cmp/generic | 2 +
tests/fe/block_mask_10.cc | 60 +
tests/fe/block_mask_10/cmp/generic | 2 +
tests/fe/block_mask_11.cc | 70 +
tests/fe/block_mask_11/cmp/generic | 2 +
tests/fe/component_mask_01.cc | 50 +
tests/fe/component_mask_01/cmp/generic | 2 +
tests/fe/component_mask_02.cc | 59 +
tests/fe/component_mask_02/cmp/generic | 3 +
tests/fe/component_mask_03.cc | 53 +
tests/fe/component_mask_03/cmp/generic | 2 +
tests/fe/component_mask_04.cc | 48 +
tests/fe/component_mask_04/cmp/generic | 2 +
tests/fe/component_mask_05.cc | 56 +
tests/fe/component_mask_05/cmp/generic | 2 +
tests/fe/component_mask_06.cc | 64 +
tests/fe/component_mask_06/cmp/generic | 4 +
tests/fe/component_mask_07.cc | 68 +
tests/fe/component_mask_07/cmp/generic | 5 +
tests/fe/component_mask_08.cc | 56 +
tests/fe/component_mask_08/cmp/generic | 2 +
tests/fe/component_mask_09.cc | 60 +
tests/fe/component_mask_09/cmp/generic | 2 +
tests/fe/component_mask_10.cc | 60 +
tests/fe/component_mask_10/cmp/generic | 2 +
tests/fe/component_mask_11.cc | 70 +
tests/fe/component_mask_11/cmp/generic | 2 +
tests/fe/component_mask_12.cc | 80 +
tests/fe/component_mask_12/cmp/generic | 2 +
tests/fe/component_mask_13.cc | 67 +
tests/fe/component_mask_13/cmp/generic | 2 +
tests/fe/component_mask_14.cc | 69 +
tests/fe/component_mask_14/cmp/generic | 34 +
tests/fe/fe_restriction_dgq/cmp/generic | 42 +-
tests/fe/fe_restriction_rt/cmp/generic | 118 +-
tests/fe/fe_restriction_sys_05/cmp/generic | 2 +-
tests/fe/fe_restriction_sys_06/cmp/generic | 8 +-
tests/fe/fe_restriction_sys_07/cmp/generic | 2 +-
tests/fe/fe_restriction_sys_08/cmp/generic | 6 +-
tests/fe/fe_restriction_sys_09/cmp/generic | 2 +-
tests/fe/fe_restriction_sys_10/cmp/generic | 6 +-
tests/fe/injection_q.cc | 3 +-
tests/fe/injection_q/cmp/generic | 70728 ++++++++--------
tests/fe/rt_5.cc | 4 +-
tests/fe/rt_5/cmp/generic | 1952 +-
133 files changed, 41295 insertions(+), 37899 deletions(-)
create mode 100644 deal.II/include/deal.II/fe/block_mask.h
create mode 100644 deal.II/include/deal.II/fe/component_mask.h
create mode 100644 deal.II/include/deal.II/fe/fe_values_extractors.h
create mode 100644 deal.II/source/fe/block_mask.cc
create mode 100644 deal.II/source/fe/component_mask.cc
create mode 100644 tests/deal.II/extract_dofs_by_component_01x.cc
create mode 100644 tests/deal.II/extract_dofs_by_component_01x/cmp/generic
create mode 100644 tests/deal.II/extract_dofs_by_component_01y.cc
create mode 100644 tests/deal.II/extract_dofs_by_component_01y/cmp/generic
create mode 100644 tests/fe/block_mask_01.cc
create mode 100644 tests/fe/block_mask_01/cmp/generic
create mode 100644 tests/fe/block_mask_02.cc
create mode 100644 tests/fe/block_mask_02/cmp/generic
create mode 100644 tests/fe/block_mask_03.cc
create mode 100644 tests/fe/block_mask_03/cmp/generic
create mode 100644 tests/fe/block_mask_04.cc
create mode 100644 tests/fe/block_mask_04/cmp/generic
create mode 100644 tests/fe/block_mask_05.cc
create mode 100644 tests/fe/block_mask_05/cmp/generic
create mode 100644 tests/fe/block_mask_06.cc
create mode 100644 tests/fe/block_mask_06/cmp/generic
create mode 100644 tests/fe/block_mask_07.cc
create mode 100644 tests/fe/block_mask_07/cmp/generic
create mode 100644 tests/fe/block_mask_08.cc
create mode 100644 tests/fe/block_mask_08/cmp/generic
create mode 100644 tests/fe/block_mask_09.cc
create mode 100644 tests/fe/block_mask_09/cmp/generic
create mode 100644 tests/fe/block_mask_10.cc
create mode 100644 tests/fe/block_mask_10/cmp/generic
create mode 100644 tests/fe/block_mask_11.cc
create mode 100644 tests/fe/block_mask_11/cmp/generic
create mode 100644 tests/fe/component_mask_01.cc
create mode 100644 tests/fe/component_mask_01/cmp/generic
create mode 100644 tests/fe/component_mask_02.cc
create mode 100644 tests/fe/component_mask_02/cmp/generic
create mode 100644 tests/fe/component_mask_03.cc
create mode 100644 tests/fe/component_mask_03/cmp/generic
create mode 100644 tests/fe/component_mask_04.cc
create mode 100644 tests/fe/component_mask_04/cmp/generic
create mode 100644 tests/fe/component_mask_05.cc
create mode 100644 tests/fe/component_mask_05/cmp/generic
create mode 100644 tests/fe/component_mask_06.cc
create mode 100644 tests/fe/component_mask_06/cmp/generic
create mode 100644 tests/fe/component_mask_07.cc
create mode 100644 tests/fe/component_mask_07/cmp/generic
create mode 100644 tests/fe/component_mask_08.cc
create mode 100644 tests/fe/component_mask_08/cmp/generic
create mode 100644 tests/fe/component_mask_09.cc
create mode 100644 tests/fe/component_mask_09/cmp/generic
create mode 100644 tests/fe/component_mask_10.cc
create mode 100644 tests/fe/component_mask_10/cmp/generic
create mode 100644 tests/fe/component_mask_11.cc
create mode 100644 tests/fe/component_mask_11/cmp/generic
create mode 100644 tests/fe/component_mask_12.cc
create mode 100644 tests/fe/component_mask_12/cmp/generic
create mode 100644 tests/fe/component_mask_13.cc
create mode 100644 tests/fe/component_mask_13/cmp/generic
create mode 100644 tests/fe/component_mask_14.cc
create mode 100644 tests/fe/component_mask_14/cmp/generic
diff --git a/deal.II/doc/doxygen/headers/glossary.h b/deal.II/doc/doxygen/headers/glossary.h
index 07894b48a3..c3dfe33c77 100644
--- a/deal.II/doc/doxygen/headers/glossary.h
+++ b/deal.II/doc/doxygen/headers/glossary.h
@@ -197,9 +197,80 @@
* More information on this topic can be found in the documentation of
* FESystem, the @ref vector_valued module and the tutorial programs
* referenced therein.
+ *
+ * Selecting blocks:
+ * Many functions allow you to restrict their operation to certain
+ * vector components or blocks. For example, this is the case for
+ * the functions that interpolate boundary values: one may want
+ * to only interpolate the boundary values for the velocity block of
+ * a finite element field but not the pressure block. The way to do
+ * this is by passing a BlockMask argument to such functions, see the
+ * @ref GlossBlockMask "block mask entry of this glossary".
*
*
*
+ * @anchor GlossBlockMask Block mask
+ *
+ *
+ * In much the same way as one can think of elements as being composed
+ * of physical vector components (see @ref GlossComponent) or logical
+ * blocks (see @ref GlossBlock), there is frequently a need to select
+ * a set of such blocks for operations that are not intended to be run
+ * on all blocks of a finite element space. Selecting which blocks
+ * to work on happens using the BlockMask class.
+ *
+ * Block masks work in much the same way as component masks, including the
+ * fact that the BlockMask class has similar semantics to the ComponentMask
+ * class. See @ref GlossComponentMask "the glossary entry on component masks"
+ * for more information.
+ *
+ * @note While components and blocks provide two alternate but equally valid
+ * viewpoints on finite elements with multiple vector components, the fact
+ * is that throughout the library there are far more places where you can
+ * pass a ComponentMask argument rather than a BlockMask argument. Fortunately,
+ * one can be converted into the other, using the syntax
+ * fe.component_mask(block_mask)
where block_mask
+ * is a variable of type BlockMask. In other words, if you have a block
+ * mask but need to call a function that only accepts a component mask, this
+ * syntax can be used to obtain the necessary component mask.
+ *
+ * Creation of block masks:
+ * Block masks are typically created by asking the finite element
+ * to generate a block mask from certain selected vector components using
+ * code such as this where we create a mask that only denotes the
+ * velocity components of a Stokes element (see @ref vector_valued):
+ * @code
+ * FESystem stokes_fe (FESystem(FE_Q(2), dim), 1, // Q2 element for the velocities
+ * FE_Q(1), 1); // Q1 element for the pressure
+ * FEValuesExtractors::Scalar pressure(dim);
+ * BlockMask pressure_mask = stokes_fe.block_mask (pressure);
+ * @endcode
+ * The result is a block mask that, in 1d as well as 2d and 3d, would have values
+ * [false, true]
. Similarly, using
+ * @code
+ * FEValuesExtractors::Vector velocities(0);
+ * BlockMask velocity_mask = stokes_fe.block_mask (velocities);
+ * @endcode
+ * would result in a mask [true, false]
in any dimension.
+ *
+ * Note, however, that if we had defined the finite element in the following
+ * way:
+ * @code
+ * FESystem stokes_fe (FE_Q(2), dim, // Q2 element for the velocities
+ * FE_Q(1), 1); // Q1 element for the pressure
+ * @endcode
+ * then the code
+ * @code
+ * FEValuesExtractors::Scalar pressure(dim);
+ * BlockMask pressure_mask = stokes_fe.block_mask (pressure);
+ * @endcode
+ * would yield a block mask that in 2d has elements
+ * [false, false, true]
because the element has
+ * dim+1
components and equally many blocks. See the
+ * discussion on what a block represents exactly in the
+ * @ref GlossBlock "block entry of this glossary".
+ *
+ *
* @anchor GlossBoundaryForm %Boundary form
*
* For a dim-dimensional triangulation in dim-dimensional space,
@@ -238,7 +309,8 @@
* equation considered in step-22 has four elements: $u=(v_x,v_y,v_z,p)^T$. We
* call the elements of the vector-valued solution components in
* deal.II. To be well-posed, for the solution to have $n$ components, there
- * need to be $n$ partial differential equations to describe them.
+ * need to be $n$ partial differential equations to describe them. This
+ * concept is discussed in great detail in the @ref vector_valued module.
*
* In finite element programs, one frequently wants to address individual
* elements (components) of this vector-valued solution, or sets of
@@ -262,13 +334,26 @@
* terms of blocks is most frequently the better strategy.
*
* For a given finite element, the number of components can be queried using
- * the FiniteElementData::n_components() function. Individual components of a
+ * the FiniteElementData::n_components() function, and you can find out
+ * which vector components are nonzero for a given finite element shape
+ * function using FiniteElement::get_nonzero_components(). The values and
+ * gradients of individual components of a
* shape function (if the element is primitive) can be queried using the
* FiniteElement::shape_value_component() and
* FiniteElement::shape_grad_component() functions on the reference cell. The
* FEValues::shape_value_component() and FEValues::shape_grad_component()
* functions do the same on a real cell. See also the documentation of the
- * FiniteElement and FEValues classes.
+ * FiniteElement and FEValues classes.
+ *
+ * Selecting components:
+ * Many functions allow you to restrict their operation to certain
+ * vector components or blocks. For example, this is the case for
+ * the functions that interpolate boundary values: one may want
+ * to only interpolate the boundary values for the velocity components of
+ * a finite element field but not the pressure component. The way to do
+ * this is by passing a ComponentMask argument to such functions, see the
+ * @ref GlossComponentMask "component mask entry of this glossary".
+ *
*
*
* @anchor GlossComponentMask Component mask
@@ -280,8 +365,8 @@
* wish to only interpolate boundary values for the velocity components
* but not the pressure. In deal.II, this is typically done by passing functions
* a component mask. Component masks are always specified as a
- * std::vector@
, i.e., an array with as many entries as the
- * finite element has components (e.g., in the Stokes case, there are
+ * ComponentMask object which one can think of as an array with
+ * as many entries as the finite element has components (e.g., in the Stokes case, there are
* dim+1
components) and where each entry is either true or false.
* In the example where we would like to interpolate boundary values only for
* the velocity components of the Stokes system, this component mask would then
@@ -296,12 +381,55 @@
* functions with these names but only some of them have a component mask
* argument.
*
- * Many of the functions that take a component mask accept a vector of length
- * zero to indicate all components, i.e., as if the vector had the
+ * Semantics of component masks:
+ * Many of the functions that take a component mask object that has been default
+ * constructed to indicate all components, i.e., as if the vector had the
* correct length and was filled with only true
values. The reason
- * is that the empty vector can be constructed in place using the code snippet
- * std::vector@()
and can thus be used as a default
+ * is that default initialized objects can be constructed in place using the code snippet
+ * ComponentMask()
and can thus be used as a default
* argument in function signatures.
+ *
+ * In other words, ComponentMask objects can be in one of two states: They can have
+ * been initialized by a vector of booleans with a nonzero length; in that case,
+ * they represent a mask of a particular length where some elements may be true
+ * and others may be false. Or, the ComponentMask may have been default initialized
+ * (using the default constructor) in which case it represents an array of indefinite
+ * length (i.e., a length appropriate to the circumstances) in which every entry
+ * is true.
+ *
+ * Creation of component masks:
+ * Component masks are typically created by asking the finite element
+ * to generate a component mask from certain selected components using
+ * code such as this where we create a mask that only denotes the
+ * velocity components of a Stokes element (see @ref vector_valued):
+ * @code
+ * FESystem stokes_fe (FE_Q(2), dim, // Q2 element for the velocities
+ * FE_Q(1), 1); // Q1 element for the pressure
+ * FEValuesExtractors::Scalar pressure(dim);
+ * ComponentMask pressure_mask = stokes_fe.component_mask (pressure);
+ * @endcode
+ * The result is a component mask that, in 2d, would have values
+ * [false, false, true]
. Similarly, using
+ * @code
+ * FEValuesExtractors::Vector velocities(0);
+ * ComponentMask velocity_mask = stokes_fe.component_mask (velocities);
+ * @endcode
+ * would result in a mask [true, true, false]
in 2d. Of
+ * course, in 3d, the result would be [true, true, true, false]
.
+ *
+ * @note Just as one can think of composed elements as being made up of
+ * @ref GlossComponent "components" or @ref GlossBlock "blocks", there are
+ * component masks (represented by the ComponentMask class) and
+ * @ref GlossBlockMask "block masks" (represented by the BlockMask class).
+ * The FiniteElement class has functions that convert between the two kinds of
+ * objects.
+ *
+ * @note Not all component masks actually make sense. For example, if you have
+ * a FE_RaviartThomas object in 2d, then it doesn't make any sense to have a
+ * component mask of the form [true, false]
because you try to
+ * select individual vector components of a finite element where each shape
+ * function has both $x$ and $y$ velocities. In essence, while you can of
+ * course create such a component mask, there is nothing you can do with it.
*
*
*
diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 6811424846..d935b2f623 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -24,12 +24,40 @@ inconvenience this causes.
+- New: In the past, deal.II used a std::vector of bools in many places
+to denote component masks (see @ref GlossComponentMask) as well as for
+block masks (see @ref GlossBlockMask). This was neither
+descriptive (the data type does not indicate what it is supposed to
+represent, nor whether the proper size for such an argument would be equal
+to the number of degrees of freedom per cell, the number of vector components
+of the finite element, or the number of degrees of freedom in total).
+
+There are now new class types ComponentMask and BlockMask that are used in these places.
+They are used both descriptively (as a return type of the function
+FiniteElement::get_nonzero_components indicating the vector components within
+which a given shape function is nonzero) as well as prescriptively (as
+input arguments to functions such as those listed in the glossary entry
+linked to above).
+
+While the descriptive places are not backward compatible (they return a
+ComponentMask which is not convertible to the std::vector of bools returned
+before), most of the prescriptive places are in fact backward compatible
+(because the std::vector of bool that was passed previously can
+implicitly be converted to an object of type ComponentMask. The sole
+exception is the function DoFTools::extract_dofs (and its multigrid
+equivalent DoFTools::extract_level_dofs) that previously
+could interpret its argument as either a component or
+a block mask, depending on a boolean flag. This function now exists
+in two different versions, one that takes a ComponentMask and
+one that takes a BlockMask. Call sites need to be adjusted.
+
+(Wolfgang Bangerth, 2012/09/22)
+
- Changed: the optional argument offset got removed from
DoFHandler and MGDoFHandler::distribute_dofs() because it was
never working correctly and it is not used.
(Timo Heister, 2012/09/03)
-
diff --git a/deal.II/include/deal.II/base/event.h b/deal.II/include/deal.II/base/event.h
index 33ac960098..c691c804d1 100644
--- a/deal.II/include/deal.II/base/event.h
+++ b/deal.II/include/deal.II/base/event.h
@@ -131,6 +131,7 @@ namespace Algorithms
/**
* The names of registered events
*/
+//TODO: This static field must be guarded by a mutex to be thread-safe!
static std::vector names;
};
diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h
index 22a9d90963..2820c3a8a9 100644
--- a/deal.II/include/deal.II/dofs/dof_tools.h
+++ b/deal.II/include/deal.II/dofs/dof_tools.h
@@ -23,6 +23,7 @@
#include
#include
#include
+#include
#include
#include
@@ -990,9 +991,8 @@ namespace DoFTools
* components of the finite element space
* shall be constrained with periodic
* boundary conditions. If it is left as
- * specified by the default value (i.e. an
- * empty array), all components are
- * constrainted. If it is different from
+ * specified by the default value all components are
+ * constrained. If it is different from
* the default value, it is assumed that
* the number of entries equals the number
* of components in the boundary functions
@@ -1007,7 +1007,7 @@ namespace DoFTools
make_periodicity_constraints (const FaceIterator &face_1,
const typename identity::type &face_2,
dealii::ConstraintMatrix &constraint_matrix,
- const std::vector &component_mask = std::vector());
+ const ComponentMask &component_mask = ComponentMask());
/**
@@ -1069,9 +1069,8 @@ namespace DoFTools
* components of the finite element space
* shall be constrained with periodic
* boundary conditions. If it is left as
- * specified by the default value (i.e. an
- * empty array), all components are
- * constrainted. If it is different from
+ * specified by the default value all components are
+ * constrained. If it is different from
* the default value, it is assumed that
* the number of entries equals the number
* of components in the boundary functions
@@ -1087,7 +1086,7 @@ namespace DoFTools
const types::boundary_id boundary_component,
const int direction,
dealii::ConstraintMatrix &constraint_matrix,
- const std::vector &component_mask = std::vector());
+ const ComponentMask &component_mask = ComponentMask());
/**
@@ -1109,7 +1108,7 @@ namespace DoFTools
dealii::Tensor<1,DH::space_dimension>
&offset,
dealii::ConstraintMatrix &constraint_matrix,
- const std::vector &component_mask = std::vector());
+ const ComponentMask &component_mask = ComponentMask());
//@}
@@ -1178,34 +1177,31 @@ namespace DoFTools
* Extract the indices of the
* degrees of freedom belonging
* to certain vector components
- * or blocks (if the last
- * argument is true) of
+ * of
* a vector-valued finite
- * element. The bit vector @p
- * select defines, which
+ * element. The @p component_mask
+ * defines which
* components or blocks of an
* FESystem are to be extracted
* from the DoFHandler @p
- * dof. The entries in @p
+ * dof. The entries in the output array @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 select must
- * equal the number of components
- * or blocks in the FiniteElement
- * used by @p dof, depending on
- * the argument
- * blocks. The size of
+ * The size of @p component_mask must
+ * be compatible with the number of components
+ * in the FiniteElement
+ * used by @p dof. The size of
* @p selected_dofs must equal
* DoFHandler::n_dofs(). Previous
* contents of this array are
* overwritten.
*
* If the finite element under
- * consideration is not primitive, that is
+ * consideration is not primitive, i.e.,
* some or all of its shape functions are
* non-zero in more than one vector
* component (which holds, for example, for
@@ -1219,37 +1215,78 @@ namespace DoFTools
* equivalent to selecting all
* vector components corresponding to this
* non-primitive base element.
+ *
+ * @note If the @p blocks argument is
+ * true,
*/
template
void
extract_dofs (const DoFHandler &dof_handler,
- const std::vector &select,
- std::vector &selected_dofs,
- const bool blocks = false);
+ const ComponentMask &component_mask,
+ std::vector &selected_dofs);
- /**
- * The same function as above,
- * but for a hp::DoFHandler.
- */
+ /**
+ * The same function as above,
+ * but for a hp::DoFHandler.
+ */
template
void
extract_dofs (const hp::DoFHandler &dof_handler,
- const std::vector &select,
- std::vector &selected_dofs,
- const bool blocks = false);
+ const ComponentMask &component_mask,
+ std::vector &selected_dofs);
- /**
- * Do the same thing as
- * extract_dofs() for one level
- * of a multi-grid DoF numbering.
- */
+ /**
+ * This function is the equivalent to the DoFTools::extract_dofs() functions above
+ * except that the selection of which degrees of freedom to extract is not done
+ * based on components (see @ref GlossComponent) but instead based on whether they
+ * are part of a particular block (see @ref GlossBlock). Consequently, the second
+ * argument is not a ComponentMask but a BlockMask object.
+ *
+ * @param dof_handler The DoFHandler object from which to extract degrees of freedom
+ * @param block_mask The block mask that describes which blocks to consider (see
+ * @ref GlossBlockMask)
+ * @param selected_dofs A vector of length DoFHandler::n_dofs() in which those
+ * entries are true that correspond to the selected blocks.
+ */
+ template
+ void
+ extract_dofs (const DoFHandler &dof_handler,
+ const BlockMask &block_mask,
+ std::vector &selected_dofs);
+
+ /**
+ * The same function as above,
+ * but for a hp::DoFHandler.
+ */
+ template
+ void
+ extract_dofs (const hp::DoFHandler &dof_handler,
+ const BlockMask &block_mask,
+ std::vector &selected_dofs);
+
+ /**
+ * Do the same thing as the corresponding
+ * extract_dofs() functino for one level
+ * of a multi-grid DoF numbering.
+ */
+ template
+ void
+ extract_level_dofs (const unsigned int level,
+ const MGDoFHandler &dof,
+ const ComponentMask &component_mask,
+ std::vector &selected_dofs);
+
+ /**
+ * Do the same thing as the corresponding
+ * extract_dofs() functino for one level
+ * of a multi-grid DoF numbering.
+ */
template
void
extract_level_dofs (const unsigned int level,
const MGDoFHandler &dof,
- const std::vector &select,
- std::vector &selected_dofs,
- const bool blocks = false);
+ const BlockMask &component_mask,
+ std::vector &selected_dofs);
/**
* Extract all degrees of freedom
@@ -1333,7 +1370,7 @@ namespace DoFTools
template
void
extract_boundary_dofs (const DH &dof_handler,
- const std::vector &component_mask,
+ const ComponentMask &component_mask,
std::vector &selected_dofs,
const std::set &boundary_indicators = std::set());
@@ -1369,8 +1406,8 @@ namespace DoFTools
template
void
extract_boundary_dofs (const DH &dof_handler,
- const std::vector &component_mask,
- IndexSet &selected_dofs,
+ const ComponentMask &component_mask,
+ IndexSet &selected_dofs,
const std::set &boundary_indicators = std::set());
/**
@@ -1405,9 +1442,9 @@ namespace DoFTools
*/
template
void
- extract_dofs_with_support_on_boundary (const DH &dof_handler,
- const std::vector &component_mask,
- std::vector &selected_dofs,
+ extract_dofs_with_support_on_boundary (const DH &dof_handler,
+ const ComponentMask &component_mask,
+ std::vector &selected_dofs,
const std::set &boundary_indicators = std::set());
/**
@@ -1732,7 +1769,7 @@ namespace DoFTools
/**
* Create a sparsity pattern which
* in each row lists the degrees of
- * freedom associated to the
+ * freedom associated with the
* corresponding cell.
*
* Ordering follows the ordering of
@@ -1957,7 +1994,7 @@ namespace DoFTools
template
void
extract_constant_modes (const DH &dof_handler,
- const std::vector &component_mask,
+ const ComponentMask &component_mask,
std::vector > &constant_modes);
/**
@@ -2597,8 +2634,8 @@ namespace DoFTools
template class DH>
void
make_zero_boundary_constraints (const DH &dof,
- ConstraintMatrix &zero_boundary_constraints,
- const std::vector &component_mask_=std::vector());
+ ConstraintMatrix &zero_boundary_constraints,
+ const ComponentMask &component_mask = ComponentMask());
/**
* Map a coupling table from the
diff --git a/deal.II/include/deal.II/fe/block_mask.h b/deal.II/include/deal.II/fe/block_mask.h
new file mode 100644
index 0000000000..5ffe92aa12
--- /dev/null
+++ b/deal.II/include/deal.II/fe/block_mask.h
@@ -0,0 +1,417 @@
+//---------------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2009, 2010, 2011, 2012 by the deal.II authors
+//
+// 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__fe_block_mask_h
+#define __deal2__fe_block_mask_h
+
+#include
+#include
+#include
+
+#include
+#ifdef HAVE_STD_IOSFWD_HEADER
+# include
+#else
+# include
+#endif
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+/**
+ * This class represents a mask that can be used to select individual
+ * vector blocks of a finite element (see also
+ * @ref GlossBlockMask "this glossary entry"). It will typically have as many
+ * elements as the finite element has blocks, and one can
+ * use operator[]
to query whether a particular block
+ * has been selected.
+ *
+ * The semantics of this class are the same as the related ComponentMask
+ * class, i.e., a default constructed mask represents all possible
+ * blocks. See there for more information about these semantics.
+ *
+ * Objects of this kind are used in many places where one wants to restrict
+ * operations to a certain subset of blocks, e.g. in DoFTools::extract_dofs.
+ * These objects can
+ * either be created by hand, or, simpler, by asking the finite element
+ * to generate a block mask from certain selected blocks using
+ * code such as this where we create a mask that only denotes the
+ * velocity block of a Stokes element (see @ref vector_valued):
+ * @code
+ * FESystem stokes_fe (FESystem(FE_Q(2), dim), 1, // Q2 element for the velocities
+ * FE_Q(1), 1); // Q1 element for the pressure
+ * FEValuesExtractors::Scalar pressure(dim);
+ * BlockMask pressure_mask = stokes_fe.block_mask (pressure);
+ * @endcode
+ * Note that by wrapping the velocity elements into a single FESystem
+ * object we make sure that the overall element has only 2 blocks.
+ * The result is a block mask that, in both 2d and 3d, would have values
+ * [false, true]
. (Compare this to the corresponding
+ * component mask discussed in the ComponentMask documentation.)
+ * Similarly, using
+ * @code
+ * FEValuesExtractors::Vector velocities(0);
+ * BlockMask velocity_mask = stokes_fe.block_mask (velocities);
+ * @endcode
+ * would result in a mask [true, false]
in both 2d and 3d.
+ *
+ * @ingroup fe
+ * @author Wolfgang Bangerth
+ * @date 2012
+ * @ingroup vector_valued
+ */
+class BlockMask
+{
+ public:
+ /**
+ * Initialize a block mask. The default is that a block
+ * mask represents a set of blocks that are all
+ * selected, i.e., calling this constructor results in
+ * a block mask that always returns true
+ * whenever asked whether a block is selected.
+ */
+ BlockMask ();
+
+ /**
+ * Initialize an object of this type with a set of selected
+ * blocks specified by the argument.
+ *
+ * @param block_mask A vector of true/false
+ * entries that determine which blocks of a finite element
+ * are selected. If the length of the given vector is zero,
+ * then this interpreted as the case where every block
+ * is selected.
+ */
+ BlockMask (const std::vector &block_mask);
+
+ /**
+ * Initialize the block mask with a number of elements that
+ * are either all true or false.
+ *
+ * @param n_blocks The number of elements of this mask
+ * @param initializer The value each of these elements is supposed to have:
+ * either true or false.
+ */
+ BlockMask (const unsigned int n_blocks,
+ const bool initializer);
+
+ /**
+ * If this block mask has been initialized with a mask of
+ * size greater than zero, then return the size of the mask
+ * represented by this object. On the other hand, if this
+ * mask has been initialized as an empty object that represents
+ * a mask that is true for every element (i.e., if this object
+ * would return true when calling represents_the_all_selected_mask())
+ * then return zero since no definite size is known.
+ */
+ unsigned int size () const;
+
+ /**
+ * Return whether a particular block is selected by this
+ * mask. If this mask represents the case of an object that
+ * selects all blocks (e.g. if it is created
+ * using the default constructor or is converted from an
+ * empty vector of type bool) then this function returns
+ * true regardless of the given argument.
+ *
+ * @param block_index The index for which the function
+ * should return whether the block is selected. If this
+ * object represents a mask in which all blocks are always
+ * selected then any index is allowed here. Otherwise, the
+ * given index needs to be between zero and the number of
+ * blocks that this mask represents.
+ */
+ bool operator[] (const unsigned int block_index) const;
+
+ /**
+ * Return whether this block mask represents a mask with
+ * exactly n
blocks. This is true if either
+ * it was initilized with a vector with exactly n
+ * entries of type bool
(in this case, @p n must
+ * equal the result of size()) or if it was initialized
+ * with an empty vector (or using the default constructor) in
+ * which case it can represent a mask with an arbitrary number
+ * of blocks and will always say that a block is
+ * selected.
+ */
+ bool
+ represents_n_blocks (const unsigned int n) const;
+
+ /**
+ * Return the number of blocks that are selected by this mask.
+ *
+ * Since empty block masks represent a block mask that
+ * would return true
for every block, this
+ * function may not know the true size of the block
+ * mask and it therefore requires an argument that denotes the
+ * overall number of blocks.
+ *
+ * If the object has been initialized with a non-empty mask (i.e.,
+ * if the size() function returns something greater than zero,
+ * or equivalently if represents_the_all_selected_mask() returns
+ * false) then the argument can be omitted and the result of size()
+ * is taken.
+ */
+ unsigned int
+ n_selected_blocks (const unsigned int overall_number_of_blocks = numbers::invalid_unsigned_int) const;
+
+ /**
+ * Return the index of the first selected block. The argument
+ * is there for the same reason it exists with the
+ * n_selected_blocks() function.
+ *
+ * The function throws an exception if no block is selected at all.
+ */
+ unsigned int
+ first_selected_block (const unsigned int overall_number_of_blocks = numbers::invalid_unsigned_int) const;
+
+ /**
+ * Return true if this mask represents a default
+ * constructed mask that corresponds to one in which
+ * all blocks are selected. If true, then the size()
+ * function will return zero.
+ */
+ bool
+ represents_the_all_selected_mask () const;
+
+ /**
+ * Return a block mask that contains the union of the
+ * blocks selected by the current object and the one
+ * passed as an argument.
+ */
+ BlockMask operator | (const BlockMask &mask) const;
+
+ /**
+ * Return a block mask that has only those elements set that
+ * are set both in the current object as well as the one
+ * passed as an argument.
+ */
+ BlockMask operator & (const BlockMask &mask) const;
+
+ /**
+ * Return whether this object and the argument are identical.
+ */
+ bool operator== (const BlockMask &mask) const;
+
+ /**
+ * Determine an estimate for the
+ * memory consumption (in bytes)
+ * of this object.
+ */
+ std::size_t
+ memory_consumption () const;
+
+ private:
+ /**
+ * The actual block mask.
+ */
+ std::vector block_mask;
+
+ // make the output operator a friend so it can access
+ // the block_mask array
+ friend
+ std::ostream & operator << (std::ostream &out,
+ const BlockMask &mask);
+};
+
+
+/**
+ * Write a block mask to an output stream. If the block
+ * mask represents one where all blocks are selected without
+ * specifying a particular size of the mask, then it
+ * writes the string [all blocks selected]
to the
+ * stream. Otherwise, it prints the block mask in a form like
+ * [true,true,true,false]
.
+ *
+ * @param out The stream to write to.
+ * @param mask The mask to write.
+ * @return A reference to the first argument.
+ */
+std::ostream & operator << (std::ostream &out,
+ const BlockMask &mask);
+
+
+// -------------------- inline functions ---------------------
+
+inline
+BlockMask::BlockMask()
+{}
+
+
+inline
+BlockMask::BlockMask(const std::vector &block_mask)
+:
+block_mask (block_mask)
+{}
+
+
+inline
+BlockMask::BlockMask(const unsigned int n_blocks,
+ const bool initializer)
+:
+block_mask (n_blocks, initializer)
+{}
+
+
+inline
+unsigned int
+BlockMask::size () const
+{
+ return block_mask.size();
+}
+
+
+inline
+bool
+BlockMask::operator [](const unsigned int block_index) const
+{
+ // if the mask represents the all-block mask
+ // then always return true
+ if (block_mask.size() == 0)
+ return true;
+ else
+ {
+ // otherwise check the validity of the index and
+ // return whatever is appropriate
+ Assert (block_index < block_mask.size(),
+ ExcIndexRange (block_index, 0, block_mask.size()));
+ return block_mask[block_index];
+ }
+}
+
+
+inline
+bool
+BlockMask::represents_n_blocks(const unsigned int n) const
+{
+ return ((block_mask.size() == 0)
+ ||
+ (block_mask.size() == n));
+}
+
+
+inline
+unsigned int
+BlockMask::n_selected_blocks(const unsigned int n) const
+{
+ if ((n != numbers::invalid_unsigned_int) && (size() > 0))
+ AssertDimension (n, size());
+
+ const unsigned int real_n = (n != numbers::invalid_unsigned_int
+ ?
+ n
+ :
+ size());
+ if (block_mask.size() == 0)
+ return real_n;
+ else
+ {
+ AssertDimension (real_n, block_mask.size());
+ unsigned int c = 0;
+ for (unsigned int i=0; i 0))
+ AssertDimension (n, size());
+
+ if (block_mask.size() == 0)
+ return 0;
+ else
+ {
+ for (unsigned int c=0; c new_mask (block_mask.size());
+ for (unsigned int i=0; i new_mask (block_mask.size());
+ for (unsigned int i=0; i
+#include
+#include
+
+#include
+#ifdef HAVE_STD_IOSFWD_HEADER
+# include
+#else
+# include
+#endif
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+/**
+ * This class represents a mask that can be used to select individual
+ * vector components of a finite element (see also
+ * @ref GlossComponentMask "this glossary entry"). It will typically have as many
+ * elements as the finite element has vector components, and one can
+ * use operator[]
to query whether a particular component
+ * has been selected.
+ *
+ * Objects of this kind are used in many places where one wants to restrict
+ * operations to a certain subset of components, e.g. in DoFTools::make_zero_boundary_values
+ * or VectorTools::interpolate_boundary_values. These objects can
+ * either be created by hand, or, simpler, by asking the finite element
+ * to generate a component mask from certain selected components using
+ * code such as this where we create a mask that only denotes the
+ * velocity components of a Stokes element (see @ref vector_valued):
+ * @code
+ * FESystem stokes_fe (FE_Q(2), dim, // Q2 element for the velocities
+ * FE_Q(1), 1); // Q1 element for the pressure
+ * FEValuesExtractors::Scalar pressure(dim);
+ * ComponentMask pressure_mask = stokes_fe.component_mask (pressure);
+ * @endcode
+ * The result is a component mask that, in 2d, would have values
+ * [false, false, true]
. Similarly, using
+ * @code
+ * FEValuesExtractors::Vector velocities(0);
+ * ComponentMask velocity_mask = stokes_fe.component_mask (velocities);
+ * @endcode
+ * would result in a mask [true, true, false]
in 2d. Of
+ * course, in 3d, the result would be [true, true, true, false]
.
+ *
+ * @ingroup fe
+ * @author Wolfgang Bangerth
+ * @date 2012
+ * @ingroup vector_valued
+ */
+class ComponentMask
+{
+ public:
+ /**
+ * Initialize a component mask. The default is that a component
+ * mask represents a set of components that are all
+ * selected, i.e., calling this constructor results in
+ * a component mask that always returns true
+ * whenever asked whether a component is selected.
+ */
+ ComponentMask ();
+
+ /**
+ * Initialize an object of this type with a set of selected
+ * components specified by the argument.
+ *
+ * @param component_mask A vector of true/false
+ * entries that determine which components of a finite element
+ * are selected. If the length of the given vector is zero,
+ * then this interpreted as the case where every component
+ * is selected.
+ */
+ ComponentMask (const std::vector &component_mask);
+
+ /**
+ * Initialize the component mask with a number of elements that
+ * are either all true or false.
+ *
+ * @param n_components The number of elements of this mask
+ * @param initializer The value each of these elements is supposed to have:
+ * either true or false.
+ */
+ ComponentMask (const unsigned int n_components,
+ const bool initializer);
+
+ /**
+ * If this component mask has been initialized with a mask of
+ * size greater than zero, then return the size of the mask
+ * represented by this object. On the other hand, if this
+ * mask has been initialized as an empty object that represents
+ * a mask that is true for every element (i.e., if this object
+ * would return true when calling represents_the_all_selected_mask())
+ * then return zero since no definite size is known.
+ */
+ unsigned int size () const;
+
+ /**
+ * Return whether a particular component is selected by this
+ * mask. If this mask represents the case of an object that
+ * selects all components (e.g. if it is created
+ * using the default constructor or is converted from an
+ * empty vector of type bool) then this function returns
+ * true regardless of the given argument.
+ *
+ * @param component_index The index for which the function
+ * should return whether the component is selected. If this
+ * object represents a mask in which all components are always
+ * selected then any index is allowed here. Otherwise, the
+ * given index needs to be between zero and the number of
+ * components that this mask represents.
+ */
+ bool operator[] (const unsigned int component_index) const;
+
+ /**
+ * Return whether this component mask represents a mask with
+ * exactly n
components. This is true if either
+ * it was initilized with a vector with exactly n
+ * entries of type bool
(in this case, @p n must
+ * equal the result of size()) or if it was initialized
+ * with an empty vector (or using the default constructor) in
+ * which case it can represent a mask with an arbitrary number
+ * of components and will always say that a component is
+ * selected.
+ */
+ bool
+ represents_n_components (const unsigned int n) const;
+
+ /**
+ * Return the number of components that are selected by this mask.
+ *
+ * Since empty component masks represent a component mask that
+ * would return true
for every component, this
+ * function may not know the true size of the component
+ * mask and it therefore requires an argument that denotes the
+ * overall number of components.
+ *
+ * If the object has been initialized with a non-empty mask (i.e.,
+ * if the size() function returns something greater than zero,
+ * or equivalently if represents_the_all_selected_mask() returns
+ * false) then the argument can be omitted and the result of size()
+ * is taken.
+ */
+ unsigned int
+ n_selected_components (const unsigned int overall_number_of_components = numbers::invalid_unsigned_int) const;
+
+ /**
+ * Return the index of the first selected component. The argument
+ * is there for the same reason it exists with the
+ * n_selected_components() function.
+ *
+ * The function throws an exception if no component is selected at all.
+ */
+ unsigned int
+ first_selected_component (const unsigned int overall_number_of_components = numbers::invalid_unsigned_int) const;
+
+ /**
+ * Return true if this mask represents a default
+ * constructed mask that corresponds to one in which
+ * all components are selected. If true, then the size()
+ * function will return zero.
+ */
+ bool
+ represents_the_all_selected_mask () const;
+
+ /**
+ * Return a component mask that contains the union of the
+ * components selected by the current object and the one
+ * passed as an argument.
+ */
+ ComponentMask operator | (const ComponentMask &mask) const;
+
+ /**
+ * Return a component mask that has only those elements set that
+ * are set both in the current object as well as the one
+ * passed as an argument.
+ */
+ ComponentMask operator & (const ComponentMask &mask) const;
+
+ /**
+ * Return whether this object and the argument are identical.
+ */
+ bool operator== (const ComponentMask &mask) const;
+
+ /**
+ * Determine an estimate for the
+ * memory consumption (in bytes)
+ * of this object.
+ */
+ std::size_t
+ memory_consumption () const;
+
+ private:
+ /**
+ * The actual component mask.
+ */
+ std::vector component_mask;
+
+ // make the output operator a friend so it can access
+ // the component_mask array
+ friend
+ std::ostream & operator << (std::ostream &out,
+ const ComponentMask &mask);
+};
+
+
+/**
+ * Write a component mask to an output stream. If the component
+ * mask represents one where all components are selected without
+ * specifying a particular size of the mask, then it
+ * writes the string [all components selected]
to the
+ * stream. Otherwise, it prints the component mask in a form like
+ * [true,true,true,false]
.
+ *
+ * @param out The stream to write to.
+ * @param mask The mask to write.
+ * @return A reference to the first argument.
+ */
+std::ostream & operator << (std::ostream &out,
+ const ComponentMask &mask);
+
+
+// -------------------- inline functions ---------------------
+
+inline
+ComponentMask::ComponentMask()
+{}
+
+
+inline
+ComponentMask::ComponentMask(const std::vector &component_mask)
+:
+component_mask (component_mask)
+{}
+
+
+inline
+ComponentMask::ComponentMask(const unsigned int n_components,
+ const bool initializer)
+:
+component_mask (n_components, initializer)
+{}
+
+
+inline
+unsigned int
+ComponentMask::size () const
+{
+ return component_mask.size();
+}
+
+
+inline
+bool
+ComponentMask::operator [](const unsigned int component_index) const
+{
+ // if the mask represents the all-component mask
+ // then always return true
+ if (component_mask.size() == 0)
+ return true;
+ else
+ {
+ // otherwise check the validity of the index and
+ // return whatever is appropriate
+ Assert (component_index < component_mask.size(),
+ ExcIndexRange (component_index, 0, component_mask.size()));
+ return component_mask[component_index];
+ }
+}
+
+
+inline
+bool
+ComponentMask::represents_n_components(const unsigned int n) const
+{
+ return ((component_mask.size() == 0)
+ ||
+ (component_mask.size() == n));
+}
+
+
+inline
+unsigned int
+ComponentMask::n_selected_components(const unsigned int n) const
+{
+ if ((n != numbers::invalid_unsigned_int) && (size() > 0))
+ AssertDimension (n, size());
+
+ const unsigned int real_n = (n != numbers::invalid_unsigned_int
+ ?
+ n
+ :
+ size());
+ if (component_mask.size() == 0)
+ return real_n;
+ else
+ {
+ AssertDimension (real_n, component_mask.size());
+ unsigned int c = 0;
+ for (unsigned int i=0; i 0))
+ AssertDimension (n, size());
+
+ if (component_mask.size() == 0)
+ return 0;
+ else
+ {
+ for (unsigned int c=0; c new_mask (component_mask.size());
+ for (unsigned int i=0; i new_mask (component_mask.size());
+ for (unsigned int i=0; i
#include
#include
+#include
+#include
+#include
DEAL_II_NAMESPACE_OPEN
@@ -394,7 +397,7 @@ class FiniteElement : public Subscriptor,
*/
FiniteElement (const FiniteElementData &fe_data,
const std::vector &restriction_is_additive_flags,
- const std::vector > &nonzero_components);
+ const std::vector &nonzero_components);
/**
* Virtual destructor. Makes sure
@@ -1352,7 +1355,7 @@ class FiniteElement : public Subscriptor,
* called non-primitive (see
* @ref GlossPrimitive).
*/
- const std::vector &
+ const ComponentMask &
get_nonzero_components (const unsigned int i) const;
/**
@@ -1573,6 +1576,159 @@ class FiniteElement : public Subscriptor,
unsigned int
component_to_block_index (const unsigned int component) const;
+ //@}
+
+ /**
+ * @name Component and block matrices
+ * @{
+ */
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * one component is true that corresponds to the given
+ * argument. See @ref GlossComponentMask "the glossary"
+ * for more information.
+ *
+ * @param scalar An object that represents a single scalar
+ * vector component of this finite element.
+ * @return A component mask that is false in all components
+ * except for the one that corresponds to the argument.
+ */
+ ComponentMask
+ component_mask (const FEValuesExtractors::Scalar &scalar) const;
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * dim
components are true that correspond to the given
+ * argument. See @ref GlossComponentMask "the glossary"
+ * for more information.
+ *
+ * @param vector An object that represents dim
+ * vector components of this finite element.
+ * @return A component mask that is false in all components
+ * except for the ones that corresponds to the argument.
+ */
+ ComponentMask
+ component_mask (const FEValuesExtractors::Vector &vector) const;
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * dim*(dim+1)/2
components are true that
+ * correspond to the given argument. See @ref GlossComponentMask "the glossary"
+ * for more information.
+ *
+ * @param sym_tensor An object that represents dim*(dim+1)/2
+ * components of this finite element that are jointly to be
+ * interpreted as forming a symmetric tensor.
+ * @return A component mask that is false in all components
+ * except for the ones that corresponds to the argument.
+ */
+ ComponentMask
+ component_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;
+
+ /**
+ * Given a block mask (see @ref GlossBlockMask "this glossary entry"),
+ * produce a component mask (see @ref GlossComponentMask "this glossary entry")
+ * that represents the components that correspond to the blocks selected in
+ * the input argument. This is essentially a conversion operator from
+ * BlockMask to ComponentMask.
+ *
+ * @param block_mask The mask that selects individual blocks of the finite
+ * element
+ * @return A mask that selects those components corresponding to the selected
+ * blocks of the input argument.
+ */
+ ComponentMask
+ component_mask (const BlockMask &block_mask) const;
+
+ /**
+ * Return a block mask with as many elements as this
+ * object has blocks and of which exactly the
+ * one component is true that corresponds to the given
+ * argument. See @ref GlossBlockMask "the glossary"
+ * for more information.
+ *
+ * @note This function will only succeed if the scalar referenced
+ * by the argument encompasses a complete block. In other words,
+ * if, for example, you pass an extractor for the single
+ * $x$ velocity and this object represents an FE_RaviartThomas
+ * object, then the single scalar object you selected is part
+ * of a larger block and consequently there is no block mask that
+ * would represent it. The function will then produce an exception.
+ *
+ * @param scalar An object that represents a single scalar
+ * vector component of this finite element.
+ * @return A component mask that is false in all components
+ * except for the one that corresponds to the argument.
+ */
+ BlockMask
+ block_mask (const FEValuesExtractors::Scalar &scalar) const;
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * dim
components are true that correspond to the given
+ * argument. See @ref GlossBlockMask "the glossary"
+ * for more information.
+ *
+ * @note The same caveat applies as to the version of the function above:
+ * The extractor object passed as argument must be so that it corresponds
+ * to full blocks and does not split blocks of this element.
+ *
+ * @param vector An object that represents dim
+ * vector components of this finite element.
+ * @return A component mask that is false in all components
+ * except for the ones that corresponds to the argument.
+ */
+ BlockMask
+ block_mask (const FEValuesExtractors::Vector &vector) const;
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * dim*(dim+1)/2
components are true that
+ * correspond to the given argument. See @ref GlossBlockMask "the glossary"
+ * for more information.
+ *
+ * @note The same caveat applies as to the version of the function above:
+ * The extractor object passed as argument must be so that it corresponds
+ * to full blocks and does not split blocks of this element.
+ *
+ * @param sym_tensor An object that represents dim*(dim+1)/2
+ * components of this finite element that are jointly to be
+ * interpreted as forming a symmetric tensor.
+ * @return A component mask that is false in all components
+ * except for the ones that corresponds to the argument.
+ */
+ BlockMask
+ block_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;
+
+ /**
+ * Given a component mask (see @ref GlossComponentMask "this glossary entry"),
+ * produce a block mask (see @ref GlossBlockMask "this glossary entry")
+ * that represents the blocks that correspond to the components selected in
+ * the input argument. This is essentially a conversion operator from
+ * ComponentMask to BlockMask.
+ *
+ * @note This function will only succeed if the components referenced
+ * by the argument encompasses complete blocks. In other words,
+ * if, for example, you pass an component mask for the single
+ * $x$ velocity and this object represents an FE_RaviartThomas
+ * object, then the single component you selected is part
+ * of a larger block and consequently there is no block mask that
+ * would represent it. The function will then produce an exception.
+ *
+ * @param component_mask The mask that selects individual components of the finite
+ * element
+ * @return A mask that selects those blocks corresponding to the selected
+ * blocks of the input argument.
+ */
+ BlockMask
+ block_mask (const ComponentMask &component_mask) const;
+
//@}
/**
@@ -2213,7 +2369,7 @@ class FiniteElement : public Subscriptor,
*/
static
std::vector
- compute_n_nonzero_components (const std::vector > &nonzero_components);
+ compute_n_nonzero_components (const std::vector &nonzero_components);
/**
* Determine the values a finite
@@ -2479,7 +2635,7 @@ class FiniteElement : public Subscriptor,
* is only one non-zero
* component.
*/
- const std::vector > nonzero_components;
+ const std::vector nonzero_components;
/**
* This array holds how many
@@ -2817,7 +2973,7 @@ FiniteElement::restriction_is_additive (const unsigned int index)
template
inline
-const std::vector &
+const ComponentMask &
FiniteElement::get_nonzero_components (const unsigned int i) const
{
Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, this->dofs_per_cell));
diff --git a/deal.II/include/deal.II/fe/fe_dg_vector.templates.h b/deal.II/include/deal.II/fe/fe_dg_vector.templates.h
index ecefcc0cca..ad46bde22f 100644
--- a/deal.II/include/deal.II/fe/fe_dg_vector.templates.h
+++ b/deal.II/include/deal.II/fe/fe_dg_vector.templates.h
@@ -28,8 +28,8 @@ FE_DGVector::FE_DGVector (
FiniteElementData(
get_dpo_vector(deg), dim, deg+1, FiniteElementData::L2, 1),
std::vector(POLY::compute_n_pols(deg), true),
- std::vector >(POLY::compute_n_pols(deg),
- std::vector(dim,true)))
+ std::vector(POLY::compute_n_pols(deg),
+ ComponentMask(dim,true)))
{
this->mapping_type = map;
const unsigned int polynomial_degree = this->tensor_degree();
diff --git a/deal.II/include/deal.II/fe/fe_poly.h b/deal.II/include/deal.II/fe/fe_poly.h
index 437dee900c..221f7b8463 100644
--- a/deal.II/include/deal.II/fe/fe_poly.h
+++ b/deal.II/include/deal.II/fe/fe_poly.h
@@ -71,7 +71,7 @@ class FE_Poly : public FiniteElement
FE_Poly (const POLY& poly_space,
const FiniteElementData &fe_data,
const std::vector &restriction_is_additive_flags,
- const std::vector > &nonzero_components);
+ const std::vector &nonzero_components);
/**
* Return the polynomial degree
diff --git a/deal.II/include/deal.II/fe/fe_poly.templates.h b/deal.II/include/deal.II/fe/fe_poly.templates.h
index cc79e573ab..70195024c3 100644
--- a/deal.II/include/deal.II/fe/fe_poly.templates.h
+++ b/deal.II/include/deal.II/fe/fe_poly.templates.h
@@ -23,7 +23,7 @@ template
FE_Poly::FE_Poly (const POLY& poly_space,
const FiniteElementData &fe_data,
const std::vector &restriction_is_additive_flags,
- const std::vector > &nonzero_components):
+ const std::vector &nonzero_components):
FiniteElement (fe_data,
restriction_is_additive_flags,
nonzero_components),
diff --git a/deal.II/include/deal.II/fe/fe_poly_face.templates.h b/deal.II/include/deal.II/fe/fe_poly_face.templates.h
index 65f65a5095..232fce4ca3 100644
--- a/deal.II/include/deal.II/fe/fe_poly_face.templates.h
+++ b/deal.II/include/deal.II/fe/fe_poly_face.templates.h
@@ -25,7 +25,7 @@ FE_PolyFace::FE_PolyFace (
const std::vector &restriction_is_additive_flags):
FiniteElement (fe_data,
restriction_is_additive_flags,
- std::vector > (1, std::vector(1,true))),
+ std::vector (1, ComponentMask(1,true))),
poly_space(poly_space)
{
AssertDimension(dim, POLY::dimension+1);
diff --git a/deal.II/include/deal.II/fe/fe_poly_tensor.h b/deal.II/include/deal.II/fe/fe_poly_tensor.h
index 14562bf88f..902b295e77 100644
--- a/deal.II/include/deal.II/fe/fe_poly_tensor.h
+++ b/deal.II/include/deal.II/fe/fe_poly_tensor.h
@@ -126,7 +126,7 @@ class FE_PolyTensor : public FiniteElement
FE_PolyTensor (const unsigned int degree,
const FiniteElementData &fe_data,
const std::vector &restriction_is_additive_flags,
- const std::vector > &nonzero_components);
+ const std::vector &nonzero_components);
/**
* Since these elements are
diff --git a/deal.II/include/deal.II/fe/fe_system.h b/deal.II/include/deal.II/fe/fe_system.h
index 2b5a38970b..dac5b6033e 100644
--- a/deal.II/include/deal.II/fe/fe_system.h
+++ b/deal.II/include/deal.II/fe/fe_system.h
@@ -1003,7 +1003,7 @@ class FESystem : public FiniteElement
* components of a composed
* finite element.
*/
- static std::vector >
+ static std::vector
compute_nonzero_components (const FiniteElement &fe1,
const unsigned int N1);
@@ -1012,7 +1012,7 @@ class FESystem : public FiniteElement
* components of a composed
* finite element.
*/
- static std::vector >
+ static std::vector
compute_nonzero_components (const FiniteElement &fe1,
const unsigned int N1,
const FiniteElement &fe2,
@@ -1023,7 +1023,7 @@ class FESystem : public FiniteElement
* components of a composed
* finite element.
*/
- static std::vector >
+ static std::vector
compute_nonzero_components (const FiniteElement &fe1,
const unsigned int N1,
const FiniteElement &fe2,
@@ -1036,7 +1036,7 @@ class FESystem : public FiniteElement
* components of a composed
* finite element.
*/
- static std::vector >
+ static std::vector
compute_nonzero_components (const FiniteElement &fe1,
const unsigned int N1,
const FiniteElement &fe2,
@@ -1049,7 +1049,7 @@ class FESystem : public FiniteElement
/**
* With 5 elements.
*/
- static std::vector >
+ static std::vector
compute_nonzero_components (const FiniteElement &fe1,
const unsigned int N1,
const FiniteElement &fe2,
@@ -1069,7 +1069,7 @@ class FESystem : public FiniteElement
* function is called from all
* the above functions.
*/
- static std::vector >
+ static std::vector
compute_nonzero_components (const std::vector*> &fes,
const std::vector &multiplicities);
diff --git a/deal.II/include/deal.II/fe/fe_values.h b/deal.II/include/deal.II/fe/fe_values.h
index fcdf99f5e9..1f43711b02 100644
--- a/deal.II/include/deal.II/fe/fe_values.h
+++ b/deal.II/include/deal.II/fe/fe_values.h
@@ -29,6 +29,7 @@
#include
#include
#include
+#include
#include
#include
#include
@@ -104,158 +105,6 @@ namespace internal
}
-/**
- * A namespace in which we declare "extractors", i.e. classes that when used
- * as subscripts in operator[] expressions on FEValues, FEFaceValues, and
- * FESubfaceValues objects extract certain components of a vector-valued
- * element. The result of applying an extractor to these objects is an object
- * with corresponding type from the namespace FEValuesViews. There are
- * extractors for single scalar components, vector components consisting of
- * dim
elements, and second order symmetric tensors consisting of
- * (dim*dim + dim)/2
components
- *
- * See the description of the @ref vector_valued module for examples how to
- * use the features of this namespace.
- *
- * @ingroup feaccess vector_valued
- */
-namespace FEValuesExtractors
-{
- /**
- * Extractor for a single scalar component
- * of a vector-valued element. The result
- * of applying an object of this type to an
- * FEValues, FEFaceValues or
- * FESubfaceValues object is of type
- * FEValuesViews::Scalar. The concept of
- * extractors is defined in the
- * documentation of the namespace
- * FEValuesExtractors and in the @ref
- * vector_valued module.
- *
- * @ingroup feaccess vector_valued
- */
- struct Scalar
- {
- /**
- * The selected scalar component of the
- * vector.
- */
- unsigned int component;
-
- /**
- * Constructor. Take the selected
- * vector component as argument.
- */
- Scalar (const unsigned int component);
- };
-
-
- /**
- * Extractor for a vector of
- * spacedim
components of a
- * vector-valued element. The value of
- * spacedim
is defined by the
- * FEValues object the extractor is applied
- * to. The result of applying an object of
- * this type to an FEValues, FEFaceValues
- * or FESubfaceValues object is of type
- * FEValuesViews::Vector.
- *
- * The concept of
- * extractors is defined in the
- * documentation of the namespace
- * FEValuesExtractors and in the @ref
- * vector_valued module.
- *
- * Note that in the current context, a
- * vector is meant in the sense physics
- * uses it: it has spacedim
- * components that behave in specific ways
- * under coordinate system
- * transformations. Examples include
- * velocity or displacement fields. This is
- * opposed to how mathematics uses the word
- * "vector" (and how we use this word in
- * other contexts in the library, for
- * example in the Vector class), where it
- * really stands for a collection of
- * numbers. An example of this latter use
- * of the word could be the set of
- * concentrations of chemical species in a
- * flame; however, these are really just a
- * collection of scalar variables, since
- * they do not change if the coordinate
- * system is rotated, unlike the components
- * of a velocity vector, and consequently,
- * this class should not be used for this
- * context.
- *
- * @ingroup feaccess vector_valued
- */
- struct Vector
- {
- /**
- * The first component of the vector
- * view.
- */
- unsigned int first_vector_component;
-
- /**
- * Constructor. Take the first
- * component of the selected vector
- * inside the FEValues object as
- * argument.
- */
- Vector (const unsigned int first_vector_component);
- };
-
-
- /**
- * Extractor for a symmetric tensor of a
- * rank specified by the template
- * argument. For a second order symmetric
- * tensor, this represents a collection of
- * (dim*dim + dim)/2
- * components of a vector-valued
- * element. The value of dim
- * is defined by the FEValues object the
- * extractor is applied to. The result of
- * applying an object of this type to an
- * FEValues, FEFaceValues or
- * FESubfaceValues object is of type
- * FEValuesViews::SymmetricTensor.
- *
- * The concept of
- * extractors is defined in the
- * documentation of the namespace
- * FEValuesExtractors and in the @ref
- * vector_valued module.
- *
- * @ingroup feaccess vector_valued
- *
- * @author Andrew McBride, 2009
- */
- template
- struct SymmetricTensor
- {
- /**
- * The first component of the tensor
- * view.
- */
- unsigned int first_tensor_component;
-
- /**
- * Constructor. Take the first
- * component of the selected tensor
- * inside the FEValues object as
- * argument.
- */
- SymmetricTensor (const unsigned int first_tensor_component);
- };
-}
-
-
/**
@@ -1575,33 +1424,63 @@ class FEValuesData
std::vector > boundary_forms;
/**
- * Indicate the first row which a
- * given shape function occupies
- * in the #shape_values,
- * #shape_gradients and
- * #shape_hessians
- * arrays. If all shape functions
- * are primitive, then this is
- * the identity mapping. If, on
- * the other hand some shape
- * functions have more than one
- * non-zero vector components,
- * then they may occupy more than
- * one row, and this array
- * indicates which is the first
- * one.
- *
- * The questions which particular
- * vector component occupies
- * which row for a given shape
- * function is answered as
- * follows: we allocate one row
- * for each non-zero component as
- * indicated by the
- * FiniteElement::get_nonzero_components()
- * function, and the rows are in
- * ascending order exactly those
- * non-zero components.
+ * When asked for the value (or
+ * gradient, or Hessian) of shape
+ * function i's c-th vector
+ * component, we need to look it
+ * up in the #shape_values,
+ * #shape_gradients and
+ * #shape_hessians arrays. The
+ * question is where in this
+ * array does the data for shape
+ * function i, component c
+ * reside. This is what this
+ * table answers.
+ *
+ * The format of the table is as
+ * follows:
+ * - It has dofs_per_cell times
+ * n_components entries.
+ * - The entry that corresponds to
+ * shape function i, component c
+ * is i * n_components + c
.
+ * - The value stored at this
+ * position indicates the row
+ * in #shape_values and the
+ * other tables where the
+ * corresponding datum is stored
+ * for all the quadrature points.
+ *
+ * In the general, vector-valued
+ * context, the number of
+ * components is larger than one,
+ * but for a given shape
+ * function, not all vector
+ * components may be nonzero
+ * (e.g., if a shape function is
+ * primitive, then exactly one
+ * vector component is non-zero,
+ * while the others are all
+ * zero). For such zero
+ * components, #shape_values and
+ * friends do not have a
+ * row. Consequently, for vector
+ * components for which shape
+ * function i is zero, the entry
+ * in the current table is
+ * numbers::invalid_unsigned_int.
+ *
+ * On the other hand, the table
+ * is guaranteed to have at least
+ * one valid index for each shape
+ * function. In particular, for a
+ * primitive finite element, each
+ * shape function has exactly one
+ * nonzero component and so for
+ * each i, there is exactly one
+ * valid index within the range
+ * [i*n_components,
+ * (i+1)*n_components)
.
*/
std::vector shape_function_to_row_table;
@@ -4612,6 +4491,8 @@ operator[] (const FEValuesExtractors::SymmetricTensor<2> &tensor) const
return fe_values_views_cache.symmetric_second_order_tensors[tensor.first_tensor_component];
}
+
+
template
inline
const double &
@@ -4630,16 +4511,19 @@ FEValuesBase::shape_value (const unsigned int i,
if (fe->is_primitive())
return this->shape_values(i,j);
else
- // otherwise, use the mapping
- // between shape function numbers
- // and rows. note that by the
- // assertions above, we know that
- // this particular shape function
- // is primitive, so there is no
- // question to which vector
- // component the call of this
- // function refers
- return this->shape_values(this->shape_function_to_row_table[i], j);
+ {
+ // otherwise, use the mapping
+ // between shape function
+ // numbers and rows. note that
+ // by the assertions above, we
+ // know that this particular
+ // shape function is primitive,
+ // so we can call
+ // system_to_component_index
+ const unsigned int
+ row = this->shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
+ return this->shape_values(row, j);
+ }
}
@@ -4658,52 +4542,18 @@ FEValuesBase::shape_value_component (const unsigned int i,
Assert (component < fe->n_components(),
ExcIndexRange(component, 0, fe->n_components()));
- // if this particular shape
- // function is primitive, then we
- // can take a short-cut by checking
- // whether the requested component
- // is the only non-zero one (note
- // that calling
- // system_to_component_table only
- // works if the shape function is
- // primitive):
- if (fe->is_primitive(i))
- {
- if (component == fe->system_to_component_index(i).first)
- return this->shape_values(this->shape_function_to_row_table[i],j);
- else
- return 0;
- }
- else
- {
- // no, this shape function is
- // not primitive. then we have
- // to loop over its components
- // to find the corresponding
- // row in the arrays of this
- // object. before that check
- // whether the shape function
- // is non-zero at all within
- // this component:
- if (fe->get_nonzero_components(i)[component] == false)
- return 0;
-
- // count how many non-zero
- // component the shape function
- // has before the one we are
- // looking for, and add this to
- // the offset of the first
- // non-zero component of this
- // shape function in the arrays
- // we index presently:
- const unsigned int
- row = (this->shape_function_to_row_table[i]
- +
- std::count (fe->get_nonzero_components(i).begin(),
- fe->get_nonzero_components(i).begin()+component,
- true));
- return this->shape_values(row, j);
- }
+ // check whether the shape function
+ // is non-zero at all within
+ // this component:
+ if (fe->get_nonzero_components(i)[component] == false)
+ return 0;
+
+ // look up the right row in the
+ // table and take the data from
+ // there
+ const unsigned int
+ row = this->shape_function_to_row_table[i * fe->n_components() + component];
+ return this->shape_values(row, j);
}
@@ -4730,16 +4580,19 @@ FEValuesBase::shape_grad (const unsigned int i,
if (fe->is_primitive())
return this->shape_gradients[i][j];
else
- // otherwise, use the mapping
- // between shape function numbers
- // and rows. note that by the
- // assertions above, we know that
- // this particular shape function
- // is primitive, so there is no
- // question to which vector
- // component the call of this
- // function refers
- return this->shape_gradients[this->shape_function_to_row_table[i]][j];
+ {
+ // otherwise, use the mapping
+ // between shape function
+ // numbers and rows. note that
+ // by the assertions above, we
+ // know that this particular
+ // shape function is primitive,
+ // so we can call
+ // system_to_component_index
+ const unsigned int
+ row = this->shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
+ return this->shape_gradients[row][j];
+ }
}
@@ -4758,52 +4611,18 @@ FEValuesBase::shape_grad_component (const unsigned int i,
Assert (component < fe->n_components(),
ExcIndexRange(component, 0, fe->n_components()));
- // if this particular shape
- // function is primitive, then we
- // can take a short-cut by checking
- // whether the requested component
- // is the only non-zero one (note
- // that calling
- // system_to_component_table only
- // works if the shape function is
- // primitive):
- if (fe->is_primitive(i))
- {
- if (component == fe->system_to_component_index(i).first)
- return this->shape_gradients[this->shape_function_to_row_table[i]][j];
- else
- return Tensor<1,spacedim>();
- }
- else
- {
- // no, this shape function is
- // not primitive. then we have
- // to loop over its components
- // to find the corresponding
- // row in the arrays of this
- // object. before that check
- // whether the shape function
- // is non-zero at all within
- // this component:
- if (fe->get_nonzero_components(i)[component] == false)
- return Tensor<1,spacedim>();
-
- // count how many non-zero
- // component the shape function
- // has before the one we are
- // looking for, and add this to
- // the offset of the first
- // non-zero component of this
- // shape function in the arrays
- // we index presently:
- const unsigned int
- row = (this->shape_function_to_row_table[i]
- +
- std::count (fe->get_nonzero_components(i).begin(),
- fe->get_nonzero_components(i).begin()+component,
- true));
- return this->shape_gradients[row][j];
- }
+ // check whether the shape function
+ // is non-zero at all within
+ // this component:
+ if (fe->get_nonzero_components(i)[component] == false)
+ return Tensor<1,spacedim>();
+
+ // look up the right row in the
+ // table and take the data from
+ // there
+ const unsigned int
+ row = this->shape_function_to_row_table[i * fe->n_components() + component];
+ return this->shape_gradients[row][j];
}
@@ -4830,16 +4649,19 @@ FEValuesBase::shape_hessian (const unsigned int i,
if (fe->is_primitive())
return this->shape_hessians[i][j];
else
- // otherwise, use the mapping
- // between shape function numbers
- // and rows. note that by the
- // assertions above, we know that
- // this particular shape function
- // is primitive, so there is no
- // question to which vector
- // component the call of this
- // function refers
- return this->shape_hessians[this->shape_function_to_row_table[i]][j];
+ {
+ // otherwise, use the mapping
+ // between shape function
+ // numbers and rows. note that
+ // by the assertions above, we
+ // know that this particular
+ // shape function is primitive,
+ // so we can call
+ // system_to_component_index
+ const unsigned int
+ row = this->shape_function_to_row_table[i * fe->n_components() + fe->system_to_component_index(i).first];
+ return this->shape_hessians[row][j];
+ }
}
@@ -4869,52 +4691,18 @@ FEValuesBase::shape_hessian_component (const unsigned int i,
Assert (component < fe->n_components(),
ExcIndexRange(component, 0, fe->n_components()));
- // if this particular shape
- // function is primitive, then we
- // can take a short-cut by checking
- // whether the requested component
- // is the only non-zero one (note
- // that calling
- // system_to_component_table only
- // works if the shape function is
- // primitive):
- if (fe->is_primitive(i))
- {
- if (component == fe->system_to_component_index(i).first)
- return this->shape_hessians[this->shape_function_to_row_table[i]][j];
- else
- return Tensor<2,spacedim>();
- }
- else
- {
- // no, this shape function is
- // not primitive. then we have
- // to loop over its components
- // to find the corresponding
- // row in the arrays of this
- // object. before that check
- // whether the shape function
- // is non-zero at all within
- // this component:
- if (fe->get_nonzero_components(i)[component] == false)
- return Tensor<2,spacedim>();
-
- // count how many non-zero
- // component the shape function
- // has before the one we are
- // looking for, and add this to
- // the offset of the first
- // non-zero component of this
- // shape function in the arrays
- // we index presently:
- const unsigned int
- row = (this->shape_function_to_row_table[i]
- +
- std::count (fe->get_nonzero_components(i).begin(),
- fe->get_nonzero_components(i).begin()+component,
- true));
- return this->shape_hessians[row][j];
- }
+ // check whether the shape function
+ // is non-zero at all within
+ // this component:
+ if (fe->get_nonzero_components(i)[component] == false)
+ return Tensor<2,spacedim>();
+
+ // look up the right row in the
+ // table and take the data from
+ // there
+ const unsigned int
+ row = this->shape_function_to_row_table[i * fe->n_components() + component];
+ return this->shape_hessians[row][j];
}
diff --git a/deal.II/include/deal.II/fe/fe_values_extractors.h b/deal.II/include/deal.II/fe/fe_values_extractors.h
new file mode 100644
index 0000000000..51b3baf675
--- /dev/null
+++ b/deal.II/include/deal.II/fe/fe_values_extractors.h
@@ -0,0 +1,175 @@
+//---------------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
+//
+// 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__fe_values_extractors_h
+#define __deal2__fe_values_extractors_h
+
+
+#include
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/**
+ * A namespace in which we declare "extractors", i.e. classes that when used
+ * as subscripts in operator[] expressions on FEValues, FEFaceValues, and
+ * FESubfaceValues objects extract certain components of a vector-valued
+ * element. The result of applying an extractor to these objects is an object
+ * with corresponding type from the namespace FEValuesViews. There are
+ * extractors for single scalar components, vector components consisting of
+ * dim
elements, and second order symmetric tensors consisting of
+ * (dim*dim + dim)/2
components
+ *
+ * See the description of the @ref vector_valued module for examples how to
+ * use the features of this namespace.
+ *
+ * @ingroup feaccess vector_valued
+ */
+namespace FEValuesExtractors
+{
+ /**
+ * Extractor for a single scalar component
+ * of a vector-valued element. The result
+ * of applying an object of this type to an
+ * FEValues, FEFaceValues or
+ * FESubfaceValues object is of type
+ * FEValuesViews::Scalar. The concept of
+ * extractors is defined in the
+ * documentation of the namespace
+ * FEValuesExtractors and in the @ref
+ * vector_valued module.
+ *
+ * @ingroup feaccess vector_valued
+ */
+ struct Scalar
+ {
+ /**
+ * The selected scalar component of the
+ * vector.
+ */
+ unsigned int component;
+
+ /**
+ * Constructor. Take the selected
+ * vector component as argument.
+ */
+ Scalar (const unsigned int component);
+ };
+
+
+ /**
+ * Extractor for a vector of
+ * spacedim
components of a
+ * vector-valued element. The value of
+ * spacedim
is defined by the
+ * FEValues object the extractor is applied
+ * to. The result of applying an object of
+ * this type to an FEValues, FEFaceValues
+ * or FESubfaceValues object is of type
+ * FEValuesViews::Vector.
+ *
+ * The concept of
+ * extractors is defined in the
+ * documentation of the namespace
+ * FEValuesExtractors and in the @ref
+ * vector_valued module.
+ *
+ * Note that in the current context, a
+ * vector is meant in the sense physics
+ * uses it: it has spacedim
+ * components that behave in specific ways
+ * under coordinate system
+ * transformations. Examples include
+ * velocity or displacement fields. This is
+ * opposed to how mathematics uses the word
+ * "vector" (and how we use this word in
+ * other contexts in the library, for
+ * example in the Vector class), where it
+ * really stands for a collection of
+ * numbers. An example of this latter use
+ * of the word could be the set of
+ * concentrations of chemical species in a
+ * flame; however, these are really just a
+ * collection of scalar variables, since
+ * they do not change if the coordinate
+ * system is rotated, unlike the components
+ * of a velocity vector, and consequently,
+ * this class should not be used for this
+ * context.
+ *
+ * @ingroup feaccess vector_valued
+ */
+ struct Vector
+ {
+ /**
+ * The first component of the vector
+ * view.
+ */
+ unsigned int first_vector_component;
+
+ /**
+ * Constructor. Take the first
+ * component of the selected vector
+ * inside the FEValues object as
+ * argument.
+ */
+ Vector (const unsigned int first_vector_component);
+ };
+
+
+ /**
+ * Extractor for a symmetric tensor of a
+ * rank specified by the template
+ * argument. For a second order symmetric
+ * tensor, this represents a collection of
+ * (dim*dim + dim)/2
+ * components of a vector-valued
+ * element. The value of dim
+ * is defined by the FEValues object the
+ * extractor is applied to. The result of
+ * applying an object of this type to an
+ * FEValues, FEFaceValues or
+ * FESubfaceValues object is of type
+ * FEValuesViews::SymmetricTensor.
+ *
+ * The concept of
+ * extractors is defined in the
+ * documentation of the namespace
+ * FEValuesExtractors and in the @ref
+ * vector_valued module.
+ *
+ * @ingroup feaccess vector_valued
+ *
+ * @author Andrew McBride, 2009
+ */
+ template
+ struct SymmetricTensor
+ {
+ /**
+ * The first component of the tensor
+ * view.
+ */
+ unsigned int first_tensor_component;
+
+ /**
+ * Constructor. Take the first
+ * component of the selected tensor
+ * inside the FEValues object as
+ * argument.
+ */
+ SymmetricTensor (const unsigned int first_tensor_component);
+ };
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/include/deal.II/hp/dof_handler.h b/deal.II/include/deal.II/hp/dof_handler.h
index a493349218..ecec9505a4 100644
--- a/deal.II/include/deal.II/hp/dof_handler.h
+++ b/deal.II/include/deal.II/hp/dof_handler.h
@@ -1285,7 +1285,10 @@ namespace hp
* not. It is used by the
* refinement listeners as a
* persistent buffer during the
- * refinement.
+ * refinement, i.e. from between
+ * when pre_refinement_action is
+ * called and when post_refinement_action
+ * runs.
*/
std::vector *> has_children;
diff --git a/deal.II/include/deal.II/hp/fe_collection.h b/deal.II/include/deal.II/hp/fe_collection.h
index ff0f550f77..a158a1616e 100644
--- a/deal.II/include/deal.II/hp/fe_collection.h
+++ b/deal.II/include/deal.II/hp/fe_collection.h
@@ -15,6 +15,8 @@
#include
#include
#include
+#include
+#include
DEAL_II_NAMESPACE_OPEN
@@ -229,6 +231,197 @@ namespace hp
*/
bool hp_constraints_are_implemented () const;
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * one component is true that corresponds to the given
+ * argument.
+ *
+ * @note This function is the equivalent of
+ * FiniteElement::component_mask() with the same arguments.
+ * It verifies that it gets the same result from every one
+ * of the elements that are stored in this FECollection. If
+ * this is not the case, it throws an exception.
+ *
+ * @param scalar An object that represents a single scalar
+ * vector component of this finite element.
+ * @return A component mask that is false in all components
+ * except for the one that corresponds to the argument.
+ */
+ ComponentMask
+ component_mask (const FEValuesExtractors::Scalar &scalar) const;
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * dim
components are true that correspond to the given
+ * argument.
+ *
+ * @note This function is the equivalent of
+ * FiniteElement::component_mask() with the same arguments.
+ * It verifies that it gets the same result from every one
+ * of the elements that are stored in this FECollection. If
+ * this is not the case, it throws an exception.
+ *
+ * @param vector An object that represents dim
+ * vector components of this finite element.
+ * @return A component mask that is false in all components
+ * except for the ones that corresponds to the argument.
+ */
+ ComponentMask
+ component_mask (const FEValuesExtractors::Vector &vector) const;
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * dim*(dim+1)/2
components are true that
+ * correspond to the given argument.
+ *
+ * @note This function is the equivalent of
+ * FiniteElement::component_mask() with the same arguments.
+ * It verifies that it gets the same result from every one
+ * of the elements that are stored in this FECollection. If
+ * this is not the case, it throws an exception.
+ *
+ * @param sym_tensor An object that represents dim*(dim+1)/2
+ * components of this finite element that are jointly to be
+ * interpreted as forming a symmetric tensor.
+ * @return A component mask that is false in all components
+ * except for the ones that corresponds to the argument.
+ */
+ ComponentMask
+ component_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;
+
+ /**
+ * Given a block mask (see @ref GlossBlockMask "this glossary entry"),
+ * produce a component mask (see @ref GlossComponentMask "this glossary entry")
+ * that represents the components that correspond to the blocks selected in
+ * the input argument. This is essentially a conversion operator from
+ * BlockMask to ComponentMask.
+ *
+ * @note This function is the equivalent of
+ * FiniteElement::component_mask() with the same arguments.
+ * It verifies that it gets the same result from every one
+ * of the elements that are stored in this FECollection. If
+ * this is not the case, it throws an exception.
+ *
+ * @param block_mask The mask that selects individual blocks of the finite
+ * element
+ * @return A mask that selects those components corresponding to the selected
+ * blocks of the input argument.
+ */
+ ComponentMask
+ component_mask (const BlockMask &block_mask) const;
+
+ /**
+ * Return a block mask with as many elements as this
+ * object has blocks and of which exactly the
+ * one component is true that corresponds to the given
+ * argument. See @ref GlossBlockMask "the glossary"
+ * for more information.
+ *
+ * @note This function will only succeed if the scalar referenced
+ * by the argument encompasses a complete block. In other words,
+ * if, for example, you pass an extractor for the single
+ * $x$ velocity and this object represents an FE_RaviartThomas
+ * object, then the single scalar object you selected is part
+ * of a larger block and consequently there is no block mask that
+ * would represent it. The function will then produce an exception.
+ *
+ * @note This function is the equivalent of
+ * FiniteElement::component_mask() with the same arguments.
+ * It verifies that it gets the same result from every one
+ * of the elements that are stored in this FECollection. If
+ * this is not the case, it throws an exception.
+ *
+ * @param scalar An object that represents a single scalar
+ * vector component of this finite element.
+ * @return A component mask that is false in all components
+ * except for the one that corresponds to the argument.
+ */
+ BlockMask
+ block_mask (const FEValuesExtractors::Scalar &scalar) const;
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * dim
components are true that correspond to the given
+ * argument. See @ref GlossBlockMask "the glossary"
+ * for more information.
+ *
+ * @note This function is the equivalent of
+ * FiniteElement::component_mask() with the same arguments.
+ * It verifies that it gets the same result from every one
+ * of the elements that are stored in this FECollection. If
+ * this is not the case, it throws an exception.
+ *
+ * @note The same caveat applies as to the version of the function above:
+ * The extractor object passed as argument must be so that it corresponds
+ * to full blocks and does not split blocks of this element.
+ *
+ * @param vector An object that represents dim
+ * vector components of this finite element.
+ * @return A component mask that is false in all components
+ * except for the ones that corresponds to the argument.
+ */
+ BlockMask
+ block_mask (const FEValuesExtractors::Vector &vector) const;
+
+ /**
+ * Return a component mask with as many elements as this
+ * object has vector components and of which exactly the
+ * dim*(dim+1)/2
components are true that
+ * correspond to the given argument. See @ref GlossBlockMask "the glossary"
+ * for more information.
+ *
+ * @note The same caveat applies as to the version of the function above:
+ * The extractor object passed as argument must be so that it corresponds
+ * to full blocks and does not split blocks of this element.
+ *
+ * @note This function is the equivalent of
+ * FiniteElement::component_mask() with the same arguments.
+ * It verifies that it gets the same result from every one
+ * of the elements that are stored in this FECollection. If
+ * this is not the case, it throws an exception.
+ *
+ * @param sym_tensor An object that represents dim*(dim+1)/2
+ * components of this finite element that are jointly to be
+ * interpreted as forming a symmetric tensor.
+ * @return A component mask that is false in all components
+ * except for the ones that corresponds to the argument.
+ */
+ BlockMask
+ block_mask (const FEValuesExtractors::SymmetricTensor<2> &sym_tensor) const;
+
+ /**
+ * Given a component mask (see @ref GlossComponentMask "this glossary entry"),
+ * produce a block mask (see @ref GlossBlockMask "this glossary entry")
+ * that represents the blocks that correspond to the components selected in
+ * the input argument. This is essentially a conversion operator from
+ * ComponentMask to BlockMask.
+ *
+ * @note This function will only succeed if the components referenced
+ * by the argument encompasses complete blocks. In other words,
+ * if, for example, you pass an component mask for the single
+ * $x$ velocity and this object represents an FE_RaviartThomas
+ * object, then the single component you selected is part
+ * of a larger block and consequently there is no block mask that
+ * would represent it. The function will then produce an exception.
+ *
+ * @note This function is the equivalent of
+ * FiniteElement::component_mask() with the same arguments.
+ * It verifies that it gets the same result from every one
+ * of the elements that are stored in this FECollection. If
+ * this is not the case, it throws an exception.
+ *
+ * @param component_mask The mask that selects individual components of the finite
+ * element
+ * @return A mask that selects those blocks corresponding to the selected
+ * blocks of the input argument.
+ */
+ BlockMask
+ block_mask (const ComponentMask &component_mask) const;
+
/**
* Exception
diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
index 970ee190aa..715527800e 100644
--- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h
+++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h
@@ -1737,7 +1737,7 @@ namespace internals
// similar as before, now with shortcut for
- // deal.II sparse matrices. this lets use
+ // deal.II sparse matrices. this lets us
// avoid using extra arrays, and does all the
// operations just in place, i.e., in the
// respective matrix row
diff --git a/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h b/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
index a5f929451c..e011893a8e 100644
--- a/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
+++ b/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
@@ -59,7 +59,7 @@ class MGConstrainedDoFs : public Subscriptor
template
void initialize(const MGDoFHandler& dof,
const typename FunctionMap::type& function_map,
- const std::vector& component_mask = std::vector());
+ const ComponentMask &component_mask = ComponentMask());
/**
* Reset the data structures.
@@ -209,7 +209,7 @@ void
MGConstrainedDoFs::initialize(
const MGDoFHandler& dof,
const typename FunctionMap::type& function_map,
- const std::vector& component_mask)
+ const ComponentMask &component_mask)
{
const unsigned int nlevels = dof.get_tria().n_levels();
boundary_indices.resize(nlevels);
diff --git a/deal.II/include/deal.II/multigrid/mg_tools.h b/deal.II/include/deal.II/multigrid/mg_tools.h
index b1695498e0..3dcf1b751e 100644
--- a/deal.II/include/deal.II/multigrid/mg_tools.h
+++ b/deal.II/include/deal.II/multigrid/mg_tools.h
@@ -234,7 +234,7 @@ namespace MGTools
make_boundary_list (const MGDoFHandler &mg_dof,
const typename FunctionMap::type &function_map,
std::vector > &boundary_indices,
- const std::vector &component_mask = std::vector());
+ const ComponentMask &component_mask = ComponentMask());
/**
* The same function as above, but return
@@ -246,7 +246,7 @@ namespace MGTools
make_boundary_list (const MGDoFHandler &mg_dof,
const typename FunctionMap::type &function_map,
std::vector &boundary_indices,
- const std::vector &component_mask = std::vector());
+ const ComponentMask &component_mask = ComponentMask());
/**
* Maybe no longer needed.
@@ -257,7 +257,7 @@ namespace MGTools
apply_boundary_values (const std::set &boundary_dofs,
SparseMatrix& matrix,
const bool preserve_symmetry,
- const bool ignore_zeros = false);
+ const bool ignore_zeros = false);
template
void
diff --git a/deal.II/include/deal.II/multigrid/mg_transfer_block.h b/deal.II/include/deal.II/multigrid/mg_transfer_block.h
index bcb3016897..d28677074c 100644
--- a/deal.II/include/deal.II/multigrid/mg_transfer_block.h
+++ b/deal.II/include/deal.II/multigrid/mg_transfer_block.h
@@ -97,6 +97,7 @@ class MGTransferBlockBase
* on the blocks having a
* true entry here.
*/
+//TODO: rename this to block_mask, in the same way as has already been done in MGTransferComponent, and give it type BlockMask
std::vector selected;
/**
@@ -116,6 +117,7 @@ class MGTransferBlockBase
* entries may be illegal
* unsigned integers.
*/
+//TODO: rename this to mg_block_mask, in the same way as has already been done in MGTransferComponent, and give it type BlockMask
std::vector mg_block;
/**
diff --git a/deal.II/include/deal.II/multigrid/mg_transfer_component.h b/deal.II/include/deal.II/multigrid/mg_transfer_component.h
index c6308d6565..73b7a98d40 100644
--- a/deal.II/include/deal.II/multigrid/mg_transfer_component.h
+++ b/deal.II/include/deal.II/multigrid/mg_transfer_component.h
@@ -92,7 +92,7 @@ class MGTransferComponentBase
* this refers to the
* renumbered components.
*/
- std::vector selected;
+ ComponentMask component_mask;
/**
* Flag of selected components.
@@ -105,7 +105,7 @@ class MGTransferComponentBase
* this refers to the
* renumbered components.
*/
- std::vector mg_selected;
+ ComponentMask mg_component_mask;
/**
* Target component of the
diff --git a/deal.II/include/deal.II/numerics/error_estimator.h b/deal.II/include/deal.II/numerics/error_estimator.h
index e856482be9..076841240c 100644
--- a/deal.II/include/deal.II/numerics/error_estimator.h
+++ b/deal.II/include/deal.II/numerics/error_estimator.h
@@ -345,7 +345,7 @@ class KellyErrorEstimator
const typename FunctionMap::type &neumann_bc,
const InputVector &solution,
Vector &error,
- const std::vector &component_mask = std::vector(),
+ const ComponentMask &component_mask = ComponentMask(),
const Function *coefficients = 0,
const unsigned int n_threads = numbers::invalid_unsigned_int,
const types::subdomain_id subdomain_id = types::invalid_subdomain_id,
@@ -362,7 +362,7 @@ class KellyErrorEstimator
const typename FunctionMap::type &neumann_bc,
const InputVector &solution,
Vector &error,
- const std::vector &component_mask = std::vector(),
+ const ComponentMask &component_mask = ComponentMask(),
const Function *coefficients = 0,
const unsigned int n_threads = multithread_info.n_default_threads,
const types::subdomain_id subdomain_id = types::invalid_subdomain_id,
@@ -402,7 +402,7 @@ class KellyErrorEstimator
const typename FunctionMap::type &neumann_bc,
const std::vector &solutions,
std::vector*> &errors,
- const std::vector &component_mask = std::vector(),
+ const ComponentMask &component_mask = ComponentMask(),
const Function *coefficients = 0,
const unsigned int n_threads = multithread_info.n_default_threads,
const types::subdomain_id subdomain_id = types::invalid_subdomain_id,
@@ -419,7 +419,7 @@ class KellyErrorEstimator
const typename FunctionMap::type &neumann_bc,
const std::vector &solutions,
std::vector*> &errors,
- const std::vector &component_mask = std::vector(),
+ const ComponentMask &component_mask = ComponentMask(),
const Function *coefficients = 0,
const unsigned int n_threads = multithread_info.n_default_threads,
const types::subdomain_id subdomain_id = types::invalid_subdomain_id,
@@ -439,7 +439,7 @@ class KellyErrorEstimator
const typename FunctionMap::type &neumann_bc,
const InputVector &solution,
Vector &error,
- const std::vector &component_mask = std::vector(),
+ const ComponentMask &component_mask = ComponentMask(),
const Function *coefficients = 0,
const unsigned int n_threads = multithread_info.n_default_threads,
const types::subdomain_id subdomain_id = types::invalid_subdomain_id,
@@ -458,7 +458,7 @@ class KellyErrorEstimator
const typename FunctionMap::type &neumann_bc,
const InputVector &solution,
Vector &error,
- const std::vector &component_mask = std::vector(),
+ const ComponentMask &component_mask = ComponentMask(),
const Function *coefficients = 0,
const unsigned int n_threads = multithread_info.n_default_threads,
const types::subdomain_id subdomain_id = types::invalid_subdomain_id,
@@ -478,7 +478,7 @@ class KellyErrorEstimator
const typename FunctionMap::type &neumann_bc,
const std::vector