]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Handle constraints between FE_Q and FE_Nothing.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 24 Feb 2011 02:43:43 +0000 (02:43 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 24 Feb 2011 02:43:43 +0000 (02:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@23444 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/fe/fe_q.cc

index bd1751d11b5756da187e8743d5addd27f3f82435..f08b5fc4dd70a64b2e6272128ccf7d6e63eea2c7 100644 (file)
@@ -973,103 +973,106 @@ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe
                                  const unsigned int        subface,
                                  FullMatrix<double>       &interpolation_matrix) const
 {
-                                  // this is only implemented, if the
-                                  // source FE is also a
-                                  // Q element
-  typedef FE_Q<dim,spacedim> FEQ;
-  typedef FiniteElement<dim,spacedim> FEL;
-  AssertThrow ((x_source_fe.get_name().find ("FE_Q<") == 0)
-               ||
-               (dynamic_cast<const FEQ*>(&x_source_fe) != 0),
-               typename FEL::
-               ExcInterpolationNotImplemented());
-
-  Assert (interpolation_matrix.n() == this->dofs_per_face,
-         ExcDimensionMismatch (interpolation_matrix.n(),
-                               this->dofs_per_face));
   Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
          ExcDimensionMismatch (interpolation_matrix.m(),
                                x_source_fe.dofs_per_face));
 
-                                  // ok, source is a Q element, so
-                                  // we will be able to do the work
-  const FE_Q<dim,spacedim> &source_fe
-    = dynamic_cast<const FE_Q<dim,spacedim>&>(x_source_fe);
-
-                                  // Make sure, that the element,
-                                   // for which the DoFs should be
-                                   // constrained is the one with
-                                   // the higher polynomial degree.
-                                   // Actually the procedure will work
-                                   // also if this assertion is not
-                                   // satisfied. But the matrices
-                                   // produced in that case might
-                                  // lead to problems in the
-                                   // hp procedures, which use this
-                                  // method.
-  Assert (this->dofs_per_face <= source_fe.dofs_per_face,
-         typename FEL::
-         ExcInterpolationNotImplemented ());
-
-                                   // generate a point on this
-                                   // cell and evaluate the
-                                   // shape functions there
-  const Quadrature<dim-1>
-    quad_face_support (source_fe.get_unit_face_support_points ());
-
-                                  // Rule of thumb for FP accuracy,
-                                  // that can be expected for a
-                                  // given polynomial degree.
-                                  // This value is used to cut
-                                  // off values close to zero.
-  double eps = 2e-13*this->degree*(dim-1);
-
-                                  // compute the interpolation
-                                  // matrix by simply taking the
-                                   // value at the support points.
+                                  // see if source is a Q element
+  if (const FE_Q<dim,spacedim> *source_fe
+      = dynamic_cast<const FE_Q<dim,spacedim> *>(&x_source_fe))
+    {
+                                      // have this test in here since
+                                      // a table of size 2x0 reports
+                                      // its size as 0x0
+      Assert (interpolation_matrix.n() == this->dofs_per_face,
+             ExcDimensionMismatch (interpolation_matrix.n(),
+                                   this->dofs_per_face));
+
+                                      // Make sure, that the element,
+                                      // for which the DoFs should be
+                                      // constrained is the one with
+                                      // the higher polynomial degree.
+                                      // Actually the procedure will work
+                                      // also if this assertion is not
+                                      // satisfied. But the matrices
+                                      // produced in that case might
+                                      // lead to problems in the
+                                      // hp procedures, which use this
+                                      // method.
+      Assert (this->dofs_per_face <= source_fe->dofs_per_face,
+             (typename FiniteElement<dim,spacedim>::
+              ExcInterpolationNotImplemented ()));
+
+                                      // generate a point on this
+                                      // cell and evaluate the
+                                      // shape functions there
+      const Quadrature<dim-1>
+       quad_face_support (source_fe->get_unit_face_support_points ());
+
+                                      // Rule of thumb for FP accuracy,
+                                      // that can be expected for a
+                                      // given polynomial degree.
+                                      // This value is used to cut
+                                      // off values close to zero.
+      double eps = 2e-13*this->degree*(dim-1);
+
+                                      // compute the interpolation
+                                      // matrix by simply taking the
+                                      // value at the support points.
 //TODO: Verify that all faces are the same with respect to
 // these support points. Furthermore, check if something has to
 // be done for the face orientation flag in 3D.
-  const Quadrature<dim> subface_quadrature
-    = QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
-  for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
-    {
-      const Point<dim> &p = subface_quadrature.point (i);
+      const Quadrature<dim> subface_quadrature
+       = QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
+      for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
+       {
+         const Point<dim> &p = subface_quadrature.point (i);
 
-      for (unsigned int j=0; j<this->dofs_per_face; ++j)
+         for (unsigned int j=0; j<this->dofs_per_face; ++j)
+           {
+             double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p);
+
+                                              // Correct the interpolated
+                                              // value. I.e. if it is close
+                                              // to 1 or 0, make it exactly
+                                              // 1 or 0. Unfortunately, this
+                                              // is required to avoid problems
+                                              // with higher order elements.
+             if (std::fabs (matrix_entry - 1.0) < eps)
+               matrix_entry = 1.0;
+             if (std::fabs (matrix_entry) < eps)
+               matrix_entry = 0.0;
+
+             interpolation_matrix(i,j) = matrix_entry;
+           }
+       }
+
+                                      // 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 j=0; j<source_fe->dofs_per_face; ++j)
        {
-         double matrix_entry = this->shape_value (this->face_to_cell_index(j, 0), p);
+         double sum = 0.;
 
-                                          // Correct the interpolated
-                                          // value. I.e. if it is close
-                                          // to 1 or 0, make it exactly
-                                          // 1 or 0. Unfortunately, this
-                                          // is required to avoid problems
-                                          // with higher order elements.
-         if (std::fabs (matrix_entry - 1.0) < eps)
-           matrix_entry = 1.0;
-         if (std::fabs (matrix_entry) < eps)
-           matrix_entry = 0.0;
+         for (unsigned int i=0; i<this->dofs_per_face; ++i)
+           sum += interpolation_matrix(j,i);
 
-         interpolation_matrix(i,j) = matrix_entry;
+         Assert (std::fabs(sum-1) < 2e-13*this->degree*this->degree*dim,
+                 ExcInternalError());
        }
     }
-
-                                  // 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 j=0; j<source_fe.dofs_per_face; ++j)
+  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != 0)
     {
-      double sum = 0.;
-
-      for (unsigned int i=0; i<this->dofs_per_face; ++i)
-        sum += interpolation_matrix(j,i);
-
-      Assert (std::fabs(sum-1) < 2e-13*this->degree*this->degree*dim,
-              ExcInternalError());
+                                      // nothing to do here, the
+                                      // FE_Nothing has no degrees of
+                                      // freedom anyway
     }
+  else
+    AssertThrow (false,
+                (typename FiniteElement<dim,spacedim>::
+                 ExcInterpolationNotImplemented()));
 }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.