// freedom on sub-objects which are
// not allocated for this
// non-active thing
- Assert (this->active() ||
+ Assert (!this->has_children() ||
(this->dof_handler->get_fe().dofs_per_cell ==
2*this->dof_handler->get_fe().dofs_per_vertex),
typename DoFAccessor<dim>::ExcNotActive());
// freedom on sub-objects which are
// not allocated for this
// non-active thing
- Assert (this->active() ||
+ Assert (!this->has_children() ||
(this->dof_handler->get_fe().dofs_per_cell ==
4*this->dof_handler->get_fe().dofs_per_vertex),
typename DoFAccessor<dim>::ExcNotActive());
// freedom on sub-objects which are
// not allocated for this
// non-active thing
- Assert (this->active() ||
+ Assert (!this->has_children() ||
(this->dof_handler->get_fe().dofs_per_cell ==
8*this->dof_handler->get_fe().dofs_per_vertex),
typename DoFAccessor<dim>::ExcNotActive());
-
//---------------------------- fe.h ---------------------------
// $Id$
// Version: $Name$
* Constructor
*/
FiniteElement (const FiniteElementData<dim> &fe_data,
- const std::vector<bool> &restriction_is_additive_flags);
+ const std::vector<bool> &restriction_is_additive_flags,
+ const std::vector<std::vector<bool> > &nonzero_components);
/**
* Virtual destructor. Makes sure
protected:
/**
- * Compute flags for initial
- * update only.
- * @ref{FEValuesBase}
+ * Given a set of flags indicating
+ * what quantities are requested
+ * from a @p{FEValues} object,
+ * return which of these can be
+ * precomputed once and for
+ * all. Often, the values of
+ * shape function at quadrature
+ * points can be precomputed, for
+ * example, in which case the
+ * return value of this function
+ * would be the logical and of
+ * the input @p{flags} and
+ * @p{update_values}.
*/
virtual UpdateFlags update_once (const UpdateFlags flags) const = 0;
/**
- * Compute flags for update on
- * each cell.
- * @ref{FEValuesBase}
+ * This is the opposite to the
+ * above function: given a set of
+ * flags indicating what we want
+ * to know, return which of these
+ * need to be computed each time
+ * we visit a new cell.
+ *
+ * If for the computation of one
+ * quantity something else is
+ * also required (for example, we
+ * often need the covariant
+ * transformation when gradients
+ * need to be computed), include
+ * this in the result as well.
*/
virtual UpdateFlags update_each (const UpdateFlags flags) const = 0;
template<int dim> class FESystem;
+
+
/**
* Dimension independent data for finite elements. See the derived
* class @ref{FiniteElementBase} class for information on its use. All
const unsigned int components;
/**
- * Default constructor. Constructs
- * an element
- * which is not so useful. Checking
- * @p{dofs_per_cell} is therefore a good way to
- * check if something went wrong.
+ * Default
+ * constructor. Constructs an
+ * element which is not so
+ * useful. Checking
+ * @p{dofs_per_cell} is therefore
+ * a good way to check if
+ * something went wrong.
*/
FiniteElementData ();
unsigned int n_components () const;
/**
- * Comparison operator. It is not clear to
- * me (WB) why we have to declare and implement
- * this one explicitly.
+ * Comparison operator.
*/
-//TODO:[WB] (compiler) remove operator and let the compiler generate it as soon as it is willing to do so
bool operator == (const FiniteElementData<dim> &) const;
/**
* introduced from the two sides are unique; it is able to handle the fact
* that the constraints for some of the dofs are entered more than once.
*
- * @author Wolfgang Bangerth, 1998, Ralf Hartmann, Guido Kanschat, 2001
+ * @author Wolfgang Bangerth, 1998, 2002, Ralf Hartmann, Guido Kanschat, 2001
*/
template <int dim>
class FiniteElementBase : public Subscriptor,
{
public:
/**
- * Basis class for internal data.
+ * Base class for internal data.
* Adds data for second derivatives to
* @ref{Mapping::InternalDataBase}
*
+ * For information about the
+ * general purpose of this class,
+ * see the documentation of the
+ * base class.
+ *
* @author Guido Kanschat, 2001
*/
class InternalDataBase : public Mapping<dim>::InternalDataBase
};
/**
- * Construct an object of this type.
- * You have to set the
- * matrices explicitly after calling
- * this base class' constructor.
+ * Construct an object of this
+ * type. You have to set some
+ * member variables, for example
+ * some matrices, explicitly
+ * after calling this base class'
+ * constructor.
*/
FiniteElementBase (const FiniteElementData<dim> &fe_data,
- const std::vector<bool> &restriction_is_additive_flags);
+ const std::vector<bool> &restriction_is_additive_flags,
+ const std::vector<std::vector<bool> > &nonzero_components);
/**
* Return the value of the
* @p{i}th shape function at the
- * point @p{p}. @p{p} is a point
- * on the reference element.
+ * point @p{p}. @p{p} is a point
+ * on the reference element. If
+ * the finite element is
+ * vector-valued, then return the
+ * value of the only non-zero
+ * component of the vector value
+ * of this shape function. If the
+ * shape function has more than
+ * one non-zero component (which
+ * we refer to with the term
+ * non-primitive), then derived
+ * classes implementing this
+ * function should throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}. In
+ * that case, use the
+ * @ref{shape_value_component}
+ * function.
*
* An
* @p{ExcUnitShapeValuesDoNotExist}
* is thrown if the shape values
* of the @p{FiniteElement} under
- * consideration depend on the
+ * consideration depends on the
* shape of the cell in real
* space.
*/
- virtual double shape_value (const unsigned int i,
- const Point<dim> &p) const;
+ virtual double shape_value (const unsigned int i,
+ const Point<dim> &p) const;
+
+ /**
+ * Just like for @p{shape_value},
+ * but this function will be
+ * called when the shape function
+ * has more than one non-zero
+ * vector component. In that
+ * case, this function should
+ * return the value of the
+ * @p{component}-th vector
+ * component of the @p{i}th shape
+ * function at point @p{p}.
+ */
+ virtual double shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
/**
* Return the gradient of the
* likewise the gradient is the
* gradient on the unit cell with
* respect to unit cell
- * coordinates.
+ * coordinates. If
+ * the finite element is
+ * vector-valued, then return the
+ * value of the only non-zero
+ * component of the vector value
+ * of this shape function. If the
+ * shape function has more than
+ * one non-zero component (which
+ * we refer to with the term
+ * non-primitive), then derived
+ * classes implementing this
+ * function should throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}. In
+ * that case, use the
+ * @ref{shape_grad_component}
+ * function.
*
* An
* @p{ExcUnitShapeValuesDoNotExist}
* is thrown if the shape values
* of the @p{FiniteElement} under
- * consideration depend on the
+ * consideration depends on the
* shape of the cell in real
* space.
*/
virtual Tensor<1,dim> shape_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Just like for @p{shape_grad},
+ * but this function will be
+ * called when the shape function
+ * has more than one non-zero
+ * vector component. In that
+ * case, this function should
+ * return the gradient of the
+ * @p{component}-th vector
+ * component of the @p{i}th shape
+ * function at point @p{p}.
+ */
+ virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the tensor of second
* derivatives of the @p{i}th
* on the unit cell. The
* derivatives are derivatives on
* the unit cell with respect to
- * unit cell coordinates.
+ * unit cell coordinates. If
+ * the finite element is
+ * vector-valued, then return the
+ * value of the only non-zero
+ * component of the vector value
+ * of this shape function. If the
+ * shape function has more than
+ * one non-zero component (which
+ * we refer to with the term
+ * non-primitive), then derived
+ * classes implementing this
+ * function should throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}. In
+ * that case, use the
+ * @ref{shape_grad_grad_component}
+ * function.
*
* An
* @p{ExcUnitShapeValuesDoNotExist}
* is thrown if the shape values
* of the @p{FiniteElement} under
- * consideration depend on the
+ * consideration depends on the
* shape of the cell in real
* space.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const;
+ const Point<dim> &p) const;
+
+ /**
+ * Just like for @p{shape_grad_grad},
+ * but this function will be
+ * called when the shape function
+ * has more than one non-zero
+ * vector component. In that
+ * case, this function should
+ * return the gradient of the
+ * @p{component}-th vector
+ * component of the @p{i}th shape
+ * function at point @p{p}.
+ */
+ virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
/**
* Projection from a fine grid
const FullMatrix<double> & prolongate (const unsigned int child) const;
/**
- * Return a readonly reference to the
- * matrix which describes the constraints
- * at the interface between a refined and
- * an unrefined cell.
+ * Return a readonly reference to
+ * the matrix which describes the
+ * constraints at the interface
+ * between a refined and an
+ * unrefined cell.
*
- * The matrix is obviously empty in only
- * one space dimension, since there are no
- * constraints then.
+ * The matrix is obviously empty
+ * in only one space dimension,
+ * since there are no constraints
+ * then.
*/
const FullMatrix<double> & constraints () const;
/**
- * Comparison operator. We also check for
- * equality of the constraint matrix,
- * which is quite an expensive operation.
- * Do therefore
- * use this function with care, if possible
- * only for debugging purposes.
+ * Comparison operator. We also
+ * check for equality of the
+ * constraint matrix, which is
+ * quite an expensive operation.
+ * Do therefore use this function
+ * with care, if possible only
+ * for debugging purposes.
*
- * Since this function is not that important,
- * we avoid an implementational question
- * about comparing arrays and do not compare
- * the matrix arrays @p{restriction} and
+ * Since this function is not
+ * that important, we avoid an
+ * implementational question
+ * about comparing arrays and do
+ * not compare the matrix arrays
+ * @p{restriction} and
* @p{prolongation}.
*/
bool operator == (const FiniteElementBase<dim> &) const;
/**
- * Compute system index from components.
+ * Given a vector component and
+ * an index of a shape function
+ * within the shape functions
+ * corresponding to this vector
+ * component, return the index of
+ * this shape function within the
+ * shape functions of this
+ * element. If this is a scalar
+ * element, then the given
+ * component may only be zero,
+ * and the given component index
+ * is also the return value.
+ *
+ * If the finite element is
+ * vector-valued and has
+ * non-primitive shape functions,
+ * i.e. some of its shape
+ * functions are non-zero in more
+ * than just one vector
+ * component, then this function
+ * cannot be used since shape
+ * functions are no more
+ * associated with individual
+ * vector components, and an
+ * exception of type
+ * @p{ExcFENotPrimitive} is
+ * thrown.
*/
unsigned int component_to_system_index (const unsigned int component,
const unsigned int component_index) const;
/**
- * Compute component and index from
- * system index.
- *
- * Return value contains first
- * component and second index in
- * component.
- */
- std::pair<unsigned int,unsigned int>
- system_to_component_index (const unsigned int index) const;
-
- /**
- * Compute system index from components on a face.
+ * Same as above, but compute the
+ * data from the index of a shape
+ * function on a face.
*/
unsigned int face_component_to_system_index (const unsigned int component,
const unsigned int component_index) const;
-
+
/**
- * Compute component and index from system
- * index for a face.
+ * Compute vector component and
+ * index of this shape function
+ * within the shape functions
+ * corresponding to this
+ * component from the index of a
+ * shape function within this
+ * finite element.
*
- * Return value contains first
- * component and second index in
- * component.
+ * If the element is scalar, then
+ * the component is always zero,
+ * and the index within this
+ * component is equal to the
+ * overall index.
+ *
+ * If the shape function
+ * referenced has more than one
+ * non-zero component, then it
+ * cannot be associated with one
+ * vector component, and an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}
+ * will be raised.
+ *
+ * Note that if the element is
+ * composed of other (base)
+ * elements, and a base element
+ * has more than one component
+ * but all its shape functions
+ * are primitive (i.e. are
+ * non-zero in only one
+ * component), then this mapping
+ * contains valid
+ * information. However, the
+ * index of a shape function of
+ * this element within one
+ * component (i.e. the second
+ * number of the respective entry
+ * of this array) does not
+ * indicate the index of the
+ * respective shape function
+ * within the base element (since
+ * that has more than one
+ * vector-component). For this
+ * information, refer to the
+ * @p{system_to_base_table}
+ * field.
+ */
+ std::pair<unsigned int,unsigned int>
+ system_to_component_index (const unsigned int index) const;
+
+ /**
+ * Same as above, but do it for
+ * shape functions and their
+ * indices on a face.
*/
std::pair<unsigned int,unsigned int>
face_system_to_component_index (const unsigned int index) const;
/**
- * The base element establishing a
- * component.
+ * Given a vector component,
+ * return an index which base
+ * element implements this
+ * component, and which vector
+ * component is this base element
+ * this is. This information is
+ * only of interest for
+ * vector-valued finite elements
+ * which are composed of several
+ * sub-elements. In that case,
+ * one may want to obtain
+ * information about the element
+ * implementing a certain vector
+ * component, which can be done
+ * using this function and the
+ * @ref{FESystem}::@p{base_element}
+ * function.
*
- * This table converts a
- * component number to the
- * @p{base_element} number. While
- * component information contains
- * multiplicity of base elements,
- * the result allows access to
- * shape functions of the base
- * element.
+ * If this is a scalar finite
+ * element, then the return value
+ * is always equalt to zero.
*/
- unsigned int component_to_base (unsigned int index) const;
+ std::pair<unsigned int,unsigned int>
+ component_to_base (unsigned int component) const;
/**
* Access the @p{restriction_is_additive_flag}
* yields a non-empty array.
*/
bool has_face_support_points () const;
+
+ /**
+ * Return in which of the vector
+ * components of this finite
+ * element the @p{i}the shape
+ * function is non-zero. The
+ * length of the returned array
+ * is equal to the number of
+ * vector components of this
+ * element.
+ *
+ * For most finite element
+ * spaces, the result of this
+ * function will be a vector with
+ * exactly one element being
+ * @p{true}, since for most
+ * spaces the individual vector
+ * components are
+ * independent. Only for those
+ * spaces that couple the
+ * components, for example to
+ * make a shape function
+ * divergence free, will there be
+ * more than one @p{true} entry.
+ */
+ const std::vector<bool> &
+ get_nonzero_components (const unsigned int i) const;
+
+ /**
+ * Return in how many vector
+ * components the @p{i}th shape
+ * function is non-zero. This
+ * value equals the number of
+ * entries equal to @p{true} in
+ * the result of the
+ * @p{get_nonzero_components}
+ * function.
+ *
+ * For most finite element
+ * spaces, the result will be
+ * equal to one. It is not equal
+ * to one only for those ansatz
+ * spaces for which vector-valued
+ * shape functions couple the
+ * individual components, for
+ * example in order to make them
+ * divergence-free.
+ */
+ unsigned int
+ n_nonzero_components (const unsigned int i) const;
+
+ /**
+ * Return whether the @p{i}th
+ * shape function is primitive in
+ * the sense that the shape
+ * function is non-zero in only
+ * one vector
+ * component. Non-primitive shape
+ * functions would then, for
+ * example, be those of
+ * divergence free ansatz spaces,
+ * in which the individual vector
+ * components are coupled.
+ *
+ * The result of the function is
+ * @p{true} if and only if the
+ * result of
+ * @p{n_nonzero_components(i)} is
+ * equal to one.
+ */
+ bool
+ is_primitive (const unsigned int i) const;
+
+ /**
+ * Return whether the entire
+ * finite element is primitive,
+ * in the sense that all its
+ * shape functions are
+ * primitive. If the finite
+ * element is scalar, then this
+ * is always the case.
+ *
+ * Since this is an extremely
+ * common operations, the result
+ * is cached in the
+ * @p{cached_primitivity}
+ * variable which is computed in
+ * the constructor.
+ */
+ bool
+ is_primitive () const;
/**
* Determine an estimate for the
*/
unsigned int memory_consumption () const;
+ /**
+ * Exception
+ */
+ DeclException1 (ExcShapeFunctionNotPrimitive,
+ int,
+ << "The shape function with index " << arg1
+ << " is not primitive, i.e. it is vector-valued and "
+ << "has more than one non-zero vector component. This "
+ << "function cannot be called for these shape functions. "
+ << "Maybe you want to use the same function with the "
+ << "_component suffix?");
+ /**
+ * Exception
+ */
+ DeclException0 (ExcFENotPrimitive);
/**
* Exception
*/
FullMatrix<double> interface_constraints;
/**
- * Map between linear dofs and
- * component dofs.
+ * Store what
+ * @p{system_to_component_index}
+ * will return.
*/
std::vector< std::pair<unsigned int, unsigned int> > system_to_component_table;
/**
- * Map between linear dofs and
- * component dofs on face.
+ * Map between linear dofs and
+ * component dofs on face. This
+ * is filled with default values
+ * in the constructor, but
+ * derived classes will have to
+ * overwrite the information if
+ * necessary.
+ *
+ * By component, we mean the
+ * vector component, not the base
+ * element. The information thus
+ * makes only sense if a shape
+ * function is non-zero in only
+ * one component.
*/
std::vector< std::pair<unsigned int, unsigned int> > face_system_to_component_table;
+ /**
+ * For each shape function, store
+ * to which base element and
+ * which instance of this base
+ * element (in case its
+ * multiplicity is greater than
+ * one) it belongs, and its index
+ * within this base element. If
+ * the element is not composed of
+ * others, then base and instance
+ * are always zero, and the index
+ * is equal to the number of the
+ * shape function. If the element
+ * is composed of single
+ * instances of other elements
+ * (i.e. all with multiplicity
+ * one) all of which are scalar,
+ * then base values and dof
+ * indices within this element
+ * are equal to the
+ * @p{system_to_component_table}. It
+ * differs only in case the
+ * element is composed of other
+ * elements and at least one of
+ * them is vector-valued itself.
+ *
+ * This array has valid values
+ * also in the case of
+ * vector-value
+ * (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;
+
+ /**
+ * Likewise for the indices on
+ * faces.
+ */
+ std::vector<std::pair<std::pair<unsigned int,unsigned int>,unsigned int> > face_system_to_base_table;
+
/**
* Map between component and
- * linear dofs.
+ * linear dofs: For each pair of
+ * vector component and index
+ * within this component, store
+ * the global dof number in the
+ * composed element. If the
+ * element is scalar, then the
+ * outer (component) index can
+ * only be zero, and the inner
+ * index is equal to the stored
+ * value.
*/
std::vector< std::vector<unsigned int> > component_to_system_table;
* a component.
*
* This table converts a
- * component number to the
- * @p{base_element} number. While
- * component information contains
+ * component number to a pair
+ * consisting of the
+ * @p{base_element} number, and
+ * the component within this base
+ * element. While component
+ * information contains
* multiplicity of base elements,
* the result allows access to
* shape functions of the base
* case, the initialization by
* the base class is sufficient.
*/
- std::vector<unsigned int> component_to_base_table;
-
+ std::vector<std::pair<unsigned int, unsigned int> > component_to_base_table;
+
/**
* Projection matrices are
* concatenated or summed up.
*/
typename std::vector<Point<dim-1> > unit_face_support_points;
+ /**
+ * For each shape function, give
+ * a vector of bools (with size
+ * equal to the number of vector
+ * components which this finite
+ * element has) indicating in
+ * which component each of these
+ * shape functions is non-zero.
+ *
+ * For primitive elements, there
+ * is only one non-zero
+ * component.
+ */
+ const std::vector<std::vector<bool> > nonzero_components;
+
+ /**
+ * This array holds how many
+ * values in the respective entry
+ * of the @p{nonzero_components}
+ * element are non-zero. The
+ * array is thus a short-cut to
+ * allow faster access to this
+ * information than if we had to
+ * count the non-zero entries
+ * upon each request for this
+ * information. The field is
+ * initialized in the constructor
+ * of this class.
+ */
+ const std::vector<unsigned int> n_nonzero_components_table;
+
+ /**
+ * Store whether all shape
+ * functions are primitive. Since
+ * finding this out is a very
+ * common operation, we cache the
+ * result, i.e. compute the value
+ * in the constructor for simpler
+ * access.
+ */
+ const bool cached_primitivity;
+
/**
* Compute second derivatives by
* finite differences of
* 1e-6.
*/
static const double fd_step_length;
-
+
+ /**
+ * Given the pattern of nonzero
+ * components for each shape
+ * function, compute for each
+ * entry how many components are
+ * non-zero for each shape
+ * function. This function is
+ * used in the constructor of
+ * this class.
+ */
+ static
+ std::vector<unsigned int>
+ compute_n_nonzero_components (const std::vector<std::vector<bool> > &nonzero_components);
/**
* Allow the FESystem class to
* that throw if the matrices are
* not already initialized.
*/
- friend class FESystem<dim>;
+ template <int dim_> friend class FESystem;
};
Assert(component_index<component_to_system_table[component].size(),
ExcIndexRange(component_index, 0,
component_to_system_table[component].size()));
+ Assert (is_primitive(),
+ typename FiniteElementBase<dim>::ExcFENotPrimitive());
+
return component_to_system_table[component][component_index];
}
std::pair<unsigned int,unsigned int>
FiniteElementBase<dim>::system_to_component_index (const unsigned int index) const
{
- Assert(index < system_to_component_table.size(),
+ Assert (index < system_to_component_table.size(),
ExcIndexRange(index, 0, system_to_component_table.size()));
+ Assert (is_primitive (index),
+ typename FiniteElementBase<dim>::ExcShapeFunctionNotPrimitive(index));
return system_to_component_table[index];
}
Assert(component_index<face_component_to_system_table[component].size(),
ExcIndexRange(component_index, 0,
face_component_to_system_table[component].size()));
+ Assert (is_primitive(),
+ typename FiniteElementBase<dim>::ExcFENotPrimitive());
+
return face_component_to_system_table[component][component_index];
}
{
Assert(index < face_system_to_component_table.size(),
ExcIndexRange(index, 0, face_system_to_component_table.size()));
+//TODO: check for primitivity of this shape function. this needs the global dof index
+// Assert (is_primitive (face_to_cell_index(index)),
+// typename FiniteElementBase<dim>::ExcShapeFunctionNotPrimitive(index));
return face_system_to_component_table[index];
}
template <int dim>
inline
-unsigned int
+std::pair<unsigned int,unsigned int>
FiniteElementBase<dim>::component_to_base (unsigned int index) const
{
- if (n_components() == 1)
- return 0;
Assert(index < component_to_base_table.size(),
ExcIndexRange(index, 0, component_to_base_table.size()));
+
return component_to_base_table[index];
}
}
+template <int dim>
+inline
+const std::vector<bool> &
+FiniteElementBase<dim>::get_nonzero_components (const unsigned int i) const
+{
+ Assert (i < dofs_per_cell, ExcIndexRange (i, 0, dofs_per_cell));
+ return nonzero_components[i];
+};
+
+
+
+template <int dim>
+inline
+unsigned int
+FiniteElementBase<dim>::n_nonzero_components (const unsigned int i) const
+{
+ Assert (i < dofs_per_cell, ExcIndexRange (i, 0, dofs_per_cell));
+ return n_nonzero_components_table[i];
+};
+
+
+
+template <int dim>
+inline
+bool
+FiniteElementBase<dim>::is_primitive (const unsigned int i) const
+{
+ Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, dofs_per_cell));
+
+ // return primitivity of a shape
+ // function by checking whether it
+ // has more than one non-zero
+ // component or not. we could cache
+ // this value in an array of bools,
+ // but accessing a bit-vector (as
+ // std::vector<bool> is) is
+ // probably more expensive than
+ // just comparing against 1
+ return (n_nonzero_components_table[i] == 1);
+};
+
+
+template <int dim>
+inline
+bool
+FiniteElementBase<dim>::is_primitive () const
+{
+ return cached_primitivity;
+};
+
+
+
+
#endif
#include <base/config.h>
#include <base/polynomial.h>
+#include <base/polynomial_space.h>
#include <fe/fe.h>
template <int dim> class PolynomialSpace;
* Constructor for tensor product
* polynomials of degree @p{k}.
*/
- FE_DGP (unsigned int k);
- /**
- * Destructor.
- */
- ~FE_DGP ();
+ FE_DGP (const unsigned int k);
/**
* Return the value of the
* @p{i}th shape function at the
- * point @p{p}. @p{p} is a point
- * on the reference element.
+ * point @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual double shape_value (const unsigned int i,
const Point<dim> &p) const;
+
+ /**
+ * Return the value of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual double shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
/**
* Return the gradient of the
* @p{i}th shape function at the
- * point @p{p}. @p{p} is a point
- * on the reference element, and
- * likewise the gradient is the
- * gradient on the unit cell with
- * respect to unit cell
- * coordinates.
+ * point @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual Tensor<1,dim> shape_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the gradient of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the tensor of second
* derivatives of the @p{i}th
* shape function at point @p{p}
- * on the unit cell. The
- * derivatives are derivatives on
- * the unit cell with respect to
- * unit cell coordinates.
+ * on the unit cell. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the second derivative
+ * of the @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the polynomial degree
* of this finite element,
/**
* Declare a nested class which
- * will has static definitions of
+ * will hold static definitions of
* various matrices such as
* constraint and embedding
* matrices. The definition of
* the various static fields are
- * in the files @p{fe_q_[123]d.cc}
+ * in the files @p{fe_dgp_[123]d.cc}
* in the source directory.
*/
struct Matrices
static std::vector<unsigned int> get_dpo_vector(unsigned int degree);
/**
- * Compute flags for initial update only.
+ * Given a set of flags indicating
+ * what quantities are requested
+ * from a @p{FEValues} object,
+ * return which of these can be
+ * precomputed once and for
+ * all. Often, the values of
+ * shape function at quadrature
+ * points can be precomputed, for
+ * example, in which case the
+ * return value of this function
+ * would be the logical and of
+ * the input @p{flags} and
+ * @p{update_values}.
+ *
+ * For the present kind of finite
+ * element, this is exactly the
+ * case.
*/
virtual UpdateFlags update_once (const UpdateFlags flags) const;
/**
- * Compute flags for update on each cell.
+ * This is the opposite to the
+ * above function: given a set of
+ * flags indicating what we want
+ * to know, return which of these
+ * need to be computed each time
+ * we visit a new cell.
+ *
+ * If for the computation of one
+ * quantity something else is
+ * also required (for example, we
+ * often need the covariant
+ * transformation when gradients
+ * need to be computed), include
+ * this in the result as well.
*/
virtual UpdateFlags update_each (const UpdateFlags flags) const;
const unsigned int degree;
/**
- * Pointer to the tensor
- * product polynomials.
+ * Pointer to an object
+ * representing the polynomial
+ * space used here.
*/
- PolynomialSpace<dim>* poly;
+ const PolynomialSpace<dim> polynomial_space;
/**
* Fields of cell-independent data.
+ *
+ * For information about the
+ * general purpose of this class,
+ * see the documentation of the
+ * base class.
*/
class InternalData : public FiniteElementBase<dim>::InternalDataBase
{
* function, containing
* values for each quadrature
* point.
+ *
+ * In this array, we store
+ * the values of the shape
+ * function in the quadrature
+ * points on the unit
+ * cell. Since these values
+ * do not change under
+ * transformation to the real
+ * cell, we only need to copy
+ * them over when visiting a
+ * concrete cell.
*/
std::vector<std::vector<double> > shape_values;
* function, containing
* values for each quadrature
* point.
+ *
+ * We store the gradients in
+ * the quadrature points on
+ * the unit cell. We then
+ * only have to apply the
+ * transformation (which is a
+ * matrix-vector
+ * multiplication) when
+ * visiting an actual cell.
*/
typename std::vector<std::vector<Tensor<1,dim> > > shape_gradients;
};
#include <base/config.h>
#include <base/polynomial.h>
+#include <base/tensor_product_polynomials.h>
#include <fe/fe.h>
template <int dim> class TensorProductPolynomials;
* Constructor for tensor product
* polynomials of degree @p{k}.
*/
- FE_DGQ (unsigned int k);
- /**
- * Destructor.
- */
- ~FE_DGQ ();
+ FE_DGQ (const unsigned int k);
/**
* Return the value of the
* @p{i}th shape function at the
- * point @p{p}. @p{p} is a point
- * on the reference element.
+ * point @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual double shape_value (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the value of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual double shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the gradient of the
* @p{i}th shape function at the
- * point @p{p}. @p{p} is a point
- * on the reference element, and
- * likewise the gradient is the
- * gradient on the unit cell with
- * respect to unit cell
- * coordinates.
+ * point @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual Tensor<1,dim> shape_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the gradient of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the tensor of second
* derivatives of the @p{i}th
* shape function at point @p{p}
- * on the unit cell. The
- * derivatives are derivatives on
- * the unit cell with respect to
- * unit cell coordinates.
+ * on the unit cell. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the second derivative
+ * of the @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the polynomial degree
* of this finite element,
/**
* Declare a nested class which
- * will has static definitions of
+ * will hold static definitions of
* various matrices such as
* constraint and embedding
* matrices. The definition of
* the various static fields are
- * in the files @p{fe_q_[123]d.cc}
+ * in the files @p{fe_dgq_[123]d.cc}
* in the source directory.
*/
struct Matrices
static std::vector<unsigned int> get_dpo_vector(unsigned int degree);
/**
- * Compute flags for initial update only.
+ * Given a set of flags indicating
+ * what quantities are requested
+ * from a @p{FEValues} object,
+ * return which of these can be
+ * precomputed once and for
+ * all. Often, the values of
+ * shape function at quadrature
+ * points can be precomputed, for
+ * example, in which case the
+ * return value of this function
+ * would be the logical and of
+ * the input @p{flags} and
+ * @p{update_values}.
+ *
+ * For the present kind of finite
+ * element, this is exactly the
+ * case.
*/
virtual UpdateFlags update_once (const UpdateFlags flags) const;
/**
- * Compute flags for update on each cell.
+ * This is the opposite to the
+ * above function: given a set of
+ * flags indicating what we want
+ * to know, return which of these
+ * need to be computed each time
+ * we visit a new cell.
+ *
+ * If for the computation of one
+ * quantity something else is
+ * also required (for example, we
+ * often need the covariant
+ * transformation when gradients
+ * need to be computed), include
+ * this in the result as well.
*/
virtual UpdateFlags update_each (const UpdateFlags flags) const;
* Pointer to the tensor
* product polynomials.
*/
- TensorProductPolynomials<dim>* poly;
+ const TensorProductPolynomials<dim> polynomial_space;
/**
* Fields of cell-independent data.
+ *
+ * For information about the
+ * general purpose of this class,
+ * see the documentation of the
+ * base class.
*/
class InternalData : public FiniteElementBase<dim>::InternalDataBase
{
* function, containing
* values for each quadrature
* point.
+ *
+ * In this array, we store
+ * the values of the shape
+ * function in the quadrature
+ * points on the unit
+ * cell. Since these values
+ * do not change under
+ * transformation to the real
+ * cell, we only need to copy
+ * them over when visiting a
+ * concrete cell.
*/
std::vector<std::vector<double> > shape_values;
* function, containing
* values for each quadrature
* point.
+ *
+ * We store the gradients in
+ * the quadrature points on
+ * the unit cell. We then
+ * only have to apply the
+ * transformation (which is a
+ * matrix-vector
+ * multiplication) when
+ * visiting an actual cell.
*/
typename std::vector<std::vector<Tensor<1,dim> > > shape_gradients;
};
* access to build_renumbering
* function.
*/
- friend class MappingQ<dim>;
+ template <int dim1> friend class MappingQ;
};
--- /dev/null
+//---------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002 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_nedelec_h
+#define __deal2__fe_nedelec_h
+
+#include <base/config.h>
+#include <base/polynomial.h>
+#include <base/tensor_product_polynomials.h>
+#include <fe/fe.h>
+
+template <int dim> class TensorProductPolynomials;
+template <int dim> class MappingQ;
+
+
+
+/**
+ * Implementation of continuous Nedelec elements for the space H_curl.
+ *
+ * The constructor of this class takes the degree @p{p} of this finite
+ * element.
+ *
+ *
+ * @sect3{Numbering of the degrees of freedom (DoFs)}
+ *
+ * Nedelec elements have their degrees of freedom on edges, with shape
+ * functions being vector valued and pointing in tangential
+ * direction. We use the standard enumeration and direction of edges
+ * in deal.II, yielding the following shape functions in 2d:
+ *
+ * @begin{verbatim}
+ * 2
+ * *--->---*
+ * | |
+ * 3^ ^1
+ * | |
+ * *--->---*
+ * 0
+ * @end{verbatim}
+ *
+ * For the 3d case, the ordering follows the same scheme: the lines
+ * are numbered as described in the documentation of the
+ * @ref{Triangulation} class, i.e.
+ * @begin{verbatim}
+ * *---6---* *---6---*
+ * /| | / /|
+ * 11 | 5 11 10 5
+ * / 7 | / / |
+ * * | | *---2---* |
+ * | *---4---* | | *
+ * | / / | 1 /
+ * 3 8 9 3 | 9
+ * |/ / | |/
+ * *---0---* *---0---*
+ * @end{verbatim}
+ * and their directions are as follows:
+ * @begin{verbatim}
+ * *--->---* *--->---*
+ * /| | / /|
+ * ^ | ^ ^ ^ ^
+ * / ^ | / / |
+ * * | | *--->---* |
+ * | *--->---* | | *
+ * | / / | ^ /
+ * ^ ^ ^ ^ | ^
+ * |/ / | |/
+ * *--->---* *--->---*
+ * @end{verbatim}
+ *
+ * The element does not make much sense in 1d, so it is not
+ * implemented there.
+ *
+ *
+ * @author Anna Schneebeli, Wolfgang Bangerth, 2002
+ */
+template <int dim>
+class FE_Nedelec : public FiniteElement<dim>
+{
+ public:
+ /**
+ * Constructor for the Nedelec
+ * element of degree @p{p}.
+ */
+ FE_Nedelec (const unsigned int p);
+
+ /**
+ * Return the value of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ */
+ virtual double shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
+ /**
+ * Return the gradient of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ */
+ virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
+ /**
+ * Return the second derivative
+ * of the @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ */
+ virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
+ /**
+ * Return the polynomial degree
+ * of this finite element,
+ * i.e. the value passed to the
+ * constructor.
+ */
+ unsigned int get_degree () const;
+
+ /**
+ * Number of base elements in a
+ * mixed discretization. Here,
+ * this is of course equal to
+ * one.
+ */
+ virtual unsigned int n_base_elements () const;
+
+ /**
+ * Access to base element
+ * objects. Since this element is
+ * atomic, @p{base_element(0)} is
+ * @p{this}, and all other
+ * indices throw an error.
+ */
+ virtual const FiniteElement<dim> & base_element (const unsigned int index) const;
+
+ /**
+ * Check for non-zero values on a face.
+ *
+ * This function returns
+ * @p{true}, if the shape
+ * function @p{shape_index} has
+ * non-zero values on the face
+ * @p{face_index}.
+ *
+ * Implementation of the
+ * interface in
+ * @ref{FiniteElement}
+ */
+ virtual bool has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const;
+
+ /**
+ * Determine an estimate for the
+ * memory consumption (in bytes)
+ * of this object.
+ *
+ * This function is made virtual,
+ * since finite element objects
+ * are usually accessed through
+ * pointers to their base class,
+ * rather than the class itself.
+ */
+ virtual unsigned int memory_consumption () const;
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNotUsefulInThisDimension);
+
+ protected:
+ /**
+ * @p{clone} function instead of
+ * a copy constructor.
+ *
+ * This function is needed by the
+ * constructors of @p{FESystem}.
+ */
+ virtual FiniteElement<dim> * clone() const;
+
+ /**
+ * Prepare internal data
+ * structures and fill in values
+ * independent of the cell.
+ */
+ virtual
+ typename Mapping<dim>::InternalDataBase *
+ get_data (const UpdateFlags,
+ const Mapping<dim>& mapping,
+ const Quadrature<dim>& quadrature) const ;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const Quadrature<dim> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_face_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const ;
+
+ /**
+ * Implementation of the same
+ * function in
+ * @ref{FiniteElement}.
+ */
+ virtual void
+ fill_fe_subface_values (const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face_no,
+ const unsigned int sub_no,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_internal,
+ typename Mapping<dim>::InternalDataBase &fe_internal,
+ FEValuesData<dim>& data) const ;
+
+ private:
+
+ /**
+ * Declare a nested class which
+ * will hold static definitions of
+ * various matrices such as
+ * constraint and embedding
+ * matrices. The definition of
+ * the various static fields are
+ * in the files @p{fe_nedelec_[23]d.cc}
+ * in the source directory.
+ */
+ struct Matrices
+ {
+ /**
+ * Embedding matrices. For
+ * each element type (the
+ * first index) there are as
+ * many embedding matrices as
+ * there are children per
+ * cell. The first index
+ * starts with linear
+ * elements and goes up in
+ * polynomial degree. The
+ * array may grow in the
+ * future with the number of
+ * elements for which these
+ * matrices have been
+ * computed. If for some
+ * element, the matrices have
+ * not been computed then you
+ * may use the element
+ * nevertheless but can not
+ * access the respective
+ * fields.
+ */
+ static const double * const
+ embedding[][GeometryInfo<dim>::children_per_cell];
+
+ /**
+ * Number of elements (first
+ * index) the above field
+ * has. Equals the highest
+ * polynomial degree for
+ * which the embedding
+ * matrices have been
+ * computed.
+ */
+ static const unsigned int n_embedding_matrices;
+
+ /**
+ * As the
+ * @p{embedding_matrices}
+ * field, but for the
+ * interface constraints. One
+ * for each element for which
+ * it has been computed.
+ */
+ static const double * const constraint_matrices[];
+
+ /**
+ * Like
+ * @p{n_embedding_matrices},
+ * but for the number of
+ * interface constraint
+ * matrices.
+ */
+ static const unsigned int n_constraint_matrices;
+ };
+
+ /**
+ * Only for internal use. Its
+ * full name is
+ * @p{get_dofs_per_object_vector}
+ * function and it creates the
+ * @p{dofs_per_object} vector that is
+ * needed within the constructor to
+ * be passed to the constructor of
+ * @p{FiniteElementData}.
+ */
+ static std::vector<unsigned int> get_dpo_vector(const unsigned int degree);
+
+ /**
+ * Initialize the
+ * @p{unit_support_points} field
+ * of the @ref{FiniteElementBase}
+ * class. Called from the
+ * constructor.
+ */
+ void initialize_unit_support_points ();
+
+ /**
+ * Initialize the
+ * @p{unit_face_support_points} field
+ * of the @ref{FiniteElementBase}
+ * class. Called from the
+ * constructor.
+ */
+ void initialize_unit_face_support_points ();
+
+ /**
+ * Given a set of flags indicating
+ * what quantities are requested
+ * from a @p{FEValues} object,
+ * return which of these can be
+ * precomputed once and for
+ * all. Often, the values of
+ * shape function at quadrature
+ * points can be precomputed, for
+ * example, in which case the
+ * return value of this function
+ * would be the logical and of
+ * the input @p{flags} and
+ * @p{update_values}.
+ *
+ * For the present kind of finite
+ * element, this is exactly the
+ * case.
+ */
+ virtual UpdateFlags update_once (const UpdateFlags flags) const;
+
+ /**
+ * This is the opposite to the
+ * above function: given a set of
+ * flags indicating what we want
+ * to know, return which of these
+ * need to be computed each time
+ * we visit a new cell.
+ *
+ * If for the computation of one
+ * quantity something else is
+ * also required (for example, we
+ * often need the covariant
+ * transformation when gradients
+ * need to be computed), include
+ * this in the result as well.
+ */
+ virtual UpdateFlags update_each (const UpdateFlags flags) const;
+
+ /**
+ * Degree of the polynomials.
+ */
+ const unsigned int degree;
+
+ /**
+ * Fields of cell-independent data.
+ *
+ * For information about the
+ * general purpose of this class,
+ * see the documentation of the
+ * base class.
+ */
+ class InternalData : public FiniteElementBase<dim>::InternalDataBase
+ {
+ public:
+ /**
+ * Array with shape function
+ * values in quadrature
+ * points. There is one
+ * vector for each shape
+ * function, containing
+ * values for each quadrature
+ * point. Since the shape
+ * functions are
+ * vector-valued (with as
+ * many components as there
+ * are space dimensions), the
+ * value is a tensor.
+ *
+ * In this array, we store
+ * the values of the shape
+ * function in the quadrature
+ * points on the unit
+ * cell. The transformation
+ * to the real space cell is
+ * then simply done by
+ * multiplication with the
+ * Jacobian of the mapping.
+ */
+ std::vector<std::vector<Tensor<1,dim> > > shape_values;
+
+ /**
+ * Array with shape function
+ * gradients in quadrature
+ * points. There is one
+ * vector for each shape
+ * function, containing
+ * values for each quadrature
+ * point.
+ *
+ * We store the gradients in
+ * the quadrature points on
+ * the unit cell. We then
+ * only have to apply the
+ * transformation (which is a
+ * matrix-vector
+ * multiplication) when
+ * visiting an actual cell.
+ */
+ typename std::vector<typename std::vector<Tensor<2,dim> > > shape_gradients;
+ };
+
+ /**
+ * Allow access from other
+ * dimensions.
+ */
+ template <int dim1> friend class FE_Nedelec;
+};
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+template <> void FE_Nedelec<1>::initialize_unit_face_support_points ();
+
+// declaration of explicit specializations of member variables, if the
+// compiler allows us to do that (the standard says we must)
+#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG
+template <>
+const double * const
+FE_Nedelec<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell];
+
+template <>
+const unsigned int FE_Nedelec<1>::Matrices::n_embedding_matrices;
+
+template <>
+const double * const FE_Nedelec<1>::Matrices::constraint_matrices[];
+
+template <>
+const unsigned int FE_Nedelec<1>::Matrices::n_constraint_matrices;
+
+template <>
+const double * const
+FE_Nedelec<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell];
+
+template <>
+const unsigned int FE_Nedelec<2>::Matrices::n_embedding_matrices;
+
+template <>
+const double * const FE_Nedelec<2>::Matrices::constraint_matrices[];
+
+template <>
+const unsigned int FE_Nedelec<2>::Matrices::n_constraint_matrices;
+
+template <>
+const double * const
+FE_Nedelec<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell];
+
+template <>
+const unsigned int FE_Nedelec<3>::Matrices::n_embedding_matrices;
+
+template <>
+const double * const FE_Nedelec<3>::Matrices::constraint_matrices[];
+
+template <>
+const unsigned int FE_Nedelec<3>::Matrices::n_constraint_matrices;
+
+#endif
+
+#endif
#include <base/config.h>
#include <base/polynomial.h>
+#include <base/tensor_product_polynomials.h>
#include <fe/fe.h>
template <int dim> class TensorProductPolynomials;
*/
FE_Q (const unsigned int p);
- /**
- * Destructor.
- */
- ~FE_Q ();
-
/**
* Return the value of the
* @p{i}th shape function at the
- * point @p{p}. @p{p} is a point
- * on the reference element.
+ * point @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual double shape_value (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the value of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual double shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the gradient of the
* @p{i}th shape function at the
- * point @p{p}. @p{p} is a point
- * on the reference element, and
- * likewise the gradient is the
- * gradient on the unit cell with
- * respect to unit cell
- * coordinates.
+ * point @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual Tensor<1,dim> shape_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the gradient of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the tensor of second
* derivatives of the @p{i}th
* shape function at point @p{p}
- * on the unit cell. The
- * derivatives are derivatives on
- * the unit cell with respect to
- * unit cell coordinates.
+ * on the unit cell. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
*/
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the second derivative
+ * of the @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is scalar,
+ * the returned value is the same
+ * as if the function without the
+ * @p{_component} suffix were
+ * called, provided that the
+ * specified component is zero.
+ */
+ virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the polynomial degree
* of this finite element,
/**
* Declare a nested class which
- * will has static definitions of
+ * will hold static definitions of
* various matrices such as
* constraint and embedding
* matrices. The definition of
void initialize_unit_face_support_points ();
/**
- * Compute flags for initial
- * update only.
+ * Given a set of flags indicating
+ * what quantities are requested
+ * from a @p{FEValues} object,
+ * return which of these can be
+ * precomputed once and for
+ * all. Often, the values of
+ * shape function at quadrature
+ * points can be precomputed, for
+ * example, in which case the
+ * return value of this function
+ * would be the logical and of
+ * the input @p{flags} and
+ * @p{update_values}.
+ *
+ * For the present kind of finite
+ * element, this is exactly the
+ * case.
*/
virtual UpdateFlags update_once (const UpdateFlags flags) const;
/**
- * Compute flags for update on
- * each cell.
+ * This is the opposite to the
+ * above function: given a set of
+ * flags indicating what we want
+ * to know, return which of these
+ * need to be computed each time
+ * we visit a new cell.
+ *
+ * If for the computation of one
+ * quantity something else is
+ * also required (for example, we
+ * often need the covariant
+ * transformation when gradients
+ * need to be computed), include
+ * this in the result as well.
*/
virtual UpdateFlags update_each (const UpdateFlags flags) const;
* Pointer to the tensor
* product polynomials.
*/
- TensorProductPolynomials<dim>* poly;
+ const TensorProductPolynomials<dim> polynomial_space;
/**
* Fields of cell-independent data.
+ *
+ * For information about the
+ * general purpose of this class,
+ * see the documentation of the
+ * base class.
*/
class InternalData : public FiniteElementBase<dim>::InternalDataBase
{
public:
/**
- * Array with shape function values
- * in quadrature points. There is one
- * vector for each shape function, containing
- * values for each quadrature point.
+ * Array with shape function
+ * values in quadrature
+ * points. There is one
+ * vector for each shape
+ * function, containing
+ * values for each quadrature
+ * point.
+ *
+ * In this array, we store
+ * the values of the shape
+ * function in the quadrature
+ * points on the unit
+ * cell. Since these values
+ * do not change under
+ * transformation to the real
+ * cell, we only need to copy
+ * them over when visiting a
+ * concrete cell.
*/
std::vector<std::vector<double> > shape_values;
/**
- * Array with shape function gradients
- * in quadrature points. There is one
- * vector for each shape function, containing
- * values for each quadrature point.
- */
+ * Array with shape function
+ * gradients in quadrature
+ * points. There is one
+ * vector for each shape
+ * function, containing
+ * values for each quadrature
+ * point.
+ *
+ * We store the gradients in
+ * the quadrature points on
+ * the unit cell. We then
+ * only have to apply the
+ * transformation (which is a
+ * matrix-vector
+ * multiplication) when
+ * visiting an actual cell.
+ */
typename std::vector<typename std::vector<Tensor<1,dim> > > shape_gradients;
};
/**
- * Allow access from other dimensions.
+ * Allow access from other
+ * dimensions. We need this since
+ * we want to call the functions
+ * @p{get_dpo_vector} and
+ * @p{lexicographic_to_hierarchic_numbering}
+ * for the faces of the finite
+ * element of dimension dim+1.
*/
template <int dim1> friend class FE_Q;
};
* into one. To the outside world, the resulting object looks just
* like a usual finite element object, which is composed of several
* other finite elements that are possibly of different type.
-
+ *
* The overall numbering of degrees of freedom is as follows: for each
* subobject (vertex, line, quad, or hex), the degrees of freedom are
* numbered such that we run over all subelements first, before
* Return the value of the
* @p{i}th shape function at the
* point @p{p}. @p{p} is a point
- * on the reference element.
+ * on the reference element. Since
+ * this finite element is always
+ * vector-valued, we return the
+ * value of the only non-zero
+ * component of the vector value
+ * of this shape function. If the
+ * shape function has more than
+ * one non-zero component (which
+ * we refer to with the term
+ * non-primitive), then throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}.
*
* An
* @p{ExcUnitShapeValuesDoNotExist}
*/
virtual double shape_value (const unsigned int i,
const Point<dim> &p) const;
+
+ /**
+ * Return the value of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is vector
+ * valued in general, it relays
+ * the computation of these
+ * values to the base elements.
+ */
+ virtual double shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
/**
* Return the gradient of the
* likewise the gradient is the
* gradient on the unit cell with
* respect to unit cell
- * coordinates.
+ * coordinates. Since
+ * this finite element is always
+ * vector-valued, we return the
+ * value of the only non-zero
+ * component of the vector value
+ * of this shape function. If the
+ * shape function has more than
+ * one non-zero component (which
+ * we refer to with the term
+ * non-primitive), then throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}.
*
* An
* @p{ExcUnitShapeValuesDoNotExist}
virtual Tensor<1,dim> shape_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the gradient of the
+ * @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is vector
+ * valued in general, it relays
+ * the computation of these
+ * values to the base elements.
+ */
+ virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Return the tensor of second
* derivatives of the @p{i}th
* on the unit cell. The
* derivatives are derivatives on
* the unit cell with respect to
- * unit cell coordinates.
+ * unit cell coordinates. Since
+ * this finite element is always
+ * vector-valued, we return the
+ * value of the only non-zero
+ * component of the vector value
+ * of this shape function. If the
+ * shape function has more than
+ * one non-zero component (which
+ * we refer to with the term
+ * non-primitive), then throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}.
*
* An
* @p{ExcUnitShapeValuesDoNotExist}
virtual Tensor<2,dim> shape_grad_grad (const unsigned int i,
const Point<dim> &p) const;
+ /**
+ * Return the second derivatives
+ * of the @p{component}th vector
+ * component of the @p{i}th shape
+ * function at the point
+ * @p{p}. See the
+ * @ref{FiniteElementBase} base
+ * class for more information
+ * about the semantics of this
+ * function.
+ *
+ * Since this element is vector
+ * valued in general, it relays
+ * the computation of these
+ * values to the base elements.
+ */
+ virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const;
+
/**
* Number of different base
* elements of this object.
* multiplicities are greater
* than one.
*/
- virtual const FiniteElement<dim> & base_element (const unsigned int index) const;
+ virtual const FiniteElement<dim> &
+ base_element (const unsigned int index) const;
/**
* Check for non-zero values on a face.
const unsigned int N2,
const FiniteElement<dim> &fe3,
const unsigned int N3);
+
+ /**
+ * Compute the non-zero vector
+ * components of a composed
+ * finite element.
+ */
+ static std::vector<std::vector<bool> >
+ compute_nonzero_components (const FiniteElement<dim> &fe1,
+ const unsigned int N1);
+
+ /**
+ * Compute the non-zero vector
+ * components of a composed
+ * finite element.
+ */
+ static std::vector<std::vector<bool> >
+ compute_nonzero_components (const FiniteElement<dim> &fe1,
+ const unsigned int N1,
+ const FiniteElement<dim> &fe2,
+ const unsigned int N2);
+
+ /**
+ * Compute the non-zero vector
+ * components of a composed
+ * finite element.
+ */
+ static std::vector<std::vector<bool> >
+ compute_nonzero_components (const FiniteElement<dim> &fe1,
+ const unsigned int N1,
+ const FiniteElement<dim> &fe2,
+ const unsigned int N2,
+ const FiniteElement<dim> &fe3,
+ const unsigned int N3);
+
+ /**
+ * Compute the nonzero components
+ * of a list of finite elements
+ * with multiplicities given in
+ * the second argument. This
+ * function is called from all
+ * the above functions.
+ */
+ static std::vector<std::vector<bool> >
+ compute_nonzero_components (const std::vector<const FiniteElement<dim>*> &fes,
+ const std::vector<unsigned int> &multiplicities);
/**
* This function is simply
/**
* Used by @p{initialize}.
*/
- void build_cell_table();
+ void build_cell_tables();
/**
* Used by @p{initialize}.
*/
- void build_face_table();
+ void build_face_tables();
/**
* Used by @p{initialize}.
* Usually: Fields of
* cell-independent data.
*
- * But for @p{FESystem} this
- * @p{InternalData} class does
+ * However, here, this class does
* not itself store the data but
* only pointers to
- * @p{InternalDatas} of the base
- * elements.
+ * @p{InternalData} objects for
+ * each of the base elements.
*/
class InternalData : public FiniteElementBase<dim>::InternalDataBase
{
* vector to
* @p{n_base_elements}.
*/
- InternalData(const unsigned int n_base_elements);
+ InternalData (const unsigned int n_base_elements);
/**
* Destructor. Deletes all
private:
/**
- * Pointers to the
- * @p{InternalDatas} of the
- * base elements. They are
+ * Pointers to
+ * @p{InternalData} objects
+ * for each of the base
+ * elements. They are
* accessed to by the
* @p{set_} and
- * @p{get_fe_data}
- * functions.
+ * @p{get_fe_data} functions.
*
* The size of this vector is
* set to @p{n_base_elements}
* by the InternalData
- * constructor. It is
- * filled by the @p{get_data}
- * function.
+ * constructor. It is filled
+ * by the @p{get_data}
+ * function. Note that since
+ * the data for each instance
+ * of a base class is
+ * necessarily the same, we
+ * only need as many of these
+ * objects as there are base
+ * elements, irrespective of
+ * their multiplicity.
*/
typename std::vector<typename FiniteElementBase<dim>::InternalDataBase *> base_fe_datas;
/**
* Pointers to the
- * @p{FEValuesDatas}
+ * @p{FEValuesData} objects
* that are given to the
* @p{fill_fe_values}
* function of the base
#include <fe/fe_update_flags.h>
#include <fe/mapping.h>
+#include <algorithm>
+
template <int dim> class Quadrature;
{
public:
/**
- * Initialize all vectors to correct size.
+ * Initialize all vectors to
+ * correct size.
*/
- void initialize (const unsigned int n_quadrature_points,
- const unsigned int n_shapes,
- const UpdateFlags flags);
+ void initialize (const unsigned int n_quadrature_points,
+ const FiniteElement<dim> &fe,
+ const UpdateFlags flags);
/**
- * Storage type for shape values.
+ * Storage type for shape
+ * values. Each row in the matrix
+ * denotes the values of a single
+ * shape function at the
+ * different points, columns are
+ * for a single point with the
+ * different shape functions.
+ *
+ * If a shape function has more
+ * than one non-zero component
+ * (in deal.II diction: it is
+ * non-primitive), then we
+ * allocate one row per non-zero
+ * component, and shift
+ * subsequent rows backward.
+ * Lookup of the correct row for
+ * a shape function is thus
+ * simple in case the entire
+ * finite element is primitive
+ * (i.e. all shape functions are
+ * primitive), since then the
+ * shape function number equals
+ * the row number. Otherwise, use
+ * the
+ * @p{shape_function_to_row_table}
+ * array to get at the first row
+ * that belongs to this
+ * particular shape function, and
+ * navigate among all the rows
+ * for this shape function using
+ * the
+ * @p{FiniteElement::get_nonzero_components}
+ * function which tells us which
+ * components are non-zero and
+ * thus have a row in the array
+ * presently under discussion.
*/
typedef FullMatrix<double> ShapeVector;
/**
- * Storage type for gradients.
+ * Storage type for
+ * gradients. The layout of data
+ * is the same as for the
+ * @ref{ShapeVector} data type.
*/
- typedef typename std::vector<typename std::vector<Tensor<1,dim> > > GradientVector;
+ typedef
+ typename std::vector<typename std::vector<Tensor<1,dim> > >
+ GradientVector;
+
+ /**
+ * Likewise for second order
+ * derivatives.
+ */
+ typedef
+ typename std::vector<typename std::vector<Tensor<2,dim> > >
+ GradGradVector;
/**
- * Store the values of the shape functions
- * at the quadrature points. Rows in the
- * matrices denote the values of a single
- * shape function at the different points,
- * columns are for a single point with the
- * different shape functions.
+ * Store the values of the shape
+ * functions at the quadrature
+ * points. See the description of
+ * the data type for the layout
+ * of the data in this field.
*/
ShapeVector shape_values;
/**
- * Store the gradients of the shape
- * functions at the quadrature points.
- * For each shape function, there is a vector
- * of the gradients in each quadrature point.
+ * Store the gradients of the
+ * shape functions at the
+ * quadrature points. See the
+ * description of the data type
+ * for the layout of the data in
+ * this field.
*/
GradientVector shape_gradients;
/**
- * Store the 2nd derivatives of the shape
- * functions at the quadrature points.
- * For each shape function, there is a vector
- * of the 2nd gradients in each quadrature point.
+ * Store the 2nd derivatives of
+ * the shape functions at the
+ * quadrature points. See the
+ * description of the data type
+ * for the layout of the data in
+ * this field.
*/
- typename std::vector<typename std::vector<Tensor<2,dim> > > shape_2nd_derivatives;
+ GradGradVector shape_2nd_derivatives;
/**
- * Store an array of weights times the
- * Jacobi determinant at the quadrature
- * points. This function is reset each time
- * @p{reinit} is called. The Jacobi determinant
- * is actually the reciprocal value of the
- * Jacobi matrices stored in this class,
- * see the general documentation of this
- * class for more information.
+ * Store an array of weights
+ * times the Jacobi determinant
+ * at the quadrature points. This
+ * function is reset each time
+ * @p{reinit} is called. The
+ * Jacobi determinant is actually
+ * the reciprocal value of the
+ * Jacobi matrices stored in this
+ * class, see the general
+ * documentation of this class
+ * for more information.
*/
std::vector<double> JxW_values;
*/
typename std::vector<Tensor<1,dim> > boundary_forms;
+ /**
+ * Indicate the first row which a
+ * given shape function occupies
+ * in the @p{ShapeVector},
+ * @p{ShapeGradient}, etc
+ * 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
+ * @p{FiniteElement::get_nonzero_components()}
+ * function, and the rows are in
+ * ascending order exactly those
+ * non-zero components.
+ */
+ std::vector<unsigned int> shape_function_to_row_table;
+
/**
* Original update flags handed
* to the constructor of
* called by the constructor and @p{reinit} functions of
* @p{FEValues*}, respectively.
*
- * \subsection{General usage}
+ * @sect3{General usage}
*
* Usually, an object of @p{FEValues*} is used in integration loops
* over all cells of a triangulation. To take full advantage of the
* like a magnifying glass being used to look at one item after the
* other. A typical piece of code looks like this:
*
- * \begin{verbatim}
+ * @begin{verbatim}
* FEValues values (mapping, finite_element, quadrature, flags);
* for (cell = dof_handler.begin_active();
* cell != dof_handler.end();
* values.reinit(cell);
* ...
* }
- * \end{verbatim}
+ * @end{verbatim}
*
*
- * \subsection{Member functions}
+ * @sect3{Member functions}
*
* The functions of this class fall into different cathegories:
- * \begin{itemize}
- * \item @p{shape_value}, @p{shape_grad}, etc: return one of the values
+ * @begin{itemize}
+ * @item @p{shape_value}, @p{shape_grad}, etc: return one of the values
* of this object at a time. These functions are inlined, so this
- * is the suggested access to all finite element values. There should be
- * no loss in performance with an optimizing compiler.
+ * is the suggested access to all finite element values. There
+ * should be no loss in performance with an optimizing compiler. If
+ * the finite element is vector valued, then these functions return
+ * the only non-zero component of the requested shape
+ * function. However, some finite elements have shape functions
+ * that have more than one non-zero component (we call them
+ * non-"primitive"), and in this case this set of functions will
+ * throw an exception since they cannot generate a useful
+ * result. Rather, use the next set of functions.
+ *
+ * @item @p{shape_value_component}, @p{shape_grad_component}, etc:
+ * This is the same set of functions as above, except that for vector
+ * valued finite elements they return only one vector component. This
+ * is useful for elements of which shape functions have more than one
+ * non-zero component, since then the above functions cannot be used,
+ * and you have to walk over all (or only the non-zero) components of
+ * the shape function using this set of functions.
*
- * \item @p{get_function_values}, @p{get_function_grads}, @p{...}:
+ * @item @p{get_function_values}, @p{get_function_grads}, @p{...}:
* Compute a finite element function or its derivative
* in quadrature points.
*
- * \item @p{reinit}: initialize the @p{FEValues} object for a certain cell.
+ * @item @p{reinit}: initialize the @p{FEValues} object for a certain cell.
* This function is not in the present class but only in the derived
* classes and has a variable call syntax.
* See the docs for the derived classes for more information.
- * \end{itemize}
+ * @end{itemize}
*
*
- * \subsection{UpdateFlags}
+ * @sect3{UpdateFlags}
*
* The @ref{UpdateFlags} object handed to the constructor is used to
* determine, which of the data fields to compute. This way, it is
const unsigned int n_quadrature_points;
/**
- * Number of shape functions
- * per cell. If we use this base class
- * to evaluate a finite element on
- * faces of cells, this is still the
- * number of degrees of freedom per
- * cell, not per face.
+ * Number of shape functions per
+ * cell. If we use this base
+ * class to evaluate a finite
+ * element on faces of cells,
+ * this is still the number of
+ * degrees of freedom per cell,
+ * not per face.
*/
const unsigned int dofs_per_cell;
~FEValuesBase ();
/**
- * Value of the @p{function_no}th shape
- * function at the @p{point_no}th quadrature
+ * Value of the @p{function_no}th
+ * shape function at the
+ * @p{point_no}th quadrature
* point on the cell, face or
* subface selected the last time
* the @p{reinit} function of the
* derived class was called.
+ *
+ * If the shape function is
+ * vector-valued, then this
+ * returns the only non-zero
+ * component. If the shape
+ * function has more than one
+ * non-zero component (i.e. it is
+ * not primitive), then throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}. In
+ * that case, use the
+ * @ref{shape_value_component}
+ * function.
*/
double shape_value (const unsigned int function_no,
const unsigned int point_no) const;
/**
- * This function is now
- * deprecated and will be removed
- * in the next release of the
- * deal.II library.
- *
- * Pointer to the matrix holding
- * all values of shape functions
- * at all integration points, on
- * the present cell, face or
- * subface selected the last time
- * the @p{reinit} function of the
- * derived class was called. For
- * the format of this matrix, see
- * the documentation for the
- * matrix itself.
- */
- const typename FEValuesData<dim>::ShapeVector & get_shape_values () const;
+ * Return one vector component of
+ * the value of a shape function
+ * at a quadrature point. If the
+ * finite element is scalar, then
+ * only component zero is allowed
+ * and the return value equals
+ * that of the @p{shape_value}
+ * function. If the finite
+ * element is vector valued but
+ * all shape functions are
+ * primitive (i.e. they are
+ * non-zero in only one
+ * component), then the value
+ * returned by @p{shape_value}
+ * equals that of this function
+ * for exactly one
+ * component. This function is
+ * therefore only of greater
+ * interest if the shape function
+ * is not primitive, but then it
+ * is necessary since the other
+ * function cannot be used.
+ */
+ double shape_value_component (const unsigned int function_no,
+ const unsigned int point_no,
+ const unsigned int component) const;
/**
* Values of the finite
* gradient's value is returned,
* there should be no major
* performance drawback.
- */
- const Tensor<1,dim> & shape_grad (const unsigned int function,
- const unsigned int quadrature_point) const;
-
- /**
- * This function is now
- * deprecated and will be removed
- * in the next release of the
- * deal.II library.
*
- * Pointer to the matrix holding
- * all gradients of shape
- * functions at all integration
- * points, on the present cell.
- * For the format of this matrix,
- * see the documentation for the
- * matrix itself.
- */
- const typename FEValuesData<dim>::GradientVector & get_shape_grads () const;
+ * If the shape function is
+ * vector-valued, then this
+ * returns the only non-zero
+ * component. If the shape
+ * function has more than one
+ * non-zero component (i.e. it is
+ * not primitive), then throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}. In
+ * that case, use the
+ * @ref{shape_grad_component}
+ * function.
+ */
+ const Tensor<1,dim> &
+ shape_grad (const unsigned int function,
+ const unsigned int quadrature_point) const;
+
+ /**
+ * Return one vector component of
+ * the gradient of a shape function
+ * at a quadrature point. If the
+ * finite element is scalar, then
+ * only component zero is allowed
+ * and the return value equals
+ * that of the @p{shape_grad}
+ * function. If the finite
+ * element is vector valued but
+ * all shape functions are
+ * primitive (i.e. they are
+ * non-zero in only one
+ * component), then the value
+ * returned by @p{shape_grad}
+ * equals that of this function
+ * for exactly one
+ * component. This function is
+ * therefore only of greater
+ * interest if the shape function
+ * is not primitive, but then it
+ * is necessary since the other
+ * function cannot be used.
+ */
+ Tensor<1,dim>
+ shape_grad_component (const unsigned int function_no,
+ const unsigned int point_no,
+ const unsigned int component) const;
/**
* Gradients of the finite
* only a reference to the
* derivative values is returned,
* there should be no major
- * performance drawback. */
- const Tensor<2,dim> & shape_2nd_derivative (const unsigned int function_no,
- const unsigned int point_no) const;
-
- /**
- * This function is now
- * deprecated and will be removed
- * in the next release of the
- * deal.II library.
+ * performance drawback.
*
- * Pointer to the
- * matrix holding all 2nd
- * derivatives of shape functions
- * at all integration points, on
- * the present cell. For the
- * format of this matrix, see the
- * documentation for the matrix
- * itself.
- */
- const typename std::vector<typename std::vector<Tensor<2,dim> > > & get_shape_2nd_derivatives () const;
+ * If the shape function is
+ * vector-valued, then this
+ * returns the only non-zero
+ * component. If the shape
+ * function has more than one
+ * non-zero component (i.e. it is
+ * not primitive), then throw an
+ * exception of type
+ * @p{ExcShapeFunctionNotPrimitive}. In
+ * that case, use the
+ * @ref{shape_grad_grad_component}
+ * function.
+ */
+ const Tensor<2,dim> &
+ shape_2nd_derivative (const unsigned int function_no,
+ const unsigned int point_no) const;
+
+
+ /**
+ * Return one vector component of
+ * the gradient of a shape
+ * function at a quadrature
+ * point. If the finite element
+ * is scalar, then only component
+ * zero is allowed and the return
+ * value equals that of the
+ * @p{shape_2nd_derivative}
+ * function. If the finite
+ * element is vector valued but
+ * all shape functions are
+ * primitive (i.e. they are
+ * non-zero in only one
+ * component), then the value
+ * returned by
+ * @p{shape_2nd_derivative}
+ * equals that of this function
+ * for exactly one
+ * component. This function is
+ * therefore only of greater
+ * interest if the shape function
+ * is not primitive, but then it
+ * is necessary since the other
+ * function cannot be used.
+ */
+ Tensor<2,dim>
+ shape_2nd_derivative_component (const unsigned int function_no,
+ const unsigned int point_no,
+ const unsigned int component) const;
/**
* Tensor of second derivatives
* Exception
*/
DeclException0 (ExcFEDontMatch);
+ /**
+ * Exception
+ */
+ DeclException1 (ExcShapeFunctionNotPrimitive,
+ int,
+ << "The shape function with index " << arg1
+ << " is not primitive, i.e. it is vector-valued and "
+ << "has more than one non-zero vector component. This "
+ << "function cannot be called for these shape functions. "
+ << "Maybe you want to use the same function with the "
+ << "_component suffix?");
+ /**
+ * Exception
+ */
+ DeclException0 (ExcFENotPrimitive);
protected:
/**
FEValuesBase<dim>::shape_value (const unsigned int i,
const unsigned int j) const
{
- Assert (update_flags & update_values, ExcAccessToUninitializedField());
- return shape_values(i,j);
+ Assert (this->update_flags & update_values,
+ ExcAccessToUninitializedField());
+ Assert (fe->is_primitive (i),
+ ExcShapeFunctionNotPrimitive(i));
+
+ // if the entire FE is primitive,
+ // then we can take a short-cut:
+ 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(shape_function_to_row_table[i], j);
}
+
+template <int dim>
+inline
+double
+FEValuesBase<dim>::shape_value_component (const unsigned int i,
+ const unsigned int j,
+ const unsigned int component) const
+{
+ Assert (this->update_flags & update_values,
+ ExcAccessToUninitializedField());
+ Assert (component < fe->n_components(),
+ ExcIndexRange(component, 0, fe->n_components()));
+
+ // if this particulat 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(i,j);
+ else
+ return 0;
+ }
+ else
+ {
+ // no, this shape function is
+ // not primitive. then we have
+ // to loop over its components
+ // and 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 = (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);
+ };
+}
+
+
+
template <int dim>
inline
const Tensor<1,dim> &
FEValuesBase<dim>::shape_grad (const unsigned int i,
const unsigned int j) const
{
- Assert (i<shape_gradients.size(),
+ Assert (i<this->shape_gradients.size(),
ExcIndexRange (i, 0, shape_gradients.size()));
Assert (j<shape_gradients[i].size(),
ExcIndexRange (j, 0, shape_gradients[i].size()));
- Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
-
- return shape_gradients[i][j];
+ Assert (this->update_flags & update_gradients,
+ ExcAccessToUninitializedField());
+ Assert (fe->is_primitive (i),
+ ExcShapeFunctionNotPrimitive(i));
+
+ // if the entire FE is primitive,
+ // then we can take a short-cut:
+ 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[shape_function_to_row_table[i]][j];
};
+template <int dim>
+inline
+Tensor<1,dim>
+FEValuesBase<dim>::shape_grad_component (const unsigned int i,
+ const unsigned int j,
+ const unsigned int component) const
+{
+ Assert (this->update_flags & update_values,
+ ExcAccessToUninitializedField());
+ Assert (component < fe->n_components(),
+ ExcIndexRange(component, 0, fe->n_components()));
+
+ // if this particulat 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[i][j];
+ else
+ return Tensor<1,dim>();
+ }
+ else
+ {
+ // no, this shape function is
+ // not primitive. then we have
+ // to loop over its components
+ // and 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,dim>();
+
+ // 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 = (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];
+ };
+}
+
+
+
template <int dim>
inline
const Tensor<2,dim> &
FEValuesBase<dim>::shape_2nd_derivative (const unsigned int i,
const unsigned int j) const
{
- Assert (i<shape_2nd_derivatives.size(),
+ Assert (i<this->shape_2nd_derivatives.size(),
ExcIndexRange (i, 0, shape_2nd_derivatives.size()));
- Assert (j<shape_2nd_derivatives[i].size(),
+ Assert (j<this->shape_2nd_derivatives[i].size(),
ExcIndexRange (j, 0, shape_2nd_derivatives[i].size()));
- Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
+ Assert (this->update_flags & update_second_derivatives,
+ ExcAccessToUninitializedField());
+ Assert (fe->is_primitive (i),
+ ExcShapeFunctionNotPrimitive(i));
+
+ // if the entire FE is primitive,
+ // then we can take a short-cut:
+ if (fe->is_primitive())
+ return this->shape_2nd_derivatives[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_2nd_derivatives[shape_function_to_row_table[i]][j];
+}
+
+
- return shape_2nd_derivatives[i][j];
+template <int dim>
+inline
+Tensor<2,dim>
+FEValuesBase<dim>::shape_2nd_derivative_component (const unsigned int i,
+ const unsigned int j,
+ const unsigned int component) const
+{
+ Assert (this->update_flags & update_values,
+ ExcAccessToUninitializedField());
+ Assert (component < fe->n_components(),
+ ExcIndexRange(component, 0, fe->n_components()));
+
+ // if this particulat 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_2nd_derivatives[i][j];
+ else
+ return Tensor<2,dim>();
+ }
+ else
+ {
+ // no, this shape function is
+ // not primitive. then we have
+ // to loop over its components
+ // and 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,dim>();
+
+ // 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 = (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_2nd_derivatives[row][j];
+ };
}
const Point<dim> &
FEFaceValuesBase<dim>::normal_vector (const unsigned int i) const
{
- Assert (i<normal_vectors.size(), ExcIndexRange(i, 0, normal_vectors.size()));
- Assert (update_flags & update_normal_vectors,
+ Assert (i<this->normal_vectors.size(), ExcIndexRange(i, 0, normal_vectors.size()));
+ Assert (this->update_flags & update_normal_vectors,
FEValuesBase<dim>::ExcAccessToUninitializedField());
- return normal_vectors[i];
+ return this->normal_vectors[i];
};
template <int dim> class FEFaceValues;
template <int dim> class FESubfaceValues;
+
/**
- * Abstract basis class for mapping classes.
+ * Abstract base class for mapping classes.
*
* The interface for filling the tables of @ref{FEValues} is provided.
* Everything else has to happen in derived classes.
*
* The following paragraph applies to the implementation of
* @ref{FEValues}. Usage of the class is as follows: first, call the
- * functionss @p{update_once} and @p{update_each} with the update
+ * functions @p{update_once} and @p{update_each} with the update
* flags you need. This includes the flags needed by the
* @ref{FiniteElement}. Then call @p{get_*_data} and with the or'd
* results. This will initialize and return some internal data
const Point<dim> &p) const=0;
/**
- * Class for internal data of finite
- * element and mapping objects.
+ * Base class for internal data
+ * of finite element and mapping
+ * objects. The internal
+ * mechanism is that upon
+ * construction of a @p{FEValues}
+ * objects, it asks the mapping
+ * and finite element classes
+ * that are to be used to
+ * allocate memory for their own
+ * purpose in which they may
+ * store data that only needs to
+ * be computed once. For example,
+ * most finite elements will
+ * store the values of the shape
+ * functions at the quadrature
+ * points in this object, since
+ * they do not change from cell
+ * to cell and only need to be
+ * computed once. Since different
+ * @p{FEValues} objects using
+ * different quadrature rules
+ * might access the same finite
+ * element object at the same
+ * time, it is necessary to
+ * create one such object per
+ * @p{FEValues} object. Ownership
+ * of this object is then
+ * transferred to the
+ * @p{FEValues} object, but a
+ * pointer to this object is
+ * passed to the finite element
+ * object every time it shall
+ * compute some data so that it
+ * has access to the precomputed
+ * values stored there.
*/
class InternalDataBase: public Subscriptor
{
// the restriction of the fine
// grid:
const unsigned int coarse_dofs_per_cell_component
- = coarse_fe.base_element(coarse_fe.component_to_base(coarse_component)).dofs_per_cell;
+ = coarse_fe.base_element(coarse_fe.component_to_base(coarse_component).first).dofs_per_cell;
// Try to find out whether the
ExcInvalidComponent (fine_component, fine_fe.n_components()));
// check whether respective finite
// elements are equal
- Assert (coarse_fe.base_element (coarse_fe.component_to_base(coarse_component))
+ Assert (coarse_fe.base_element (coarse_fe.component_to_base(coarse_component).first)
==
- fine_fe.base_element (fine_fe.component_to_base(fine_component)),
+ fine_fe.base_element (fine_fe.component_to_base(fine_component).first),
ExcFiniteElementsDontMatch());
#ifdef DEBUG
#include <dofs/dof_accessor.h>
#include <grid/tria_boundary.h>
+#include <algorithm>
+#include <functional>
+
+
// if necessary try to work around a bug in the IBM xlC compiler
#ifdef XLC_WORK_AROUND_STD_BUG
using namespace std;
template <int dim>
FiniteElementBase<dim>::FiniteElementBase (const FiniteElementData<dim> &fe_data,
- const std::vector<bool> &restriction_is_additive_flags)
+ const std::vector<bool> &restriction_is_additive_flags,
+ const std::vector<std::vector<bool> > &nonzero_components)
:
FiniteElementData<dim> (fe_data),
system_to_component_table(dofs_per_cell),
face_system_to_component_table(dofs_per_face),
+ system_to_base_table(dofs_per_cell),
+ face_system_to_base_table(dofs_per_face),
component_to_system_table(components, std::vector<unsigned>(dofs_per_cell)),
- face_component_to_system_table(components, std::vector<unsigned>(dofs_per_face)),
- component_to_base_table (components, 0),
- restriction_is_additive_flags(restriction_is_additive_flags)
+ face_component_to_system_table(components, std::vector<unsigned>(dofs_per_face)),
+ component_to_base_table (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,
- ExcDimensionMismatch(restriction_is_additive_flags.size(),fe_data.components));
-
+ Assert (restriction_is_additive_flags.size()==fe_data.components,
+ ExcDimensionMismatch(restriction_is_additive_flags.size(),
+ fe_data.components));
+ Assert (nonzero_components.size() == dofs_per_cell,
+ ExcInternalError());
+ for (unsigned int i=0; i<nonzero_components.size(); ++i)
+ {
+ Assert (nonzero_components[i].size() == n_components(),
+ ExcInternalError());
+ Assert (std::count (nonzero_components[i].begin(),
+ nonzero_components[i].end(),
+ true)
+ >= 1,
+ ExcInternalError());
+ Assert (n_nonzero_components_table[i] >= 1,
+ ExcInternalError());
+ Assert (n_nonzero_components_table[i] <= n_components(),
+ ExcInternalError());
+ };
+
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
{
restriction[i].reinit (dofs_per_cell, dofs_per_cell);
for (unsigned int j=0 ; j<dofs_per_cell ; ++j)
{
system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
+ system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);
component_to_system_table[0][j] = j;
}
for (unsigned int j=0 ; j<dofs_per_face ; ++j)
{
face_system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
+ face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);
face_component_to_system_table[0][j] = j;
}
};
+template <int dim>
+double
+FiniteElementBase<dim>::shape_value_component (const unsigned int,
+ const Point<dim> &,
+ const unsigned int) const
+{
+ AssertThrow(false, ExcUnitShapeValuesDoNotExist());
+ return 0.;
+}
+
+
+
template <int dim>
Tensor<1,dim>
FiniteElementBase<dim>::shape_grad (const unsigned int,
+template <int dim>
+Tensor<1,dim>
+FiniteElementBase<dim>::shape_grad_component (const unsigned int,
+ const Point<dim> &,
+ const unsigned int) const
+{
+ AssertThrow(false, ExcUnitShapeValuesDoNotExist());
+ return Tensor<1,dim> ();
+}
+
+
+
template <int dim>
Tensor<2,dim>
FiniteElementBase<dim>::shape_grad_grad (const unsigned int,
}
+
+template <int dim>
+Tensor<2,dim>
+FiniteElementBase<dim>::shape_grad_grad_component (const unsigned int,
+ const Point<dim> &,
+ const unsigned int) const
+{
+ AssertThrow(false, ExcUnitShapeValuesDoNotExist());
+ return Tensor<2,dim> ();
+}
+
+
template <int dim>
const FullMatrix<double> &
FiniteElementBase<dim>::restrict (const unsigned int child) const
MemoryConsumption::memory_consumption (interface_constraints) +
MemoryConsumption::memory_consumption (system_to_component_table) +
MemoryConsumption::memory_consumption (face_system_to_component_table) +
+ MemoryConsumption::memory_consumption (system_to_base_table) +
+ MemoryConsumption::memory_consumption (face_system_to_base_table) +
MemoryConsumption::memory_consumption (component_to_system_table) +
MemoryConsumption::memory_consumption (face_component_to_system_table) +
MemoryConsumption::memory_consumption (component_to_base_table) +
- MemoryConsumption::memory_consumption (restriction_is_additive_flags));
+ MemoryConsumption::memory_consumption (restriction_is_additive_flags) +
+ MemoryConsumption::memory_consumption (nonzero_components) +
+ MemoryConsumption::memory_consumption (n_nonzero_components_table));
};
}
+
+template <int dim>
+std::vector<unsigned int>
+FiniteElementBase<dim>::
+compute_n_nonzero_components (const std::vector<std::vector<bool> > &nonzero_components)
+{
+ std::vector<unsigned int> retval (nonzero_components.size());
+ for (unsigned int i=0; i<nonzero_components.size(); ++i)
+ retval[i] = std::count (nonzero_components[i].begin(),
+ nonzero_components[i].end(),
+ true);
+ return retval;
+};
+
+
+
/*------------------------------- FiniteElement ----------------------*/
template <int dim>
FiniteElement<dim>::FiniteElement (const FiniteElementData<dim> &fe_data,
- const std::vector<bool> &restriction_is_additive_flags) :
+ const std::vector<bool> &restriction_is_additive_flags,
+ const std::vector<std::vector<bool> > &nonzero_components) :
FiniteElementBase<dim> (fe_data,
- restriction_is_additive_flags)
+ restriction_is_additive_flags,
+ nonzero_components)
{}
//----------------------------------------------------------------
#include <base/quadrature.h>
-#include <base/polynomial.h>
-#include <base/polynomial_space.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
FE_DGP<dim>::FE_DGP (unsigned int degree)
:
FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
- std::vector<bool>(1,true)),
+ std::vector<bool>(1,true),
+ std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+ std::vector<bool>(1,true))),
degree(degree),
- poly(0)
+ polynomial_space (Legendre<double>::generate_complete_basis(degree))
{
- // create array of Legendre polynomials
- std::vector<Legendre<double> > v;
- for (unsigned int i=0;i<=degree;++i)
- v.push_back(Legendre<double>(i));
- poly = new PolynomialSpace<dim> (v);
-
// if defined, copy over matrices
// from precomputed arrays
if ((degree < Matrices::n_embedding_matrices) &&
else
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
prolongation[i].reinit(0,0);
+
// // same as above: copy over matrix
// // from predefined values and
// // generate all others by rotation
template <int dim>
-FE_DGP<dim>::~FE_DGP ()
+FiniteElement<dim> *
+FE_DGP<dim>::clone() const
{
- // delete poly member and set it to
- // zero to prevent accidental use
- delete poly;
- poly = 0;
+ return new FE_DGP<dim>(degree);
}
template <int dim>
-FiniteElement<dim> *
-FE_DGP<dim>::clone() const
+double
+FE_DGP<dim>::shape_value (const unsigned int i,
+ const Point<dim> &p) const
{
- return new FE_DGP<dim>(degree);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_value(i, p);
}
template <int dim>
double
-FE_DGP<dim>::shape_value (const unsigned int i,
- const Point<dim> &p) const
+FE_DGP<dim>::shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
{
- return poly->compute_value(i, p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_value(i, p);
}
FE_DGP<dim>::shape_grad (const unsigned int i,
const Point<dim> &p) const
{
- return poly->compute_grad(i, p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_grad(i, p);
+}
+
+
+template <int dim>
+Tensor<1,dim>
+FE_DGP<dim>::shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_grad(i, p);
}
FE_DGP<dim>::shape_grad_grad (const unsigned int i,
const Point<dim> &p) const
{
- return poly->compute_grad_grad(i, p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_grad_grad(i, p);
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_DGP<dim>::shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_grad_grad(i, p);
}
UpdateFlags
FE_DGP<dim>::update_once (const UpdateFlags flags) const
{
- UpdateFlags out = update_default;
-
- if (flags & update_values)
- out |= update_values;
-
- return out;
+ // for this kind of elements, only
+ // the values can be precomputed
+ // once and for all. set this flag
+ // if the values are requested at
+ // all
+ return (update_default | (flags & update_values));
}
const Mapping<dim> &mapping,
const Quadrature<dim> &quadrature) const
{
+ // generate a new data object
InternalData* data = new InternalData;
- std::vector<double> values(0);
- std::vector<Tensor<1,dim> > grads(0);
- std::vector<Tensor<2,dim> > grad_grads(0);
-
// check what needs to be
// initialized only once and what
// on every cell/face/subface we
const UpdateFlags flags(data->update_flags);
const unsigned int n_q_points = quadrature.n_quadrature_points;
+
+ // have some scratch arrays
+ std::vector<double> values(0);
+ std::vector<Tensor<1,dim> > grads(0);
+ std::vector<Tensor<2,dim> > grad_grads(0);
+ // initialize fields only if really
+ // necessary. otherwise, don't
+ // allocate memory
if (flags & update_values)
{
values.resize (dofs_per_cell);
if (flags & update_second_derivatives)
data->initialize_2nd (this, mapping, quadrature);
-
+ // next already fill those fields
+ // of which we have information by
+ // now. note that the shape
+ // gradients are only those on the
+ // unit cell, and need to be
+ // transformed when visiting an
+ // actual cell
if (flags & (update_values | update_gradients))
for (unsigned int i=0; i<n_q_points; ++i)
{
- poly->compute(quadrature.point(i), values, grads, grad_grads);
+ polynomial_space.compute(quadrature.point(i),
+ values, grads, grad_grads);
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
if (flags & update_values)
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
- for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
- if (flags & update_values)
+ if (flags & update_values)
+ for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
data.shape_values(k,i) = fe_data.shape_values[k][i];
if (flags & update_gradients)
FE_DGQ<dim>::FE_DGQ (unsigned int degree)
:
FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
- std::vector<bool>(1,true)),
+ std::vector<bool>(1,true),
+ std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+ std::vector<bool>(1,true))),
degree(degree),
- poly(0)
+ polynomial_space (LagrangeEquidistant::generate_complete_basis(degree))
{
- if (degree==0)
- {
- // create constant polynomial
- std::vector<Polynomial<double> >
- v(1, Polynomial<double> (std::vector<double> (1,1.)));
- poly = new TensorProductPolynomials<dim> (v);
- }
- else
- {
- // create array of Lagrange polynomials
- std::vector<LagrangeEquidistant> v;
- for (unsigned int i=0;i<=degree;++i)
- v.push_back(LagrangeEquidistant(degree,i));
- poly = new TensorProductPolynomials<dim> (v);
- }
-
// generate permutation/rotation
// index sets to generate some
// matrices from others
template <int dim>
-FE_DGQ<dim>::~FE_DGQ ()
+FiniteElement<dim> *
+FE_DGQ<dim>::clone() const
{
- // delete poly member and set it to
- // zero to prevent accidental use
- delete poly;
- poly = 0;
+ return new FE_DGQ<dim>(degree);
}
template <int dim>
-FiniteElement<dim> *
-FE_DGQ<dim>::clone() const
+double
+FE_DGQ<dim>::shape_value (const unsigned int i,
+ const Point<dim> &p) const
{
- return new FE_DGQ<dim>(degree);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_value(i, p);
}
template <int dim>
double
-FE_DGQ<dim>::shape_value (const unsigned int i,
- const Point<dim> &p) const
+FE_DGQ<dim>::shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
{
- return poly->compute_value(i, p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_value(i, p);
}
FE_DGQ<dim>::shape_grad (const unsigned int i,
const Point<dim> &p) const
{
- return poly->compute_grad(i, p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_grad(i, p);
+}
+
+
+template <int dim>
+Tensor<1,dim>
+FE_DGQ<dim>::shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_grad(i, p);
}
FE_DGQ<dim>::shape_grad_grad (const unsigned int i,
const Point<dim> &p) const
{
- return poly->compute_grad_grad(i, p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_grad_grad(i, p);
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_DGQ<dim>::shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_grad_grad(i, p);
}
UpdateFlags
FE_DGQ<dim>::update_once (const UpdateFlags flags) const
{
- UpdateFlags out = update_default;
-
- if (flags & update_values)
- out |= update_values;
-
- return out;
+ // for this kind of elements, only
+ // the values can be precomputed
+ // once and for all. set this flag
+ // if the values are requested at
+ // all
+ return (update_default | (flags & update_values));
}
const Mapping<dim> &mapping,
const Quadrature<dim> &quadrature) const
{
+ // generate a new data object
InternalData* data = new InternalData;
- std::vector<double> values(0);
- std::vector<Tensor<1,dim> > grads(0);
- std::vector<Tensor<2,dim> > grad_grads(0);
// check what needs to be
// initialized only once and what
const UpdateFlags flags(data->update_flags);
const unsigned int n_q_points = quadrature.n_quadrature_points;
-
+
+ // have some scratch arrays
+ std::vector<double> values(0);
+ std::vector<Tensor<1,dim> > grads(0);
+ std::vector<Tensor<2,dim> > grad_grads(0);
+
+ // initialize fields only if really
+ // necessary. otherwise, don't
+ // allocate memory
if (flags & update_values)
{
values.resize (dofs_per_cell);
if (flags & update_second_derivatives)
data->initialize_2nd (this, mapping, quadrature);
-
+ // next already fill those fields
+ // of which we have information by
+ // now. note that the shape
+ // gradients are only those on the
+ // unit cell, and need to be
+ // transformed when visiting an
+ // actual cell
if (flags & (update_values | update_gradients))
for (unsigned int i=0; i<n_q_points; ++i)
{
- poly->compute(quadrature.point(i), values, grads, grad_grads);
+ polynomial_space.compute(quadrature.point(i),
+ values, grads, grad_grads);
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
if (flags & update_values)
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
- for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
- if (flags & update_values)
+ if (flags & update_values)
+ for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
data.shape_values(k,i) = fe_data.shape_values[k][i];
if (flags & update_gradients)
FE_DGQ<dim>::has_support_on_face (const unsigned int shape_index,
const unsigned int face_index) const
{
+ Assert (shape_index < dofs_per_cell,
+ ExcIndexRange (shape_index, 0, dofs_per_cell));
+ Assert (face_index < GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
unsigned int n = degree+1;
unsigned int n2 = n*n;
--- /dev/null
+//----------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002 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.
+//
+//----------------------------------------------------------------
+
+#include <base/quadrature.h>
+#include <base/polynomial.h>
+#include <base/tensor_product_polynomials.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+#include <fe/mapping.h>
+#include <fe/fe_nedelec.h>
+#include <fe/fe_values.h>
+
+
+template <int dim>
+FE_Nedelec<dim>::FE_Nedelec (const unsigned int degree)
+ :
+ FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),
+ dim),
+//TODO: I'd think this element is actually additive in the restriction
+ std::vector<bool> (dim,false),
+ std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),dim).dofs_per_cell,
+ std::vector<bool>(dim,true))),
+ degree(degree)
+{
+ Assert (dim >= 2, ExcNotUsefulInThisDimension());
+
+ // copy constraint matrices if they
+ // are defined. otherwise set them
+ // to invalid size
+ if (degree<Matrices::n_constraint_matrices+1)
+ interface_constraints.fill (Matrices::constraint_matrices[degree-1]);
+ else
+ interface_constraints.reinit(0,0);
+
+ // next copy over embedding
+ // matrices if they are defined
+ if ((degree < Matrices::n_embedding_matrices+1) &&
+ (Matrices::embedding[degree-1][0] != 0))
+ for (unsigned int c=0; c<GeometryInfo<dim>::children_per_cell; ++c)
+ prolongation[c].fill (Matrices::embedding[degree-1][c]);
+ else
+ for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
+ prolongation[i].reinit(0,0);
+
+ // then fill restriction
+ // matrices. they are hardcoded for
+ // the first few elements
+ switch (dim)
+ {
+ case 2: // 2d
+ {
+ switch (degree)
+ {
+ case 1:
+ {
+ // DoF on bottom line
+ // of coarse cell will
+ // be mean value of
+ // bottom DoFs on the
+ // two adjacent child
+ // cells
+ restriction[0](0,0) = 0.5;
+ restriction[1](0,0) = 0.5;
+ // same for other DoFs
+ restriction[1](1,1) = 0.5;
+ restriction[2](1,1) = 0.5;
+
+ restriction[2](2,2) = 0.5;
+ restriction[3](2,2) = 0.5;
+
+ restriction[3](3,3) = 0.5;
+ restriction[0](3,3) = 0.5;
+
+ break;
+ };
+
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ restriction[i].reinit(0,0);
+ };
+ };
+
+ break;
+ };
+
+
+ case 3: // 3d
+ {
+ switch (degree)
+ {
+ case 1:
+ {
+ // same principle as in
+ // 2d
+ restriction[0](0,0) = 0.5;
+ restriction[1](0,0) = 0.5;
+
+ restriction[1](1,1) = 0.5;
+ restriction[2](1,1) = 0.5;
+
+ restriction[2](2,2) = 0.5;
+ restriction[3](2,2) = 0.5;
+
+ restriction[3](3,3) = 0.5;
+ restriction[0](3,3) = 0.5;
+
+ restriction[4](4,4) = 0.5;
+ restriction[5](4,4) = 0.5;
+
+ restriction[5](5,5) = 0.5;
+ restriction[6](5,5) = 0.5;
+
+ restriction[6](6,6) = 0.5;
+ restriction[7](6,6) = 0.5;
+
+ restriction[7](7,7) = 0.5;
+ restriction[4](7,7) = 0.5;
+
+
+ restriction[1](8,8) = 0.5;
+ restriction[5](8,8) = 0.5;
+
+ restriction[2](9,9) = 0.5;
+ restriction[6](9,9) = 0.5;
+
+ restriction[3](10,10) = 0.5;
+ restriction[7](10,10) = 0.5;
+
+ restriction[0](11,11) = 0.5;
+ restriction[5](11,11) = 0.5;
+
+ break;
+ };
+
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ restriction[i].reinit(0,0);
+ };
+ };
+
+ break;
+ };
+
+ default:
+ Assert (false,ExcNotImplemented());
+ }
+
+ // finally fill in support points
+ // on cell and face
+ initialize_unit_support_points ();
+ initialize_unit_face_support_points ();
+};
+
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_Nedelec<dim>::clone() const
+{
+ return new FE_Nedelec<dim>(degree);
+}
+
+
+
+template <int dim>
+double
+FE_Nedelec<dim>::shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component < dim, ExcIndexRange (component, 0, dim));
+
+ switch (dim)
+ {
+ case 2: // 2D
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
+ {
+ switch (i)
+ {
+ // (1-y, 0)
+ case 0: return (component == 0 ? 1-p(1) : 0);
+ // (0,x)
+ case 1: return (component == 0 ? 0 : p(0));
+ // (y, 0)
+ case 2: return (component == 0 ? p(1) : 0);
+ // (0, 1-x)
+ case 3: return (component == 0 ? 0 : 1-p(0));
+
+ // there are
+ // only four
+ // shape
+ // functions!?
+ default:
+ Assert (false, ExcInternalError());
+ return 0;
+ };
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+
+ case 3: // 3D
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
+ {
+ // note that the
+ // degrees of freedom
+ // on opposite faces
+ // have a common vector
+ // direction, so simply
+ // that a little. these
+ // directions are:
+ //
+ // for lines 0, 2, 4, 6:
+ // (1,0,0)
+ // for lines 1, 3, 5, 7:
+ // (0,0,1)
+ // for lines 8, 9, 10, 11:
+ // (0,1,0)
+ //
+ // thus, sort out all
+ // those cases where
+ // the component is
+ // zero anyway, and
+ // only otherwise
+ // compute the
+ // spatially dependent
+ // part which is then
+ // also the return
+ // value
+ if (((i<8) && (((i%2==0) && (component!=0)) ||
+ ((i%2==1) && (component!=2)))) ||
+ ((i>=8) && (component != 1)))
+ return 0;
+
+ // now we know that the
+ // only non-zero
+ // component is
+ // requested:
+//TODO[Anna]: check
+ const double x = p(0),
+ y = p(1),
+ z = p(2);
+ switch (i)
+ {
+ case 0: return (1-y)*(1-z);
+ case 2: return (1-y)*z;
+ case 1: return x*(1-y);
+ case 3: return (1-x)*(1-y);
+
+ case 4: return y*(1-z);
+ case 6: return y*z;
+ case 5: return x*y;
+ case 7: return (1-x)*y;
+
+ case 8: return (1-x)*(1-z);
+ case 9: return x*(1-z);
+ case 10: return x*z;
+ case 11: return (1-x)*z;
+ default:
+ Assert (false, ExcInternalError());
+ return 0;
+ };
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+
+ // presently no other space
+ // dimension implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ return 0;
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FE_Nedelec<dim>::shape_grad_component (const unsigned int i,
+ const Point<dim> &/*p*/,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component < dim, ExcIndexRange (component, 0, dim));
+
+ switch (dim)
+ {
+ case 2: // 2D
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
+ {
+ // on the unit cell,
+ // the gradients of
+ // these shape
+ // functions are
+ // constant, so we pack
+ // them into a table
+ // for simpler lookup
+ //
+ // the format is: first
+ // index=shape function
+ // number; second
+ // index=vector
+ // component, thrid
+ // index=component
+ // within gradient
+//TODO[Anna]: check
+ static const double unit_gradients[4][2][2]
+ = { { {0.,-1.}, {0.,0.} },
+ { {0.,0.}, {1.,0.} },
+ { {0.,+1.}, {0.,0.} },
+ { {0.,0.}, {-1.,0.} } };
+ return Tensor<1,dim>(unit_gradients[i][component]);
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+
+ // presently no other space
+ // dimension implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ return Tensor<1,dim>();
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_Nedelec<dim>::shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &/*p*/,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component < dim, ExcIndexRange (component, 0, dim));
+
+ switch (dim)
+ {
+ case 2: // 2D
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements. their second
+ // derivatives on the
+ // unit cell are zero
+ case 1:
+ {
+ return Tensor<2,dim>();
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+
+ case 3: // 3D
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements. their second
+ // derivatives on the
+ // unit cell are zero
+ case 1:
+ {
+ return Tensor<2,dim>();
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
+
+
+ // presently no other space
+ // dimension implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+
+ return Tensor<2,dim>();
+}
+
+
+//----------------------------------------------------------------------
+// Auxiliary functions
+//----------------------------------------------------------------------
+
+
+
+template <int dim>
+void FE_Nedelec<dim>::initialize_unit_support_points ()
+{
+//TODO: fix for higher orders. correct now for lowest order, all dimensions
+// is this correct? all DoFs on lines, none on faces or bubbles?
+
+ // all degrees of freedom are on
+ // edges, and their order is the
+ // same as the edges themselves
+ unit_support_points.resize(GeometryInfo<dim>::lines_per_cell * degree);
+ unsigned int index = 0;
+ for (unsigned int line=0; line<GeometryInfo<dim>::lines_per_cell; ++line)
+ {
+ const unsigned int
+ vertex_index_0 = GeometryInfo<dim>::vertices_adjacent_to_line(line,0),
+ vertex_index_1 = GeometryInfo<dim>::vertices_adjacent_to_line(line,1);
+
+ const Point<dim>
+ vertex_0 = GeometryInfo<dim>::unit_cell_vertex(vertex_index_0),
+ vertex_1 = GeometryInfo<dim>::unit_cell_vertex(vertex_index_1);
+
+ // place dofs equispaced
+ // between the vertices of each
+ // line
+ for (unsigned int d=0; d<degree; ++d, ++index)
+ unit_support_points[index]
+ = (vertex_0*(d+1) + vertex_1*(degree-d)) / (degree+1);
+ };
+};
+
+
+#if deal_II_dimension == 1
+
+template <>
+void FE_Nedelec<1>::initialize_unit_face_support_points ()
+{
+ // no faces in 1d, so nothing to do
+};
+
+#endif
+
+
+template <int dim>
+void FE_Nedelec<dim>::initialize_unit_face_support_points ()
+{
+//TODO: fix for higher orders. correct now for lowest order, all dimensions
+// is this correct? all DoFs on lines, none on faces or bubbles?
+ // do this the same as above, but
+ // for one dimension less
+ unit_face_support_points.resize(GeometryInfo<dim-1>::lines_per_cell * degree);
+ unsigned int index = 0;
+ for (unsigned int line=0; line<GeometryInfo<dim-1>::lines_per_cell; ++line)
+ {
+ const unsigned int
+ vertex_index_0 = GeometryInfo<dim-1>::vertices_adjacent_to_line(line,0),
+ vertex_index_1 = GeometryInfo<dim-1>::vertices_adjacent_to_line(line,1);
+
+ const Point<dim-1>
+ vertex_0 = GeometryInfo<dim-1>::unit_cell_vertex(vertex_index_0),
+ vertex_1 = GeometryInfo<dim-1>::unit_cell_vertex(vertex_index_1);
+
+ // place dofs equispaced
+ // between the vertices of each
+ // line
+ for (unsigned int d=0; d<degree; ++d, ++index)
+ unit_face_support_points[index]
+ = (vertex_0*(d+1) + vertex_1*(degree-d)) / (degree+1);
+ };
+};
+
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_Nedelec<dim>::get_dpo_vector(const unsigned int degree)
+{
+//TODO: fix for higher orders. correct now for lowest order, all dimensions
+ std::vector<unsigned int> dpo(dim+1, 0);
+// can this be done in a dimension independent and degree independent way?
+// if DoFs are located only on lines, the the following is the correct way
+
+ // put all degrees of freedom on
+ // the lines, and in particular
+ // @p{degree} DoFs per line:
+ dpo[1] = degree;
+ return dpo;
+}
+
+
+
+template <int dim>
+UpdateFlags
+FE_Nedelec<dim>::update_once (const UpdateFlags flags) const
+{
+//TODO: think about what this actually means here???
+ // for this kind of elements, only
+ // the values can be precomputed
+ // once and for all. set this flag
+ // if the values are requested at
+ // all
+ return (update_default | (flags & update_values));
+}
+
+
+
+template <int dim>
+UpdateFlags
+FE_Nedelec<dim>::update_each (const UpdateFlags flags) const
+{
+//TODO: think about what this actually means here???
+
+ UpdateFlags out = update_default;
+
+ if (flags & update_gradients)
+ out |= update_gradients | update_covariant_transformation;
+ if (flags & update_second_derivatives)
+ out |= update_second_derivatives | update_covariant_transformation;
+
+ return out;
+}
+
+
+
+//----------------------------------------------------------------------
+// Data field initialization
+//----------------------------------------------------------------------
+
+template <int dim>
+typename Mapping<dim>::InternalDataBase *
+FE_Nedelec<dim>::get_data (const UpdateFlags /*update_flags*/,
+ const Mapping<dim> &/*mapping*/,
+ const Quadrature<dim> &/*quadrature*/) const
+{
+ return 0;
+//TODO
+// // generate a new data object and
+// // initialize some fields
+// InternalData* data = new InternalData;
+
+// // check what needs to be
+// // initialized only once and what
+// // on every cell/face/subface we
+// // visit
+// data->update_once = update_once(update_flags);
+// data->update_each = update_each(update_flags);
+// data->update_flags = data->update_once | data->update_each;
+
+// const UpdateFlags flags(data->update_flags);
+// const unsigned int n_q_points = quadrature.n_quadrature_points;
+
+// // some scratch arrays
+// std::vector<double> values(0);
+// std::vector<Tensor<1,dim> > grads(0);
+// std::vector<Tensor<2,dim> > grad_grads(0);
+
+// // initialize fields only if really
+// // necessary. otherwise, don't
+// // allocate memory
+// if (flags & update_values)
+// {
+// values.resize (dofs_per_cell);
+// data->shape_values.resize(dofs_per_cell,
+// std::vector<double>(n_q_points));
+// }
+
+// if (flags & update_gradients)
+// {
+// grads.resize (dofs_per_cell);
+// data->shape_gradients.resize(dofs_per_cell,
+// std::vector<Tensor<1,dim> >(n_q_points));
+// }
+
+// // if second derivatives through
+// // finite differencing is required,
+// // then initialize some objects for
+// // that
+// if (flags & update_second_derivatives)
+// data->initialize_2nd (this, mapping, quadrature);
+
+// // next already fill those fields
+// // of which we have information by
+// // now. note that the shape
+// // gradients are only those on the
+// // unit cell, and need to be
+// // transformed when visiting an
+// // actual cell
+// if (flags & (update_values | update_gradients))
+// for (unsigned int i=0; i<n_q_points; ++i)
+// {
+// polynomial_space.compute(quadrature.point(i),
+// values, grads, grad_grads);
+
+// if (flags & update_values)
+// for (unsigned int k=0; k<dofs_per_cell; ++k)
+// data->shape_values[renumber[k]][i] = values[k];
+
+// if (flags & update_gradients)
+// for (unsigned int k=0; k<dofs_per_cell; ++k)
+// data->shape_gradients[renumber[k]][i] = grads[k];
+// }
+// return data;
+}
+
+
+
+
+//----------------------------------------------------------------------
+// Fill data of FEValues
+//----------------------------------------------------------------------
+
+template <int dim>
+void
+FE_Nedelec<dim>::fill_fe_values (const Mapping<dim> &/*mapping*/,
+ const typename DoFHandler<dim>::cell_iterator &/*cell*/,
+ const Quadrature<dim> &/*quadrature*/,
+ typename Mapping<dim>::InternalDataBase &/*mapping_data*/,
+ typename Mapping<dim>::InternalDataBase &/*fedata*/,
+ FEValuesData<dim> &/*data*/) const
+{
+//TODO!!
+// // convert data object to internal
+// // data for this class. fails with
+// // an exception if that is not
+// // possible
+// InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+// const UpdateFlags flags(fe_data.current_update_flags());
+
+// for (unsigned int k=0; k<dofs_per_cell; ++k)
+// {
+// if (flags & update_values)
+// for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
+// data.shape_values(k,i) = fe_data.shape_values[k][i];
+
+// if (flags & update_gradients)
+// mapping.transform_covariant(data.shape_gradients[k],
+// fe_data.shape_gradients[k],
+// mapping_data, 0);
+// }
+
+// if (flags & update_second_derivatives)
+// compute_2nd (mapping, cell, 0, mapping_data, fe_data, data);
+
+// fe_data.first_cell = false;
+}
+
+
+
+template <int dim>
+void
+FE_Nedelec<dim>::fill_fe_face_values (const Mapping<dim> &/*mapping*/,
+ const typename DoFHandler<dim>::cell_iterator &/*cell*/,
+ const unsigned int /*face*/,
+ const Quadrature<dim-1> &/*quadrature*/,
+ typename Mapping<dim>::InternalDataBase &/*mapping_data*/,
+ typename Mapping<dim>::InternalDataBase &/*fedata*/,
+ FEValuesData<dim> &/*data*/) const
+{
+//TODO!!
+// // convert data object to internal
+// // data for this class. fails with
+// // an exception if that is not
+// // possible
+// InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+// // offset determines which data set
+// // to take (all data sets for all
+// // faces are stored contiguously)
+// const unsigned int offset = face * quadrature.n_quadrature_points;
+
+// const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+
+// for (unsigned int k=0; k<dofs_per_cell; ++k)
+// {
+// for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
+// if (flags & update_values)
+// data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
+
+// if (flags & update_gradients)
+// mapping.transform_covariant(data.shape_gradients[k],
+// fe_data.shape_gradients[k],
+// mapping_data, offset);
+// }
+
+// if (flags & update_second_derivatives)
+// compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+
+// fe_data.first_cell = false;
+}
+
+
+
+template <int dim>
+void
+FE_Nedelec<dim>::fill_fe_subface_values (const Mapping<dim> &/*mapping*/,
+ const typename DoFHandler<dim>::cell_iterator &/*cell*/,
+ const unsigned int /*face*/,
+ const unsigned int /*subface*/,
+ const Quadrature<dim-1> &/*quadrature*/,
+ typename Mapping<dim>::InternalDataBase &/*mapping_data*/,
+ typename Mapping<dim>::InternalDataBase &/*fedata*/,
+ FEValuesData<dim> &/*data*/) const
+{
+//TODO!!
+// // convert data object to internal
+// // data for this class. fails with
+// // an exception if that is not
+// // possible
+// InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+// // offset determines which data set
+// // to take (all data sets for all
+// // sub-faces are stored contiguously)
+// const unsigned int offset = (face * GeometryInfo<dim>::subfaces_per_face + subface)
+// * quadrature.n_quadrature_points;
+
+// const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+
+// for (unsigned int k=0; k<dofs_per_cell; ++k)
+// {
+// for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
+// if (flags & update_values)
+// data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
+
+// if (flags & update_gradients)
+// mapping.transform_covariant(data.shape_gradients[k],
+// fe_data.shape_gradients[k],
+// mapping_data, offset);
+// }
+
+// if (flags & update_second_derivatives)
+// compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+
+// fe_data.first_cell = false;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_Nedelec<dim>::n_base_elements () const
+{
+ return 1;
+};
+
+
+
+template <int dim>
+const FiniteElement<dim> &
+FE_Nedelec<dim>::base_element (const unsigned int index) const
+{
+ Assert (index==0, ExcIndexRange(index, 0, 1));
+ return *this;
+};
+
+
+
+template <int dim>
+bool
+FE_Nedelec<dim>::has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ Assert (shape_index < dofs_per_cell,
+ ExcIndexRange (shape_index, 0, dofs_per_cell));
+ Assert (face_index < GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
+
+//TODO: fix for higher orders. correct now for lowest order, all dimensions
+//TODO!!
+// can this be done in a way that is dimension and degree independent?
+
+ // all degrees of freedom are on
+ // lines, so also on a face. the
+ // question is whether it has
+ // support on this particular face
+ Assert (false, ExcNotImplemented());
+ return true;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_Nedelec<dim>::memory_consumption () const
+{
+ Assert (false, ExcNotImplemented ());
+ return 0;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_Nedelec<dim>::get_degree () const
+{
+ return degree;
+};
+
+
+
+template class FE_Nedelec<deal_II_dimension>;
--- /dev/null
+//----------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002 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.
+//
+//----------------------------------------------------------------
+
+
+
+// only compile this file if in 1d. note that Nedelec elemets do not
+// make much sense in 1d, so this file only contains dummy
+// implementations to avoid linker errors due to missing symbols
+#if deal_II_dimension == 1
+
+
+#include <fe/fe_nedelec.h>
+
+
+template <>
+const double * const
+FE_Nedelec<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell] =
+{};
+
+
+template <>
+const unsigned int
+FE_Nedelec<1>::Matrices::n_embedding_matrices = 0;
+
+
+
+// No constraints in 1d
+template <>
+const unsigned int
+FE_Nedelec<1>::Matrices::n_constraint_matrices = 0;
+
+
+template <>
+const double * const
+FE_Nedelec<1>::Matrices::constraint_matrices[] = {};
+
+
+#else // #if deal_II_dimension
+// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
+// files, so provide some dummy code
+namespace { void dummy () {}; };
+#endif // #if deal_II_dimension == 1
+
--- /dev/null
+//----------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2001, 2002 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.
+//
+//----------------------------------------------------------------
+
+
+// only compile this file if in 2d
+#if deal_II_dimension == 2
+
+
+#include <fe/fe_nedelec.h>
+
+// Transfer matrices for finite elements: have one matrix for each of
+// the four child cells which tells us how the degrees of freedom on
+// the child cell are obtained from the degrees of freedom on the
+// mother cell
+namespace FE_Nedelec_2d
+{
+ static const double q1_into_q1_refined_0[] =
+ {
+ 1., 0, 0, 0,
+ 0, 0.5,0, 0.5,
+ 0.5, 0, 0.5,0,
+ 0, 0, 0, 1
+ };
+
+ static const double q1_into_q1_refined_1[] =
+ {
+//TODO[Anna]
+// 13.5/27., 13.5/27., 0., 0.,
+// 0., 1., 0., 0.,
+// 0., 13.5/27., 13.5/27., 0.,
+// 6.75/27., 6.75/27., 6.75/27., 6.75/27.,
+ };
+
+ static const double q1_into_q1_refined_2[] =
+ {
+// 6.75/27., 6.75/27., 6.75/27., 6.75/27.,
+// 0., 13.5/27., 13.5/27., 0.,
+// 0., 0., 1., 0.,
+// 0., 0., 13.5/27., 13.5/27.,
+ };
+
+ static const double q1_into_q1_refined_3[] =
+ {
+// 13.5/27., 0., 0., 13.5/27.,
+// 6.75/27., 6.75/27., 6.75/27., 6.75/27.,
+// 0., 0., 13.5/27., 13.5/27.,
+// 0., 0., 0., 1.,
+ };
+}; // namespace FE_Nedelec_2d
+
+
+// embedding matrices
+
+template <>
+const double * const
+FE_Nedelec<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell] =
+{
+ { FE_Nedelec_2d::q1_into_q1_refined_0, FE_Nedelec_2d::q1_into_q1_refined_1,
+ FE_Nedelec_2d::q1_into_q1_refined_2, FE_Nedelec_2d::q1_into_q1_refined_3 }
+};
+
+
+template <>
+const unsigned int
+FE_Nedelec<2>::Matrices::n_embedding_matrices
+ = sizeof(FE_Nedelec<2>::Matrices::embedding) /
+ sizeof(FE_Nedelec<2>::Matrices::embedding[0]);
+
+
+// Constraint matrices: how do the new value on child faces depend on
+// the values on the mother face if that face has a hanging node
+namespace FE_Nedelec_2d
+{
+ static const double constraint_q1[] =
+ {
+ // the function is constant
+ // along each edge, so each
+ // degree of freedom on the
+ // refined edge has the same
+ // value as that on the
+ // coarse edge
+ 1., 1.
+ };
+
+};
+
+
+template <>
+const double * const
+FE_Nedelec<2>::Matrices::constraint_matrices[] =
+{
+ FE_Nedelec_2d::constraint_q1
+};
+
+
+template <>
+const unsigned int
+FE_Nedelec<2>::Matrices::n_constraint_matrices
+ = sizeof(FE_Nedelec<2>::Matrices::constraint_matrices) /
+ sizeof(FE_Nedelec<2>::Matrices::constraint_matrices[0]);
+
+
+
+#else // #if deal_II_dimension
+// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
+// files, so provide some dummy code
+namespace { void dummy () {}; };
+#endif // #if deal_II_dimension == 2
--- /dev/null
+//----------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2001, 2002 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.
+//
+//----------------------------------------------------------------
+
+// Transfer matrices for finite elements
+
+
+// only compile this file if in 3d
+#if deal_II_dimension == 3
+
+#include <fe/fe_nedelec.h>
+
+namespace FE_Nedelec_3d
+{
+ static const double q1_into_q1_refined_0[] =
+ {
+// 1., 0., 0., 0., 0., 0., 0., 0.,
+// 13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0.,
+// 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0.,
+// 13.5/27., 0., 0., 13.5/27., 0., 0., 0., 0.,
+// 13.5/27., 0., 0., 0., 13.5/27., 0., 0., 0.,
+// 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0.,
+// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
+// 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27.,
+ };
+
+ static const double q1_into_q1_refined_1[] =
+ {
+// 13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0.,
+// 0., 1., 0., 0., 0., 0., 0., 0.,
+// 0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0.,
+// 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0.,
+// 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0.,
+// 0., 13.5/27., 0., 0., 0., 13.5/27., 0., 0.,
+// 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0.,
+// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
+ };
+
+ static const double q1_into_q1_refined_2[] =
+ {
+// 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0.,
+// 0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0.,
+// 0., 0., 1., 0., 0., 0., 0., 0.,
+// 0., 0., 13.5/27., 13.5/27., 0., 0., 0., 0.,
+// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
+// 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0.,
+// 0., 0., 13.5/27., 0., 0., 0., 13.5/27., 0.,
+// 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27.,
+ };
+
+ static const double q1_into_q1_refined_3[] =
+ {
+// 13.5/27., 0., 0., 13.5/27., 0., 0., 0., 0.,
+// 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0.,
+// 0., 0., 13.5/27., 13.5/27., 0., 0., 0., 0.,
+// 0., 0., 0., 1., 0., 0., 0., 0.,
+// 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27.,
+// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
+// 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27.,
+// 0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27.,
+ };
+
+ static const double q1_into_q1_refined_4[] =
+ {
+// 13.5/27., 0., 0., 0., 13.5/27., 0., 0., 0.,
+// 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0.,
+// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
+// 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27.,
+// 0., 0., 0., 0., 1., 0., 0., 0.,
+// 0., 0., 0., 0., 13.5/27., 13.5/27., 0., 0.,
+// 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27.,
+// 0., 0., 0., 0., 13.5/27., 0., 0., 13.5/27.,
+ };
+
+ static const double q1_into_q1_refined_5[] =
+ {
+// 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0.,
+// 0., 13.5/27., 0., 0., 0., 13.5/27., 0., 0.,
+// 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0.,
+// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
+// 0., 0., 0., 0., 13.5/27., 13.5/27., 0., 0.,
+// 0., 0., 0., 0., 0., 1., 0., 0.,
+// 0., 0., 0., 0., 0., 13.5/27., 13.5/27., 0.,
+// 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27.,
+ };
+
+ static const double q1_into_q1_refined_6[] =
+ {
+// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
+// 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0.,
+// 0., 0., 13.5/27., 0., 0., 0., 13.5/27., 0.,
+// 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27.,
+// 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27.,
+// 0., 0., 0., 0., 0., 13.5/27., 13.5/27., 0.,
+// 0., 0., 0., 0., 0., 0., 1., 0.,
+// 0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27.,
+ };
+
+ static const double q1_into_q1_refined_7[] =
+ {
+// 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27.,
+// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27.,
+// 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27.,
+// 0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27.,
+// 0., 0., 0., 0., 13.5/27., 0., 0., 13.5/27.,
+// 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27.,
+// 0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27.,
+// 0., 0., 0., 0., 0., 0., 0., 1.,
+ };
+
+}; // namespace FE_Nedelec_3d
+
+
+// embedding matrices
+
+template <>
+const double * const
+FE_Nedelec<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell] =
+{
+ { FE_Nedelec_3d::q1_into_q1_refined_0, FE_Nedelec_3d::q1_into_q1_refined_1,
+ FE_Nedelec_3d::q1_into_q1_refined_2, FE_Nedelec_3d::q1_into_q1_refined_3,
+ FE_Nedelec_3d::q1_into_q1_refined_4, FE_Nedelec_3d::q1_into_q1_refined_5,
+ FE_Nedelec_3d::q1_into_q1_refined_6, FE_Nedelec_3d::q1_into_q1_refined_7 }
+};
+
+
+template <>
+const unsigned int
+FE_Nedelec<3>::Matrices::n_embedding_matrices
+ = sizeof(FE_Nedelec<3>::Matrices::embedding) /
+ sizeof(FE_Nedelec<3>::Matrices::embedding[0]);
+
+
+
+// Constraint matrices
+
+namespace FE_Nedelec_3d
+{
+ static const double constraint_q1[] =
+ {
+//TODO[WB]
+// .25,.25,.25,.25,
+// .5,.5,0.,0.,
+// 0.,.5,.5,0.,
+// 0.,0.,.5,.5,
+// .5,0.,0.,.5
+ };
+
+};
+
+
+
+template <>
+const double * const
+FE_Nedelec<3>::Matrices::constraint_matrices[] =
+{
+ FE_Nedelec_3d::constraint_q1
+};
+
+
+
+template <>
+const unsigned int
+FE_Nedelec<3>::Matrices::n_constraint_matrices
+ = sizeof(FE_Nedelec<3>::Matrices::constraint_matrices) /
+ sizeof(FE_Nedelec<3>::Matrices::constraint_matrices[0]);
+
+
+
+#else // #if deal_II_dimension
+// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
+// files, so provide some dummy code
+namespace { void dummy () {}; };
+#endif // #if deal_II_dimension == 3
#include <fe/fe_values.h>
-
-
template <int dim>
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> (1,false),
+ std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
+ std::vector<bool>(1,true))),
degree(degree),
renumber(dofs_per_cell, 0),
renumber_inverse(dofs_per_cell, 0),
face_renumber(dofs_per_face, 0),
- poly(0)
+ polynomial_space(LagrangeEquidistant::generate_complete_basis(degree))
{
- // Q0 elements cannot be
- // continuous, use FE_DGQ<dim>(0)
- // instead
- Assert (degree>0, ExcNotImplemented());
- std::vector<LagrangeEquidistant> v;
- for (unsigned int i=0;i<=degree;++i)
- v.push_back(LagrangeEquidistant(degree,i));
-
- poly = new TensorProductPolynomials<dim> (v);
-
// do some internal book-keeping on
// cells and faces. if in 1d, the
// face function is empty
// then fill restriction
// matrices. they are hardcoded for
- // the first few elements
+ // the first few elements. in
+ // contrast to the other matrices,
+ // these are not stored in the
+ // files fe_q_[123]d.cc, since they
+ // contain only a rather small
+ // number of zeros, and storing
+ // them element-wise is more
+ // expensive than just setting the
+ // nonzero elements as done here
switch (dim)
{
- case 1:
- switch (degree)
- {
- case 1:
- restriction[0](0,0) = 1;
- restriction[1](1,1) = 1;
- break;
- case 2:
- restriction[0](0,0) = 1;
- restriction[0](2,1) = 1;
- restriction[1](1,1) = 1;
- restriction[1](1,1) = 1;
- break;
- case 3:
- restriction[0](0,0) = 1;
- restriction[0](2,3) = 1;
- restriction[1](1,1) = 1;
- restriction[1](3,2) = 1;
- break;
- case 4:
- restriction[0](0,0) = 1;
- restriction[0](2,3) = 1;
- restriction[0](3,1) = 1;
- restriction[1](1,1) = 1;
- restriction[1](3,0) = 1;
- restriction[1](4,3) = 1;
- break;
- default:
- for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
- restriction[i].reinit(0,0);
- }
- break;
- case 2:
- switch (degree)
- {
- case 1:
- restriction[0](0,0) = 1;
- restriction[1](1,1) = 1;
- restriction[2](2,2) = 1;
- restriction[3](3,3) = 1;
- break;
- case 2:
- restriction[0](0,0) = 1;
- restriction[0](4,1) = 1;
- restriction[0](7,3) = 1;
- restriction[0](8,2) = 1;
- restriction[1](1,1) = 1;
- restriction[1](4,0) = 1;
- restriction[1](5,2) = 1;
- restriction[1](8,3) = 1;
- restriction[2](2,2) = 1;
- restriction[2](5,1) = 1;
- restriction[2](6,3) = 1;
- restriction[2](8,0) = 1;
- restriction[3](3,3) = 1;
- restriction[3](6,2) = 1;
- restriction[3](7,0) = 1;
- restriction[3](8,1) = 1;
- break;
- case 3:
- restriction[0](0,0) = 1;
- restriction[0](4,5) = 1;
- restriction[0](10,11) = 1;
- restriction[0](12,15) = 1;
- restriction[1](1,1) = 1;
- restriction[1](5,4) = 1;
- restriction[1](6,7) = 1;
- restriction[1](13,14) = 1;
- restriction[2](2,2) = 1;
- restriction[2](7,6) = 1;
- restriction[2](9,8) = 1;
- restriction[2](15,12) = 1;
- restriction[3](3,3) = 1;
- restriction[3](8,9) = 1;
- restriction[3](11,10) = 1;
- restriction[3](14,13) = 1;
- break;
- case 4:
- restriction[0](0,0) = 1;
- restriction[0](4,5) = 1;
- restriction[0](5,1) = 1;
- restriction[0](13,14) = 1;
- restriction[0](14,3) = 1;
- restriction[0](16,20) = 1;
- restriction[0](17,8) = 1;
- restriction[0](19,11) = 1;
- restriction[0](20,2) = 1;
- restriction[1](1,1) = 1;
- restriction[1](5,0) = 1;
- restriction[1](6,5) = 1;
- restriction[1](7,8) = 1;
- restriction[1](8,2) = 1;
- restriction[1](17,14) = 1;
- restriction[1](18,20) = 1;
- restriction[1](20,3) = 1;
- restriction[1](21,11) = 1;
- restriction[2](2,2) = 1;
- restriction[2](8,1) = 1;
- restriction[2](9,8) = 1;
- restriction[2](11,3) = 1;
- restriction[2](12,11) = 1;
- restriction[2](20,0) = 1;
- restriction[2](21,5) = 1;
- restriction[2](23,14) = 1;
- restriction[2](24,20) = 1;
- restriction[3](3,3) = 1;
- restriction[3](10,11) = 1;
- restriction[3](11,2) = 1;
- restriction[3](14,0) = 1;
- restriction[3](15,14) = 1;
- restriction[3](19,5) = 1;
- restriction[3](20,1) = 1;
- restriction[3](22,20) = 1;
- restriction[3](23,8) = 1;
- break;
- default:
- for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
- restriction[i].reinit(0,0);
- }
- break;
- case 3:
- switch (degree)
- {
- case 1:
- restriction[0](0,0) = 1;
- restriction[1](1,1) = 1;
- restriction[2](2,2) = 1;
- restriction[3](3,3) = 1;
- restriction[4](4,4) = 1;
- restriction[5](5,5) = 1;
- restriction[6](6,6) = 1;
- restriction[7](7,7) = 1;
- break;
- case 2:
- restriction[0](0,0) = 1;
- restriction[0](8,1) = 1;
- restriction[0](11,3) = 1;
- restriction[0](16,4) = 1;
- restriction[0](20,2) = 1;
- restriction[0](22,5) = 1;
- restriction[0](25,7) = 1;
- restriction[0](26,6) = 1;
- restriction[1](1,1) = 1;
- restriction[1](8,0) = 1;
- restriction[1](9,2) = 1;
- restriction[1](17,5) = 1;
- restriction[1](20,3) = 1;
- restriction[1](22,4) = 1;
- restriction[1](23,6) = 1;
- restriction[1](26,7) = 1;
- restriction[2](2,2) = 1;
- restriction[2](9,1) = 1;
- restriction[2](10,3) = 1;
- restriction[2](18,6) = 1;
- restriction[2](20,0) = 1;
- restriction[2](23,5) = 1;
- restriction[2](24,7) = 1;
- restriction[2](26,4) = 1;
- restriction[3](3,3) = 1;
- restriction[3](10,2) = 1;
- restriction[3](11,0) = 1;
- restriction[3](19,7) = 1;
- restriction[3](20,1) = 1;
- restriction[3](24,6) = 1;
- restriction[3](25,4) = 1;
- restriction[3](26,5) = 1;
- restriction[4](4,4) = 1;
- restriction[4](12,5) = 1;
- restriction[4](15,7) = 1;
- restriction[4](16,0) = 1;
- restriction[4](21,6) = 1;
- restriction[4](22,1) = 1;
- restriction[4](25,3) = 1;
- restriction[4](26,2) = 1;
- restriction[5](5,5) = 1;
- restriction[5](12,4) = 1;
- restriction[5](13,6) = 1;
- restriction[5](17,1) = 1;
- restriction[5](21,7) = 1;
- restriction[5](22,0) = 1;
- restriction[5](23,2) = 1;
- restriction[5](26,3) = 1;
- restriction[6](6,6) = 1;
- restriction[6](13,5) = 1;
- restriction[6](14,7) = 1;
- restriction[6](18,2) = 1;
- restriction[6](21,4) = 1;
- restriction[6](23,1) = 1;
- restriction[6](24,3) = 1;
- restriction[6](26,0) = 1;
- restriction[7](7,7) = 1;
- restriction[7](14,6) = 1;
- restriction[7](15,4) = 1;
- restriction[7](19,3) = 1;
- restriction[7](21,5) = 1;
- restriction[7](24,2) = 1;
- restriction[7](25,0) = 1;
- restriction[7](26,1) = 1;
- break;
- case 3:
- restriction[0](0,0) = 1;
- restriction[0](8,9) = 1;
- restriction[0](14,15) = 1;
- restriction[0](24,25) = 1;
- restriction[0](32,35) = 1;
- restriction[0](40,43) = 1;
- restriction[0](52,55) = 1;
- restriction[0](56,63) = 1;
- restriction[1](1,1) = 1;
- restriction[1](9,8) = 1;
- restriction[1](10,11) = 1;
- restriction[1](26,27) = 1;
- restriction[1](33,34) = 1;
- restriction[1](41,42) = 1;
- restriction[1](44,47) = 1;
- restriction[1](57,62) = 1;
- restriction[2](2,2) = 1;
- restriction[2](11,10) = 1;
- restriction[2](13,12) = 1;
- restriction[2](28,29) = 1;
- restriction[2](35,32) = 1;
- restriction[2](46,45) = 1;
- restriction[2](49,50) = 1;
- restriction[2](61,58) = 1;
- restriction[3](3,3) = 1;
- restriction[3](12,13) = 1;
- restriction[3](15,14) = 1;
- restriction[3](30,31) = 1;
- restriction[3](34,33) = 1;
- restriction[3](48,51) = 1;
- restriction[3](54,53) = 1;
- restriction[3](60,59) = 1;
- restriction[4](4,4) = 1;
- restriction[4](16,17) = 1;
- restriction[4](22,23) = 1;
- restriction[4](25,24) = 1;
- restriction[4](36,39) = 1;
- restriction[4](42,41) = 1;
- restriction[4](53,54) = 1;
- restriction[4](58,61) = 1;
- restriction[5](5,5) = 1;
- restriction[5](17,16) = 1;
- restriction[5](18,19) = 1;
- restriction[5](27,26) = 1;
- restriction[5](37,38) = 1;
- restriction[5](43,40) = 1;
- restriction[5](45,46) = 1;
- restriction[5](59,60) = 1;
- restriction[6](6,6) = 1;
- restriction[6](19,18) = 1;
- restriction[6](21,20) = 1;
- restriction[6](29,28) = 1;
- restriction[6](39,36) = 1;
- restriction[6](47,44) = 1;
- restriction[6](51,48) = 1;
- restriction[6](63,56) = 1;
- restriction[7](7,7) = 1;
- restriction[7](20,21) = 1;
- restriction[7](23,22) = 1;
- restriction[7](31,30) = 1;
- restriction[7](38,37) = 1;
- restriction[7](50,49) = 1;
- restriction[7](55,52) = 1;
- restriction[7](62,57) = 1;
- break;
- case 4:
- restriction[0](0,0) = 1;
- restriction[0](8,9) = 1;
- restriction[0](9,1) = 1;
- restriction[0](17,18) = 1;
- restriction[0](18,3) = 1;
- restriction[0](32,33) = 1;
- restriction[0](33,4) = 1;
- restriction[0](44,48) = 1;
- restriction[0](45,12) = 1;
- restriction[0](47,15) = 1;
- restriction[0](48,2) = 1;
- restriction[0](62,66) = 1;
- restriction[0](63,36) = 1;
- restriction[0](65,21) = 1;
- restriction[0](66,5) = 1;
- restriction[0](89,93) = 1;
- restriction[0](90,30) = 1;
- restriction[0](92,42) = 1;
- restriction[0](93,7) = 1;
- restriction[0](98,111) = 1;
- restriction[0](99,75) = 1;
- restriction[0](101,57) = 1;
- restriction[0](102,24) = 1;
- restriction[0](107,84) = 1;
- restriction[0](108,39) = 1;
- restriction[0](110,27) = 1;
- restriction[0](111,6) = 1;
- restriction[1](1,1) = 1;
- restriction[1](9,0) = 1;
- restriction[1](10,9) = 1;
- restriction[1](11,12) = 1;
- restriction[1](12,2) = 1;
- restriction[1](35,36) = 1;
- restriction[1](36,5) = 1;
- restriction[1](45,18) = 1;
- restriction[1](46,48) = 1;
- restriction[1](48,3) = 1;
- restriction[1](49,15) = 1;
- restriction[1](63,33) = 1;
- restriction[1](64,66) = 1;
- restriction[1](66,4) = 1;
- restriction[1](67,21) = 1;
- restriction[1](71,75) = 1;
- restriction[1](72,24) = 1;
- restriction[1](74,39) = 1;
- restriction[1](75,6) = 1;
- restriction[1](99,93) = 1;
- restriction[1](100,111) = 1;
- restriction[1](102,30) = 1;
- restriction[1](103,57) = 1;
- restriction[1](108,42) = 1;
- restriction[1](109,84) = 1;
- restriction[1](111,7) = 1;
- restriction[1](112,27) = 1;
- restriction[2](2,2) = 1;
- restriction[2](12,1) = 1;
- restriction[2](13,12) = 1;
- restriction[2](15,3) = 1;
- restriction[2](16,15) = 1;
- restriction[2](38,39) = 1;
- restriction[2](39,6) = 1;
- restriction[2](48,0) = 1;
- restriction[2](49,9) = 1;
- restriction[2](51,18) = 1;
- restriction[2](52,48) = 1;
- restriction[2](74,36) = 1;
- restriction[2](75,5) = 1;
- restriction[2](77,75) = 1;
- restriction[2](78,24) = 1;
- restriction[2](81,42) = 1;
- restriction[2](82,84) = 1;
- restriction[2](84,7) = 1;
- restriction[2](85,27) = 1;
- restriction[2](108,33) = 1;
- restriction[2](109,66) = 1;
- restriction[2](111,4) = 1;
- restriction[2](112,21) = 1;
- restriction[2](117,93) = 1;
- restriction[2](118,111) = 1;
- restriction[2](120,30) = 1;
- restriction[2](121,57) = 1;
- restriction[3](3,3) = 1;
- restriction[3](14,15) = 1;
- restriction[3](15,2) = 1;
- restriction[3](18,0) = 1;
- restriction[3](19,18) = 1;
- restriction[3](41,42) = 1;
- restriction[3](42,7) = 1;
- restriction[3](47,9) = 1;
- restriction[3](48,1) = 1;
- restriction[3](50,48) = 1;
- restriction[3](51,12) = 1;
- restriction[3](80,84) = 1;
- restriction[3](81,39) = 1;
- restriction[3](83,27) = 1;
- restriction[3](84,6) = 1;
- restriction[3](92,33) = 1;
- restriction[3](93,4) = 1;
- restriction[3](95,93) = 1;
- restriction[3](96,30) = 1;
- restriction[3](107,66) = 1;
- restriction[3](108,36) = 1;
- restriction[3](110,21) = 1;
- restriction[3](111,5) = 1;
- restriction[3](116,111) = 1;
- restriction[3](117,75) = 1;
- restriction[3](119,57) = 1;
- restriction[3](120,24) = 1;
- restriction[4](4,4) = 1;
- restriction[4](20,21) = 1;
- restriction[4](21,5) = 1;
- restriction[4](29,30) = 1;
- restriction[4](30,7) = 1;
- restriction[4](33,0) = 1;
- restriction[4](34,33) = 1;
- restriction[4](53,57) = 1;
- restriction[4](54,24) = 1;
- restriction[4](56,27) = 1;
- restriction[4](57,6) = 1;
- restriction[4](65,9) = 1;
- restriction[4](66,1) = 1;
- restriction[4](68,66) = 1;
- restriction[4](69,36) = 1;
- restriction[4](90,18) = 1;
- restriction[4](91,93) = 1;
- restriction[4](93,3) = 1;
- restriction[4](94,42) = 1;
- restriction[4](101,48) = 1;
- restriction[4](102,12) = 1;
- restriction[4](104,111) = 1;
- restriction[4](105,75) = 1;
- restriction[4](110,15) = 1;
- restriction[4](111,2) = 1;
- restriction[4](113,84) = 1;
- restriction[4](114,39) = 1;
- restriction[5](5,5) = 1;
- restriction[5](21,4) = 1;
- restriction[5](22,21) = 1;
- restriction[5](23,24) = 1;
- restriction[5](24,6) = 1;
- restriction[5](36,1) = 1;
- restriction[5](37,36) = 1;
- restriction[5](54,30) = 1;
- restriction[5](55,57) = 1;
- restriction[5](57,7) = 1;
- restriction[5](58,27) = 1;
- restriction[5](66,0) = 1;
- restriction[5](67,9) = 1;
- restriction[5](69,33) = 1;
- restriction[5](70,66) = 1;
- restriction[5](72,12) = 1;
- restriction[5](73,75) = 1;
- restriction[5](75,2) = 1;
- restriction[5](76,39) = 1;
- restriction[5](102,18) = 1;
- restriction[5](103,48) = 1;
- restriction[5](105,93) = 1;
- restriction[5](106,111) = 1;
- restriction[5](111,3) = 1;
- restriction[5](112,15) = 1;
- restriction[5](114,42) = 1;
- restriction[5](115,84) = 1;
- restriction[6](6,6) = 1;
- restriction[6](24,5) = 1;
- restriction[6](25,24) = 1;
- restriction[6](27,7) = 1;
- restriction[6](28,27) = 1;
- restriction[6](39,2) = 1;
- restriction[6](40,39) = 1;
- restriction[6](57,4) = 1;
- restriction[6](58,21) = 1;
- restriction[6](60,30) = 1;
- restriction[6](61,57) = 1;
- restriction[6](75,1) = 1;
- restriction[6](76,36) = 1;
- restriction[6](78,12) = 1;
- restriction[6](79,75) = 1;
- restriction[6](84,3) = 1;
- restriction[6](85,15) = 1;
- restriction[6](87,42) = 1;
- restriction[6](88,84) = 1;
- restriction[6](111,0) = 1;
- restriction[6](112,9) = 1;
- restriction[6](114,33) = 1;
- restriction[6](115,66) = 1;
- restriction[6](120,18) = 1;
- restriction[6](121,48) = 1;
- restriction[6](123,93) = 1;
- restriction[6](124,111) = 1;
- restriction[7](7,7) = 1;
- restriction[7](26,27) = 1;
- restriction[7](27,6) = 1;
- restriction[7](30,4) = 1;
- restriction[7](31,30) = 1;
- restriction[7](42,3) = 1;
- restriction[7](43,42) = 1;
- restriction[7](56,21) = 1;
- restriction[7](57,5) = 1;
- restriction[7](59,57) = 1;
- restriction[7](60,24) = 1;
- restriction[7](83,15) = 1;
- restriction[7](84,2) = 1;
- restriction[7](86,84) = 1;
- restriction[7](87,39) = 1;
- restriction[7](93,0) = 1;
- restriction[7](94,33) = 1;
- restriction[7](96,18) = 1;
- restriction[7](97,93) = 1;
- restriction[7](110,9) = 1;
- restriction[7](111,1) = 1;
- restriction[7](113,66) = 1;
- restriction[7](114,36) = 1;
- restriction[7](119,48) = 1;
- restriction[7](120,12) = 1;
- restriction[7](122,111) = 1;
- restriction[7](123,75) = 1;
- break;
- default:
- for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell;++i)
- restriction[i].reinit(0,0);
- }
- break;
+ case 1: // 1d
+ {
+ switch (degree)
+ {
+ case 1:
+ restriction[0](0,0) = 1;
+ restriction[1](1,1) = 1;
+ break;
+ case 2:
+ restriction[0](0,0) = 1;
+ restriction[0](2,1) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](1,1) = 1;
+ break;
+ case 3:
+ restriction[0](0,0) = 1;
+ restriction[0](2,3) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](3,2) = 1;
+ break;
+ case 4:
+ restriction[0](0,0) = 1;
+ restriction[0](2,3) = 1;
+ restriction[0](3,1) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](3,0) = 1;
+ restriction[1](4,3) = 1;
+ break;
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ restriction[i].reinit(0,0);
+ };
+ }
+ break;
+ };
+
+ case 2: // 2d
+ {
+ switch (degree)
+ {
+ case 1:
+ restriction[0](0,0) = 1;
+ restriction[1](1,1) = 1;
+ restriction[2](2,2) = 1;
+ restriction[3](3,3) = 1;
+ break;
+ case 2:
+ restriction[0](0,0) = 1;
+ restriction[0](4,1) = 1;
+ restriction[0](7,3) = 1;
+ restriction[0](8,2) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](4,0) = 1;
+ restriction[1](5,2) = 1;
+ restriction[1](8,3) = 1;
+ restriction[2](2,2) = 1;
+ restriction[2](5,1) = 1;
+ restriction[2](6,3) = 1;
+ restriction[2](8,0) = 1;
+ restriction[3](3,3) = 1;
+ restriction[3](6,2) = 1;
+ restriction[3](7,0) = 1;
+ restriction[3](8,1) = 1;
+ break;
+ case 3:
+ restriction[0](0,0) = 1;
+ restriction[0](4,5) = 1;
+ restriction[0](10,11) = 1;
+ restriction[0](12,15) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](5,4) = 1;
+ restriction[1](6,7) = 1;
+ restriction[1](13,14) = 1;
+ restriction[2](2,2) = 1;
+ restriction[2](7,6) = 1;
+ restriction[2](9,8) = 1;
+ restriction[2](15,12) = 1;
+ restriction[3](3,3) = 1;
+ restriction[3](8,9) = 1;
+ restriction[3](11,10) = 1;
+ restriction[3](14,13) = 1;
+ break;
+ case 4:
+ restriction[0](0,0) = 1;
+ restriction[0](4,5) = 1;
+ restriction[0](5,1) = 1;
+ restriction[0](13,14) = 1;
+ restriction[0](14,3) = 1;
+ restriction[0](16,20) = 1;
+ restriction[0](17,8) = 1;
+ restriction[0](19,11) = 1;
+ restriction[0](20,2) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](5,0) = 1;
+ restriction[1](6,5) = 1;
+ restriction[1](7,8) = 1;
+ restriction[1](8,2) = 1;
+ restriction[1](17,14) = 1;
+ restriction[1](18,20) = 1;
+ restriction[1](20,3) = 1;
+ restriction[1](21,11) = 1;
+ restriction[2](2,2) = 1;
+ restriction[2](8,1) = 1;
+ restriction[2](9,8) = 1;
+ restriction[2](11,3) = 1;
+ restriction[2](12,11) = 1;
+ restriction[2](20,0) = 1;
+ restriction[2](21,5) = 1;
+ restriction[2](23,14) = 1;
+ restriction[2](24,20) = 1;
+ restriction[3](3,3) = 1;
+ restriction[3](10,11) = 1;
+ restriction[3](11,2) = 1;
+ restriction[3](14,0) = 1;
+ restriction[3](15,14) = 1;
+ restriction[3](19,5) = 1;
+ restriction[3](20,1) = 1;
+ restriction[3](22,20) = 1;
+ restriction[3](23,8) = 1;
+ break;
+
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ restriction[i].reinit(0,0);
+ };
+ }
+ break;
+ };
+
+ case 3: // 3d
+ {
+ switch (degree)
+ {
+ case 1:
+ restriction[0](0,0) = 1;
+ restriction[1](1,1) = 1;
+ restriction[2](2,2) = 1;
+ restriction[3](3,3) = 1;
+ restriction[4](4,4) = 1;
+ restriction[5](5,5) = 1;
+ restriction[6](6,6) = 1;
+ restriction[7](7,7) = 1;
+ break;
+ case 2:
+ restriction[0](0,0) = 1;
+ restriction[0](8,1) = 1;
+ restriction[0](11,3) = 1;
+ restriction[0](16,4) = 1;
+ restriction[0](20,2) = 1;
+ restriction[0](22,5) = 1;
+ restriction[0](25,7) = 1;
+ restriction[0](26,6) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](8,0) = 1;
+ restriction[1](9,2) = 1;
+ restriction[1](17,5) = 1;
+ restriction[1](20,3) = 1;
+ restriction[1](22,4) = 1;
+ restriction[1](23,6) = 1;
+ restriction[1](26,7) = 1;
+ restriction[2](2,2) = 1;
+ restriction[2](9,1) = 1;
+ restriction[2](10,3) = 1;
+ restriction[2](18,6) = 1;
+ restriction[2](20,0) = 1;
+ restriction[2](23,5) = 1;
+ restriction[2](24,7) = 1;
+ restriction[2](26,4) = 1;
+ restriction[3](3,3) = 1;
+ restriction[3](10,2) = 1;
+ restriction[3](11,0) = 1;
+ restriction[3](19,7) = 1;
+ restriction[3](20,1) = 1;
+ restriction[3](24,6) = 1;
+ restriction[3](25,4) = 1;
+ restriction[3](26,5) = 1;
+ restriction[4](4,4) = 1;
+ restriction[4](12,5) = 1;
+ restriction[4](15,7) = 1;
+ restriction[4](16,0) = 1;
+ restriction[4](21,6) = 1;
+ restriction[4](22,1) = 1;
+ restriction[4](25,3) = 1;
+ restriction[4](26,2) = 1;
+ restriction[5](5,5) = 1;
+ restriction[5](12,4) = 1;
+ restriction[5](13,6) = 1;
+ restriction[5](17,1) = 1;
+ restriction[5](21,7) = 1;
+ restriction[5](22,0) = 1;
+ restriction[5](23,2) = 1;
+ restriction[5](26,3) = 1;
+ restriction[6](6,6) = 1;
+ restriction[6](13,5) = 1;
+ restriction[6](14,7) = 1;
+ restriction[6](18,2) = 1;
+ restriction[6](21,4) = 1;
+ restriction[6](23,1) = 1;
+ restriction[6](24,3) = 1;
+ restriction[6](26,0) = 1;
+ restriction[7](7,7) = 1;
+ restriction[7](14,6) = 1;
+ restriction[7](15,4) = 1;
+ restriction[7](19,3) = 1;
+ restriction[7](21,5) = 1;
+ restriction[7](24,2) = 1;
+ restriction[7](25,0) = 1;
+ restriction[7](26,1) = 1;
+ break;
+ case 3:
+ restriction[0](0,0) = 1;
+ restriction[0](8,9) = 1;
+ restriction[0](14,15) = 1;
+ restriction[0](24,25) = 1;
+ restriction[0](32,35) = 1;
+ restriction[0](40,43) = 1;
+ restriction[0](52,55) = 1;
+ restriction[0](56,63) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](9,8) = 1;
+ restriction[1](10,11) = 1;
+ restriction[1](26,27) = 1;
+ restriction[1](33,34) = 1;
+ restriction[1](41,42) = 1;
+ restriction[1](44,47) = 1;
+ restriction[1](57,62) = 1;
+ restriction[2](2,2) = 1;
+ restriction[2](11,10) = 1;
+ restriction[2](13,12) = 1;
+ restriction[2](28,29) = 1;
+ restriction[2](35,32) = 1;
+ restriction[2](46,45) = 1;
+ restriction[2](49,50) = 1;
+ restriction[2](61,58) = 1;
+ restriction[3](3,3) = 1;
+ restriction[3](12,13) = 1;
+ restriction[3](15,14) = 1;
+ restriction[3](30,31) = 1;
+ restriction[3](34,33) = 1;
+ restriction[3](48,51) = 1;
+ restriction[3](54,53) = 1;
+ restriction[3](60,59) = 1;
+ restriction[4](4,4) = 1;
+ restriction[4](16,17) = 1;
+ restriction[4](22,23) = 1;
+ restriction[4](25,24) = 1;
+ restriction[4](36,39) = 1;
+ restriction[4](42,41) = 1;
+ restriction[4](53,54) = 1;
+ restriction[4](58,61) = 1;
+ restriction[5](5,5) = 1;
+ restriction[5](17,16) = 1;
+ restriction[5](18,19) = 1;
+ restriction[5](27,26) = 1;
+ restriction[5](37,38) = 1;
+ restriction[5](43,40) = 1;
+ restriction[5](45,46) = 1;
+ restriction[5](59,60) = 1;
+ restriction[6](6,6) = 1;
+ restriction[6](19,18) = 1;
+ restriction[6](21,20) = 1;
+ restriction[6](29,28) = 1;
+ restriction[6](39,36) = 1;
+ restriction[6](47,44) = 1;
+ restriction[6](51,48) = 1;
+ restriction[6](63,56) = 1;
+ restriction[7](7,7) = 1;
+ restriction[7](20,21) = 1;
+ restriction[7](23,22) = 1;
+ restriction[7](31,30) = 1;
+ restriction[7](38,37) = 1;
+ restriction[7](50,49) = 1;
+ restriction[7](55,52) = 1;
+ restriction[7](62,57) = 1;
+ break;
+ case 4:
+ restriction[0](0,0) = 1;
+ restriction[0](8,9) = 1;
+ restriction[0](9,1) = 1;
+ restriction[0](17,18) = 1;
+ restriction[0](18,3) = 1;
+ restriction[0](32,33) = 1;
+ restriction[0](33,4) = 1;
+ restriction[0](44,48) = 1;
+ restriction[0](45,12) = 1;
+ restriction[0](47,15) = 1;
+ restriction[0](48,2) = 1;
+ restriction[0](62,66) = 1;
+ restriction[0](63,36) = 1;
+ restriction[0](65,21) = 1;
+ restriction[0](66,5) = 1;
+ restriction[0](89,93) = 1;
+ restriction[0](90,30) = 1;
+ restriction[0](92,42) = 1;
+ restriction[0](93,7) = 1;
+ restriction[0](98,111) = 1;
+ restriction[0](99,75) = 1;
+ restriction[0](101,57) = 1;
+ restriction[0](102,24) = 1;
+ restriction[0](107,84) = 1;
+ restriction[0](108,39) = 1;
+ restriction[0](110,27) = 1;
+ restriction[0](111,6) = 1;
+ restriction[1](1,1) = 1;
+ restriction[1](9,0) = 1;
+ restriction[1](10,9) = 1;
+ restriction[1](11,12) = 1;
+ restriction[1](12,2) = 1;
+ restriction[1](35,36) = 1;
+ restriction[1](36,5) = 1;
+ restriction[1](45,18) = 1;
+ restriction[1](46,48) = 1;
+ restriction[1](48,3) = 1;
+ restriction[1](49,15) = 1;
+ restriction[1](63,33) = 1;
+ restriction[1](64,66) = 1;
+ restriction[1](66,4) = 1;
+ restriction[1](67,21) = 1;
+ restriction[1](71,75) = 1;
+ restriction[1](72,24) = 1;
+ restriction[1](74,39) = 1;
+ restriction[1](75,6) = 1;
+ restriction[1](99,93) = 1;
+ restriction[1](100,111) = 1;
+ restriction[1](102,30) = 1;
+ restriction[1](103,57) = 1;
+ restriction[1](108,42) = 1;
+ restriction[1](109,84) = 1;
+ restriction[1](111,7) = 1;
+ restriction[1](112,27) = 1;
+ restriction[2](2,2) = 1;
+ restriction[2](12,1) = 1;
+ restriction[2](13,12) = 1;
+ restriction[2](15,3) = 1;
+ restriction[2](16,15) = 1;
+ restriction[2](38,39) = 1;
+ restriction[2](39,6) = 1;
+ restriction[2](48,0) = 1;
+ restriction[2](49,9) = 1;
+ restriction[2](51,18) = 1;
+ restriction[2](52,48) = 1;
+ restriction[2](74,36) = 1;
+ restriction[2](75,5) = 1;
+ restriction[2](77,75) = 1;
+ restriction[2](78,24) = 1;
+ restriction[2](81,42) = 1;
+ restriction[2](82,84) = 1;
+ restriction[2](84,7) = 1;
+ restriction[2](85,27) = 1;
+ restriction[2](108,33) = 1;
+ restriction[2](109,66) = 1;
+ restriction[2](111,4) = 1;
+ restriction[2](112,21) = 1;
+ restriction[2](117,93) = 1;
+ restriction[2](118,111) = 1;
+ restriction[2](120,30) = 1;
+ restriction[2](121,57) = 1;
+ restriction[3](3,3) = 1;
+ restriction[3](14,15) = 1;
+ restriction[3](15,2) = 1;
+ restriction[3](18,0) = 1;
+ restriction[3](19,18) = 1;
+ restriction[3](41,42) = 1;
+ restriction[3](42,7) = 1;
+ restriction[3](47,9) = 1;
+ restriction[3](48,1) = 1;
+ restriction[3](50,48) = 1;
+ restriction[3](51,12) = 1;
+ restriction[3](80,84) = 1;
+ restriction[3](81,39) = 1;
+ restriction[3](83,27) = 1;
+ restriction[3](84,6) = 1;
+ restriction[3](92,33) = 1;
+ restriction[3](93,4) = 1;
+ restriction[3](95,93) = 1;
+ restriction[3](96,30) = 1;
+ restriction[3](107,66) = 1;
+ restriction[3](108,36) = 1;
+ restriction[3](110,21) = 1;
+ restriction[3](111,5) = 1;
+ restriction[3](116,111) = 1;
+ restriction[3](117,75) = 1;
+ restriction[3](119,57) = 1;
+ restriction[3](120,24) = 1;
+ restriction[4](4,4) = 1;
+ restriction[4](20,21) = 1;
+ restriction[4](21,5) = 1;
+ restriction[4](29,30) = 1;
+ restriction[4](30,7) = 1;
+ restriction[4](33,0) = 1;
+ restriction[4](34,33) = 1;
+ restriction[4](53,57) = 1;
+ restriction[4](54,24) = 1;
+ restriction[4](56,27) = 1;
+ restriction[4](57,6) = 1;
+ restriction[4](65,9) = 1;
+ restriction[4](66,1) = 1;
+ restriction[4](68,66) = 1;
+ restriction[4](69,36) = 1;
+ restriction[4](90,18) = 1;
+ restriction[4](91,93) = 1;
+ restriction[4](93,3) = 1;
+ restriction[4](94,42) = 1;
+ restriction[4](101,48) = 1;
+ restriction[4](102,12) = 1;
+ restriction[4](104,111) = 1;
+ restriction[4](105,75) = 1;
+ restriction[4](110,15) = 1;
+ restriction[4](111,2) = 1;
+ restriction[4](113,84) = 1;
+ restriction[4](114,39) = 1;
+ restriction[5](5,5) = 1;
+ restriction[5](21,4) = 1;
+ restriction[5](22,21) = 1;
+ restriction[5](23,24) = 1;
+ restriction[5](24,6) = 1;
+ restriction[5](36,1) = 1;
+ restriction[5](37,36) = 1;
+ restriction[5](54,30) = 1;
+ restriction[5](55,57) = 1;
+ restriction[5](57,7) = 1;
+ restriction[5](58,27) = 1;
+ restriction[5](66,0) = 1;
+ restriction[5](67,9) = 1;
+ restriction[5](69,33) = 1;
+ restriction[5](70,66) = 1;
+ restriction[5](72,12) = 1;
+ restriction[5](73,75) = 1;
+ restriction[5](75,2) = 1;
+ restriction[5](76,39) = 1;
+ restriction[5](102,18) = 1;
+ restriction[5](103,48) = 1;
+ restriction[5](105,93) = 1;
+ restriction[5](106,111) = 1;
+ restriction[5](111,3) = 1;
+ restriction[5](112,15) = 1;
+ restriction[5](114,42) = 1;
+ restriction[5](115,84) = 1;
+ restriction[6](6,6) = 1;
+ restriction[6](24,5) = 1;
+ restriction[6](25,24) = 1;
+ restriction[6](27,7) = 1;
+ restriction[6](28,27) = 1;
+ restriction[6](39,2) = 1;
+ restriction[6](40,39) = 1;
+ restriction[6](57,4) = 1;
+ restriction[6](58,21) = 1;
+ restriction[6](60,30) = 1;
+ restriction[6](61,57) = 1;
+ restriction[6](75,1) = 1;
+ restriction[6](76,36) = 1;
+ restriction[6](78,12) = 1;
+ restriction[6](79,75) = 1;
+ restriction[6](84,3) = 1;
+ restriction[6](85,15) = 1;
+ restriction[6](87,42) = 1;
+ restriction[6](88,84) = 1;
+ restriction[6](111,0) = 1;
+ restriction[6](112,9) = 1;
+ restriction[6](114,33) = 1;
+ restriction[6](115,66) = 1;
+ restriction[6](120,18) = 1;
+ restriction[6](121,48) = 1;
+ restriction[6](123,93) = 1;
+ restriction[6](124,111) = 1;
+ restriction[7](7,7) = 1;
+ restriction[7](26,27) = 1;
+ restriction[7](27,6) = 1;
+ restriction[7](30,4) = 1;
+ restriction[7](31,30) = 1;
+ restriction[7](42,3) = 1;
+ restriction[7](43,42) = 1;
+ restriction[7](56,21) = 1;
+ restriction[7](57,5) = 1;
+ restriction[7](59,57) = 1;
+ restriction[7](60,24) = 1;
+ restriction[7](83,15) = 1;
+ restriction[7](84,2) = 1;
+ restriction[7](86,84) = 1;
+ restriction[7](87,39) = 1;
+ restriction[7](93,0) = 1;
+ restriction[7](94,33) = 1;
+ restriction[7](96,18) = 1;
+ restriction[7](97,93) = 1;
+ restriction[7](110,9) = 1;
+ restriction[7](111,1) = 1;
+ restriction[7](113,66) = 1;
+ restriction[7](114,36) = 1;
+ restriction[7](119,48) = 1;
+ restriction[7](120,12) = 1;
+ restriction[7](122,111) = 1;
+ restriction[7](123,75) = 1;
+ break;
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ restriction[i].reinit(0,0);
+ };
+ }
+ break;
+ };
+
default:
Assert (false,ExcNotImplemented());
}
-template <int dim>
-FE_Q<dim>::~FE_Q ()
-{
- // delete poly member and set it to
- // zero to prevent accidental use
- delete poly;
- poly = 0;
-}
-
-
-
template <int dim>
FiniteElement<dim> *
FE_Q<dim>::clone() const
FE_Q<dim>::shape_value (const unsigned int i,
const Point<dim> &p) const
{
- return poly->compute_value(renumber_inverse[i], p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_value(renumber_inverse[i], p);
+}
+
+
+template <int dim>
+double
+FE_Q<dim>::shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_value(renumber_inverse[i], p);
}
FE_Q<dim>::shape_grad (const unsigned int i,
const Point<dim> &p) const
{
- return poly->compute_grad(renumber_inverse[i], p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_grad(renumber_inverse[i], p);
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FE_Q<dim>::shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_grad(renumber_inverse[i], p);
}
FE_Q<dim>::shape_grad_grad (const unsigned int i,
const Point<dim> &p) const
{
- return poly->compute_grad_grad(renumber_inverse[i], p);
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_Q<dim>::shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i,0,dofs_per_cell));
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
}
UpdateFlags
FE_Q<dim>::update_once (const UpdateFlags flags) const
{
- UpdateFlags out = update_default;
-
- if (flags & update_values)
- out |= update_values;
-
- return out;
+ // for this kind of elements, only
+ // the values can be precomputed
+ // once and for all. set this flag
+ // if the values are requested at
+ // all
+ return (update_default | (flags & update_values));
}
const Mapping<dim> &mapping,
const Quadrature<dim> &quadrature) const
{
+ // generate a new data object and
+ // initialize some fields
InternalData* data = new InternalData;
- std::vector<double> values(0);
- std::vector<Tensor<1,dim> > grads(0);
- std::vector<Tensor<2,dim> > grad_grads(0);
// check what needs to be
// initialized only once and what
const UpdateFlags flags(data->update_flags);
const unsigned int n_q_points = quadrature.n_quadrature_points;
-
+
+ // some scratch arrays
+ std::vector<double> values(0);
+ std::vector<Tensor<1,dim> > grads(0);
+ std::vector<Tensor<2,dim> > grad_grads(0);
+
+ // initialize fields only if really
+ // necessary. otherwise, don't
+ // allocate memory
if (flags & update_values)
{
values.resize (dofs_per_cell);
// that
if (flags & update_second_derivatives)
data->initialize_2nd (this, mapping, quadrature);
-
+
+ // next already fill those fields
+ // of which we have information by
+ // now. note that the shape
+ // gradients are only those on the
+ // unit cell, and need to be
+ // transformed when visiting an
+ // actual cell
if (flags & (update_values | update_gradients))
for (unsigned int i=0; i<n_q_points; ++i)
{
- poly->compute(quadrature.point(i), values, grads, grad_grads);
- for (unsigned int k=0; k<dofs_per_cell; ++k)
- {
- if (flags & update_values)
- data->shape_values[renumber[k]][i] = values[k];
- if (flags & update_gradients)
- data->shape_gradients[renumber[k]][i] = grads[k];
- }
+ polynomial_space.compute(quadrature.point(i),
+ values, grads, grad_grads);
+
+ if (flags & update_values)
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ data->shape_values[renumber[k]][i] = values[k];
+
+ if (flags & update_gradients)
+ for (unsigned int k=0; k<dofs_per_cell; ++k)
+ data->shape_gradients[renumber[k]][i] = grads[k];
}
return data;
}
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
- for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
- if (flags & update_values)
+ if (flags & update_values)
+ for (unsigned int i=0; i<quadrature.n_quadrature_points; ++i)
data.shape_values(k,i) = fe_data.shape_values[k][i];
if (flags & update_gradients)
FE_Q<dim>::has_support_on_face (const unsigned int shape_index_,
const unsigned int face_index) const
{
+ Assert (shape_index_ < dofs_per_cell,
+ ExcIndexRange (shape_index_, 0, dofs_per_cell));
+ Assert (face_index < GeometryInfo<dim>::faces_per_cell,
+ ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
+
unsigned int shape_index = shape_index_;
if (dim==1)
template <int dim>
-FESystem<dim>::FESystem (const FiniteElement<dim> &fe, const unsigned int n_elements) :
+FESystem<dim>::FESystem (const FiniteElement<dim> &fe,
+ const unsigned int n_elements) :
FiniteElement<dim> (multiply_dof_numbers(fe, n_elements),
- compute_restriction_is_additive_flags (fe, n_elements)),
+ compute_restriction_is_additive_flags (fe, n_elements),
+ compute_nonzero_components(fe, n_elements)),
base_elements(1)
{
base_elements[0] = ElementPair(fe.clone(), n_elements);
template <int dim>
-FESystem<dim>::FESystem (const FiniteElement<dim> &fe1, const unsigned int n1,
- const FiniteElement<dim> &fe2, const unsigned int n2) :
+FESystem<dim>::FESystem (const FiniteElement<dim> &fe1,
+ const unsigned int n1,
+ const FiniteElement<dim> &fe2,
+ const unsigned int n2) :
FiniteElement<dim> (multiply_dof_numbers(fe1, n1, fe2, n2),
compute_restriction_is_additive_flags (fe1, n1,
- fe2, n2)),
+ fe2, n2),
+ compute_nonzero_components(fe1, n1,
+ fe2, n2)),
base_elements(2)
{
base_elements[0] = ElementPair(fe1.clone(), n1);
template <int dim>
-FESystem<dim>::FESystem (const FiniteElement<dim> &fe1, const unsigned int n1,
- const FiniteElement<dim> &fe2, const unsigned int n2,
- const FiniteElement<dim> &fe3, const unsigned int n3) :
+FESystem<dim>::FESystem (const FiniteElement<dim> &fe1,
+ const unsigned int n1,
+ const FiniteElement<dim> &fe2,
+ const unsigned int n2,
+ const FiniteElement<dim> &fe3,
+ const unsigned int n3) :
FiniteElement<dim> (multiply_dof_numbers(fe1, n1,
fe2, n2,
fe3, n3),
compute_restriction_is_additive_flags (fe1, n1,
fe2, n2,
- fe3, n3)),
+ fe3, n3),
+ compute_nonzero_components(fe1, n1,
+ fe2, n2,
+ fe3, n3)),
base_elements(3)
{
base_elements[0] = ElementPair(fe1.clone(), n1);
FESystem<dim>::shape_value (const unsigned int i,
const Point<dim> &p) const
{
- Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (is_primitive(i),
+ typename FiniteElementBase<dim>::ExcShapeFunctionNotPrimitive(i));
- std::pair<unsigned,unsigned> comp = system_to_component_index(i);
+ return (base_element(system_to_base_table[i].first.first)
+ .shape_value(system_to_base_table[i].second, p));
+}
+
+
+
+template <int dim>
+double
+FESystem<dim>::shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (component < n_components(),
+ ExcIndexRange (component, 0, n_components()));
- return base_element(component_to_base_table[comp.first])
- .shape_value(comp.second, p);
+ // first find out to which of the
+ // base elements this desired
+ // component belongs, and which
+ // component within this base
+ // element it is
+ const unsigned int base = component_to_base(i).first;
+ const unsigned int component_in_base = component_to_base(i).second;
+
+ // then get value from base
+ // element. note that that will
+ // throw an error should the
+ // respective shape function not be
+ // primitive; thus, there is no
+ // need to check this here
+ return (base_element(base).
+ shape_value_component(system_to_base_table[i].second,
+ p,
+ component_in_base));
}
FESystem<dim>::shape_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (is_primitive(i),
+ typename FiniteElementBase<dim>::ExcShapeFunctionNotPrimitive(i));
- std::pair<unsigned,unsigned> comp = system_to_component_index(i);
+ return (base_element(system_to_base_table[i].first.first)
+ .shape_grad(system_to_base_table[i].second, p));
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FESystem<dim>::shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (component < n_components(),
+ ExcIndexRange (component, 0, n_components()));
- return base_element(component_to_base_table[comp.first])
- .shape_grad(comp.second, p);
+ // first find out to which of the
+ // base elements this desired
+ // component belongs, and which
+ // component within this base
+ // element it is
+ const unsigned int base = component_to_base(i).first;
+ const unsigned int component_in_base = component_to_base(i).second;
+
+ // then get value from base
+ // element. note that that will
+ // throw an error should the
+ // respective shape function not be
+ // primitive; thus, there is no
+ // need to check this here
+ return (base_element(base).
+ shape_grad_component(system_to_base_table[i].second,
+ p,
+ component_in_base));
}
FESystem<dim>::shape_grad_grad (const unsigned int i,
const Point<dim> &p) const
{
- Assert((i<dofs_per_cell), ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (is_primitive(i),
+ typename FiniteElementBase<dim>::ExcShapeFunctionNotPrimitive(i));
+
+ return (base_element(system_to_base_table[i].first.first)
+ .shape_grad_grad(system_to_base_table[i].second, p));
+}
+
- std::pair<unsigned,unsigned> comp = system_to_component_index(i);
+template <int dim>
+Tensor<2,dim>
+FESystem<dim>::shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<dofs_per_cell, ExcIndexRange(i, 0, dofs_per_cell));
+ Assert (component < n_components(),
+ ExcIndexRange (component, 0, n_components()));
- return base_element(component_to_base_table[comp.first])
- .shape_grad_grad(comp.second, p);
+ // first find out to which of the
+ // base elements this desired
+ // component belongs, and which
+ // component within this base
+ // element it is
+ const unsigned int base = component_to_base(i).first;
+ const unsigned int component_in_base = component_to_base(i).second;
+
+ // then get value from base
+ // element. note that that will
+ // throw an error should the
+ // respective shape function not be
+ // primitive; thus, there is no
+ // need to check this here
+ return (base_element(base).
+ shape_grad_grad_component(system_to_base_table[i].second,
+ p,
+ component_in_base));
}
UpdateFlags sub_flags = flags;
// if second derivatives through
- // finite differencing is required,
+ // finite differencing are required,
// then initialize some objects for
// that
data->compute_second_derivatives = flags & update_second_derivatives;
FEValuesData<dim> *base_data = new FEValuesData<dim>();
data->set_fe_values_data(base_no, base_data);
}
- data->update_flags=data->update_once | data->update_each;
+ data->update_flags = data->update_once |
+ data->update_each;
return data;
}
FEValuesData<dim> &base_data=fe_data.get_fe_values_data(base_no);
const FiniteElement<dim> &base_fe=base_element(base_no);
base_data.initialize(quadrature.n_quadrature_points,
- base_fe.dofs_per_cell,
+ base_fe,
base_update_flags);
}
}
template <int dim>
void
-FESystem<dim>::build_cell_table()
+FESystem<dim>::build_cell_tables()
{
unsigned total_index = 0;
- for (unsigned int base=0 ; base < n_base_elements() ; ++base)
+ for (unsigned int base=0; base < n_base_elements(); ++base)
for (unsigned int m = 0; m < element_multiplicity(base); ++m)
for (unsigned int k=0; k<base_element(base).n_components(); ++k)
- component_to_base_table[total_index++] = base;
+ component_to_base_table[total_index++] = std::make_pair(base,k);
Assert (total_index == component_to_base_table.size(),
ExcInternalError());
- // Initialize index table
- // Multi-component base elements have
- // to be thought of.
+ // Initialize index tables.
+ // Multi-component base elements
+ // have to be thought of. For
+ // non-primitive shape functions,
+ // have a special invalid index.
+ const std::pair<unsigned int, unsigned int>
+ non_primitive_index (static_cast<unsigned int>(-1),
+ static_cast<unsigned int>(-1));
- // 1. Vertices
+ // First enumerate vertex indices,
+ // where we first enumerate all
+ // indices on the first vertex in
+ // the order of the base elements,
+ // then of the second vertex, etc
total_index = 0;
- for (unsigned int vertex_number= 0 ;
- vertex_number < GeometryInfo<dim>::vertices_per_cell ;
+ for (unsigned int vertex_number=0;
+ vertex_number<GeometryInfo<dim>::vertices_per_cell;
++vertex_number)
{
unsigned comp_start = 0;
- for(unsigned int base = 0; base < n_base_elements() ;
- ++base)
- {
- for (unsigned int m = 0; m < element_multiplicity(base); ++m)
+ for(unsigned int base=0; base<n_base_elements(); ++base)
+ for (unsigned int m=0; m<element_multiplicity(base);
+ ++m, comp_start+=base_element(base).n_components())
+ for (unsigned int local_index = 0;
+ local_index < base_element(base).dofs_per_vertex;
+ ++local_index, ++total_index)
{
- for (unsigned int local_index = 0 ;
- local_index < base_element(base).dofs_per_vertex ;
- ++local_index)
+ const unsigned int index_in_base
+ = (base_element(base).dofs_per_vertex*vertex_number +
+ local_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base, m), index_in_base);
+
+ if (base_element(base).is_primitive(index_in_base))
{
- system_to_component_table[total_index++]
- = std::make_pair (comp_start+m,
- vertex_number*base_element(base).dofs_per_vertex
- +local_index);
+ const unsigned int comp_in_base
+ = base_element(base).system_to_component_index(index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int index_in_comp
+ = base_element(base).system_to_component_index(index_in_base).second;
+ system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
}
+ else
+ system_to_component_table[total_index] = non_primitive_index;
}
- comp_start += element_multiplicity(base);
- }
}
// 2. Lines
- for (unsigned int line_number= 0 ; ((line_number != GeometryInfo<dim>::lines_per_cell) &&
- (GeometryInfo<dim>::lines_per_cell > 0));
- ++line_number)
- {
- unsigned comp_start = 0;
- for(unsigned int base = 0; base < n_base_elements() ;
- ++base)
- {
- for (unsigned int m = 0; m < element_multiplicity(base); ++m)
- {
- for (unsigned int local_index = 0 ;
- local_index < base_element(base).dofs_per_line ;
- ++local_index)
- {
- system_to_component_table[total_index++]
- = std::pair<unsigned,unsigned>
- (comp_start+m,
- line_number*base_element(base).dofs_per_line
- +local_index+base_element(base).first_line_index);
- }
- }
- comp_start += element_multiplicity(base);
- }
- }
+ if (GeometryInfo<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ unsigned comp_start = 0;
+ for (unsigned int base=0; base<n_base_elements(); ++base)
+ for (unsigned int m=0; m<element_multiplicity(base);
+ ++m, comp_start+=base_element(base).n_components())
+ for (unsigned int local_index = 0;
+ local_index < base_element(base).dofs_per_line;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (base_element(base).dofs_per_line*line_number +
+ local_index +
+ base_element(base).first_line_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
+
+ if (base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = base_element(base).system_to_component_index(index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int index_in_comp
+ = base_element(base).system_to_component_index(index_in_base).second;
+ system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
+ }
+ else
+ system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
// 3. Quads
- for (unsigned int quad_number= 0 ;
- ((quad_number != GeometryInfo<dim>::quads_per_cell) &&
- (GeometryInfo<dim>::quads_per_cell > 0));
- ++quad_number)
- {
- unsigned int comp_start = 0;
- for(unsigned int base = 0; base < n_base_elements() ;
- ++base)
- {
- for (unsigned int m = 0; m < element_multiplicity(base); ++m)
- {
- for (unsigned int local_index = 0 ;
- local_index < base_element(base).dofs_per_quad ;
- ++local_index)
- {
- system_to_component_table[total_index++]
- = std::make_pair (comp_start+m,
- quad_number*base_element(base).dofs_per_quad
- +local_index+base_element(base).first_quad_index);
- }
- }
- comp_start += element_multiplicity(base);
- }
- }
+ if (GeometryInfo<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<n_base_elements(); ++base)
+ for (unsigned int m=0; m<element_multiplicity(base);
+ ++m, comp_start += base_element(base).n_components())
+ for (unsigned int local_index = 0;
+ local_index < base_element(base).dofs_per_quad;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (base_element(base).dofs_per_quad*quad_number +
+ local_index +
+ base_element(base).first_quad_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
+
+ if (base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = base_element(base).system_to_component_index(index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int index_in_comp
+ = base_element(base).system_to_component_index(index_in_base).second;
+ system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
+ }
+ else
+ system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
- // 4. Hex
- for (unsigned int hex_number= 0 ;
- ((hex_number != GeometryInfo<dim>::hexes_per_cell) &&
- (GeometryInfo<dim>::hexes_per_cell > 0));
- ++hex_number)
+ // 4. Hexes
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ unsigned int comp_start = 0;
+ for(unsigned int base=0; base<n_base_elements(); ++base)
+ for (unsigned int m=0; m<element_multiplicity(base);
+ ++m, comp_start+=base_element(base).n_components())
+ for (unsigned int local_index = 0;
+ local_index < base_element(base).dofs_per_hex;
+ ++local_index, ++total_index)
+ {
+ const unsigned int index_in_base
+ = (base_element(base).dofs_per_hex*hex_number +
+ local_index +
+ base_element(base).first_hex_index);
+
+ system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), index_in_base);
+
+ if (base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = base_element(base).system_to_component_index(index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int index_in_comp
+ = base_element(base).system_to_component_index(index_in_base).second;
+ system_to_component_table[total_index]
+ = std::make_pair (comp, index_in_comp);
+ }
+ else
+ system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+
+ // Initialize mapping from
+ // components to linear
+ // index. Fortunately, this is the
+ // inverse of what we just did.
+ std::vector<unsigned int> dofs_per_component (n_components(), 0);
+ for (unsigned int sys=0; sys<dofs_per_cell; ++sys)
+ ++dofs_per_component[system_to_component_index(sys).first];
+ for (unsigned int component=0; component<n_components(); ++component)
+ component_to_system_table[component].resize(dofs_per_component[component]);
+
+ // then go the reverse way to fill the array
+ for (unsigned int sys=0; sys<dofs_per_cell; ++sys)
{
- unsigned int comp_start = 0;
- for(unsigned int base = 0; base < n_base_elements() ;
- ++base)
- {
- for (unsigned int m = 0; m < element_multiplicity(base); ++m)
- {
- for (unsigned int local_index = 0 ;
- local_index < base_element(base).dofs_per_hex ;
- ++local_index)
- {
- system_to_component_table[total_index++]
- = std::make_pair (comp_start+m,
- hex_number*base_element(base).dofs_per_hex
- +local_index+base_element(base).first_hex_index);
- }
- }
- comp_start += element_multiplicity(base);
-
- }
- }
- // Initialize mapping from component
- // to base element
- // Initialize mapping from components to
- // linear index. Fortunately, this is
- // the inverse of what we just did.
- for (unsigned int comp=0 ; comp<n_components() ; ++comp)
- component_to_system_table[comp]
- .resize(base_element(component_to_base_table[comp]).dofs_per_cell);
-
- for (unsigned int sys=0 ; sys < dofs_per_cell ; ++sys)
- component_to_system_table[system_to_component_table[sys].first]
- [system_to_component_table[sys].second] = sys;
-}
+ const unsigned int
+ comp = system_to_component_index(sys).first,
+ index_in_comp = system_to_component_index(sys).second;
+
+ Assert (comp < component_to_system_table.size(),
+ ExcInternalError());
+ Assert (index_in_comp < component_to_system_table[comp].size(),
+ ExcInternalError());
+ component_to_system_table[comp][index_in_comp] = sys;
+ };
+};
template <int dim>
void
-FESystem<dim>::build_face_table()
+FESystem<dim>::build_face_tables()
{
- // Initialize index table
- // Multi-component base elements have
- // to be thought of.
+ // Initialize index tables. do this
+ // in the same way as done for the
+ // cell tables, except that we now
+ // loop over the objects of faces
+
+ // For non-primitive shape
+ // functions, have a special
+ // invalid index
+ const std::pair<unsigned int, unsigned int>
+ non_primitive_index (static_cast<unsigned int>(-1),
+ static_cast<unsigned int>(-1));
// 1. Vertices
unsigned int total_index = 0;
- for (unsigned int vertex_number= 0 ; vertex_number < GeometryInfo<dim>::vertices_per_face ;
+ for (unsigned int vertex_number=0;
+ vertex_number<GeometryInfo<dim>::vertices_per_face;
++vertex_number)
{
unsigned int comp_start = 0;
- for(unsigned int base = 0; base < n_base_elements() ;
- ++base)
- {
- for (unsigned int m = 0; m < element_multiplicity(base); ++m)
+ for(unsigned int base=0; base<n_base_elements(); ++base)
+ for (unsigned int m=0; m<element_multiplicity(base);
+ ++m, comp_start += base_element(base).n_components())
+ for (unsigned int local_index = 0;
+ local_index < base_element(base).dofs_per_vertex;
+ ++local_index, ++total_index)
{
- for (unsigned int local_index = 0 ;
- local_index < base_element(base).dofs_per_vertex ;
- ++local_index)
+ // get (cell) index of
+ // this shape function
+ // inside the base
+ // element to see
+ // whether the shape
+ // function is
+ // primitive (assume
+ // that all shape
+ // functions on
+ // vertices share the
+ // same primitivity
+ // property; assume
+ // likewise for all
+ // shape functions
+ // located on lines,
+ // quads, etc. this
+ // way, we can ask for
+ // primitivity of only
+ // _one_ shape
+ // function, which is
+ // taken as
+ // representative for
+ // all others located
+ // on the same type of
+ // object):
+ const unsigned int index_in_base
+ = (base_element(base).dofs_per_vertex*vertex_number +
+ local_index);
+
+ const unsigned int face_index_in_base
+ = (base_element(base).dofs_per_vertex*vertex_number +
+ local_index);
+
+ face_system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), face_index_in_base);
+
+ if (base_element(base).is_primitive(index_in_base))
{
- face_system_to_component_table[total_index++]
- = std::pair<unsigned,unsigned>
- (comp_start+m,
- vertex_number*base_element(base).dofs_per_vertex+local_index);
+ const unsigned int comp_in_base
+ = base_element(base).face_system_to_component_index(face_index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int face_index_in_comp
+ = base_element(base).face_system_to_component_index(face_index_in_base).second;
+ face_system_to_component_table[total_index]
+ = std::make_pair (comp, face_index_in_comp);
}
+ else
+ face_system_to_component_table[total_index] = non_primitive_index;
}
- comp_start += element_multiplicity(base);
- }
}
- Assert (total_index <= face_system_to_component_table.size(),
- ExcInternalError());
// 2. Lines
- for (unsigned line_number= 0 ; ((line_number != GeometryInfo<dim>::lines_per_face) &&
- (GeometryInfo<dim>::lines_per_cell > 0));
- ++line_number)
- {
- unsigned comp_start = 0;
- for(unsigned base = 0; base < n_base_elements() ;
- ++base)
- {
- for (unsigned m = 0; m < element_multiplicity(base); ++m)
- {
- for (unsigned local_index = 0 ;
- local_index < base_element(base).dofs_per_line ;
- ++local_index)
- {
- face_system_to_component_table[total_index++]
- = std::pair<unsigned,unsigned>
- (comp_start+m,
- line_number*base_element(base).dofs_per_line
- +local_index+base_element(base).first_face_line_index);
- }
- }
- comp_start += element_multiplicity(base);
- }
- }
- Assert (total_index <= face_system_to_component_table.size(),
- ExcInternalError());
+ if (GeometryInfo<dim>::lines_per_face > 0)
+ for (unsigned line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_face;
+ ++line_number)
+ {
+ unsigned comp_start = 0;
+ for(unsigned base = 0; base < n_base_elements(); ++base)
+ for (unsigned m=0; m<element_multiplicity(base);
+ ++m, comp_start += base_element(base).n_components())
+ for (unsigned local_index = 0;
+ local_index < base_element(base).dofs_per_line;
+ ++local_index, ++total_index)
+ {
+ // do everything
+ // alike for this
+ // type of object
+ const unsigned int index_in_base
+ = (base_element(base).dofs_per_line*line_number +
+ local_index +
+ base_element(base).first_line_index);
+
+ const unsigned int face_index_in_base
+ = (base_element(base).first_face_line_index +
+ base_element(base).dofs_per_line * line_number +
+ local_index);
+
+ face_system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), face_index_in_base);
+
+ if (base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = base_element(base).face_system_to_component_index(face_index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int face_index_in_comp
+ = base_element(base).face_system_to_component_index(face_index_in_base).second;
+ face_system_to_component_table[total_index]
+ = std::make_pair (comp, face_index_in_comp);
+ }
+ else
+ face_system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
// 3. Quads
- for (unsigned quad_number= 0 ; ((quad_number != GeometryInfo<dim>::quads_per_face) &&
- (GeometryInfo<dim>::quads_per_cell > 0));
- ++quad_number)
- {
- unsigned comp_start = 0;
- for(unsigned base = 0; base < n_base_elements() ;
- ++base)
- {
- for (unsigned m = 0; m < element_multiplicity(base); ++m)
- {
- for (unsigned local_index = 0 ;
- local_index < base_element(base).dofs_per_quad ;
- ++local_index)
- {
- face_system_to_component_table[total_index++]
- = std::pair<unsigned,unsigned>
- (comp_start+m,
- quad_number*base_element(base).dofs_per_quad
- +local_index+base_element(base).first_face_quad_index);
- }
- }
- comp_start += element_multiplicity(base);
- }
- }
- Assert (total_index <= face_system_to_component_table.size(),
+ if (GeometryInfo<dim>::quads_per_face > 0)
+ for (unsigned quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_face;
+ ++quad_number)
+ {
+ unsigned comp_start = 0;
+ for(unsigned base=0; base<n_base_elements(); ++base)
+ for (unsigned m=0; m<element_multiplicity(base);
+ ++m, comp_start += base_element(base).n_components())
+ for (unsigned local_index = 0;
+ local_index < base_element(base).dofs_per_quad;
+ ++local_index, ++total_index)
+ {
+ // do everything
+ // alike for this
+ // type of object
+ const unsigned int index_in_base
+ = (base_element(base).dofs_per_quad*quad_number +
+ local_index +
+ base_element(base).first_quad_index);
+
+ const unsigned int face_index_in_base
+ = (base_element(base).first_face_quad_index +
+ base_element(base).dofs_per_quad * quad_number +
+ local_index);
+
+ face_system_to_base_table[total_index]
+ = std::make_pair (std::make_pair(base,m), face_index_in_base);
+
+ if (base_element(base).is_primitive(index_in_base))
+ {
+ const unsigned int comp_in_base
+ = base_element(base).face_system_to_component_index(face_index_in_base).first;
+ const unsigned int comp
+ = comp_start + comp_in_base;
+ const unsigned int face_index_in_comp
+ = base_element(base).face_system_to_component_index(face_index_in_base).second;
+ face_system_to_component_table[total_index]
+ = std::make_pair (comp, face_index_in_comp);
+ }
+ else
+ face_system_to_component_table[total_index] = non_primitive_index;
+ }
+ }
+ Assert (total_index == dofs_per_face, ExcInternalError());
+ Assert (total_index == face_system_to_component_table.size(),
ExcInternalError());
-
- // Initialize mapping from component
- // to base element
- // Initialize mapping from components to
- // linear index. Fortunately, this is
- // the inverse of what we just did.
- for (unsigned comp=0 ; comp<n_components() ; ++comp)
- face_component_to_system_table[comp]
- .resize(base_element(component_to_base_table[comp]).dofs_per_cell);
-
- for (unsigned sys=0 ; sys < dofs_per_face ; ++sys)
- face_component_to_system_table[face_system_to_component_table[sys].first]
- [face_system_to_component_table[sys].second] = sys;
-}
+ Assert (total_index == face_system_to_base_table.size(),
+ ExcInternalError());
+
+ // finally, initialize reverse mapping
+ std::vector<unsigned int> dofs_per_component (n_components(), 0);
+ for (unsigned int sys=0; sys<dofs_per_face; ++sys)
+ ++dofs_per_component[face_system_to_component_index(sys).first];
+ for (unsigned int component=0; component<n_components(); ++component)
+ face_component_to_system_table[component].resize(dofs_per_component[component]);
+
+ // then go the reverse way to fill
+ // the array
+ for (unsigned int sys=0; sys<dofs_per_face; ++sys)
+ {
+ const unsigned int
+ comp = face_system_to_component_index(sys).first,
+ index_in_comp = face_system_to_component_index(sys).second;
+
+ Assert (comp < face_component_to_system_table.size(),
+ ExcInternalError());
+ Assert (index_in_comp < face_component_to_system_table[comp].size(),
+ ExcInternalError());
+ face_component_to_system_table[comp][index_in_comp] = sys;
+ };
+};
template <int dim>
void FESystem<dim>::build_interface_constraints ()
{
- // the layout of the constraints matrix is
- // described in the FiniteElement class. you
- // may want to look there first before trying
- // to understand the following, especially
- // the mapping of the @p{n} index.
+//TODO: Review for base elements with more than one component. do it in the same way as for building the cell restriction and interpolation tables, i.e. top-down
+
+ // the layout of the constraints
+ // matrix is described in the
+ // FiniteElement class. you may
+ // want to look there first before
+ // trying to understand the
+ // following, especially the
+ // mapping of the @p{m} index.
//
- // in order to map it to the fe-system class,
- // we have to know which base element a
- // degree of freedom within a vertex, line,
- // etc belongs to. this can be accomplished
- // by the system_to_component_index
+ // in order to map it to the
+ // fe-system class, we have to know
+ // which base element a degree of
+ // freedom within a vertex, line,
+ // etc belongs to. this can be
+ // accomplished by the
+ // system_to_component_index
// function in conjunction with the
- // numbers first_{line,quad,...}_index
+ // numbers
+ // first_{line,quad,...}_index
for (unsigned int n=0; n<interface_constraints.n(); ++n)
for (unsigned int m=0; m<interface_constraints.m(); ++m)
{
- // for the pair (n,m) find out which
- // component they belong to and
- // the number therein
+ // for the pair (n,m) find
+ // out which base element
+ // they belong to and the
+ // number therein
//
- // first value in pair is component,
- // second is index
- const std::pair<unsigned int, unsigned int> n_index
- = face_system_to_component_index (n);
-
- std::pair<unsigned int, unsigned int> m_index;
+ // first for the n
+ // index. this is simple
+ // since the n indices are in
+ // the same order as they are
+ // usually on a face. note
+ // that for the data type,
+ // first value in pair is
+ // (base element,instance of
+ // base element), second is
+ // index within this instance
+ const std::pair<std::pair<unsigned int,unsigned int>, unsigned int> n_index
+ = face_system_to_base_table[n];
+
+ // likewise for the m
+ // index. this is more
+ // complicated due to the
+ // strange ordering we have
+ // for the dofs on the
+ // refined faces.
+ std::pair<std::pair<unsigned int,unsigned int>, unsigned int> m_index;
switch (dim)
{
case 1:
// as for the first vertex
// of the whole cell, so we
// can use the
- // system_to_component_index
- // function (using the
- // face_s_t_c_i function would
+ // system_to_base_table
+ // variable (using the
+ // face_s_t_base_t function would
// yield the same)
if (m < dofs_per_vertex)
- m_index = system_to_component_index (m);
+ m_index = system_to_base_table[m];
else
// then come the two sets of
// line indices
const unsigned int sub_line
= (m-dofs_per_vertex) / dofs_per_line;
Assert (sub_line < 2, ExcInternalError());
-
- // get the component by
- // asking s_t_c_index and
- // tweaking the index a bit
- m_index.first = system_to_component_index
- (GeometryInfo<2>::vertices_per_cell * dofs_per_vertex
- + index_in_line).first;
- // first find out the how-many'th
- // line index of that component
- // this was
- m_index.second = (system_to_component_index
- (GeometryInfo<2>::vertices_per_cell * dofs_per_vertex
- + index_in_line).second
- - base_element (component_to_base_table[m_index.first]).first_line_index)
- // then add the number of dofs
- // per vertex to get the index
- // on the first line
- + base_element(component_to_base_table[m_index.first]).dofs_per_vertex
- // if on the second line: add
- // some more
- + base_element(component_to_base_table[m_index.first]).dofs_per_line * sub_line;
+
+ // from this
+ // information, try
+ // to get base
+ // element and
+ // instance of base
+ // element. we do
+ // so by
+ // constructing the
+ // corresponding
+ // face index of m
+ // in the present
+ // element, then
+ // use
+ // face_system_to_base_table
+ const unsigned int tmp1 = 2*dofs_per_vertex+index_in_line;
+ m_index.first = face_system_to_base_table[tmp1].first;
+
+ // what we are
+ // still missing is
+ // the index of m
+ // within the base
+ // elements
+ // interface_constraints
+ // table
+ //
+ // here, the second
+ // value of
+ // face_system_to_base_table
+ // can help: it
+ // denotes the face
+ // index of that
+ // shape function
+ // within the base
+ // element. since
+ // we know that it
+ // is a line dof,
+ // we can construct
+ // the rest: tmp2
+ // will denote the
+ // index of this
+ // shape function
+ // among the line
+ // shape functions:
+ Assert (face_system_to_base_table[tmp1].second >=
+ 2*base_element(m_index.first.first).dofs_per_vertex,
+ ExcInternalError());
+ const unsigned int tmp2 = face_system_to_base_table[tmp1].second -
+ 2*base_element(m_index.first.first).dofs_per_vertex;
+ Assert (tmp2 < base_element(m_index.first.first).dofs_per_line,
+ ExcInternalError());
+ m_index.second = base_element(m_index.first.first).dofs_per_vertex +
+ base_element(m_index.first.first).dofs_per_line*sub_line +
+ tmp2;
};
break;
};
case 3:
{
+ // same way as above,
+ // although a little
+ // more complicated...
+
// the indices
// m=0..5*d_v-1 are
// from the center and
// so we can use the
// simple arithmetic
if (m < 5*dofs_per_vertex)
- {
- m_index.first = system_to_component_index(m % dofs_per_vertex).first;
- m_index.second = m / dofs_per_vertex;
- }
+ m_index = system_to_base_table[m];
else
// then come the 12 sets of
// line indices
if (m < 5*dofs_per_vertex + 12*dofs_per_line)
- {
+ {
+ // for the
+ // meaning of all
+ // this, see the
+ // 2d part
const unsigned int index_in_line
= (m-5*dofs_per_vertex) % dofs_per_line;
const unsigned int sub_line
= (m-5*dofs_per_vertex) / dofs_per_line;
Assert (sub_line < 12, ExcInternalError());
-
- // get the component by
- // asking s_t_c_index and
- // tweaking the index a bit
- m_index.first = system_to_component_index
- (GeometryInfo<3>::vertices_per_cell * dofs_per_vertex
- + index_in_line).first;
-
- // first find out the how-many'th
- // line index of that component
- // this was
- m_index.second = (system_to_component_index
- (GeometryInfo<3>::vertices_per_cell * dofs_per_vertex
- + index_in_line).second
- - base_element (component_to_base_table[m_index.first]).first_line_index)
- // then add the number of dofs
- // for the five vertices to get
- // the index on the first line
- + 5*base_element(component_to_base_table[m_index.first]).dofs_per_vertex
- // and correct for the
- // how-many'th line
- + base_element(component_to_base_table[m_index.first]).dofs_per_line * sub_line;
+
+ const unsigned int tmp1 = 4*dofs_per_vertex+index_in_line;
+ m_index.first = face_system_to_base_table[tmp1].first;
+
+ Assert (face_system_to_base_table[tmp1].second >=
+ 4*base_element(m_index.first.first).dofs_per_vertex,
+ ExcInternalError());
+ const unsigned int tmp2 = face_system_to_base_table[tmp1].second -
+ 4*base_element(m_index.first.first).dofs_per_vertex;
+ Assert (tmp2 < base_element(m_index.first.first).dofs_per_line,
+ ExcInternalError());
+ m_index.second = base_element(m_index.first.first).dofs_per_vertex +
+ base_element(m_index.first.first).dofs_per_line*sub_line +
+ tmp2;
}
else
// on one of the four sub-quads
{
+ // for the
+ // meaning of all
+ // this, see the
+ // 2d part
const unsigned int index_in_quad
= (m-5*dofs_per_vertex-12*dofs_per_line) % dofs_per_line;
const unsigned int sub_quad
= (m-5*dofs_per_vertex-12*dofs_per_line) / dofs_per_line;
Assert (sub_quad < 4, ExcInternalError());
-
- // get the component by
- // asking s_t_c_index and
- // tweaking the index a bit
- m_index.first = system_to_component_index
- (GeometryInfo<3>::vertices_per_cell * dofs_per_vertex
- + GeometryInfo<3>::lines_per_cell * dofs_per_line
- + index_in_quad).first;
-
- // first find out the how-many'th
- // quad index of that component
- // this was
- m_index.second = (system_to_component_index
- (GeometryInfo<3>::vertices_per_cell * dofs_per_vertex
- + GeometryInfo<3>::lines_per_cell * dofs_per_line
- + index_in_quad).second
- - base_element (component_to_base_table[m_index.first]).first_quad_index)
- // then add the number of dofs
- // for the five vertices and 12 lines
- // to get the index on the first quad
- + 5*base_element(component_to_base_table[m_index.first]).dofs_per_vertex
- + 12*base_element(component_to_base_table[m_index.first]).dofs_per_line
- // and correct for the
- // how-many'th line
- + base_element(component_to_base_table[m_index.first]).dofs_per_quad * sub_quad;
+
+ const unsigned int tmp1 = 4*dofs_per_vertex+4*dofs_per_line+index_in_quad;
+ m_index.first = face_system_to_base_table[tmp1].first;
+
+ Assert (face_system_to_base_table[tmp1].second >=
+ 4*base_element(m_index.first.first).dofs_per_vertex +
+ 4*base_element(m_index.first.first).dofs_per_line,
+ ExcInternalError());
+ const unsigned int tmp2 = face_system_to_base_table[tmp1].second -
+ 4*base_element(m_index.first.first).dofs_per_vertex -
+ 4*base_element(m_index.first.first).dofs_per_line;
+ Assert (tmp2 < base_element(m_index.first.first).dofs_per_quad,
+ ExcInternalError());
+ m_index.second = base_element(m_index.first.first).dofs_per_vertex +
+ 12*base_element(m_index.first.first).dofs_per_line +
+ base_element(m_index.first.first).dofs_per_quad*sub_quad +
+ tmp2;
};
break;
Assert (false, ExcNotImplemented());
};
- // now that we gathered all information:
- // use it to build the matrix. note
- // that if n and m belong to different
- // components, there definitely will be
- // no coupling
+ // now that we gathered all
+ // information: use it to
+ // build the matrix. note
+ // that if n and m belong to
+ // different base elements or
+ // instances, then there
+ // definitely will be no
+ // coupling
if (n_index.first == m_index.first)
interface_constraints(m,n)
- = (base_element(component_to_base_table[n_index.first])
- .constraints()(m_index.second,
- n_index.second));
+ = (base_element(n_index.first.first).constraints()(m_index.second,
+ n_index.second));
};
};
template <int dim>
void FESystem<dim>::initialize ()
{
- build_cell_table();
- build_face_table();
+ build_cell_tables();
+ build_face_tables();
// Check if some of the matrices of
// the base elements are void.
for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
prolongation[i].reinit(0,0);
- // distribute the matrices of the base
- // finite elements to the matrices of
- // this object
- for (unsigned int component=0; component<n_components(); ++component)
- // transform restriction and
- // prolongation matrices
- for (unsigned int i=0; i<base_element(component_to_base_table[component]).dofs_per_cell; ++i)
- for (unsigned int j=0; j<base_element(component_to_base_table[component]).dofs_per_cell; ++j)
- // only fill block diagonals, no
- // intermixing of subelements
+ // distribute the matrices of the
+ // base finite elements to the
+ // matrices of this object. for
+ // this, loop over all degrees of
+ // freedom and take the respective
+ // entry of the underlying base
+ // element.
+ //
+ // note that we by definition of a
+ // base element, they are
+ // independent, i.e. do not
+ // couple. only DoFs that belong to
+ // the same instance of a base
+ // element may couple
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ // first find out to which
+ // base element indices i and
+ // j belong, and which
+ // instance thereof in case
+ // the base element has a
+ // multiplicity greater than
+ // one. if they should not
+ // happen to belong to the
+ // same instance of a base
+ // element, then they cannot
+ // couple, so go on with the
+ // next index
+ if (system_to_base_table[i].first !=
+ system_to_base_table[j].first)
+ continue;
+
+ // so get the common base
+ // element and the indices
+ // therein:
+ const unsigned int
+ base = system_to_base_table[i].first.first;
+
+ const unsigned int
+ base_index_i = system_to_base_table[i].second,
+ base_index_j = system_to_base_table[j].second;
+
+ // if we are sure that DoFs i
+ // and j may couple, then
+ // copy entries of the
+ // matrices:
for (unsigned int child=0; child<GeometryInfo<dim>::children_per_cell; ++child)
{
if (do_restriction)
- restriction[child] (component_to_system_index (component,i),
- component_to_system_index (component, j))
- = base_element(component_to_base_table[component]).restrict(child)(i,j);
+ restriction[child] (i,j)
+ = base_element(base).restrict(child)(base_index_i,base_index_j);
+
if (do_prolongation)
- prolongation[child] (component_to_system_index (component,i),
- component_to_system_index (component, j))
- = base_element(component_to_base_table[component]).prolongate(child)(i,j);
+ prolongation[child] (i,j)
+ = base_element(base).prolongate(child)(base_index_i,base_index_j);
};
-
+ };
// now set up the interface constraints.
// this is kind'o hairy, so don't try
// an empty array to
// demonstrate that
// fact
- for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
+ for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
if (!base_element(base_el).has_support_points())
{
unit_support_points.resize(0);
// from unit support points of sub
// elements
unit_support_points.resize(dofs_per_cell);
-
- unsigned int comp = 0;
- for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
{
- // we know that there are
- // support points on the cell,
- // so collect them
const unsigned int
- base_element_dofs_per_cell = base_element(base_el).dofs_per_cell;
-
- const typename std::vector<Point<dim> >
- & base_unit_support_points = base_element(base_el).get_unit_support_points ();
-
- // otherwise distribute the
- // support points of this base
- // element to all degrees of
- // freedom contributed by this
- // base element
- Assert(base_unit_support_points.size()==base_element_dofs_per_cell,
- ExcInternalError());
- for (unsigned int n=0; n<element_multiplicity(base_el); ++n, ++comp)
- for (unsigned int i=0; i<base_element_dofs_per_cell; ++i)
- unit_support_points[component_to_system_index(comp,i)]
- = base_unit_support_points[i];
- }
+ base = system_to_base_table[i].first.first,
+ base_index = system_to_base_table[i].second;
+ Assert (base<n_base_elements(), ExcInternalError());
+ Assert (base_index<base_element(base).unit_support_points.size(),
+ ExcInternalError());
+ unit_support_points[i] = base_element(base).unit_support_points[base_index];
+ };
}
FESystem<dim>::
initialize_unit_face_support_points ()
{
- // if one of the base elements
- // has no support points, then
- // it makes no sense to define
- // support points for the
- // composed element, so return
- // an empty array to
- // demonstrate that fact (note
- // that we ask whether the base
- // element has no support
- // points at all, not only none
- // on the face!)
- for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
+ // if one of the base elements has
+ // no support points, then it makes
+ // no sense to define support
+ // points for the composed element,
+ // so return an empty array to
+ // demonstrate that fact (note that
+ // we ask whether the base element
+ // has no support points at all,
+ // not only none on the face!)
+ for (unsigned int base_el=0; base_el<n_base_elements(); ++base_el)
if (!base_element(base_el).has_support_points())
{
unit_face_support_points.resize(0);
// from unit support points of sub
// elements
unit_face_support_points.resize(dofs_per_face);
-
- unsigned int comp = 0;
- for (unsigned int base_el=0 ; base_el<n_base_elements(); ++base_el)
+
+ for (unsigned int i=0; i<dofs_per_face; ++i)
{
- // in some cases, finite
- // elements have support points
- // (we have made sure that they
- // have above) but don't have
- // any on the face (e.g. DG
- // elements). in that case,
- // don't even bother with this
- // base element and directly go
- // to the next one:
- if (!base_element(base_el).has_face_support_points())
- {
- comp += element_multiplicity(base_el);
- continue;
- };
+ const unsigned int base_i = face_system_to_base_table[i].first.first;
+ const unsigned int index_in_base = face_system_to_base_table[i].second;
- // otherwise, we know that
- // there are support points on
- // the face, so collect them
- const unsigned int
- base_element_dofs_per_face = base_element(base_el).dofs_per_face;
-
- const typename std::vector<Point<dim-1> > &
- base_unit_support_points = base_element(base_el).get_unit_face_support_points ();
+ Assert (index_in_base < base_element(base_i).unit_face_support_points.size(),
+ ExcInternalError());
- // distribute the support
- // points of this base element
- // to all degrees of freedom
- // contributed by this base
- // element
- Assert(base_unit_support_points.size()==base_element_dofs_per_face,
- ExcNotImplemented());
- for (unsigned int n=0; n<element_multiplicity(base_el); ++n, ++comp)
- for (unsigned int i=0; i<base_element_dofs_per_face; ++i)
- unit_face_support_points[face_component_to_system_index(comp,i)]
- = base_unit_support_points[i];
- }
+ unit_face_support_points[i]
+ = base_element(base_i).unit_face_support_points[index_in_base];
+ };
}
+template <int dim>
+std::vector<std::vector<bool> >
+FESystem<dim>::compute_nonzero_components (const FiniteElement<dim> &fe1,
+ const unsigned int N1)
+{
+ std::vector<const FiniteElement<dim>*> fe_list;
+ std::vector<unsigned int> multiplicities;
+
+ fe_list.push_back (&fe1);
+ multiplicities.push_back (N1);
+
+ return compute_nonzero_components (fe_list, multiplicities);
+};
+
+
+
+template <int dim>
+std::vector<std::vector<bool> >
+FESystem<dim>::compute_nonzero_components (const FiniteElement<dim> &fe1,
+ const unsigned int N1,
+ const FiniteElement<dim> &fe2,
+ const unsigned int N2)
+{
+ 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_nonzero_components (fe_list, multiplicities);
+};
+
+
+
+template <int dim>
+std::vector<std::vector<bool> >
+FESystem<dim>::compute_nonzero_components (const FiniteElement<dim> &fe1,
+ const unsigned int N1,
+ const FiniteElement<dim> &fe2,
+ const unsigned int N2,
+ const FiniteElement<dim> &fe3,
+ const unsigned int N3)
+{
+ 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_nonzero_components (fe_list, multiplicities);
+};
+
+
+
+template <int dim>
+std::vector<std::vector<bool> >
+FESystem<dim>::
+compute_nonzero_components (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];
+
+ unsigned int n_components = 0;
+ for (unsigned int i=0; i<fes.size(); ++i)
+ n_components += fes[i]->n_components() * multiplicities[i];
+
+ // generate the array that will
+ // hold the output
+ std::vector<std::vector<bool> >
+ retval (n_shape_functions, std::vector<bool> (n_components, 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 non-zero flags from the base
+ // element to this one, taking into
+ // account multiplicities, multiple
+ // components in base elements, and
+ // other complications
+ unsigned int total_index = 0;
+ for (unsigned int vertex_number=0;
+ vertex_number<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex_number)
+ {
+ unsigned comp_start = 0;
+ for(unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_components())
+ 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 (comp_start+fes[base]->n_components() <=
+ retval[total_index].size(),
+ ExcInternalError());
+ for (unsigned int c=0; c<fes[base]->n_components(); ++c)
+ {
+ Assert (index_in_base < fes[base]->nonzero_components.size(),
+ ExcInternalError());
+ Assert (c < fes[base]->nonzero_components[index_in_base].size(),
+ ExcInternalError());
+ retval[total_index][comp_start+c]
+ = fes[base]->nonzero_components[index_in_base][c];
+ };
+ }
+ }
+
+ // 2. Lines
+ if (GeometryInfo<dim>::lines_per_cell > 0)
+ for (unsigned int line_number= 0;
+ line_number != GeometryInfo<dim>::lines_per_cell;
+ ++line_number)
+ {
+ unsigned comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_components())
+ 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 (comp_start+fes[base]->n_components() <=
+ retval[total_index].size(),
+ ExcInternalError());
+ for (unsigned int c=0; c<fes[base]->n_components(); ++c)
+ {
+ Assert (index_in_base < fes[base]->nonzero_components.size(),
+ ExcInternalError());
+ Assert (c < fes[base]->nonzero_components[index_in_base].size(),
+ ExcInternalError());
+ retval[total_index][comp_start+c]
+ = fes[base]->nonzero_components[index_in_base][c];
+ };
+ }
+ }
+
+ // 3. Quads
+ if (GeometryInfo<dim>::quads_per_cell > 0)
+ for (unsigned int quad_number= 0;
+ quad_number != GeometryInfo<dim>::quads_per_cell;
+ ++quad_number)
+ {
+ unsigned int comp_start = 0;
+ for (unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_components())
+ 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 (comp_start+fes[base]->n_components() <=
+ retval[total_index].size(),
+ ExcInternalError());
+ for (unsigned int c=0; c<fes[base]->n_components(); ++c)
+ {
+ Assert (index_in_base < fes[base]->nonzero_components.size(),
+ ExcInternalError());
+ Assert (c < fes[base]->nonzero_components[index_in_base].size(),
+ ExcInternalError());
+ retval[total_index][comp_start+c]
+ = fes[base]->nonzero_components[index_in_base][c];
+ };
+ }
+ }
+
+ // 4. Hexes
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (unsigned int hex_number= 0;
+ hex_number != GeometryInfo<dim>::hexes_per_cell;
+ ++hex_number)
+ {
+ unsigned int comp_start = 0;
+ for(unsigned int base=0; base<fes.size(); ++base)
+ for (unsigned int m=0; m<multiplicities[base];
+ ++m, comp_start+=fes[base]->n_components())
+ 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 (comp_start+fes[base]->n_components() <=
+ retval[total_index].size(),
+ ExcInternalError());
+ for (unsigned int c=0; c<fes[base]->n_components(); ++c)
+ {
+ Assert (index_in_base < fes[base]->nonzero_components.size(),
+ ExcInternalError());
+ Assert (c < fes[base]->nonzero_components[index_in_base].size(),
+ ExcInternalError());
+ retval[total_index][comp_start+c]
+ = fes[base]->nonzero_components[index_in_base][c];
+ };
+ }
+ }
+
+ Assert (total_index == retval.size(), ExcInternalError());
+
+ return retval;
+};
+
+
+
template <int dim>
const FiniteElement<dim> &
FESystem<dim>::base_element (const unsigned int index) const
{
const std::pair<unsigned int, unsigned int> component
= system_to_component_index(shape_index);
- const unsigned int base = component_to_base(component.first);
+ const unsigned int base = component_to_base(component.first).first;
return base_element(base).has_support_on_face(component.second,
face_index);
}
template <int dim>
void
-FEValuesData<dim>::initialize (const unsigned int n_quadrature_points,
- const unsigned int n_shapes,
- const UpdateFlags flags)
+FEValuesData<dim>::initialize (const unsigned int n_quadrature_points,
+ const FiniteElement<dim> &fe,
+ const UpdateFlags flags)
{
update_flags = flags;
-
+
+ // initialize the table mapping
+ // from shape function number to
+ // the rows in the tables denoting
+ // its first non-zero
+ // component. with this also count
+ // the total number of non-zero
+ // components accumulated over all
+ // shape functions
+ shape_function_to_row_table.resize (fe.dofs_per_cell);
+ unsigned int row = 0;
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ shape_function_to_row_table[i] = row;
+ row += fe.n_nonzero_components (i);
+ };
+
+ const unsigned int n_nonzero_shape_components = row;
+ Assert (n_nonzero_shape_components >= fe.dofs_per_cell,
+ ExcInternalError());
+
+ // with the number of rows now
+ // known, initialize those fields
+ // that we will need to their
+ // correct size
if (flags & update_values)
- shape_values.reinit(n_shapes, n_quadrature_points);
+ shape_values.reinit(n_nonzero_shape_components, n_quadrature_points);
if (flags & update_gradients)
{
- shape_gradients.resize(n_shapes);
- for (unsigned int i=0;i<n_shapes;++i)
+ shape_gradients.resize(n_nonzero_shape_components);
+ for (unsigned int i=0; i<n_nonzero_shape_components; ++i)
shape_gradients[i].resize(n_quadrature_points);
}
if (flags & update_second_derivatives)
{
- shape_2nd_derivatives.resize(n_shapes);
- for (unsigned int i=0;i<n_shapes;++i)
+ shape_2nd_derivatives.resize(n_nonzero_shape_components);
+ for (unsigned int i=0; i<n_nonzero_shape_components; ++i)
shape_2nd_derivatives[i].resize(n_quadrature_points);
}
template <int dim>
FEValuesBase<dim>::~FEValuesBase ()
{
- if (fe_data)
+ // delete those fields that were
+ // created by the mapping and
+ // finite element objects,
+ // respectively, but of which we
+ // have assumed ownership
+ if (fe_data != 0)
{
typename Mapping<dim>::InternalDataBase *tmp1=fe_data;
fe_data=0;
delete tmp1;
}
- if (mapping_data)
+ if (mapping_data != 0)
{
typename Mapping<dim>::InternalDataBase *tmp1=mapping_data;
mapping_data=0;
std::fill_n (values.begin(), n_quadrature_points, 0);
// add up contributions of trial
+ // functions. note that here we
+ // deal with scalar finite
+ // elements, so no need to check
+ // for non-primitivity of shape
// functions
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
values[point] += (dof_values(shape_func) *
- shape_values(shape_func, point));
+ shape_value(shape_func, point));
};
{
Assert (n_quadrature_points == values.size(),
ExcWrongVectorSize(values.size(), n_quadrature_points));
+
+ const unsigned int n_components = fe->n_components();
for (unsigned i=0;i<values.size();++i)
- Assert (values[i].size() == fe->n_components(),
- ExcWrongNoOfComponents());
+ Assert (values[i].size() == n_components, ExcWrongNoOfComponents());
+
Assert (update_flags & update_values, ExcAccessToUninitializedField());
Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs()));
std::fill_n (values[i].begin(), values[i].size(), 0);
// add up contributions of trial
- // functions
+ // functions. now check whether the
+ // shape function is primitive or
+ // not. if it is, then set its only
+ // non-zero component, otherwise
+ // loop over components
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
- values[point](fe->system_to_component_index(shape_func).first)
- += (dof_values(shape_func) * shape_values(shape_func, point));
-};
-
-
-
-template <int dim>
-const typename FEValuesData<dim>::ShapeVector &
-FEValuesBase<dim>::get_shape_values () const
-{
- Assert (update_flags & update_values, ExcAccessToUninitializedField());
- return shape_values;
-};
-
-
-
-template <int dim>
-const typename FEValuesData<dim>::GradientVector &
-FEValuesBase<dim>::get_shape_grads () const
-{
- Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
- return shape_gradients;
-};
-
-
-
-
-template <int dim>
-const typename std::vector<typename std::vector<Tensor<2,dim> > > &
-FEValuesBase<dim>::get_shape_2nd_derivatives () const
-{
- Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
- return shape_2nd_derivatives;
+ if (fe->is_primitive(shape_func))
+ values[point](fe->system_to_component_index(shape_func).first)
+ += (dof_values(shape_func) * shape_value(shape_func, point));
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ values[point](c) += (dof_values(shape_func) *
+ shape_value_component(shape_func, point, c));
};
std::fill_n (gradients.begin(), n_quadrature_points, Tensor<1,dim>());
// add up contributions of trial
+ // functions. note that here we
+ // deal with scalar finite
+ // elements, so no need to check
+ // for non-primitivity of shape
// functions
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
- Tensor<1,dim> tmp(shape_gradients[shape_func][point]);
+ Tensor<1,dim> tmp = shape_grad(shape_func,point);
tmp *= dof_values(shape_func);
gradients[point] += tmp;
};
{
Assert (n_quadrature_points == gradients.size(),
ExcWrongNoOfComponents());
- for (unsigned i=0;i<gradients.size();++i)
- Assert (gradients[i].size() == fe->n_components(),
- ExcWrongVectorSize(gradients[i].size(), fe->n_components()));
+
+ const unsigned int n_components = fe->n_components();
+ for (unsigned i=0; i<gradients.size(); ++i)
+ Assert (gradients[i].size() == n_components, ExcWrongNoOfComponents());
+
Assert (update_flags & update_gradients, ExcAccessToUninitializedField());
Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs()));
std::fill_n (gradients[i].begin(), gradients[i].size(), Tensor<1,dim>());
// add up contributions of trial
- // functions
+ // functions. now check whether the
+ // shape function is primitive or
+ // not. if it is, then set its only
+ // non-zero component, otherwise
+ // loop over components
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
- {
- Tensor<1,dim> tmp(shape_gradients[shape_func][point]);
- tmp *= dof_values(shape_func);
- gradients[point][fe->system_to_component_index(shape_func).first]
- += tmp;
- };
+ if (fe->is_primitive (shape_func))
+ {
+ Tensor<1,dim> tmp = shape_grad(shape_func,point);
+ tmp *= dof_values(shape_func);
+ gradients[point][fe->system_to_component_index(shape_func).first]
+ += tmp;
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ Tensor<1,dim> tmp = shape_grad_component(shape_func,point,c);
+ tmp *= dof_values(shape_func);
+ gradients[point][c] += tmp;
+ };
};
std::fill_n (second_derivatives.begin(), n_quadrature_points, Tensor<2,dim>());
// add up contributions of trial
+ // functions. note that here we
+ // deal with scalar finite
+ // elements, so no need to check
+ // for non-primitivity of shape
// functions
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
{
- Tensor<2,dim> tmp(shape_2nd_derivatives[shape_func][point]);
+ Tensor<2,dim> tmp = shape_2nd_derivative(shape_func,point);
tmp *= dof_values(shape_func);
second_derivatives[point] += tmp;
};
{
Assert (n_quadrature_points == second_derivs.size(),
ExcWrongNoOfComponents());
+
+ const unsigned int n_components = fe->n_components();
for (unsigned i=0;i<second_derivs.size();++i)
- Assert (second_derivs[i].size() == fe->n_components(),
- ExcWrongVectorSize(second_derivs[i].size(), fe->n_components()));
+ Assert (second_derivs[i].size() == n_components, ExcWrongNoOfComponents());
+
Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField());
Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(),
ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs()));
// functions
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<dofs_per_cell; ++shape_func)
- {
- Tensor<2,dim> tmp(shape_2nd_derivatives[shape_func][point]);
- tmp *= dof_values(shape_func);
- second_derivs[point][fe->system_to_component_index(shape_func).first]
- += tmp;
- };
+ if (fe->is_primitive(shape_func))
+ {
+ Tensor<2,dim> tmp(shape_2nd_derivative(shape_func,point));
+ tmp *= dof_values(shape_func);
+ second_derivs[point][fe->system_to_component_index(shape_func).first]
+ += tmp;
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ Tensor<2,dim> tmp = shape_2nd_derivative_component(shape_func,point,c);
+ tmp *= dof_values(shape_func);
+ second_derivs[point][c] += tmp;
+ };
};
MemoryConsumption::memory_consumption (mapping_data) +
MemoryConsumption::memory_consumption (*mapping_data) +
MemoryConsumption::memory_consumption (fe_data) +
- MemoryConsumption::memory_consumption (*fe_data));
+ MemoryConsumption::memory_consumption (*fe_data) +
+ MemoryConsumption::memory_consumption (shape_function_to_row_table));
};
fe_data = fe->get_data(flags, *mapping, quadrature);
// set up objects within this class
- FEValuesData<dim>::initialize(n_quadrature_points,
- dofs_per_cell,
- flags);
+ FEValuesData<dim>::initialize (n_quadrature_points, *fe, flags);
};
fe_data = fe->get_face_data(flags, *mapping, quadrature);
// set up objects within this class
- FEValuesData<dim>::initialize(n_quadrature_points,
- dofs_per_cell,
- flags);
+ FEValuesData<dim>::initialize(n_quadrature_points, *fe, flags);
};
fe_data = fe->get_subface_data(flags, *mapping, quadrature);
// set up objects within this class
- FEValuesData<dim>::initialize(n_quadrature_points,
- dofs_per_cell,
- flags);
+ FEValuesData<dim>::initialize(n_quadrature_points, *fe, flags);
};
<h3>General</h3>
<ol>
+ <li> <p>
+ Removed: The functions <code
+ class="member">FEValuesBase::get_shape_values</code>, <code
+ class="member">FEValuesBase::get_shape_grads</code>, and <code
+ class="member">FEValuesBase::get_shape_2nd_derivatives</code> are
+ now removed as they expose too much of the internal data
+ structure of their class, and interfere with plans for the
+ extension of this and related classes. These functions, which
+ had been deprecated in the previous version, are thus not
+ part of the release after version 3.3 any more. Use the <code
+ class="member">FEValuesBase::shape_value</code> and alike
+ functions as a replacement.
+ <br>
+ For more information, read
+ <a href="http://www.dealii.org/mail/msg00638.html" target="body">this mail</a>.
+ <br>
+ (WB 2002/06/10)
+ </p>
- <li> Improved: make knows a target <tt>veryclean</tt>
- now. <tt>clean</tt> leaves the libraries now, removing everything
- that is not needed to use
- <acronym>deal.II</acronym>. <tt>veryclean</tt> removes even the
- libraries, leaving the directory more or less in the state like
- after <tt>configure</tt>.
- <br>
- (GK 2002/06/07)
+ <li> <p>
+ New: <acronym>deal.II</acronym> now also supports vector-valued
+ finite elements with shape functions for which more than just
+ one vector component is non-zero. Such elements are, for
+ example, the Nedelec and Raviart-Thomas families. Previously,
+ vector-valued elements were only supported insofar as they
+ could be composed of scalar ones; in that case, each
+ (vector-valued) shape function had only one non-zero vector
+ component.
+ <br>
+ (WB 2002/06/10)
+ </p>
+
+ <li> <p>
+ Improved: make knows a target <tt>veryclean</tt>
+ now. <tt>clean</tt> leaves the libraries now, removing
+ everything that is not needed to use
+ <acronym>deal.II</acronym>. <tt>veryclean</tt> removes even the
+ libraries, leaving the directory more or less in the state like
+ after <tt>configure</tt>.
+ <br>
+ (GK 2002/06/07)
</ol>
<dir>
<dl>
+ <dt>
+ <strong>2002/06/10: Supporting framework for vector-valued
+ finite elements committed</strong>
+ </dt>
+ <dd>
+ <acronym>deal.II</acronym> supports vector-valued finite
+ elements for a long time, but only by composing them of
+ scalar elements. This excludes Nedelec or Raviart-Thomas
+ type elements, for which every shape function is
+ vector-valued and non-zero in all components. Support for
+ such elements is now available, and will be part of the next
+ version.
+ </p>
+ </dd>
+
<dt>
<strong>Changes between version 3.4 and the
main branch</strong>
derivatives.exe : derivatives.go $(libraries)
fe_data_test.exe : fe_data_test.go $(libraries)
+fe_traits.exe : fe_traits.go $(libraries)
+internals.exe : internals.go $(libraries)
mapping.exe : mapping.go $(libraries)
mapping_c1.exe : mapping_c1.go $(libraries)
-shapes.exe : shapes.go $(libraries)
-numbering.exe : numbering.go $(libraries)
mapping_q1_eulerian.exe : mapping_q1_eulerian.go $(libraries)
-internals.exe : internals.go $(libraries)
+non_primitive_1.exe : non_primitive_1.go $(libraries)
+numbering.exe : numbering.go $(libraries)
+shapes.exe : shapes.go $(libraries)
-tests = derivatives fe_data_test mapping mapping_c1 shapes \
- numbering mapping_q1_eulerian internals
+tests = derivatives fe_data_test fe_traits internals \
+ mapping mapping_c1 mapping_q1_eulerian \
+ non_primitive_1 numbering shapes
############################################################
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
#include <fe/fe_system.h>
FE_Q<dim> (3), 1));
fe_datas.push_back(new FESystem<dim>(FE_Q<dim> (4), 2));
}
+
+ // have systems of systems, and
+ // construct hierarchies of
+ // subsequently weirder elements by
+ // taking each of them in turn as
+ // basis of others
+ fe_datas.push_back (new FESystem<dim> (FESystem<dim> (FE_Q<dim>(1),2),2));
+ fe_datas.push_back (new FESystem<dim> (FESystem<dim> (FE_Q<dim>(1),2),1,
+ FESystem<dim> (FE_DGQ<dim>(1),2),1));
+ fe_datas.push_back (new FESystem<dim> (FESystem<dim> (FE_Q<dim>(1),1,
+ FE_Q<dim>(2),1),1,
+ FESystem<dim> (FE_Q<dim>(2),2),1,
+ FESystem<dim> (FE_DGQ<dim>(2),2),1));
+ fe_datas.push_back (new FESystem<dim> (*fe_datas[fe_datas.size()-3], 2,
+ *fe_datas[fe_datas.size()-2], 1,
+ *fe_datas[fe_datas.size()-1], 2));
deallog << std::endl << "dim=" << dim << std::endl;
for (unsigned int n=0; n<fe_datas.size(); ++n)
DEAL::dofs_per_face=2
DEAL::dofs_per_cell=10
DEAL::components=2
+DEAL::fe_data[7]:
+DEAL::dofs_per_vertex=4
+DEAL::dofs_per_line=0
+DEAL::dofs_per_quad=0
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=8
+DEAL::first_quad_index=8
+DEAL::first_hex_index=8
+DEAL::first_face_line_index=4
+DEAL::first_face_quad_index=8
+DEAL::dofs_per_face=4
+DEAL::dofs_per_cell=8
+DEAL::components=4
+DEAL::fe_data[8]:
+DEAL::dofs_per_vertex=2
+DEAL::dofs_per_line=4
+DEAL::dofs_per_quad=0
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=4
+DEAL::first_quad_index=8
+DEAL::first_hex_index=8
+DEAL::first_face_line_index=2
+DEAL::first_face_quad_index=4
+DEAL::dofs_per_face=2
+DEAL::dofs_per_cell=8
+DEAL::components=4
+DEAL::fe_data[9]:
+DEAL::dofs_per_vertex=4
+DEAL::dofs_per_line=9
+DEAL::dofs_per_quad=0
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=8
+DEAL::first_quad_index=17
+DEAL::first_hex_index=17
+DEAL::first_face_line_index=4
+DEAL::first_face_quad_index=8
+DEAL::dofs_per_face=4
+DEAL::dofs_per_cell=17
+DEAL::components=6
+DEAL::fe_data[10]:
+DEAL::dofs_per_vertex=18
+DEAL::dofs_per_line=22
+DEAL::dofs_per_quad=0
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=36
+DEAL::first_quad_index=58
+DEAL::first_hex_index=58
+DEAL::first_face_line_index=18
+DEAL::first_face_quad_index=36
+DEAL::dofs_per_face=18
+DEAL::dofs_per_cell=58
+DEAL::components=24
DEAL::
DEAL::dim=2
DEAL::fe_data[0]:
DEAL::dofs_per_face=10
DEAL::dofs_per_cell=50
DEAL::components=2
+DEAL::fe_data[7]:
+DEAL::dofs_per_vertex=4
+DEAL::dofs_per_line=0
+DEAL::dofs_per_quad=0
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=16
+DEAL::first_quad_index=16
+DEAL::first_hex_index=16
+DEAL::first_face_line_index=8
+DEAL::first_face_quad_index=16
+DEAL::dofs_per_face=8
+DEAL::dofs_per_cell=16
+DEAL::components=4
+DEAL::fe_data[8]:
+DEAL::dofs_per_vertex=2
+DEAL::dofs_per_line=0
+DEAL::dofs_per_quad=8
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=8
+DEAL::first_quad_index=8
+DEAL::first_hex_index=16
+DEAL::first_face_line_index=4
+DEAL::first_face_quad_index=8
+DEAL::dofs_per_face=4
+DEAL::dofs_per_cell=16
+DEAL::components=4
+DEAL::fe_data[9]:
+DEAL::dofs_per_vertex=4
+DEAL::dofs_per_line=3
+DEAL::dofs_per_quad=21
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=16
+DEAL::first_quad_index=28
+DEAL::first_hex_index=49
+DEAL::first_face_line_index=8
+DEAL::first_face_quad_index=19
+DEAL::dofs_per_face=11
+DEAL::dofs_per_cell=49
+DEAL::components=6
+DEAL::fe_data[10]:
+DEAL::dofs_per_vertex=18
+DEAL::dofs_per_line=6
+DEAL::dofs_per_quad=50
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=72
+DEAL::first_quad_index=96
+DEAL::first_hex_index=146
+DEAL::first_face_line_index=36
+DEAL::first_face_quad_index=78
+DEAL::dofs_per_face=42
+DEAL::dofs_per_cell=146
+DEAL::components=24
DEAL::
DEAL::dim=3
DEAL::fe_data[0]:
DEAL::dofs_per_face=17
DEAL::dofs_per_cell=43
DEAL::components=3
+DEAL::fe_data[4]:
+DEAL::dofs_per_vertex=4
+DEAL::dofs_per_line=0
+DEAL::dofs_per_quad=0
+DEAL::dofs_per_hex=0
+DEAL::first_line_index=32
+DEAL::first_quad_index=32
+DEAL::first_hex_index=32
+DEAL::first_face_line_index=16
+DEAL::first_face_quad_index=16
+DEAL::dofs_per_face=16
+DEAL::dofs_per_cell=32
+DEAL::components=4
+DEAL::fe_data[5]:
+DEAL::dofs_per_vertex=2
+DEAL::dofs_per_line=0
+DEAL::dofs_per_quad=0
+DEAL::dofs_per_hex=16
+DEAL::first_line_index=16
+DEAL::first_quad_index=16
+DEAL::first_hex_index=16
+DEAL::first_face_line_index=8
+DEAL::first_face_quad_index=8
+DEAL::dofs_per_face=8
+DEAL::dofs_per_cell=32
+DEAL::components=4
+DEAL::fe_data[6]:
+DEAL::dofs_per_vertex=4
+DEAL::dofs_per_line=3
+DEAL::dofs_per_quad=3
+DEAL::dofs_per_hex=57
+DEAL::first_line_index=32
+DEAL::first_quad_index=68
+DEAL::first_hex_index=86
+DEAL::first_face_line_index=16
+DEAL::first_face_quad_index=28
+DEAL::dofs_per_face=31
+DEAL::dofs_per_cell=143
+DEAL::components=6
+DEAL::fe_data[7]:
+DEAL::dofs_per_vertex=18
+DEAL::dofs_per_line=6
+DEAL::dofs_per_quad=6
+DEAL::dofs_per_hex=130
+DEAL::first_line_index=144
+DEAL::first_quad_index=216
+DEAL::first_hex_index=252
+DEAL::first_face_line_index=72
+DEAL::first_face_quad_index=96
+DEAL::dofs_per_face=102
+DEAL::dofs_per_cell=382
+DEAL::components=24
--- /dev/null
+//---------------------------- fe_traits_test.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002 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.
+//
+//---------------------------- fe_traits_test.cc ---------------------------
+
+
+#include <iostream>
+#include <fstream>
+
+#include <base/logstream.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_dgp.h>
+#include <fe/fe_system.h>
+
+
+template <int dim>
+void check (const FiniteElement<dim> &fe)
+{
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+ // first check whether shape
+ // functions are primitive:
+ deallog << " Primitivity: ";
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ deallog << (fe.is_primitive(i) ? 1 : 0);
+ deallog << std::endl;
+
+ deallog << " Overall primitivity: " << fe.is_primitive() << std::endl;
+
+ // then check n_nonzero_components
+ deallog << " n_nonzero_components: ";
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ deallog << fe.n_nonzero_components(i);
+ deallog << std::endl;
+
+ // finally check component pattern
+ // for each shape function
+ deallog << " component pattern for each shape function:" << std::endl;
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ deallog << " shape function i=" << i
+ << " [" << fe.system_to_component_index(i).first
+ << ','
+ << fe.system_to_component_index(i).second
+ << "]: ";
+
+ for (unsigned int j=0; j<fe.n_components(); ++j)
+ {
+ Assert (fe.get_nonzero_components(i).size() == fe.n_components(),
+ ExcInternalError());
+ deallog << (fe.get_nonzero_components(i)[j] ? 1 : 0);
+ };
+ deallog << std::endl;
+ };
+};
+
+
+
+template <int dim>
+void check ()
+{
+ // check usual Lagrange elements
+ for (unsigned int p=1; p<3; ++p)
+ {
+ deallog << "Checking FE_Q<" << dim << ">(" << p << "): "
+ << std::endl;
+ check (FE_Q<dim>(p));
+ };
+
+ // check DG Lagrange elements
+ for (unsigned int p=0; p<3; ++p)
+ {
+ deallog << "Checking FE_DGQ<" << dim << ">(" << p << "): "
+ << std::endl;
+ check (FE_DGQ<dim>(p));
+ };
+
+ // check DG-P elements
+ for (unsigned int p=0; p<3; ++p)
+ {
+ deallog << "Checking FE_DGP<" << dim << ">(" << p << "): "
+ << std::endl;
+ check (FE_DGP<dim>(p));
+ };
+
+ // check systems of Q-elements
+ for (unsigned int p=1; p<3; ++p)
+ {
+ deallog << "Checking FE_Q<" << dim << ">(" << p << ")^2: "
+ << std::endl;
+ check (FESystem<dim> (FE_Q<dim>(p),2));
+ };
+
+ // check systems of systems of
+ // Q-elements
+ for (unsigned int p=1; p<3; ++p)
+ {
+ deallog << "Checking FE_Q<" << dim << ">(" << p << ")^2^2: "
+ << std::endl;
+ check (FESystem<dim> (FESystem<dim> (FE_Q<dim>(p),2), 2));
+ };
+};
+
+
+
+
+int main ()
+{
+ std::ofstream logfile("fe_traits.output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ check<1> ();
+ check<2> ();
+ check<3> ();
+};
+
QTrapez<1> q_trapez;
QIterated<dim> q(q_trapez, div);
- FEValues<dim> fe(mapping, finel, q, UpdateFlags(update_values));
+ FEValues<dim> fe(mapping, finel, q, UpdateFlags(update_values |
+ update_gradients |
+ update_second_derivatives));
fe.reinit(c);
sprintf(fname, "Shapes%dd-%s", dim, name);
for (unsigned int i=0;i<finel.dofs_per_cell;++i)
{
deallog << " " << fe.shape_value(i,k) + 1.;
+
+ // some additional
+ // checks
+ for (unsigned int c=0; c<fe.get_fe().n_components(); ++c)
+ {
+ if (fe.get_fe().system_to_component_index(i).first == c)
+ Assert ((fe.shape_value(i,k) ==
+ fe.shape_value_component(i,k,c))
+ &&
+ (fe.shape_grad(i,k) ==
+ fe.shape_grad_component(i,k,c))
+ &&
+ (fe.shape_2nd_derivative(i,k) ==
+ fe.shape_2nd_derivative_component(i,k,c)),
+ ExcInternalError())
+ else
+ Assert ((fe.shape_value_component(i,k,c) == 0) &&
+ (fe.shape_grad_component(i,k,c) == Tensor<1,dim>()) &&
+ (fe.shape_2nd_derivative_component(i,k,c) == Tensor<2,dim>()),
+ ExcInternalError());
+ };
}
deallog << std::endl;
k++;
QTrapez<1> q_trapez;
QIterated<dim-1> q(q_trapez, div);
FEFaceValues<dim> fe(mapping, finel, q, UpdateFlags(update_values
- | update_q_points));
+ | update_gradients
+ | update_second_derivatives
+ | update_q_points));
FESubfaceValues<dim> sub(mapping, finel, q, UpdateFlags(update_values
- | update_q_points));
+ | update_gradients
+ | update_second_derivatives
+ | update_q_points));
sprintf(fname, "ShapesFace%dd-%s", dim, name);
deallog.push(fname);
for (unsigned int i=0;i<finel.dofs_per_cell;++i)
{
deallog << " " << fe.shape_value(i,k) + 1.;
+
+ // some additional
+ // checks
+ for (unsigned int c=0; c<fe.get_fe().n_components(); ++c)
+ {
+ if (fe.get_fe().system_to_component_index(i).first == c)
+ Assert ((fe.shape_value(i,k) ==
+ fe.shape_value_component(i,k,c))
+ &&
+ (fe.shape_grad(i,k) ==
+ fe.shape_grad_component(i,k,c))
+ &&
+ (fe.shape_2nd_derivative(i,k) ==
+ fe.shape_2nd_derivative_component(i,k,c)),
+ ExcInternalError())
+ else
+ Assert ((fe.shape_value_component(i,k,c) == 0) &&
+ (fe.shape_grad_component(i,k,c) == Tensor<1,dim>()) &&
+ (fe.shape_2nd_derivative_component(i,k,c) == Tensor<2,dim>()),
+ ExcInternalError());
+ };
}
deallog << std::endl;
k++;
for (unsigned int i=0;i<finel.dofs_per_cell;++i)
{
deallog << " " << sub.shape_value(i,k) + 1.;
+
+ // some additional
+ // checks
+ for (unsigned int c=0; c<fe.get_fe().n_components(); ++c)
+ {
+ if (fe.get_fe().system_to_component_index(i).first == c)
+ Assert ((sub.shape_value(i,k) ==
+ sub.shape_value_component(i,k,c))
+ &&
+ (sub.shape_grad(i,k) ==
+ sub.shape_grad_component(i,k,c))
+ &&
+ (sub.shape_2nd_derivative(i,k) ==
+ sub.shape_2nd_derivative_component(i,k,c)),
+ ExcInternalError())
+ else
+ Assert ((sub.shape_value_component(i,k,c) == 0) &&
+ (sub.shape_grad_component(i,k,c) == Tensor<1,dim>()) &&
+ (sub.shape_2nd_derivative_component(i,k,c) == Tensor<2,dim>()),
+ ExcInternalError());
+ };
}
deallog << std::endl;
k++;
bool coincide=true;
for (unsigned int x=0; x<q.n_quadrature_points; ++x)
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- if (fabs(fe_values.shape_value(i,x)-fe.shape_value(i,q.point(x)))>1e-14)
- coincide=false;
+ {
+ if (fabs(fe_values.shape_value(i,x)-fe.shape_value(i,q.point(x)))>1e-14)
+ coincide=false;
+ for (unsigned int c=0; c<fe.n_components(); ++c)
+ Assert (((c == fe.system_to_component_index(i).first) &&
+ (fe_values.shape_value(i,x) == fe_values.shape_value_component(i,x,c)))
+ ||
+ ((c != fe.system_to_component_index(i).first) &&
+ (fe_values.shape_value_component(i,x,c) == 0)),
+ ExcInternalError());
+ };
+
if (!coincide)
deallog << "Error in fe.shape_value for " << name << std::endl;
tmp-=fe.shape_grad(i,q.point(x));
if (sqrt(tmp*tmp)>1e-14)
coincide=false;
+
+ for (unsigned int c=0; c<fe.n_components(); ++c)
+ Assert (((c == fe.system_to_component_index(i).first) &&
+ (fe_values.shape_grad(i,x) == fe_values.shape_grad_component(i,x,c)))
+ ||
+ ((c != fe.system_to_component_index(i).first) &&
+ (fe_values.shape_grad_component(i,x,c) == Tensor<1,dim>())),
+ ExcInternalError());
}
if (!coincide)
if (diff>max_diff) max_diff=diff;
if (fabs(tmp[j][k])>1e-6)
coincide=false;
- }
-
+ }
+
+ for (unsigned int c=0; c<fe.n_components(); ++c)
+ Assert (((c == fe.system_to_component_index(i).first) &&
+ (fe_values.shape_2nd_derivative(i,x) ==
+ fe_values.shape_2nd_derivative_component(i,x,c)))
+ ||
+ ((c != fe.system_to_component_index(i).first) &&
+ (fe_values.shape_2nd_derivative_component(i,x,c) == Tensor<2,dim>())),
+ ExcInternalError());
}
if (!coincide)