/**
* Construct an object of this
- * type. You have to set some
+ * type. You have to set some
* member variables, for example
* some matrices, explicitly
* after calling this base class'
- * constructor.
+ * constructor. For this see the
+ * existing finite element
+ * classes. For the second and
+ * third parameter of this
+ * constructor, see the document
+ * of the respective member
+ * variables.
*/
FiniteElementBase (const FiniteElementData<dim> &fe_data,
const std::vector<bool> &restriction_is_additive_flags,
component_to_base (unsigned int component) const;
/**
- * Access the @p{restriction_is_additive_flag}
- * field. See there for more information on
- * its contents.
+ * Access the
+ * @p{restriction_is_additive_flag}
+ * field. See there for more
+ * information on its contents.
+ *
+ * The index must be between zero
+ * and the number of shape
+ * functions of this element.
*/
- bool restriction_is_additive (const unsigned int component) const;
+ bool restriction_is_additive (const unsigned int index) const;
/**
* Return the support points of
*
* This array has valid values
* also in the case of
- * vector-value
+ * vector-valued
* (i.e. non-primitive) shape
* functions, in contrast to the
* @p{system_to_component_table}.
*/
- std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> > system_to_base_table;
+ std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
+ system_to_base_table;
/**
* Likewise for the indices on
* faces.
*/
- std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> > face_system_to_base_table;
+ std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> >
+ face_system_to_base_table;
/**
* Map between component and
* build the complete matrix. The
* flag should be @p{true}.
*
- * There is one flag per
- * component in vector valued
- * elements.
+ * For examples of use of these
+ * flags, see the places in the
+ * library where it is queried.
+ *
+ * There is one flag per shape
+ * function, indicating whether
+ * it belongs to the class of
+ * shape functions that are
+ * additive in the restriction or
+ * not.
+ *
+ * Note that in previous versions
+ * of the library, there was one
+ * flag per vector component of
+ * the element. This is based on
+ * the fact that all the shape
+ * functions that belong to the
+ * same vector component must
+ * necessarily behave in the same
+ * way, to make things
+ * reasonable. However, the
+ * problem is that it is
+ * sometimes impossible to query
+ * this flag in the vector-valued
+ * case: this used to be done
+ * with the
+ * @p{system_to_component_index}
+ * function that returns which
+ * vector component a shape
+ * function is associated
+ * with. The point is that since
+ * we now support shape functions
+ * that are associated with more
+ * than one vector component (for
+ * example the shape functions of
+ * Raviart-Thomas, or Nedelec
+ * elements), that function can
+ * no more be used, so it can be
+ * difficult to find out which
+ * for vector component we would
+ * like to query the
+ * restriction-is-additive flags.
*/
const std::vector<bool> restriction_is_additive_flags;
template <int dim>
inline
bool
-FiniteElementBase<dim>::restriction_is_additive (const unsigned int component) const
+FiniteElementBase<dim>::restriction_is_additive (const unsigned int index) const
{
- Assert(component<this->n_components(),
- ExcIndexRange(component, 0, this->n_components()));
- return restriction_is_additive_flags[component];
+ Assert(index < this->dofs_per_cell,
+ ExcIndexRange(index, 0, this->dofs_per_cell));
+ return restriction_is_additive_flags[index];
}
/**
- * Helper function used in the constructor:
- * takes a @p{FiniteElement} object
- * and returns an boolean vector including
- * the @p{restriction_is_additive_flags} of
- * the mixed element consisting of @p{N}
- * elements of the sub-element @p{fe}.
+ * Helper function used in the
+ * constructor: takes a
+ * @p{FiniteElement} object and
+ * returns an boolean vector
+ * including the
+ * @p{restriction_is_additive_flags}
+ * of the mixed element
+ * consisting of @p{N} elements
+ * of the sub-element @p{fe}.
*/
static std::vector<bool>
compute_restriction_is_additive_flags (const FiniteElement<dim> &fe,
const FiniteElement<dim> &fe3,
const unsigned int N3);
+ /**
+ * Compute the named flags for a
+ * list of finite elements with
+ * multiplicities given in the
+ * second argument. This function
+ * is called from all the above
+ * functions.
+ */
+ static std::vector<bool>
+ compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim>*> &fes,
+ const std::vector<unsigned int> &multiplicities);
+
+
/**
* Compute the non-zero vector
* components of a composed
template <int dim>
-FiniteElementBase<dim>::FiniteElementBase (const FiniteElementData<dim> &fe_data,
- const std::vector<bool> &restriction_is_additive_flags,
- const std::vector<std::vector<bool> > &nonzero_components)
+FiniteElementBase<dim>::
+FiniteElementBase (const FiniteElementData<dim> &fe_data,
+ const std::vector<bool> &restriction_is_additive_flags,
+ const std::vector<std::vector<bool> > &nonzero_components)
:
FiniteElementData<dim> (fe_data),
- system_to_component_table(this->dofs_per_cell),
- face_system_to_component_table(this->dofs_per_face),
- system_to_base_table(this->dofs_per_cell),
- face_system_to_base_table(this->dofs_per_face),
- component_to_system_table(this->components,
- std::vector<unsigned>(this->dofs_per_cell)),
- face_component_to_system_table(this->components,
- std::vector<unsigned>(this->dofs_per_face)),
- component_to_base_table (this->components,
- std::make_pair(0U, 0U)),
- restriction_is_additive_flags(restriction_is_additive_flags),
- nonzero_components (nonzero_components),
- n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)),
- cached_primitivity (std::find_if (n_nonzero_components_table.begin(),
- n_nonzero_components_table.end(),
- std::bind2nd(std::not_equal_to<unsigned int>(),
- 1U))
- ==
- n_nonzero_components_table.end())
+ system_to_component_table (this->dofs_per_cell),
+ face_system_to_component_table(this->dofs_per_face),
+ system_to_base_table(this->dofs_per_cell),
+ face_system_to_base_table(this->dofs_per_face),
+ component_to_system_table(this->components,
+ std::vector<unsigned>(this->dofs_per_cell)),
+ face_component_to_system_table(this->components,
+ std::vector<unsigned>(this->dofs_per_face)),
+ component_to_base_table (this->components,
+ std::make_pair(0U, 0U)),
+ restriction_is_additive_flags(restriction_is_additive_flags),
+ nonzero_components (nonzero_components),
+ n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)),
+ cached_primitivity (std::find_if (n_nonzero_components_table.begin(),
+ n_nonzero_components_table.end(),
+ std::bind2nd(std::not_equal_to<unsigned int>(),
+ 1U))
+ ==
+ n_nonzero_components_table.end())
{
- Assert (restriction_is_additive_flags.size()==fe_data.components,
+ Assert (restriction_is_additive_flags.size() == this->dofs_per_cell,
ExcDimensionMismatch(restriction_is_additive_flags.size(),
- fe_data.components));
+ this->dofs_per_cell));
Assert (nonzero_components.size() == this->dofs_per_cell,
ExcInternalError());
for (unsigned int i=0; i<nonzero_components.size(); ++i)
FE_DGP<dim>::FE_DGP (const unsigned int degree)
:
FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
- std::vector<bool>(1,true),
+ std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,true),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
std::vector<bool>(1,true))),
- degree(degree),
- polynomial_space (Legendre<double>::generate_complete_basis(degree))
+ degree(degree),
+ polynomial_space (Legendre<double>::generate_complete_basis(degree))
{
// if defined, copy over matrices
// from precomputed arrays
else
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
this->prolongation[i].reinit(0,0);
-
+
+ // restriction can be defined
+ // through projection for
+ // discontinuous elements, but is
+ // presently not implemented for DGP
+ // elements.
+ //
+ // if it were, then the following
+ // snippet would be the right code
// if ((degree < Matrices::n_projection_matrices) &&
// (Matrices::projection_matrices[degree] != 0))
// {
// else
// // matrix undefined, set size to zero
// for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
-// restriction[i].reinit(0);
+// restriction[i].reinit(0, 0);
+ // since not implemented, set to
+ // "empty"
+ for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
+ restriction[i].reinit(0, 0);
+
+ // note further, that these
+ // elements have neither support
+ // nor face-support points, so
+ // leave these fields empty
};
FE_DGQ<dim>::FE_DGQ (const unsigned int degree)
:
FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
- std::vector<bool>(1,true),
+ std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+ true),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
std::vector<bool>(1,true))),
degree(degree),
// from precomputed arrays and
// generate all other matrices by
// permutations
+ //
+ // (note: the matrix is defined if
+ // something was entered into the
+ // respective table, and what was
+ // entered is not a NULL pointer --
+ // this would allow for "holes")
if ((degree < Matrices::n_embedding_matrices) &&
(Matrices::embedding[degree] != 0))
{
FE_Q<dim>::FE_Q (const unsigned int degree)
:
FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
- std::vector<bool> (1,false),
+ std::vector<bool> (FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+ false),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
std::vector<bool>(1,true))),
degree(degree),
FESystem<dim>::compute_restriction_is_additive_flags (const FiniteElement<dim> &fe,
const unsigned int n_elements)
{
- std::vector<bool> tmp;
- for (unsigned int i=0; i<n_elements; ++i)
- for (unsigned int component=0; component<fe.n_components(); ++component)
- tmp.push_back (fe.restriction_is_additive (component));
- return tmp;
+ std::vector<const FiniteElement<dim>*> fe_list;
+ std::vector<unsigned int> multiplicities;
+
+ fe_list.push_back (&fe);
+ multiplicities.push_back (n_elements);
+
+ return compute_restriction_is_additive_flags (fe_list, multiplicities);
};
const FiniteElement<dim> &fe2,
const unsigned int N2)
{
- std::vector<bool> tmp;
- for (unsigned int i=0; i<N1; ++i)
- for (unsigned int component=0; component<fe1.n_components(); ++component)
- tmp.push_back (fe1.restriction_is_additive (component));
- for (unsigned int i=0; i<N2; ++i)
- for (unsigned int component=0; component<fe2.n_components(); ++component)
- tmp.push_back (fe2.restriction_is_additive (component));
- return tmp;
+ std::vector<const FiniteElement<dim>*> fe_list;
+ std::vector<unsigned int> multiplicities;
+
+ fe_list.push_back (&fe1);
+ multiplicities.push_back (N1);
+
+ fe_list.push_back (&fe2);
+ multiplicities.push_back (N2);
+
+ return compute_restriction_is_additive_flags (fe_list, multiplicities);
};
const FiniteElement<dim> &fe3,
const unsigned int N3)
{
- std::vector<bool> tmp;
- for (unsigned int i=0; i<N1; ++i)
- for (unsigned int component=0; component<fe1.n_components(); ++component)
- tmp.push_back (fe1.restriction_is_additive (component));
- for (unsigned int i=0; i<N2; ++i)
- for (unsigned int component=0; component<fe2.n_components(); ++component)
- tmp.push_back (fe2.restriction_is_additive (component));
- for (unsigned int i=0; i<N3; ++i)
- for (unsigned int component=0; component<fe3.n_components(); ++component)
- tmp.push_back (fe3.restriction_is_additive (component));
- return tmp;
+ std::vector<const FiniteElement<dim>*> fe_list;
+ std::vector<unsigned int> multiplicities;
+
+ fe_list.push_back (&fe1);
+ multiplicities.push_back (N1);
+
+ fe_list.push_back (&fe2);
+ multiplicities.push_back (N2);
+
+ fe_list.push_back (&fe3);
+ multiplicities.push_back (N3);
+
+ return compute_restriction_is_additive_flags (fe_list, multiplicities);
+};
+
+
+
+template <int dim>
+std::vector<bool>
+FESystem<dim>::
+compute_restriction_is_additive_flags (const std::vector<const FiniteElement<dim>*> &fes,
+ const std::vector<unsigned int> &multiplicities)
+{
+ Assert (fes.size() == multiplicities.size(), ExcInternalError());
+
+ // first count the number of dofs
+ // and components that will emerge
+ // from the given FEs
+ unsigned int n_shape_functions = 0;
+ for (unsigned int i=0; i<fes.size(); ++i)
+ n_shape_functions += fes[i]->dofs_per_cell * multiplicities[i];
+
+ // generate the array that will
+ // hold the output
+ std::vector<bool> retval (n_shape_functions, false);
+
+ // finally go through all the shape
+ // functions of the base elements,
+ // and copy their flags. this
+ // somehow copies the code in
+ // build_cell_table, which is not
+ // nice as it uses too much
+ // implicit knowledge about the
+ // layout of the individual bases
+ // in the composed FE, but there
+ // seems no way around...
+ //
+ // for each shape function, copy
+ // the flags from the base element
+ // to this one, taking into account
+ // multiplicities, and other
+ // complications
+ unsigned int total_index = 0;
+ for (unsigned int vertex_number=0;
+ vertex_number<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex_number)
+ {
+ for(unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ for (unsigned int local_index = 0;
+ local_index < fes[base]->dofs_per_vertex;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fes[base]->dofs_per_vertex*vertex_number +
+ local_index);
+
+ Assert (index_in_base < fes[base]->dofs_per_cell,
+ ExcInternalError());
+ retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+ }
+ }
+
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ for (unsigned int local_index = 0;
+ local_index < fes[base]->dofs_per_line;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fes[base]->dofs_per_line*line_number +
+ local_index +
+ fes[base]->first_line_index);
+
+ Assert (index_in_base < fes[base]->dofs_per_cell,
+ ExcInternalError());
+ retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+ }
+ }
+
+ // 3. Quads
+ if (GeometryInfo<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ for (unsigned int local_index = 0;
+ local_index < fes[base]->dofs_per_quad;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fes[base]->dofs_per_quad*quad_number +
+ local_index +
+ fes[base]->first_quad_index);
+
+ Assert (index_in_base < fes[base]->dofs_per_cell,
+ ExcInternalError());
+ retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+ }
+ }
+
+ // 4. Hexes
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ for(unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base]; ++m)
+ for (unsigned int local_index = 0;
+ local_index < fes[base]->dofs_per_hex;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (fes[base]->dofs_per_hex*hex_number +
+ local_index +
+ fes[base]->first_hex_index);
+
+ Assert (index_in_base < fes[base]->dofs_per_cell,
+ ExcInternalError());
+ retval[total_index] = fes[base]->restriction_is_additive(index_in_base);
+ }
+ }
+
+ Assert (total_index == n_shape_functions, ExcInternalError());
+
+ return retval;
};
}
}
- Assert (total_index == retval.size(), ExcInternalError());
+ Assert (total_index == n_shape_functions, ExcInternalError());
return retval;
};
+
template <int dim>
const FiniteElement<dim> &
FESystem<dim>::base_element (const unsigned int index) const