//---------------------------- fe_tools.cc ---------------------------
-#include <base/quadrature.h>
+#include <base/quadrature_lib.h>
#include <lac/full_matrix.h>
#include <lac/vector.h>
#include <lac/block_vector.h>
template <int dim, typename number>
void FETools::get_projection_matrix (const FiniteElement<dim> &fe1,
const FiniteElement<dim> &fe2,
- FullMatrix<number> &interpolation_matrix)
+ FullMatrix<number> &matrix)
{
+ Assert (fe1.n_components() == 1, ExcNotImplemented());
Assert (fe1.n_components() == fe2.n_components(),
ExcDimensionMismatch(fe1.n_components(), fe2.n_components()));
- Assert(interpolation_matrix.m()==fe2.dofs_per_cell &&
- interpolation_matrix.n()==fe1.dofs_per_cell,
- ExcMatrixDimensionMismatch(interpolation_matrix.m(),
- interpolation_matrix.n(),
+ Assert(matrix.m()==fe2.dofs_per_cell && matrix.n()==fe1.dofs_per_cell,
+ ExcMatrixDimensionMismatch(matrix.m(), matrix.n(),
fe2.dofs_per_cell,
fe1.dofs_per_cell));
+ matrix.clear();
- unsigned int m = fe1.dofs_per_cell;
- unsigned int n = fe2.dofs_per_cell;
+ unsigned int n1 = fe1.dofs_per_cell;
+ unsigned int n2 = fe2.dofs_per_cell;
// First, create a mass matrix. We
// cannot use MatrixTools, since we
// want to have a FullMatrix.
//
- // This happens in the target space
- FullMatrix<double> mass (n, n);
Triangulation<dim> tr;
GridGenerator::hyper_cube(tr);
- DoFHandler<dim> dof;
- dof.distribute_dofs(tr, fe2);
- typename DoFHandler<dim>::active_cell_iterator cell = dof.begin();
+
+ DoFHandler<dim> dof1(tr);
+ dof1.distribute_dofs(fe1);
+ typename DoFHandler<dim>::active_cell_iterator cell1 = dof1.begin();
+ DoFHandler<dim> dof2(tr);
+ dof2.distribute_dofs(fe2);
+ typename DoFHandler<dim>::active_cell_iterator cell2 = dof2.begin();
+
+ // Choose a quadrature rule
+ // Gauss is exact up to degree 2n-1
+ const unsigned int degree = std::max(fe1.tensor_degree(), fe2.tensor_degree());
+ Assert (degree != deal_II_numbers::invalid_unsigned_int,
+ ExcNotImplemented());
+
+ QGauss<dim> quadrature(degree+1);
+ // Set up FEValues.
+ const UpdateFlags flags = update_values | update_q_points | update_JxW_values;
+ FEValues<dim> val1 (fe1, quadrature, update_values);
+ val1.reinit (cell1);
+ FEValues<dim> val2 (fe2, quadrature, flags);
+ val2.reinit (cell2);
+
+ // Integrate and invert mass matrix
+ // This happens in the target space
+ FullMatrix<double> mass (n2, n2);
- // Set up FEValues. Here, we
+ for (unsigned int k=0;k<quadrature.n_quadrature_points;++k)
+ {
+ const double w = val2.JxW(k);
+ for (unsigned int i=0;i<n2;++i)
+ {
+ const double v = val2.shape_value(i,k);
+ for (unsigned int j=0;j<n2;++j)
+ mass(i,j) += w*v * val2.shape_value(j,k);
+ }
+ }
+ // Gauss-Jordan should be
+ // sufficient since we expect the
+ // mass matrix to be
+ // well-conditioned
+ mass.gauss_jordan();
+
+ // Now, test every function of fe1
+ // with test functions of fe2 and
+ // compute the projection of each
+ // unit vector.
+ Vector<double> b(n2);
+ Vector<double> x(n2);
+
+ for (unsigned int j=0;j<n1;++j)
+ {
+ b = 0.;
+ for (unsigned int i=0;i<n2;++i)
+ {
+ for (unsigned int k=0;k<quadrature.n_quadrature_points;++k)
+ {
+ const double w = val2.JxW(k);
+ const double u = val1.shape_value(j,k);
+ const double v = val2.shape_value(i,k);
+ b(i) += u*v*w;
+ }
+ }
+ // Multiply by the inverse
+ mass.vmult(x,b);
+ for (unsigned int i=0;i<n2;++i)
+ matrix(i,j) = x(i);
+ }
}
const FiniteElement<deal_II_dimension> &,
FullMatrix<float> &);
+template
+void FETools::get_projection_matrix<deal_II_dimension>
+(const FiniteElement<deal_II_dimension> &,
+ const FiniteElement<deal_II_dimension> &,
+ FullMatrix<double> &);
template
void FETools::interpolate<deal_II_dimension>
--- /dev/null
+//--------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//--------------------------------------------------------------------------
+
+// Test the cell matrices generated in FETools.
+
+#include <iostream>
+#include <fstream>
+
+#include <base/logstream.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_dgp.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_system.h>
+
+#include <fe/fe_tools.h>
+
+
+template<int dim>
+void test_projection (const FiniteElement<dim>& fe1,
+ const FiniteElement<dim>& fe2,
+ std::ostream& out)
+{
+ out << fe1.get_name() << " -> "
+ << fe2.get_name() << std::endl;
+
+ const unsigned int n1 = fe1.dofs_per_cell;
+ const unsigned int n2 = fe2.dofs_per_cell;
+
+ FullMatrix<double> P(n2,n1);
+
+ FETools::get_projection_matrix(fe1, fe2, P);
+ P.print_formatted(out, 2, false, 5);
+}
+
+
+template<int dim>
+void test_projection (std::ostream& out)
+{
+ FE_DGQ<dim> q0(0);
+ FE_DGQ<dim> q1(1);
+ FE_DGQ<dim> q2(2);
+ FE_DGQ<dim> q3(3);
+ FE_DGQ<dim> q4(4);
+
+ FE_DGP<dim> p0(0);
+ FE_DGP<dim> p1(1);
+ FE_DGP<dim> p2(2);
+ FE_DGP<dim> p3(3);
+ FE_DGP<dim> p4(4);
+
+ test_projection(p1,p0, out);
+ test_projection(p0,p1, out);
+ test_projection(p2,p1, out);
+ test_projection(p1,p2, out);
+ test_projection(p2,p0, out);
+ test_projection(p0,p2, out);
+ test_projection(p3,p2, out);
+ test_projection(p2,p3, out);
+ test_projection(p4,p3, out);
+ test_projection(p3,p4, out);
+
+ test_projection(q1,q0, out);
+ test_projection(q2,q0, out);
+ test_projection(q3,q0, out);
+ test_projection(q4,q0, out);
+ test_projection(q2,q1, out);
+ test_projection(q1,q2, out);
+ test_projection(q3,q2, out);
+ test_projection(q4,q3, out);
+}
+
+
+int main()
+{
+ std::ofstream logfile("fe_tools.output");
+
+ test_projection<1>(logfile);
+ test_projection<2>(logfile);
+ test_projection<3>(logfile);
+}