template <int dim> class DoFHandler;
template <typename number> class Vector;
template <int dim> class FE_Q;
+class ConstraintMatrix;
+
#include <base/config.h>
#include <base/exceptions.h>
* at the DoF values on the
* discontinuities.
*
- * Note that the resulting output
+ * Note that for continuous
+ * elements on grids with hanging
+ * nodes (i.e. locally refined
+ * grids) this function does not
+ * give the expected output.
+ * Indeed, the resulting output
* vector does not necessarily
* respect continuity
* requirements at hanging nodes:
* triangulation, although it is
* of course Q1 on each cell.
*
- * To make the field conforming
- * again, you have to have a
- * hanging node constraints
- * object, and call its
- * @p{distribute} function on the
- * result of this function.
+ * For this case (continuous
+ * elements on grids with hanging
+ * nodes), please use the
+ * @p{interpolate} function with
+ * an additional
+ * @p{ConstraintMatrix} argument,
+ * see below, or make the field
+ * conforming yourself by calling
+ * the @p{distribute} function of
+ * your hanging node constraints
+ * object.
*/
template <int dim, typename number>
static void interpolate (const DoFHandler<dim> &dof1,
const DoFHandler<dim> &dof2,
Vector<number> &u2);
+ /**
+ * Gives the interpolation of a
+ * the @p{dof1}-function @p{u1}
+ * to a @p{dof2}-function
+ * @p{u2}. @p{dof1} and @p{dof2}
+ * need to be @ref{DoFHandler}s
+ * based on the same
+ * triangulation.
+ * @p{constraints} is a hanging
+ * node constraints object
+ * corresponding to
+ * @p{dof2}. This object is
+ * particular important when
+ * interpolating onto continuous
+ * elements on grids with hanging
+ * nodes (locally refined grids).
+ *
+ * If the elements @p{fe1} and @p{fe2}
+ * are either both continuous or
+ * both discontinuous then this
+ * interpolation is the usual point
+ * interpolation. The same is true
+ * if @p{fe1} is a continuous and
+ * @p{fe2} is a discontinuous finite
+ * element. For the case that @p{fe1}
+ * is a discontinuous and @p{fe2} is
+ * a continuous finite element
+ * there is no point interpolation
+ * defined at the discontinuities.
+ * Therefore the meanvalue is taken
+ * at the DoF values on the
+ * discontinuities.
+ */
+ template <int dim, typename number>
+ static void interpolate (const DoFHandler<dim> &dof1,
+ const Vector<number> &u1,
+ const DoFHandler<dim> &dof2,
+ const ConstraintMatrix &constraints,
+ Vector<number> &u2);
+
+
/**
* Gives the interpolation of the
* @p{fe1}-function @p{u1} to a
* @p{fe1}-function named
* @p{u1_interpolated}.
*
- * Note, that this function only
- * makes sense if the finite
- * element space due to @p{fe1}
- * is not a subset of the finite
- * element space due to @p{fe2},
- * as if it were a subset then
- * @p{u1_interpolated} would be
- * equal to @p{u1}.
- *
* Note, that this function does
* not work on continuous
- * elements at hanging nodes.
+ * elements at hanging nodes. For
+ * that case use the
+ * @p{back_interpolate} function,
+ * below, that takes an
+ * additional
+ * @p{ConstraintMatrix} object.
+ *
+ * Furthermore note, that for the
+ * specific case when the finite
+ * element space corresponding to
+ * @p{fe1} is a subset of the
+ * finite element space
+ * corresponding to @p{fe2}, this
+ * function is simply an identity
+ * mapping.
*/
template <int dim, typename number>
static void back_interpolate (const DoFHandler<dim> &dof1,
const Vector<number> &u1,
const FiniteElement<dim> &fe2,
+ Vector<number> &u1_interpolated);
+
+ /**
+ * Gives the interpolation of the
+ * @p{dof1}-function @p{u1} to a
+ * @p{dof2}-function, and
+ * interpolates this to a second
+ * @p{dof1}-function named
+ * @p{u1_interpolated}.
+ * @p{constraints1} and
+ * @p{constraints2} are the
+ * hanging node constraints
+ * corresponding to @p{dof1} and
+ * @p{dof2}, respectively. These
+ * objects are particular
+ * important when continuous
+ * elements on grids with hanging
+ * nodes (locally refined grids)
+ * are involved.
+ *
+ * Furthermore note, that for the
+ * specific case when the finite
+ * element space corresponding to
+ * @p{dof1} is a subset of the
+ * finite element space
+ * corresponding to @p{dof2}, this
+ * function is simply an identity
+ * mapping.
+ */
+ template <int dim, typename number>
+ static void back_interpolate (const DoFHandler<dim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const Vector<number> &u1,
+ const DoFHandler<dim> &dof2,
+ const ConstraintMatrix &constraints2,
Vector<number> &u1_interpolated);
/**
*
* Note, that this function does
* not work on continuous
- * elements at hanging nodes.
+ * elements at hanging nodes. For
+ * that case use the
+ * @p{interpolation_difference}
+ * function, below, that takes an
+ * additional
+ * @p{ConstraintMatrix} object.
*/
template <int dim, typename number>
static void interpolation_difference(const DoFHandler<dim> &dof1,
const FiniteElement<dim> &fe2,
Vector<number> &z1_difference);
+ /**
+ * Gives $(Id-I_h)z1$ for a given
+ * @p{dof1}-function @p{z1}, where $I_h$
+ * is the interpolation from @p{fe1}
+ * to @p{fe2}. $(Id-I_h)z1$ is
+ * denoted by @p{z1_difference}.
+ * @p{constraints1} and
+ * @p{constraints2} are the
+ * hanging node constraints
+ * corresponding to @p{dof1} and
+ * @p{dof2}, respectively. These
+ * objects are particular
+ * important when continuous
+ * elements on grids with hanging
+ * nodes (locally refined grids)
+ * are involved.
+ */
+ template <int dim, typename number>
+ static void interpolation_difference(const DoFHandler<dim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const Vector<number> &z1,
+ const DoFHandler<dim> &dof2,
+ const ConstraintMatrix &constraints2,
+ Vector<number> &z1_difference);
+
/**
* Gives the patchwise
* extrapolation of a @p{dof1}
* Exception
*/
DeclException0 (ExcTriangulationMismatch);
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcHangingNodesNotAllowed);
/**
* Exception
#include <fe/fe_q.h>
#include <fe/fe_values.h>
#include <fe/mapping_q1.h>
+#include <dofs/dof_constraints.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
}
+template <int dim, typename number>
+void FETools::interpolate(const DoFHandler<dim> &dof1,
+ const Vector<number> &u1,
+ const DoFHandler<dim> &dof2,
+ Vector<number> &u2)
+{
+ interpolate(dof1, u1, dof2, ConstraintMatrix(), u2);
+}
+
template <int dim, typename number>
void FETools::interpolate(const DoFHandler<dim> &dof1,
const Vector<number> &u1,
const DoFHandler<dim> &dof2,
+ const ConstraintMatrix &constraints,
Vector<number> &u2)
{
Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch());
Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
Assert(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs()));
+ // for continuous elements on grids
+ // with hanging nodes we need
+ // hanging node
+ // constraints. Consequentely, if
+ // there are no constraints then
+ // hanging nodes are not allowed.
+ const bool hanging_nodes_not_allowed=
+ (dof2.get_fe().dofs_per_vertex != 0)
+ && (constraints.n_constraints()==0);
+
const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell;
const unsigned int dofs_per_cell2=dof2.get_fe().dofs_per_cell;
std::vector<unsigned int> index_multiplicity(dof2.n_dofs(),0);
std::vector<unsigned int> dofs (dofs_per_cell2);
u2.clear ();
-
+
for (; cell1!=endc1, cell2!=endc2; ++cell1, ++cell2)
{
+ if (hanging_nodes_not_allowed)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ Assert (cell1->at_boundary(face) ||
+ cell1->neighbor(face)->level() == cell1->level(),
+ ExcHangingNodesNotAllowed());
+
cell1->get_dof_values(u1, u1_local);
interpolation_matrix.vmult(u2_local, u1_local);
cell2->get_dof_indices(dofs);
++index_multiplicity[dofs[i]];
}
}
+
+ // when a discontinuous element is
+ // interpolated to a continuous
+ // one, we take the mean values.
for (unsigned int i=0; i<dof2.n_dofs(); ++i)
{
Assert(index_multiplicity[i]!=0, ExcInternalError());
u2(i) /= index_multiplicity[i];
}
-}
+ // Apply hanging node constraints.
+ constraints.distribute(u2);
+}
template <int dim, typename number>
ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components()));
Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
Assert(u1_interpolated.size()==dof1.n_dofs(),
- ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));
+ ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));
+
+ // For continuous elements on grids
+ // with hanging nodes we need
+ // hnaging node
+ // constraints. Consequently, when
+ // the elements are continuous no
+ // hanging node constraints are
+ // allowed.
+ const bool hanging_nodes_not_allowed=
+ (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0);
const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell;
Vector<number> u1_local(dofs_per_cell1);
Vector<number> u1_int_local(dofs_per_cell1);
+ typename DoFHandler<dim>::active_cell_iterator cell = dof1.begin_active(),
+ endc = dof1.end();
+
FullMatrix<number> interpolation_matrix(dofs_per_cell1, dofs_per_cell1);
FETools::get_back_interpolation_matrix(dof1.get_fe(), fe2,
interpolation_matrix);
-
- typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
- endc1 = dof1.end();
-
- for (; cell1!=endc1; ++cell1)
+ for (; cell!=endc; ++cell)
{
- cell1->get_dof_values(u1, u1_local);
+ if (hanging_nodes_not_allowed)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ Assert (cell->at_boundary(face) ||
+ cell->neighbor(face)->level() == cell->level(),
+ ExcHangingNodesNotAllowed());
+
+ cell->get_dof_values(u1, u1_local);
interpolation_matrix.vmult(u1_int_local, u1_local);
- cell1->set_dof_values(u1_int_local, u1_interpolated);
+ cell->set_dof_values(u1_int_local, u1_interpolated);
}
}
+
+template <int dim, typename number>
+void FETools::back_interpolate(const DoFHandler<dim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const Vector<number> &u1,
+ const DoFHandler<dim> &dof2,
+ const ConstraintMatrix &constraints2,
+ Vector<number> &u1_interpolated)
+{
+ // For discontinuous elements
+ // without constraints take the
+ // simpler version of the
+ // back_interpolate function.
+ if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0
+ && constraints1.n_constraints()==0 && constraints2.n_constraints()==0)
+ back_interpolate(dof1, u1, dof2.get_fe(), u1_interpolated);
+ else
+ {
+ Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(),
+ ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components()));
+ Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
+ Assert(u1_interpolated.size()==dof1.n_dofs(),
+ ExcDimensionMismatch(u1_interpolated.size(), dof1.n_dofs()));
+
+ // For continuous elements
+ // first interpolate to dof2,
+ // taking into account
+ // constraints2, and then
+ // interpolate back to dof1
+ // taking into account
+ // constraints1
+ Vector<double> u2(dof2.n_dofs());
+ interpolate(dof1, u1, dof2, constraints2, u2);
+ interpolate(dof2, u2, dof1, constraints1, u1_interpolated);
+ }
+}
+
+
template <int dim, typename number>
void FETools::interpolation_difference(const DoFHandler<dim> &dof1,
const Vector<number> &u1,
ExcDimensionMismatch(dof1.get_fe().n_components(), fe2.n_components()));
Assert(u1.size()==dof1.n_dofs(), ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
Assert(u1_difference.size()==dof1.n_dofs(),
- ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs()));
-
- const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell;
-
- Vector<number> u1_local(dofs_per_cell1);
- Vector<number> u1_diff_local(dofs_per_cell1);
+ ExcDimensionMismatch(u1_difference.size(), dof1.n_dofs()));
+
+ // For continuous elements on grids
+ // with hanging nodes we need
+ // hnaging node
+ // constraints. Consequently, when
+ // the elements are continuous no
+ // hanging node constraints are
+ // allowed.
+ const bool hanging_nodes_not_allowed=
+ (dof1.get_fe().dofs_per_vertex != 0) || (fe2.dofs_per_vertex != 0);
+
+ const unsigned int dofs_per_cell=dof1.get_fe().dofs_per_cell;
+
+ Vector<number> u1_local(dofs_per_cell);
+ Vector<number> u1_diff_local(dofs_per_cell);
- FullMatrix<number> difference_matrix(dofs_per_cell1, dofs_per_cell1);
+ FullMatrix<number> difference_matrix(dofs_per_cell, dofs_per_cell);
FETools::get_interpolation_difference_matrix(dof1.get_fe(), fe2,
difference_matrix);
- typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active(),
- endc1 = dof1.end();
+ typename DoFHandler<dim>::active_cell_iterator cell = dof1.begin_active(),
+ endc = dof1.end();
- for (; cell1!=endc1; ++cell1)
+ for (; cell!=endc; ++cell)
{
- cell1->get_dof_values(u1, u1_local);
+ if (hanging_nodes_not_allowed)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ Assert (cell->at_boundary(face) ||
+ cell->neighbor(face)->level() == cell->level(),
+ ExcHangingNodesNotAllowed());
+
+ cell->get_dof_values(u1, u1_local);
difference_matrix.vmult(u1_diff_local, u1_local);
- cell1->set_dof_values(u1_diff_local, u1_difference);
+ cell->set_dof_values(u1_diff_local, u1_difference);
}
}
+template <int dim, typename number>
+void FETools::interpolation_difference(const DoFHandler<dim> &dof1,
+ const ConstraintMatrix &constraints1,
+ const Vector<number> &u1,
+ const DoFHandler<dim> &dof2,
+ const ConstraintMatrix &constraints2,
+ Vector<number> &u1_difference)
+{
+ // For discontinuous elements
+ // without constraints take the
+ // cheaper version of the
+ // interpolation_difference function.
+ if (dof1.get_fe().dofs_per_vertex==0 && dof2.get_fe().dofs_per_vertex==0
+ && constraints1.n_constraints()==0 && constraints2.n_constraints()==0)
+ interpolation_difference(dof1, u1, dof2.get_fe(), u1_difference);
+ else
+ {
+ back_interpolate(dof1, constraints1, u1, dof2, constraints2, u1_difference);
+ u1_difference.sadd(-1, u1);
+ }
+}
+
+
template <int dim, typename number>
void FETools::extrapolate(const DoFHandler<dim> &dof1,
const Vector<number> &u1,
const DoFHandler<deal_II_dimension> &,
Vector<double> &);
template
+void FETools::interpolate<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &,
+ const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ Vector<double> &);
+template
void FETools::back_interpolate<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &,
const Vector<double> &,
const FiniteElement<deal_II_dimension> &,
Vector<double> &);
template
+void FETools::back_interpolate<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ Vector<double> &);
+template
void FETools::interpolation_difference<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &,
const Vector<double> &,
const FiniteElement<deal_II_dimension> &,
Vector<double> &);
template
+void FETools::interpolation_difference<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &,
+ const ConstraintMatrix &,
+ Vector<double> &);
+template
void FETools::extrapolate<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &,
const Vector<double> &,