const FiniteElement<dim> &fe2,
FullMatrix<number> &difference_matrix);
+ /**
+ * Compute the local
+ * $L^2$-projection matrix from
+ * fe1 to fe2.
+ */
template <int dim, typename number>
static void get_projection_matrix(const FiniteElement<dim> &fe1,
const FiniteElement<dim> &fe2,
const InVector& z1,
const DoFHandler<dim>& dof2,
const ConstraintMatrix& constraints2,
- OutVector& z1_difference);
-
+ OutVector& z1_difference);
+
+ /**
+ * $L^2$ projection for
+ * discontinuous
+ * elements. Operates the same
+ * direction as interpolate.
+ *
+ * The global projection can be
+ * computed by local matrices if
+ * the finite element spaces are
+ * discontinuous. With continuous
+ * elements, this is impossible,
+ * since a global mass matrix
+ * must be inverted.
+ */
+ template <int dim, class InVector, class OutVector>
+ static void project_dg (const DoFHandler<dim>& dof1,
+ const InVector& u1,
+ const DoFHandler<dim>& dof2,
+ OutVector& u2);
+
/**
* Gives the patchwise
* extrapolation of a @p{dof1}
}
+template <int dim, class InVector, class OutVector>
+void FETools::project_dg(const DoFHandler<dim> &dof1,
+ const InVector &u1,
+ const DoFHandler<dim> &dof2,
+ OutVector &u2)
+{
+ Assert(&dof1.get_tria()==&dof2.get_tria(), ExcTriangulationMismatch());
+ 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(u2.size()==dof2.n_dofs(), ExcDimensionMismatch(u2.size(), dof2.n_dofs()));
+
+ typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin_active();
+ typename DoFHandler<dim>::active_cell_iterator cell2 = dof2.begin_active();
+ typename DoFHandler<dim>::active_cell_iterator end = dof2.end();
+
+ const unsigned int n1 = dof1.get_fe().dofs_per_cell;
+ const unsigned int n2 = dof2.get_fe().dofs_per_cell;
+
+ Vector<double> u1_local(n1);
+ Vector<double> u2_local(n2);
+ std::vector<unsigned int> dofs(n2);
+
+ FullMatrix<double> matrix(n2,n1);
+ get_projection_matrix(dof1.get_fe(), dof2.get_fe(), matrix);
+
+ while (cell2 != end)
+ {
+ cell1->get_dof_values(u1, u1_local);
+ matrix.vmult(u2_local, u1_local);
+ cell2->get_dof_indices(dofs);
+ for (unsigned int i=0; i<n2; ++i)
+ {
+ u2(dofs[i])+=u2_local(i);
+ }
+
+ ++cell1;
+ ++cell2;
+ }
+}
+
+
template <int dim, class InVector, class OutVector>
void FETools::extrapolate(const DoFHandler<dim> &dof1,
const InVector &u1,
const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
Vector<double> &);
template
+void FETools::project_dg<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
+ const DoFHandler<deal_II_dimension> &, Vector<double> &);
+template
void FETools::extrapolate<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &, const Vector<double> &,
const DoFHandler<deal_II_dimension> &, Vector<double> &);
const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
Vector<float> &);
template
+void FETools::project_dg<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
+ const DoFHandler<deal_II_dimension> &, Vector<float> &);
+template
void FETools::extrapolate<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &, const Vector<float> &,
const DoFHandler<deal_II_dimension> &, Vector<float> &);
const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
BlockVector<double> &);
template
+void FETools::project_dg<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
+ const DoFHandler<deal_II_dimension> &, BlockVector<double> &);
+template
void FETools::extrapolate<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &, const BlockVector<double> &,
const DoFHandler<deal_II_dimension> &, BlockVector<double> &);
const DoFHandler<deal_II_dimension> &, const ConstraintMatrix &,
BlockVector<float> &);
template
+void FETools::project_dg<deal_II_dimension>
+(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
+ const DoFHandler<deal_II_dimension> &, BlockVector<float> &);
+template
void FETools::extrapolate<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &, const BlockVector<float> &,
const DoFHandler<deal_II_dimension> &, BlockVector<float> &);