From 0c73d45889a8daabe0bd072f0bc1477683c0c933 Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Fri, 26 May 2000 16:25:03 +0000 Subject: [PATCH] add a new class FETools that provides interpolation, back_interpolation, and interpolation_difference functions and the respective local matrices. Remove the obsolete create_interpolation_matrix function in MatrixCreator. git-svn-id: https://svn.dealii.org/trunk@2959 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_tools.h | 206 +++++++++++ deal.II/deal.II/include/numerics/matrices.h | 21 -- deal.II/deal.II/source/fe/fe_tools.cc | 384 ++++++++++++++++++++ deal.II/deal.II/source/numerics/matrices.cc | 36 -- 4 files changed, 590 insertions(+), 57 deletions(-) create mode 100644 deal.II/deal.II/include/fe/fe_tools.h create mode 100644 deal.II/deal.II/source/fe/fe_tools.cc diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h new file mode 100644 index 0000000000..06ec845035 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -0,0 +1,206 @@ +//---------------------------- fe_tools.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000 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. +// +//---------------------------- fe_tools.h --------------------------- +#ifndef __deal2__fe_tools_H +#define __deal2__fe_tools_H + + + +template class FullMatrix; +template class FiniteElement; +template class DoFHandler; +template class Vector; + +#include + +/** + * This class performs interpolations and extrapolations of discrete + * functions of one #FiniteElement# #fe1# to another #FiniteElement# + * #fe2#. + * + * It also provides the local interpolation matrices that interpolate + * on each cell. Furthermore it provides the difference matrix + * $id-I_h$ that is needed for evaluating $(id-I_h)z$ for e.g. the + * dual solution $z$. + * + * @author Ralf Hartmann, 2000 + */ +class FETools +{ + public: + /** + * Gives the interpolation matrix + * that interpolates a #fe1#- + * function to a #fe2#-function on + * each cell. The interpolation_matrix + * needs to be of size + * #(fe2.dofs_per_cell, fe1.dofs_per_cell)#. + * + * Note, that if the finite element + * space #fe1# is a subset of + * the finite element space + * #fe2# than the #interpolation_matrix# + * is an embedding matrix. + */ + template + static void get_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix); + + /** + * Gives the interpolation matrix + * that interpolates a #fe1#- + * function to a #fe2#-function, and + * interpolates this to a second + * #fe1#-function on + * each cell. The interpolation_matrix + * needs to be of size + * #(fe1.dofs_per_cell, fe1.dofs_per_cell)#. + * + * Note, that this function only + * makes sense if the finite element + * space due to #fe1# is not a subset of + * the finite element space due to + * #fe2#, as if it were a subset then + * the #interpolation_matrix# would be + * only the unit matrix. + */ + template + static void get_back_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix); + + /** + * Gives the unit matrix minus the + * back interpolation matrix. + * The #difference_matrix# + * needs to be of size + * #(fe1.dofs_per_cell, fe1.dofs_per_cell)#. + * + * This function gives + * the matrix that transforms a + * #fe1# function $z$ to $z-I_hz$ + * where $I_h$ denotes the interpolation + * operator from the #fe1# space to + * the #fe2# space. This matrix hence + * is useful to evaluate + * error-representations where $z$ + * denotes the dual solution. + */ + template + static void get_interpolation_difference_matrix( + const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &difference_matrix); + + /** + * Gives the interpolation of a the + * #dof1#-function #u1# to a + * #dof2#-function #u2#. #dof1# and + * #dof2# need to be #DoFHandler#s + * based on the same triangulation. + * + * If the elements #fe1# and #fe2# + * are either both continuous or + * both discontinuous then this + * interpolation is the usual point + * interpolation. The same is true + * if #fe1# is a continuous and + * #fe2# is a discontinuous finite + * element. For the case that #fe1# + * is a discontinuous and #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 + static void interpolate(const DoFHandler &dof1, + const Vector &u1, + const DoFHandler &dof2, + Vector &u2); + + /** + * Gives the interpolation of the #fe1#- + * function #u1# to a #fe2#-function, and + * interpolates this to a second + * #fe1#-function named #u1_interpolated#. + * + * Note, that this function only + * makes sense if the finite element + * space due to #fe1# is not a subset of + * the finite element space due to + * #fe2#, as if it were a subset then + * #u1_interpolated# would be equal to #u1#. + */ + template + static void back_interpolate(const DoFHandler &dof1, + const Vector &u1, + const FiniteElement &fe2, + Vector &u1_interpolated); + + /** + * Gives $(Id-I_h)z2$ for a given + * #fe2#-function #z2#, where $I_h$ + * is the interpolation from #fe2# + * to #fe1#. $(Id-I_h)z2$ is + * denoted by #z2_difference#. + */ + template + static void interpolation_difference(const DoFHandler &dof1, + const Vector &z1, + const FiniteElement &fe2, + Vector &z1_difference); + + /** + * Gives the patchwise + * extrapolation of a #dof1# + * function #z1# to a #dof2# + * function #z2#. #dof1# and + * #dof2# need to be #DoFHandler# + * based on the same triangulation. + * + * This function is interesting for + * e.g. extrapolating patchwise a + * piecewise linear dual solution + * to a piecewise quadratic dual + * solution. + */ + template + static void extrapolate(const DoFHandler &dof1, + const Vector &z1, + const DoFHandler &dof2, + Vector &z2); + + + /** + * Exception + */ + DeclException0 (ExcTriangulationMismatch); + + /** + * Exception + */ + DeclException4 (ExcMatrixDimensionMismatch, + int, int, int, int, + << "This is a " << arg1 << "x" << arg2 << " matrix, " + << "but should be a " << arg1 << "x" << arg2 << " matrix."); +}; + + + +/*---------------------------- fe_tools.h ---------------------------*/ +/* end of #ifndef __deal2__fe_tools_H */ +#endif +/*---------------------------- fe_tools.h ---------------------------*/ diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index c12b39d506..b8c183e989 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -329,27 +329,6 @@ class MatrixCreator Vector &rhs_vector, const Function *a = 0); - /** - * Lagrange interpolation - * matrix for different - * elements. - * - * This function builds a matrix - * $A$ such that a function - * $u_{high}$ is interpolated to - * a function of lower order - * $u_{low}$ by cell-wise - * multiplication - * $u_{low} = A u_{high}$. - */ - static void create_interpolation_matrix(const FiniteElement &high, - const FiniteElement &low, - FullMatrix& result); - - /** - * Exception - */ - DeclException0 (ExcInvalidFE); /** * Exception */ diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc new file mode 100644 index 0000000000..f698b4f5e4 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -0,0 +1,384 @@ +//---------------------------- fe_tools.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000 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. +// +//---------------------------- fe_tools.cc --------------------------- + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +template +void FETools::get_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) +{ + 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(), + fe2.dofs_per_cell, + fe1.dofs_per_cell)); + + // Initialize FEValues for fe1 at + // the unit support points of the + // fe2 element. + vector phantom_weights(fe2.dofs_per_cell,1.); + vector > fe2_support_points (fe2.dofs_per_cell); + fe2.get_unit_support_points (fe2_support_points); + Quadrature fe2_support_points_quadrature(fe2_support_points, + phantom_weights); + + FEValues fe_values( + fe1, fe2_support_points_quadrature, update_values); + + + for (unsigned int i=0; i > unit_support_points (fe2.dofs_per_cell); + fe2.get_unit_support_points (unit_support_points); + + for (unsigned int i=0; i +void FETools::get_back_interpolation_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &interpolation_matrix) +{ + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(interpolation_matrix.m()==fe1.dofs_per_cell && + interpolation_matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(interpolation_matrix.m(), + interpolation_matrix.n(), + fe1.dofs_per_cell, + fe1.dofs_per_cell)); + + FullMatrix first_matrix (fe2.dofs_per_cell, fe1.dofs_per_cell); + FullMatrix second_matrix(fe1.dofs_per_cell, fe2.dofs_per_cell); + + get_interpolation_matrix(fe1, fe2, first_matrix); + get_interpolation_matrix(fe2, fe1, second_matrix); + + // int_matrix=second_matrix*first_matrix + second_matrix.mmult(interpolation_matrix, first_matrix); +} + + + +template +void FETools::get_interpolation_difference_matrix(const FiniteElement &fe1, + const FiniteElement &fe2, + FullMatrix &difference_matrix) +{ + Assert (fe1.n_components() == fe2.n_components(), + ExcDimensionMismatch(fe1.n_components(), fe2.n_components())); + Assert(difference_matrix.m()==fe1.dofs_per_cell && + difference_matrix.n()==fe1.dofs_per_cell, + ExcMatrixDimensionMismatch(difference_matrix.m(), + difference_matrix.n(), + fe1.dofs_per_cell, + fe1.dofs_per_cell)); + + FullMatrix interpolation_matrix(fe1.dofs_per_cell); + get_back_interpolation_matrix(fe1, fe2, interpolation_matrix); + + for (unsigned int i=0; i +void FETools::interpolate(const DoFHandler &dof1, + const Vector &u1, + const DoFHandler &dof2, + Vector &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())); + + const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell; + const unsigned int dofs_per_cell2=dof2.get_fe().dofs_per_cell; + + Vector u1_local(dofs_per_cell1); + Vector u2_local(dofs_per_cell2); + + FullMatrix interpolation_matrix(dofs_per_cell2, + dofs_per_cell1); + FETools::get_interpolation_matrix(dof1.get_fe(), dof2.get_fe(), + interpolation_matrix); + + DoFHandler::active_cell_iterator cell1 = dof1.begin_active(), + endc1 = dof1.end(), + cell2 = dof2.begin_active(), + endc2 = dof2.end(); + + vector index_multiplicity(dof2.n_dofs(),0); + vector dofs (dofs_per_cell2); + u2=0; + + for (; cell1!=endc1, cell2!=endc2; ++cell1, ++cell2) + { + cell1->get_dof_values(u1, u1_local); + interpolation_matrix.vmult(u2_local, u1_local); + cell2->get_dof_indices(dofs); + for (unsigned int i=0; i(index_multiplicity[i]); + } +} + + + +template +void FETools::back_interpolate(const DoFHandler &dof1, + const Vector &u1, + const FiniteElement &fe2, + Vector &u1_interpolated) +{ + Assert(dof1.get_fe().n_components() == fe2.n_components(), + 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())); + + const unsigned int dofs_per_cell1=dof1.get_fe().dofs_per_cell; + + Vector u1_local(dofs_per_cell1); + Vector u1_int_local(dofs_per_cell1); + + FullMatrix interpolation_matrix(dofs_per_cell1, dofs_per_cell1); + FETools::get_back_interpolation_matrix(dof1.get_fe(), fe2, + interpolation_matrix); + + DoFHandler::active_cell_iterator cell1 = dof1.begin_active(), + endc1 = dof1.end(); + + for (; cell1!=endc1; ++cell1) + { + cell1->get_dof_values(u1, u1_local); + interpolation_matrix.vmult(u1_int_local, u1_local); + cell1->set_dof_values(u1_int_local, u1_interpolated); + } +} + + + +template +void FETools::interpolation_difference(const DoFHandler &dof1, + const Vector &u1, + const FiniteElement &fe2, + Vector &u1_difference) +{ + Assert(dof1.get_fe().n_components() == fe2.n_components(), + 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 u1_local(dofs_per_cell1); + Vector u1_diff_local(dofs_per_cell1); + + FullMatrix difference_matrix(dofs_per_cell1, dofs_per_cell1); + FETools::get_interpolation_difference_matrix(dof1.get_fe(), fe2, + difference_matrix); + + DoFHandler::active_cell_iterator cell1 = dof1.begin_active(), + endc1 = dof1.end(); + + for (; cell1!=endc1; ++cell1) + { + cell1->get_dof_values(u1, u1_local); + difference_matrix.vmult(u1_diff_local, u1_local); + cell1->set_dof_values(u1_diff_local, u1_difference); + } +} + + +template +void FETools::extrapolate(const DoFHandler &dof1, + const Vector &u1, + const DoFHandler &dof2, + Vector &u2) +{ + Assert(dof1.get_fe().n_components() == dof2.get_fe().n_components(), + ExcDimensionMismatch(dof1.get_fe().n_components(), dof2.get_fe().n_components())); + 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())); + + Vector u3(dof2.n_dofs()); + interpolate(dof1, u1, dof2, u3); + + const unsigned int dofs_per_cell = dof2.get_fe().dofs_per_cell; + const Triangulation &tria=dof1.get_tria(); + Vector dof_values(dofs_per_cell); + + for (unsigned int level=0; level::cell_iterator cell=dof2.begin(level), + endc=dof2.end(level); + + for (; cell!=endc; ++cell) + { + if (!cell->active()) + { + bool active_children=false; + for (unsigned int child_n=0; + child_n::children_per_cell; ++child_n) + if (cell->child(child_n)->active()) + { + active_children=true; + break; + } + + if (active_children) + { + cell->get_interpolated_dof_values(u3, dof_values); + cell->set_dof_values_by_interpolation(dof_values, u2); + } + } + } + } +} + + + + +/*-------------- Explicit Instantiations -------------------------------*/ + +template +void FETools::get_interpolation_matrix(const FiniteElement &, + const FiniteElement &, + FullMatrix &); +template +void FETools::get_back_interpolation_matrix(const FiniteElement &, + const FiniteElement &, + FullMatrix &); +template +void FETools::get_interpolation_difference_matrix(const FiniteElement &, + const FiniteElement &, + FullMatrix &); +template +void FETools::interpolate(const DoFHandler &, + const Vector &, + const DoFHandler &, + Vector &); +template +void FETools::back_interpolate(const DoFHandler &, + const Vector &, + const FiniteElement &, + Vector &); +template +void FETools::interpolation_difference(const DoFHandler &, + const Vector &, + const FiniteElement &, + Vector &); +template +void FETools::extrapolate(const DoFHandler &, + const Vector &, + const DoFHandler &, + Vector &); + + +template +void FETools::get_interpolation_matrix(const FiniteElement &, + const FiniteElement &, + FullMatrix &); +template +void FETools::get_back_interpolation_matrix(const FiniteElement &, + const FiniteElement &, + FullMatrix &); +template +void FETools::get_interpolation_difference_matrix(const FiniteElement &, + const FiniteElement &, + FullMatrix &); +template +void FETools::interpolate(const DoFHandler &, + const Vector &, + const DoFHandler &, + Vector &); +template +void FETools::back_interpolate(const DoFHandler &, + const Vector &, + const FiniteElement &, + Vector &); +template +void FETools::interpolation_difference(const DoFHandler &, + const Vector &, + const FiniteElement &, + Vector &); +template +void FETools::extrapolate(const DoFHandler &, + const Vector &, + const DoFHandler &, + Vector &); + +/*---------------------------- fe_tools.cc ---------------------------*/ diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index aeebcbefac..6a9648598b 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -1368,42 +1368,6 @@ void LaplaceMatrix::assemble (Vector &rhs, -template -void -MatrixCreator::create_interpolation_matrix(const FiniteElement &high, - const FiniteElement &low, - FullMatrix& result) -{ - Assert (high.n_components() == low.n_components(), - ExcInvalidFE()); - - Assert (result.m() == low.dofs_per_cell, - ExcDimensionMismatch(result.m(), low.dofs_per_cell)); - Assert (result.n() == high.dofs_per_cell, - ExcDimensionMismatch(result.n(), high.dofs_per_cell)); - - - // Initialize FEValues at the support points - // of the low element. - vector phantom_weights(low.dofs_per_cell,1.); - vector > support_points(low.dofs_per_cell); - low.get_unit_support_points(support_points); - Quadrature low_points(support_points, - phantom_weights); - - FEValues fe(high, low_points, update_values); - - for (unsigned int i=0; i; -- 2.39.5