*/
const unsigned int n_transform_functions;
+ /**
+ * Number of components.
+ */
+ const unsigned int n_components;
+
/**
* Default constructor. Constructs an element
* which is not so useful. Checking
*/
FiniteElementData (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
- const unsigned int n_transform_functions);
+ const unsigned int n_transform_functions,
+ const unsigned int n_components = 1);
/**
* Another frequently used useful
*/
const unsigned int n_transform_functions;
+ /**
+ * Number of components.
+ */
+ const unsigned int n_components;
+
/**
* Default constructor. Constructs an element
* which is not so useful. Checking
FiniteElementData (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
const unsigned int dofs_per_quad,
- const unsigned int n_transform_functions);
+ const unsigned int n_transform_functions,
+ const unsigned int n_components = 1);
/**
* Another frequently used useful
virtual void get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &cell,
const Boundary<dim> &boundary,
dFMatrix &local_mass_matrix) const =0;
+
+ /**
+ * Compute system index from components.
+ */
+ virtual unsigned component_to_system_index (unsigned component,
+ unsigned component_index) const;
+
+ /**
+ * Compute component and index from system index.
+ */
+ virtual pair<unsigned,unsigned> system_to_component_index (unsigned index) const;
/**
* Exception
#include <fe/fe.h>
#include <base/quadrature.h>
-
+#include <math.h>
/**
* This class implements a rather unusual macro element, the so-called
#include <fe/fe_linear_mapping.h>
-
/**
* Define a (bi-, tri-, etc)linear finite element in #dim# space dimensions,
* along with (bi-, tri-)linear (therefore isoparametric) transforms from the
/*---------------------------- fe_linear_mapping.h ---------------------------*/
+#include <math.h>
#include <fe/fe.h>
* shape functions. These #N# shape functions for each degree of freedom
* of the basic finite element are numbered consecutively, i.e. for
* the common case of a velocity #(u,v,w)#, the sequence of basis functions
- * would be #u1, v1, w1, u2, v2, w2, ..., uN, vN, wN#. The other possibility
- * would have been #u1, ..., uN, v1, ..., vN, w1, ...wN#, but we chose the
- * first way.
+ * will be #u1, v1, w1, u2, v2, w2, ..., uN, vN, wN# compared to
+ * #u1, ..., uN, v1, ..., vN, w1, ...wN#.
*
* Using this scheme, the overall numbering of degrees of freedom is as
* follows: for each subobject (vertex, line, quad, or hex), the degrees
* \end{itemize}
*
* In the most cases, the composed element behaves as if it were a usual element
- * with more degrees of freedom. Howeverm the underlying structure is visible in
+ * with more degrees of freedom. However the underlying structure is visible in
* the restriction, prolongation and interface constraint matrices, which do not
* couple the degrees of freedom of the subobject. E.g. the continuity requirement
* is imposed for the shape functions of the subobjects separately; no requirement
* @author Wolfgang Bangerth, 1999
*/
template <int dim>
-class FESystem : public FiniteElement<dim> {
+class FESystem : public FiniteElement<dim>
+{
+ /**
+ * Copy constructor prohibited.
+ */
+ FESystem(const FESystem<dim>&);
+
+
public:
/**
const unsigned int subface_no,
const vector<Point<dim-1> > &unit_points,
vector<Point<dim> > &normal_vectors) const;
-
- private:
- /**
- * Pointer to an object of the underlying
- * finite element class. This object is
- * created by the constructor.
- */
- const FiniteElement<dim> *const base_element;
/**
* Number of subelements of this object.
- * Since these objects may have
+ * Since these objects can have
* subobjects themselves, this may be
* smaller than the total number of finite
* elements composed into this structure.
+ * This is definitely not what you'd
+ * usally intend, so don't do it!
*/
const unsigned int n_sub_elements;
+
+ /**
+ * Access to the single valued element.
+ *
+ * If you assemble your system
+ * matrix, you usually will not
+ * want to have an FEValues object
+ * with a lot of equal entries. Ok,
+ * so initialize your FEValues with
+ * the #base_element# yuo get by
+ * this function.
+ *
+ */
+ const FiniteElement<dim>& get_base_element() const;
+
+ /**
+ * Calculate the actual position.
+ *
+ * For a given #component#
+ * (e.g. u,v,w in the example
+ * above) of the #base# function of
+ * the #base_element#, return the
+ * actual index in the local
+ * degrees of freedom vector of the
+ * system.
+ *
+ */
+ unsigned index(unsigned component, unsigned base) const;
+
+ private:
+
+ /**
+ * Pointer to an object of the underlying
+ * finite element class. This object is
+ * created by the constructor.
+ */
+ const FiniteElement<dim> *const base_element;
+
/**
* Helper function used in the constructor:
/* ------------------------- template functions ------------------------- */
+template <int dim>
+inline const FiniteElement<dim>&
+FESystem<dim>::get_base_element() const
+{
+ return *base_element;
+}
+
+template <int dim>
+inline unsigned
+FESystem<dim>::index(unsigned component, unsigned base) const
+{
+ return n_sub_elements * base + component;
+}
+
template <int dim>
template <typename FE>
FESystem<dim>::FESystem (const FE &fe, const unsigned int n_elements) :
FiniteElement (multiply_dof_numbers(fe, n_elements)),
- base_element (new FE()),
- n_sub_elements (n_elements)
+ n_sub_elements (n_elements),
+ base_element (new FE())
{
base_element->subscribe ();
initialize_matrices ();
* of the derived class was called,
* at the quadrature points.
*
+ * To get values of
+ * multi-component elements,
+ * there is another
+ * #get_function_values#
+ * returning a vector of vectors of
+ * results.
+ *
* The function assumes that the
* #values# object already has the
- * right size.
+ * right size.
*/
void get_function_values (const dVector &fe_function,
vector<double> &values) const;
+ /**
+ * Access to vector valued finite element functions.
+ *
+ * This function does the same as
+ * the other #get_function_values#,
+ * but applied to multi-component
+ * elements.
+ */
+ void get_function_values (const dVector &fe_function,
+ vector<vector<double> > &values) const;
+
/**
* Return the gradient of the #i#th shape
* function at the #j# quadrature point.
/**
* Exception
*/
+ DeclException0 (ExcWrongNoOfComponents);
+ /**
+ * Exception.
+ */
DeclException2 (ExcWrongVectorSize,
int, int,
<< "Vector has wrong size " << arg1
#include <base/point.h>
#include <grid/geometry_info.h>
+#include <math.h>
dofs_per_line(0),
dofs_per_face(0),
total_dofs(0),
- n_transform_functions(0) {
+ n_transform_functions(0),
+ n_components(0)
+{
Assert (false, ExcInternalError());
};
FiniteElementData<1>::FiniteElementData (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
- const unsigned int n_transform_functions) :
+ const unsigned int n_transform_functions,
+ const unsigned int n_components) :
dofs_per_vertex(dofs_per_vertex),
dofs_per_line(dofs_per_line),
dofs_per_face(dofs_per_vertex),
total_dofs (2*dofs_per_vertex+dofs_per_line),
- n_transform_functions(n_transform_functions)
+ n_transform_functions(n_transform_functions),
+ n_components(n_components)
{};
dofs_per_line(fe_data.dofs_per_line),
dofs_per_face(fe_data.dofs_per_face),
total_dofs (fe_data.total_dofs),
- n_transform_functions(fe_data.n_transform_functions)
+ n_transform_functions(fe_data.n_transform_functions),
+ n_components(fe_data.n_components)
{};
bool FiniteElementData<1>::operator== (const FiniteElementData<1> &f) const {
return ((dofs_per_vertex == f.dofs_per_vertex) &&
(dofs_per_line == f.dofs_per_line) &&
- (total_dofs == f.total_dofs));
+ (total_dofs == f.total_dofs) &&
+ (n_components == f.n_components));
};
#endif
dofs_per_quad(0),
dofs_per_face(0),
total_dofs(0),
- n_transform_functions(0) {
+ n_transform_functions(0),
+ n_components(0)
+{
Assert (false, ExcInternalError());
};
FiniteElementData<2>::FiniteElementData (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
const unsigned int dofs_per_quad,
- const unsigned int n_transform_functions) :
+ const unsigned int n_transform_functions,
+ const unsigned int n_components) :
dofs_per_vertex(dofs_per_vertex),
dofs_per_line(dofs_per_line),
dofs_per_quad(dofs_per_quad),
total_dofs (4*dofs_per_vertex+
4*dofs_per_line+
dofs_per_quad),
- n_transform_functions (n_transform_functions)
+ n_transform_functions (n_transform_functions),
+ n_components(n_components)
{};
dofs_per_quad(fe_data.dofs_per_quad),
dofs_per_face(fe_data.dofs_per_face),
total_dofs (fe_data.total_dofs),
- n_transform_functions (fe_data.n_transform_functions)
+ n_transform_functions (fe_data.n_transform_functions),
+ n_components(fe_data.n_components)
{};
return ((dofs_per_vertex == f.dofs_per_vertex) &&
(dofs_per_line == f.dofs_per_line) &&
(dofs_per_quad == f.dofs_per_quad) &&
- (total_dofs == f.total_dofs));
+ (total_dofs == f.total_dofs) &&
+ (n_components == f.n_components));
};
#endif
template <int dim>
-FESystem<dim>::~FESystem () {
+FESystem<dim>::~FESystem ()
+{
base_element->unsubscribe ();
delete base_element;
};
template <int dim>
-void FESystem<dim>::initialize_matrices () {
+void FESystem<dim>::initialize_matrices ()
+{
// distribute the matrices of the base
// finite element to the matrices of
// this object
- for (unsigned int i=0; i<base_element->total_dofs; ++i)
+/* for (unsigned int i=0; i<base_element->total_dofs; ++i)
for (unsigned int j=0; j<base_element->total_dofs; ++j)
for (unsigned int n=0; n<n_sub_elements; ++n)
// only fill diagonals of the blocks
j*n_sub_elements + n)
= base_element->constraints()(i,j);
};
-};
+*/};
const unsigned int N) {
return FiniteElementData<1> (fe_data.dofs_per_vertex * N,
fe_data.dofs_per_line * N,
- fe_data.n_transform_functions);
+ fe_data.n_transform_functions,
+ fe_data.n_components * N);
};
#endif
return FiniteElementData<2> (fe_data.dofs_per_vertex * N,
fe_data.dofs_per_line * N,
fe_data.dofs_per_quad * N,
- fe_data.n_transform_functions);
+ fe_data.n_transform_functions,
+ fe_data.n_components * N);
};
#endif
template <int dim>
void FEValuesBase<dim>::get_function_values (const dVector &fe_function,
- vector<double> &values) const {
+ vector<double> &values) const
+{
+ Assert (fe.n_components == 1,
+ ExcWrongNoOfComponents());
Assert (selected_dataset<shape_values.size(),
ExcInvalidIndex (selected_dataset, shape_values.size()));
Assert (values.size() == n_quadrature_points,
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.templates.h>
#include <grid/geometry_info.h>
+#include <math.h>
/* Note: explicit instantiations at the end of the different sections! */
#include <numeric>
#include <algorithm>
+#include <math.h>