From c9c10241f75e8286cdeac70d033959c94ed5da7d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 19 Aug 2013 02:46:37 +0000 Subject: [PATCH] Implement interpolation between FE_Nothing and FE_Q. git-svn-id: https://svn.dealii.org/trunk@30341 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 +++ deal.II/source/fe/fe_q_base.cc | 109 ++++++++++++++++++--------------- 2 files changed, 67 insertions(+), 49 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 93d16b5ed8..50d1a2dd86 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -56,6 +56,13 @@ inconvenience this causes.

Specific improvements

    +
  1. + Fixed: SolutionTransfer used to crash whenever one transfered in the hp + context between cells that use FE_Nothing and FE_Q. This is now fixed. +
    + (Krzyszof Bzowski, Wolfgang Bangerth, 2013/08/18) +
    2. +
  2. Fixed: Under some circumstances (see http://code.google.com/p/dealii/issues/detail?id=82) the DoFTools::make_periodicity_constraints() function could create cycles in diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc index a43e3be528..25f5a9f134 100644 --- a/deal.II/source/fe/fe_q_base.cc +++ b/deal.II/source/fe/fe_q_base.cc @@ -459,10 +459,6 @@ FE_Q_Base:: get_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &interpolation_matrix) const { - // this is only implemented, if the source FE is also a Q element - AssertThrow ((dynamic_cast *>(&x_source_fe) != 0), - (typename FiniteElement::ExcInterpolationNotImplemented())); - Assert (interpolation_matrix.m() == this->dofs_per_cell, ExcDimensionMismatch (interpolation_matrix.m(), this->dofs_per_cell)); @@ -470,60 +466,75 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, ExcDimensionMismatch (interpolation_matrix.m(), x_source_fe.dofs_per_cell)); - // ok, source is a Q element, so we will be able to do the work - const FE_Q_Base &source_fe - = dynamic_cast&>(x_source_fe); + // go through the list of elements we can interpolate from + if (const FE_Q_Base *source_fe + = dynamic_cast*>(&x_source_fe)) + { + // ok, source is a Q element, so we will be able to do the work - // only evaluate Q dofs - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); - const unsigned int source_q_dofs_per_cell = Utilities::fixed_power(source_fe.degree+1); + // only evaluate Q dofs + const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int source_q_dofs_per_cell = Utilities::fixed_power(source_fe->degree+1); - // evaluation is simply done by evaluating the other FE's basis functions on - // the unit support points (FE_Q has the property that the cell - // interpolation matrix is a unit matrix, so no need to evaluate it and - // invert it) - for (unsigned int j=0; j p = this->unit_support_points[j]; + // evaluation is simply done by evaluating the other FE's basis functions on + // the unit support points (FE_Q has the property that the cell + // interpolation matrix is a unit matrix, so no need to evaluate it and + // invert it) + for (unsigned int j=0; j p = this->unit_support_points[j]; - // FE_Q element evaluates to 1 in unit support point and to zero in all - // other points by construction - Assert(std::abs(this->poly_space.compute_value (j, p)-1.)<1e-13, - ExcInternalError()); + // FE_Q element evaluates to 1 in unit support point and to zero in all + // other points by construction + Assert(std::abs(this->poly_space.compute_value (j, p)-1.)<1e-13, + ExcInternalError()); - for (unsigned int i=0; ipoly_space.compute_value (i, p); + } - // for FE_Q_DG0, add one last row of identity - if (q_dofs_per_cell < this->dofs_per_cell) - { - AssertDimension(source_q_dofs_per_cell+1, source_fe.dofs_per_cell); - for (unsigned int i=0; idofs_per_cell) + { + AssertDimension(source_q_dofs_per_cell+1, source_fe->dofs_per_cell); + for (unsigned int i=0; idegree*dim; - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdegree*dim; + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + if (std::fabs(interpolation_matrix(i,j)) < eps) + interpolation_matrix(i,j) = 0.; - // make sure that the row sum of each of the matrices is 1 at this - // point. this must be so since the shape functions sum up to 1 - for (unsigned int i=0; idofs_per_cell; ++i) - { - double sum = 0.; - for (unsigned int j=0; jdofs_per_cell; ++i) + { + double sum = 0.; + for (unsigned int j=0; jdofs_per_cell; ++j) + sum += interpolation_matrix(i,j); - Assert (std::fabs(sum-1) < eps, ExcInternalError()); + Assert (std::fabs(sum-1) < eps, ExcInternalError()); + } } + else if (const FE_Nothing *source_fe + = dynamic_cast*>(&x_source_fe)) + { + // the element we want to interpolate from is an FE_Nothing. this + // element represents a function that is constant zero and has no + // degrees of freedom, so the interpolation is simply a multiplication + // with a n_dofs x 0 matrix. there is nothing to do here + } + else + AssertThrow (false, + (typename FiniteElement::ExcInterpolationNotImplemented())); + } -- 2.39.5