From: Guido Kanschat Date: Tue, 8 Feb 2005 11:10:39 +0000 (+0000) Subject: new computation routine for embedding matrices X-Git-Tag: v8.0.0~14589 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bdc738058848072e5ea8726385388aa56b266c17;p=dealii.git new computation routine for embedding matrices git-svn-id: https://svn.dealii.org/trunk@9901 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index da4b702484..f1281424b2 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors +// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -160,6 +160,36 @@ class FETools const FiniteElement &fe2, FullMatrix &matrix); + /** + * Compute the embedding matrices + * from a coarse cell to + * 2dim child cells. + * + * This function computes the + * coarse grid function in a + * sufficiently large number of + * quadrature points and fits the + * fine grid functions using + * least squares + * approximation. Therefore, the + * use of this function is + * restricted to the case that + * the finite element spaces are + * actually nested. + * + * @arg fe The finite element + * class for which we compute the + * embedding matrices. + * @arg matrices A pointer to + * 2dim FullMatrix + * objects. This is the format + * used in FiniteElementBase, + * where we want to use ths + * function mostly. + */ + template + static void compute_embedding_matrices(const FiniteElement &fe, + FullMatrix* matrices); /** * Gives the interpolation of a the @@ -666,6 +696,14 @@ class FETools static std::pair *, unsigned int> get_fe_from_name_aux (const std::string &name); + + /** + * Exception thrown if an + * embedding matrix was computed + * inaccurately. + */ + DeclException1(ExcLeastSquaresError, double, + << "Least squares fit leaves a gap of " << arg1); }; /*@}*/ diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 12ed9265aa..0d2bf7d48d 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors +// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -13,6 +13,7 @@ #include +#include #include #include #include @@ -30,6 +31,7 @@ #include #include #include +#include #include #include #include @@ -469,6 +471,117 @@ void FETools::get_projection_matrix (const FiniteElement &fe1, } +template +void +FETools::compute_embedding_matrices(const FiniteElement& fe, + FullMatrix* matrices) +{ + const unsigned int nc = GeometryInfo::children_per_cell; + const unsigned int n = fe.dofs_per_cell; + const unsigned int degree = fe.degree; + + for (unsigned int i=0;i tr_coarse; + Triangulation tr_fine; + GridGenerator::hyper_cube (tr_coarse, 0, 1); + GridGenerator::hyper_cube (tr_fine, 0, 1); + tr_fine.refine_global(1); + DoFHandler dof_coarse(tr_coarse); + dof_coarse.distribute_dofs(fe); + DoFHandler dof_fine(tr_fine); + dof_fine.distribute_dofs(fe); + + MappingCartesian mapping; + QGauss q_fine(degree+1); + const unsigned int nq = q_fine.n_quadrature_points; + + FEValues fine (mapping, fe, q_fine, + update_q_points | update_JxW_values | update_values); + + typename DoFHandler::active_cell_iterator coarse_cell + = dof_coarse.begin_active(); + typename DoFHandler::active_cell_iterator fine_cell; + + /** + * We search for the polynomial on + * the small cell, being equal to + * the coarse polynomial in all + * quadrature points. + * + * First build the matrix for this + * least squares problem. This + * contains the values of the fine + * cell polynomials in the fine + * cell grid points. + * + * This matrix is the same for all + * children. + */ + fine.reinit(dof_fine.begin_active()); + FullMatrix A(nq, n); + for (unsigned int i=0;i v_coarse(nq); + Vector v_fine(n); + + unsigned int cell_number = 0; + for (fine_cell = dof_fine.begin_active(); + fine_cell != dof_fine.end(); + ++fine_cell, ++cell_number) + { + FullMatrix& matrix = matrices[cell_number]; + fine.reinit(fine_cell); + Quadrature q_coarse (fine.get_quadrature_points(), + fine.get_JxW_values()); + FEValues coarse (mapping, fe, q_coarse, update_values); + coarse.reinit(coarse_cell); + + // Compute this once for each + // coarse grid basis function + for (unsigned int i=0;i B = A; + const double result = B.least_squares(v_fine, v_coarse); + Assert (result < 1.e-10, ExcLeastSquaresError(result)); + + // Finally copy into the + // result matrix. Since the + // matrix maps a coarse + // grid function to a fine + // grid function, the + // columns are fine grid. + for (unsigned int j=0;j @@ -1659,6 +1772,10 @@ void FETools::get_projection_matrix const FiniteElement &, FullMatrix &); +template +void FETools::compute_embedding_matrices +(const FiniteElement &, FullMatrix*); + template void FETools::interpolate (const DoFHandler &, const Vector &, diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 6d2014b9ee..39c7a29841 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -56,6 +56,16 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

deal.II

    +
  1. + New: Function FETools::compute_embedding_matrices allows for + automatic computation of embedding matrices under the sole + assumption that the coarse grid space is embedded into the fine + grid space. In particular, no interpolation points are required. +
    + (GK, 2005/02/08) +

    +
  2. Fixed: Several wrong assertions in the Raviart-Thomas finite element class have been fixed, that were triggered when integrating face terms.