/**
* Number of components.
*/
- const unsigned int n_components;
+ const unsigned int n_components;
/**
* Default constructor. Constructs an element
* @author Wolfgang Bangerth, 1998
*/
template <int dim>
-struct FiniteElementBase :
+class FiniteElementBase :
public Subscriptor,
public FiniteElementData<dim>
/**
* Compute component and index from system index.
+ *
+ * Return value contains first
+ * component and second index in
+ * component.
*/
pair<unsigned,unsigned> system_to_component_index (unsigned index) const;
* space dimension.
*/
dFMatrix interface_constraints;
+
+ /**
+ * Map between linear dofs and component dofs.
+ */
+ vector< pair<unsigned, unsigned> > system_to_component_table;
};
DeclException0 (ExcJacobiDeterminantHasWrongSign);
};
+template <int dim>
+inline unsigned
+FiniteElementBase<dim>::component_to_system_index (unsigned component,
+ unsigned component_index) const
+{
+ Assert(false, ExcNotImplemented());
+ return component_index;
+}
-
+template <int dim>
+inline pair<unsigned,unsigned>
+FiniteElementBase<dim>::system_to_component_index (unsigned index) const
+{
+ Assert(index < system_to_component_table.size(), ExcInvalidIndex(index));
+ return system_to_component_table[index];
+}
/*---------------------------- fe.h ---------------------------*/
/* end of #ifndef __fe_H */
/*---------------------------- fe_lib.system.h ---------------------------*/
/* $Id$ */
-#ifndef __fe_system_H
-#define __fe_system_H
+#ifndef __deal_fe_system_H
+#define __deal_fe_system_H
/*---------------------------- fe_lib.system.h ---------------------------*/
#include <fe/fe.h>
-
+#include <vector>
+#include <pair.h>
/**
* This class provides an interface to group several equal elements together
template <int dim>
class FESystem : public FiniteElement<dim>
{
- /**
- * Copy constructor prohibited.
- */
- FESystem(const FESystem<dim>&);
-
+ /**
+ * Copy constructor prohibited.
+ */
+ FESystem(const FESystem<dim>&);
public:
const unsigned int face_no,
const unsigned int subface_no,
const vector<Point<dim-1> > &unit_points,
- vector<Point<dim> > &normal_vectors) const;
+ vector<Point<dim> > &normal_vectors) const;
- /**
- * Number of subelements of this object.
+ /**
+ *Number of different composing
+ * elements of this object.
+ *
* 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!
+ * multiplicity and subobjects
+ * themselves, this may be
+ * smaller than the total number
+ * of finite elements composed
+ * into this structure.
*/
- 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;
-
+ unsigned n_component_elements() const;
+
+ /**
+ * How often is a composing element used.
+ *
+ */
+ unsigned element_multiplicity(unsigned index) const;
+
+ /**
+ * Access to a composing 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# you get by
+ * this function. In a mixed
+ * discretization, you can choose
+ * the different base element types
+ * by index.
+ *
+ */
+ const FiniteElement<dim>& base_element(unsigned index) const;
+
private:
/**
- * Pointer to an object of the underlying
- * finite element class. This object is
- * created by the constructor.
+ * Pairs of multiplicity and element type.
+ */
+ typedef pair<const FiniteElement<dim> *, unsigned > ElementPair;
+
+ /**
+ * Pointer to underlying finite
+ * element classes.
+ *
+ * This object contains a pointer
+ * to each contributing element
+ * of a mixed discretization and
+ * its multiplicity. It is
+ * created by the constructor and
+ * constant afterwards.
*/
- const FiniteElement<dim> *const base_element;
+ vector< ElementPair > base_elements;
/**
* transformation from unit to real
* cell.
*/
- static FiniteElementData<dim> multiply_dof_numbers (const FiniteElementData<dim> &fe_data,
- const unsigned int N);
+ static FiniteElementData<dim>
+ multiply_dof_numbers (const FiniteElementData<dim> &fe_data,
+ const unsigned int N);
/**
* This function is simply singled out of
- * the constructor; it sets up the
+ * the constructor. It sets up the
+ * index table for the system as well as
* #restriction# and #prolongation#
* matrices. Since the operation of this
* function can be done without explicit
* the general template definition in
* the #.h# file.
*/
- void initialize_matrices ();
+ void initialize();
};
/* ------------------------- template functions ------------------------- */
-template <int dim>
-inline const FiniteElement<dim>&
-FESystem<dim>::get_base_element() const
+template<int dim>
+inline unsigned
+FESystem<dim>::n_component_elements() const
{
- return *base_element;
+ return base_elements.size();
}
-template <int dim>
+
+template<int dim>
inline unsigned
-FESystem<dim>::index(unsigned component, unsigned base) const
+FESystem<dim>::element_multiplicity(unsigned index) const
{
- return n_sub_elements * base + component;
+ return base_elements[index].second;
}
+
template <int dim>
-template <typename FE>
+inline const FiniteElement<dim>&
+FESystem<dim>::base_element(unsigned index) const
+{
+ return *base_elements[index].first;
+}
+
+
+template <int dim>
+template <class FE>
FESystem<dim>::FESystem (const FE &fe, const unsigned int n_elements) :
FiniteElement (multiply_dof_numbers(fe, n_elements)),
- n_sub_elements (n_elements),
- base_element (new FE())
+ base_elements(1)
{
- base_element->subscribe ();
- initialize_matrices ();
+ base_elements[0] = ElementPair(&fe, n_elements);
+ base_elements[0].first -> subscribe ();
+ initialize ();
};
-/*---------------------------- fe_lib.system.h ---------------------------*/
+/*---------------------------- fe_lib.system.h ---------------------------*/
/* end of #ifndef __fe_system_H */
#endif
-/*---------------------------- fe_lib.system.h ---------------------------*/
+/*---------------------------- fe_lib.system.h ---------------------------*/
/*------------------------------- FiniteElementBase ----------------------*/
-
#if deal_II_dimension == 1
template <>
};
interface_constraints.reinit (1,1);
interface_constraints(0,0)=1.;
+
+ for (unsigned int j=0 ; j<total_dofs ; ++j)
+ system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
};
#endif
};
interface_constraints.reinit (dofs_per_vertex+2*dofs_per_line,
2*dofs_per_vertex+dofs_per_line);
+ for (unsigned int j=0 ; j<total_dofs ; ++j)
+ system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
};
#endif
FiniteElement<dim>::FiniteElement (const FiniteElementData<dim> &fe_data) :
FiniteElementBase<dim> (fe_data) {};
-
-
-template <int dim>
-unsigned
-FiniteElementBase<dim>::component_to_system_index (unsigned component,
- unsigned component_index) const
-{
- Assert(false, ExcInvalidIndex(component));
- return component_index;
-}
-
-
-template <int dim>
-pair<unsigned,unsigned>
-FiniteElementBase<dim>::system_to_component_index (unsigned index) const
-{
- Assert(false, ExcInvalidIndex(index));
- return pair<unsigned,unsigned>(0,0);
-}
-
-
-
#if deal_II_dimension == 1
// declare this function to be explicitely specialized before first use
template <int dim>
FESystem<dim>::~FESystem ()
{
- base_element->unsubscribe ();
- delete base_element;
+ for (unsigned i=0;i<base_elements.size();++i)
+ {
+ base_elements[i].first -> unsubscribe ();
+ delete base_elements[i].first;
+ }
};
template <int dim>
-void FESystem<dim>::initialize_matrices ()
+void FESystem<dim>::initialize ()
{
+ // Initialize index table
+ // Multi-component base elements have to be thought of.
+ // 1. Vertices
+ unsigned total_index = 0;
+ for(unsigned comp = 0; comp < n_component_elements() ; ++comp)
+ {
+ for (unsigned local_index = 0 ;
+ local_index < base_element(comp).dofs_per_vertex ;
+ ++local_index)
+ {
+ system_to_component_table[total_index++] = pair<unsigned,unsigned>
+ (comp, 0);
+ }
+ }
+ // 2. Lines
+ for(unsigned comp = 0; comp < n_component_elements() ; ++comp)
+ {
+ for (unsigned local_index = 0 ;
+ local_index < base_element(comp).dofs_per_vertex ;
+ ++local_index)
+ {
+ system_to_component_table[total_index++] = pair<unsigned,unsigned>
+ (comp, 0);
+ }
+ }
+ // 3. Quads
+ for(unsigned comp = 0; comp < n_component_elements() ; ++comp)
+ {
+ for (unsigned local_index = 0 ;
+ local_index < base_element(comp).dofs_per_vertex ;
+ ++local_index)
+ {
+ system_to_component_table[total_index++] = pair<unsigned,unsigned>
+ (comp, 0);
+ }
+ }
+
+
// distribute the matrices of the base
// finite element to the matrices of
// this object
template <>
FiniteElementData<2>
FESystem<2>::multiply_dof_numbers (const FiniteElementData<2> &fe_data,
- const unsigned int N) {
+ const unsigned int N)
+{
return FiniteElementData<2> (fe_data.dofs_per_vertex * N,
fe_data.dofs_per_line * N,
fe_data.dofs_per_quad * N,
+/*
template <int dim>
double FESystem<dim>::shape_value (const unsigned int i,
- const Point<dim> &p) const {
+ const Point<dim> &p) const
+{
Assert((i<total_dofs), ExcInvalidIndex(i));
- return base_element->shape_value (i / n_sub_elements, p);
+// return base_element->shape_value (i / n_sub_elements, p);
};
const Point<dim> &p) const {
Assert((i<total_dofs), ExcInvalidIndex(i));
- return base_element->shape_grad (i / n_sub_elements, p);
+// return base_element->shape_grad (i / n_sub_elements, p);
};
const Point<dim> &p) const {
Assert((i<total_dofs), ExcInvalidIndex(i));
- return base_element->shape_grad_grad (i / n_sub_elements, p);
+// return base_element->shape_grad_grad (i / n_sub_elements, p);
};
base_element->get_normal_vectors (cell, face_no, subface_no, unit_points, normal_vectors);
};
-
+*/