Implement interpolation between FE_Nothing and FE_Q.

git-svn-id: https://svn.dealii.org/trunk@30341 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 93d16b5ed883ce946a1b95168c766419a9a66a09..50d1a2dd8678f5f10365bf10fd0b28e5ac76852c 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -56,6 +56,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>
+  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.
+  <br>
+  (Krzyszof Bzowski, Wolfgang Bangerth, 2013/08/18)
+  </li>
+
   <li>
   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 a43e3be528b3c054cf27988f54a2efea188ce3ec..25f5a9f13493b9e851f301984aefea5a9ccceeba 100644 (file)
--- 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<POLY,dim,spacedim>::
 get_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
                           FullMatrix<double>       &interpolation_matrix) const
 {
-  // this is only implemented, if the source FE is also a Q element
-  AssertThrow ((dynamic_cast<const FE_Q_Base<POLY,dim,spacedim> *>(&x_source_fe) != 0),
-               (typename FiniteElement<dim,spacedim>::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<dim,spacedim> &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<POLY,dim,spacedim> &source_fe
-    = dynamic_cast<const FE_Q_Base<POLY,dim,spacedim>&>(x_source_fe);
+  // go through the list of elements we can interpolate from
+  if (const FE_Q_Base<POLY,dim,spacedim> *source_fe
+      = dynamic_cast<const FE_Q_Base<POLY,dim,spacedim>*>(&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<dim>(this->degree+1);
-  const unsigned int source_q_dofs_per_cell = Utilities::fixed_power<dim>(source_fe.degree+1);
+      // only evaluate Q dofs
+      const unsigned int q_dofs_per_cell = Utilities::fixed_power<dim>(this->degree+1);
+      const unsigned int source_q_dofs_per_cell = Utilities::fixed_power<dim>(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<q_dofs_per_cell; ++j)
-    {
-      // read in a point on this cell and evaluate the shape functions there
-      const Point<dim> 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<q_dofs_per_cell; ++j)
+        {
+          // read in a point on this cell and evaluate the shape functions there
+          const Point<dim> 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; i<source_q_dofs_per_cell; ++i)
-        interpolation_matrix(j,i) = source_fe.poly_space.compute_value (i, p);
-    }
+          for (unsigned int i=0; i<source_q_dofs_per_cell; ++i)
+            interpolation_matrix(j,i) = source_fe->poly_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; i<source_q_dofs_per_cell; ++i)
-        interpolation_matrix(q_dofs_per_cell, i) = 0.;
-      for (unsigned int j=0; j<q_dofs_per_cell; ++j)
-        interpolation_matrix(j, source_q_dofs_per_cell) = 0.;
-      interpolation_matrix(q_dofs_per_cell, source_q_dofs_per_cell) = 1.;
-    }
+      // 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; i<source_q_dofs_per_cell; ++i)
+            interpolation_matrix(q_dofs_per_cell, i) = 0.;
+          for (unsigned int j=0; j<q_dofs_per_cell; ++j)
+            interpolation_matrix(j, source_q_dofs_per_cell) = 0.;
+          interpolation_matrix(q_dofs_per_cell, source_q_dofs_per_cell) = 1.;
+        }
 
-  // cut off very small values
-  const double eps = 2e-13*this->degree*dim;
-  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-    for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
-      if (std::fabs(interpolation_matrix(i,j)) < eps)
-        interpolation_matrix(i,j) = 0.;
+      // cut off very small values
+      const double eps = 2e-13*this->degree*dim;
+      for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+        for (unsigned int j=0; j<source_fe->dofs_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; i<this->dofs_per_cell; ++i)
-    {
-      double sum = 0.;
-      for (unsigned int j=0; j<source_fe.dofs_per_cell; ++j)
-        sum += interpolation_matrix(i,j);
+      // 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; i<this->dofs_per_cell; ++i)
+        {
+          double sum = 0.;
+          for (unsigned int j=0; j<source_fe->dofs_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<dim> *source_fe
+      = dynamic_cast<const FE_Nothing<dim>*>(&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<dim,spacedim>::ExcInterpolationNotImplemented()));
+
 }