/**
- * Number of components and dimension of the image space.
+ * Number of components and dimension of
+ * the image space.
*/
const unsigned int n_components;
+ /**
+ * This flag determines how the restriction
+ * of data from child cells to its mother
+ * is to be done. In this, it also
+ * determines in which way the restriction
+ * matrices of the derived class are to
+ * be used.
+ *
+ * For most elements, the mode is the
+ * following. Consider a 1d linear element,
+ * with two children and nodal values
+ * 1 and 2 on the first child, and 2 and 4
+ * on the second child. The restriction
+ * to the mother child then yields the
+ * values 1 and four, i.e. the values on
+ * the mother cell can be obtained by
+ * pointwise interpolation, where for each
+ * nodal value on the mother child one
+ * point on exactly one child exists.
+ * However, already on the quadratic
+ * element, the midpoint on the mother
+ * element can be obtained from any of
+ * the two children, which however would
+ * both yield the same value due to
+ * continuity. What we do in practice
+ * is to compute them from both sides
+ * and set them, rather than add them up.
+ * This makes some things much easier. In
+ * practice, if a degree of freedom on
+ * one of the child cells yields a
+ * nonzero contribution to one of the
+ * degrees of freedom on the mother
+ * cell, we overwrite the value on
+ * the mother cell. This way, when setting
+ * up the restriction matrices, we do not
+ * have to track which child is responsible
+ * for setting a given value on the mother
+ * cell. We call this the non-additive
+ * mode.
+ *
+ * The other possibility would be to
+ * add up the contributions from the
+ * different children. This would mean
+ * that both of the inner endpoint of
+ * the quadratic child elements above
+ * would have a weight of 1/2 with
+ * respect to the midpoint value on
+ * the mother cell. However, this also
+ * means that we have to first compute
+ * the restriction to the mother cell
+ * by addition from the child cells, and
+ * afterwards set them to the global
+ * vector. The same process, adding
+ * up the local contributions to the
+ * global vector is not possible since
+ * we do not know how many coarse cells
+ * contribute to nodes on the boundary.
+ *
+ * In contrast to the non-additive mode
+ * described above, which is the simplest
+ * way for elements can be interpolated
+ * from its children, interpolation is
+ * not possible for piecewise constant
+ * elements, to name only one example.
+ * Here, the value on the mother cell
+ * has to be taken the average of the
+ * values on the children, i.e. all
+ * children contribute alike to the
+ * one degree of freedom. Here, we have
+ * to sum up the contributions of all
+ * child cells with the same weight,
+ * and the non-additive mode of above
+ * would only set the value on the mother
+ * cell to the value of one of the child
+ * cell, irrespective of the values on the
+ * other cells.
+ *
+ * Similarly, for discontinuous linear
+ * elements, it might be better to not
+ * interpolate the values at the corners
+ * from the child cells, but to take a
+ * better average, for example
+ * interpolating at the centers of the
+ * child cells; in that case, the
+ * contributions of the child cells
+ * have to be additive as well.
+ *
+ * Given these notes, the flag under
+ * consideration has to be set to #false#
+ * for the usual continuous Lagrange
+ * elements, and #true# for the other
+ * cases mentioned above. The main function
+ * where it is used is
+ * #DoFAccessor::get_interpolated_dof_values#.
+ */
+ const bool restriction_is_additive;
+
/**
- * Default constructor. Constructs an element
+ * Default constructor. Constructs
+ * an element
* which is not so useful. Checking
* #total_dofs# is therefore a good way to
- * check if something went wrong. */
+ * check if something went wrong.
+ */
FiniteElementData ();
/**
FiniteElementData (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
const unsigned int n_transform_functions,
- const unsigned int n_components);
+ const unsigned int n_components,
+ const bool restriction_is_additive);
/**
* Constructor for a 2-dimensional
const unsigned int dofs_per_line,
const unsigned int dofs_per_quad,
const unsigned int n_transform_functions,
- const unsigned int n_components);
+ const unsigned int n_components,
+ const bool restriction_is_additive);
/**
* Constructor for a 3-dimensional
const unsigned int dofs_per_quad,
const unsigned int dofs_per_hex,
const unsigned int n_transform_functions,
- const unsigned int n_components);
+ const unsigned int n_components,
+ const bool restriction_is_additive);
/**
* Declare this destructor virtual in
unsigned int component_index) const;
/**
- * Compute component and index from system index.
+ * Compute component and index from
+ * system index.
*
* Return value contains first
* component and second index in
* are implemented here and do not have to be taken care of later.
*/
template <int dim>
-class FEQ1Mapping
- :
- public FiniteElement<dim>
+class FEQ1Mapping : public FiniteElement<dim>
{
public:
/**
* shall be zero.
*/
FEQ1Mapping (const unsigned int dofs_per_vertex,
- const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad=0,
- const unsigned int dofs_per_hex =0,
- const unsigned int n_components =1);
+ const unsigned int dofs_per_line,
+ const unsigned int dofs_per_quad =0,
+ const unsigned int dofs_per_hex =0,
+ const unsigned int n_components =1,
+ const bool restriction_is_additive=false);
/**
* Return the value of the #i#th shape
const unsigned int dofs_per_quad,
const unsigned int dofs_per_hex,
const unsigned int n_transform_functions,
- const unsigned int n_components) :
+ const unsigned int n_components,
+ const bool restriction_is_additive) :
dofs_per_vertex(dofs_per_vertex),
dofs_per_line(dofs_per_line),
dofs_per_quad(dofs_per_quad),
* dofs_per_line),
total_dofs (first_hex_index+dofs_per_hex),
n_transform_functions (n_transform_functions),
- n_components(n_components)
+ n_components(n_components),
+ restriction_is_additive(restriction_is_additive)
{
Assert(dim==3, ExcDimensionMismatch(3,dim));
};
const unsigned int dofs_per_line,
const unsigned int dofs_per_quad,
const unsigned int n_transform_functions,
- const unsigned int n_components) :
+ const unsigned int n_components,
+ const bool restriction_is_additive) :
dofs_per_vertex(dofs_per_vertex),
dofs_per_line(dofs_per_line),
dofs_per_quad(dofs_per_quad),
* dofs_per_line),
total_dofs (first_quad_index+dofs_per_quad),
n_transform_functions (n_transform_functions),
- n_components(n_components)
+ n_components(n_components),
+ restriction_is_additive(restriction_is_additive)
{
Assert(dim==2, ExcDimensionMismatch(2,dim));
};
FiniteElementData<dim>::FiniteElementData (const unsigned int dofs_per_vertex,
const unsigned int dofs_per_line,
const unsigned int n_transform_functions,
- const unsigned int n_components) :
+ const unsigned int n_components,
+ const bool restriction_is_additive) :
dofs_per_vertex(dofs_per_vertex),
dofs_per_line(dofs_per_line),
dofs_per_quad(0),
* dofs_per_line),
total_dofs (first_line_index+dofs_per_line),
n_transform_functions (n_transform_functions),
- n_components(n_components)
+ n_components(n_components),
+ restriction_is_additive(restriction_is_additive)
{
Assert(dim==1, ExcDimensionMismatch(1,dim));
};
(dofs_per_quad == f.dofs_per_quad) &&
(dofs_per_hex == f.dofs_per_hex) &&
(n_transform_functions == f.n_transform_functions) &&
- (n_components == f.n_components));
+ (n_components == f.n_components) &&
+ (restriction_is_additive == f.restriction_is_additive));
};
-#if deal_II_dimension == 1
-template <>
-FEDG_Q0<1>::FEDG_Q0 () :
- FEQ1Mapping<1> (0, 1)
+template <int dim>
+FEDG_Q0<dim>::FEDG_Q0 () :
+ FEQ1Mapping<dim> (0,
+ (dim==1 ? 1 : 0),
+ (dim==2 ? 1 : 0),
+ (dim==3 ? 1 : 0),
+ 1,
+ true)
{
- // for restriction and prolongation matrices:
- // note that we do not add up all the
- // contributions since then we would get
- // two summands per vertex in 1d (four
- // in 2d, etc), but only one per line dof.
- // We could accomplish for that by dividing
- // the vertex dof values by 2 (4, etc), but
- // would get into trouble at the boundary
- // of the domain since there only one
- // cell contributes to a vertex. Rather,
- // we do not add up the contributions but
- // set them right into the matrices!
-
- // The restriction matrices got crazy values
- // as it is yet not clear how they should work
- // in the DG(0) case. In general
- // the use of the restriction matrices
- // is not yet finally decided about, too.
- restriction[0](0,0) = 1e8;
- restriction[1](0,0) = 1e8;
-
- prolongation[0](0,0) = 1.0;
- prolongation[1](0,0) = 1.0;
+ for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
+ {
+ restriction[i](0,0) = 1./GeometryInfo<dim>::children_per_cell;
+ prolongation[i](0,0) = 1.0;
+ }
};
+#if deal_II_dimension == 1
+
+
template <>
void
FEDG_Q0<1>::get_face_support_points (const DoFHandler<1>::face_iterator &,
#endif
-
-
-#if deal_II_dimension == 2
-
-template <>
-FEDG_Q0<2>::FEDG_Q0 () :
- FEQ1Mapping<2> (0, 0, 1)
-{
- // The restriction matrices got crazy values
- // as it is yet not clear how they should work
- // in the DG(0) case. In general
- // the use of the restriction matrices
- // is not yet finally decided about, too.
- restriction[0](0,0) = 1e8;
- restriction[1](0,0) = 1e8;
- restriction[2](0,0) = 1e8;
- restriction[3](0,0) = 1e8;
-
- prolongation[0](0,0) = 1.0;
-
- prolongation[1](0,0) = 1.0;
-
- prolongation[2](0,0) = 1.0;
-
- prolongation[3](0,0) = 1.0;
-};
-
-
-
-#endif
-
-
-
-
template <int dim>
inline
double
template <>
FEQ1Mapping<1>::FEQ1Mapping (const unsigned int dofs_per_vertex,
- const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad,
- const unsigned int dofs_per_hex,
- const unsigned int n_components) :
+ const unsigned int dofs_per_line,
+ const unsigned int dofs_per_quad,
+ const unsigned int dofs_per_hex,
+ const unsigned int n_components,
+ const bool restriction_is_additive) :
FiniteElement<1> (FiniteElementData<1> (dofs_per_vertex,
dofs_per_line,
GeometryInfo<1>::vertices_per_cell,
- n_components))
+ n_components,
+ restriction_is_additive))
{
Assert (dofs_per_quad==0, ExcInvalidData());
Assert (dofs_per_hex==0, ExcInvalidData());
template <>
FEQ1Mapping<2>::FEQ1Mapping (const unsigned int dofs_per_vertex,
- const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad,
- const unsigned int dofs_per_hex,
- const unsigned int n_components) :
+ const unsigned int dofs_per_line,
+ const unsigned int dofs_per_quad,
+ const unsigned int dofs_per_hex,
+ const unsigned int n_components,
+ const bool restriction_is_additive) :
FiniteElement<2> (FiniteElementData<2> (dofs_per_vertex,
dofs_per_line,
dofs_per_quad,
GeometryInfo<2>::vertices_per_cell,
- n_components))
+ n_components,
+ restriction_is_additive))
{
Assert (dofs_per_hex == 0, ExcInvalidData());
};
template <>
FEQ1Mapping<3>::FEQ1Mapping (const unsigned int dofs_per_vertex,
- const unsigned int dofs_per_line,
- const unsigned int dofs_per_quad,
- const unsigned int dofs_per_hex,
- const unsigned int n_components) :
+ const unsigned int dofs_per_line,
+ const unsigned int dofs_per_quad,
+ const unsigned int dofs_per_hex,
+ const unsigned int n_components,
+ const bool restriction_is_additive) :
FiniteElement<3> (FiniteElementData<3> (dofs_per_vertex,
dofs_per_line,
dofs_per_quad,
dofs_per_hex,
GeometryInfo<3>::vertices_per_cell,
- n_components))
+ n_components,
+ restriction_is_additive))
{};