From fcd6600fe1989c5b8f0c4eda6ae2c6e503da6f69 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Fri, 28 Nov 2003 15:03:27 +0000 Subject: [PATCH] local projection matrices git-svn-id: https://svn.dealii.org/trunk@8212 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_tools.cc | 93 +++++++++++++++++++++++---- tests/fe/Makefile | 3 +- tests/fe/fe_tools.cc | 92 ++++++++++++++++++++++++++ 3 files changed, 173 insertions(+), 15 deletions(-) create mode 100644 tests/fe/fe_tools.cc diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 23b76f20e0..2ac088be83 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -12,7 +12,7 @@ //---------------------------- fe_tools.cc --------------------------- -#include +#include #include #include #include @@ -314,32 +314,92 @@ void FETools::get_interpolation_difference_matrix (const FiniteElement &fe1 template void FETools::get_projection_matrix (const FiniteElement &fe1, const FiniteElement &fe2, - FullMatrix &interpolation_matrix) + FullMatrix &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 mass (n, n); Triangulation tr; GridGenerator::hyper_cube(tr); - DoFHandler dof; - dof.distribute_dofs(tr, fe2); - typename DoFHandler::active_cell_iterator cell = dof.begin(); + + DoFHandler dof1(tr); + dof1.distribute_dofs(fe1); + typename DoFHandler::active_cell_iterator cell1 = dof1.begin(); + DoFHandler dof2(tr); + dof2.distribute_dofs(fe2); + typename DoFHandler::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 quadrature(degree+1); + // Set up FEValues. + const UpdateFlags flags = update_values | update_q_points | update_JxW_values; + FEValues val1 (fe1, quadrature, update_values); + val1.reinit (cell1); + FEValues val2 (fe2, quadrature, flags); + val2.reinit (cell2); + + // Integrate and invert mass matrix + // This happens in the target space + FullMatrix mass (n2, n2); - // Set up FEValues. Here, we + for (unsigned int k=0;k b(n2); + Vector x(n2); + + for (unsigned int j=0;j const FiniteElement &, FullMatrix &); +template +void FETools::get_projection_matrix +(const FiniteElement &, + const FiniteElement &, + FullMatrix &); template void FETools::interpolate diff --git a/tests/fe/Makefile b/tests/fe/Makefile index 3fadc7e152..37499af31c 100644 --- a/tests/fe/Makefile +++ b/tests/fe/Makefile @@ -24,6 +24,7 @@ default: run-tests derivatives.exe : derivatives.g.$(OBJEXT) $(libraries) fe_data_test.exe : fe_data_test.g.$(OBJEXT) $(libraries) +fe_tools.exe : fe_tools.g.$(OBJEXT) $(libraries) fe_tools_test.exe : fe_tools_test.g.$(OBJEXT) $(libraries) fe_traits.exe : fe_traits.g.$(OBJEXT) $(libraries) internals.exe : internals.g.$(OBJEXT) $(libraries) @@ -52,7 +53,7 @@ shapes.exe : shapes.g.$(OBJEXT) $(libraries) transfer.exe : transfer.g.$(OBJEXT) $(libraries) up_and_down.exe : up_and_down.g.$(OBJEXT) $(libraries) -tests = fe_data_test fe_traits fe_tools_test mapping \ +tests = fe_data_test fe_traits fe_tools fe_tools_test mapping \ mapping_c1 shapes derivatives numbering mapping_q1_eulerian \ transfer internals \ dgq_1 \ diff --git a/tests/fe/fe_tools.cc b/tests/fe/fe_tools.cc new file mode 100644 index 0000000000..7fed85a08f --- /dev/null +++ b/tests/fe/fe_tools.cc @@ -0,0 +1,92 @@ +//-------------------------------------------------------------------------- +// $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 +#include + +#include + +#include +#include +#include +#include +#include + +#include + + +template +void test_projection (const FiniteElement& fe1, + const FiniteElement& 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 P(n2,n1); + + FETools::get_projection_matrix(fe1, fe2, P); + P.print_formatted(out, 2, false, 5); +} + + +template +void test_projection (std::ostream& out) +{ + FE_DGQ q0(0); + FE_DGQ q1(1); + FE_DGQ q2(2); + FE_DGQ q3(3); + FE_DGQ q4(4); + + FE_DGP p0(0); + FE_DGP p1(1); + FE_DGP p2(2); + FE_DGP p3(3); + FE_DGP 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); +} -- 2.39.5